1       subroutine sisvat(SnoMod,BloMod,jjtime)
 
  312       common/sisvat_loc_abc/seplab,fillab     
 
  316       common/sisvat_loc_num/nwunit            
 
  319       integer   ikl   ,isn   ,isl   ,ist      
 
  320       integer   ist__s,ist__w                 
 
  324       integer   IceMsk,IcIndx(
klonv)          
 
  327       real      drr_Ca,rrCa_n,drip            
 
  328       real      dsn_Ca,snCa_n,FallOK(
klonv)   
 
  329       real      roSMin,roSn_1,roSn_2,roSn_3   
 
  330       real      Dendr1,Dendr2,Dendr3          
 
  331       real      Spher1,Spher2,Spher3,Spher4   
 
  333       real      PorSno,Por_BS,Salt_f,PorRef   
 
  416       data      t__min / 200.00/              
 
  417       data      tapole / 263.15/              
 
  429       data      emisol /   0.99999999/        
 
  430       data      emiveg /   0.99999999/        
 
  431       data      emiwat /   0.99999999/        
 
  432       data      emisno /   0.99999999/        
 
  436       data      z0mbs0 /   0.5e-6/            
 
  439       data      z0m_s0/    0.00005/           
 
  462  600    
format(/,
'### SISVAT ERROR, Soil IR Upward  not defined ###',
 
  463      .         /,
'###               Initialize and Store IRs_SV ###')
 
  465           write(*,*)
'ikl',ikl,
'IR',
irs_sv(ikl)
 
  502      .       /,              a18,
'|           Grid Point ',2i4,
 
  504      .         
'    | Z0m =',f12.6,
' | Albedo = ',f6.3,
' |',
 
  505      .       /,
' -------+',7(
'---------+'),2(
'--------+'))
 
  583             rrca_n      = 
rrcasv(ikl)  +drr_ca       
 
  585             drip        = rrca_n       -
rrmxsv(ikl)  
 
  586             drip        =      max(
zer0,drip)        
 
  587             rrca_n      = rrca_n       -drip         
 
  599             snca_n      = 
sncasv(ikl)  +dsn_ca       
 
  600             drip        = snca_n       -
rrmxsv(ikl)  
 
  601             drip        =      max(
zer0,drip)        
 
  602             snca_n      = snca_n       -drip         
 
  681             bufs_n      =      
bufssv(ikl) +d_bufs      
 
  688  6601       
format(/,
'Buffer *: ',3e15.6)
 
  696             buf_ro      =  max( rosmin,                 
 
  698      .            +rosn_3*sqrt( 
vv__sv(ikl)))           
 
  704             bros_n      = (1. - polair) *   buf_ro      
 
  723             bros_n      = bros_n     * (1.0-fallok(ikl))
 
  724      .                  + 300.       *      fallok(ikl) 
 
  728             brossv(ikl) =(bros_n     *      d_bufs      
 
  738  6602       
format(
'rho    *: ',3e15.6,
'    dsnbSV: ',e15.6)
 
  743      .               min(dendr1*
vv__sv(ikl)-dendr2,     
 
  745             buf_g2      =  min( spher4,                 
 
  746      .               max(spher1*
vv__sv(ikl)+spher2,     
 
  748             buf_g1      = (1. - polair) *   buf_g1      
 
  750             buf_g2      = (1. - polair) *   buf_g2      
 
  760  6603       
format(
'G1,G2  *: ',3e15.6)
 
  821             bg1__n      =((1. - fallok(ikl))*   g1      
 
  822      .                        + fallok(ikl) *   99.)    
 
  823      .                  *       d_bufs/max(
eps6,d_bufs) 
 
  824             bg2__n      =((1. - fallok(ikl))*   g2      
 
  825      .                        + fallok(ikl) *   30.)    
 
  826      .                  *       d_bufs/max(
eps6,d_bufs) 
 
  839  6604       
format(
'G1,G2 F*: ',3e15.6,
'    T__Veg: ',e15.6)
 
  854             zronew  =(    typ__1  *d_bufs                 
 
  855      .              + (1.-typ__1) *
bufssv(ikl))           
 
  857             g1_new  =     typ__1  *buf_g1                 
 
  858      .              + (1.-typ__1) *
bg1ssv(ikl)            
 
  859             g2_new  =     typ__1  *buf_g2                 
 
  860      .              + (1.-typ__1) *
bg2ssv(ikl)            
 
  861             zroold  =((1.-typ__1) *d_bufs                 
 
  864             g1_old  = (1.-typ__1) *buf_g1                 
 
  866             g2_old  = (1.-typ__1) *buf_g2                 
 
  875             siz_av  =   ( zronew  *siznew+zroold*sizold)  
 
  876             sph_av = min( zronew  *sphnew+zroold*sphold   
 
  878             den_av = min((siz_av  - (    sph_av *
dscdsv    
  894             g1diff  =(   -dendok *den_av
 
  895      .               +(1.-dendok)*sph_av) *
g1_dsv 
  896             g2diff  =     dendok *sph_av  *
g1_dsv 
  897      .               +(1.-dendok)*siz_av
 
  899      .               +(1.-sameok)*g1diff
 
  901      .               +(1.-sameok)*g2diff
 
  904      .                  *       bufs_n/max(
eps6,bufs_n) 
 
  906      .                  *       bufs_n/max(
eps6,bufs_n) 
 
  915  6605       
format(
'B1,Typ  : ',2e15.6,11
x,
'OK,Den,Sph,Siz: ',4e15.6
 
  916      .          ,/,
'          ',30
x   ,11
x,
'sam,dif,G1    : ',3e15.6)
 
  941  6606       
format(
'G1,G2 N*: ',2e15.6,i15,e27.6)
 
  962  6004   
format(i3,
'  dsn+Buf=',f6.2,6
x,
'z dz *ro =',10f6.2,
 
  973      .                      + bdzssv(ikl)     *    
nlaysv(ikl) 
 
 1001             z_snsv(ikl)     = z_snsv(ikl) + 
dzsnsv(ikl,isn)
 
 1002             zzsnsv(ikl,isn) = z_snsv(ikl)
 
 1038  6006   
format(i3,
'  dsn+Buf=',f6.2,6
x,
'* dz *ro =',10f6.2,
 
 1104             lsnmsk     =     min( 1,
isnosv(ikl))
 
 1106             evg_sv(ikl)=  emiveg*(1-lsnmsk)+emisno*lsnmsk  
 
 1107             eso_sv(ikl)=  emisol*(1-lsnmsk)+emisno*lsnmsk  
 
 1111      .    + emiwat                 *(1-
lsmask(ikl)))*(1-lsnmsk)
 
 1132           snowat      =                     min(
isnosv(ikl),0)
 
 1150           icindx(ikl) = max(icindx(ikl),
 
 1153      .                                   int(
ro__sv(ikl,isn)-900.))))
 
 1158           lismsk    =     min(
iicesv(ikl),1     )
 
 1159           lismsk    =     max(
lsmask(ikl),lismsk)
 
 1160           icemsk    =     max(0,sign(1   ,icindx(ikl)-1)  )
 
 1191           growth      =min(max(0,7-
ivgtsv(ikl)),1)
 
 1204           z0mlnd      =max( z0mlnd   ,
 
 1206      .                     +z0_ice   *   icemsk )
 
 1209      .                     -
zzsnsv(ikl,max(icindx(ikl),0)))/7. 
 
 1210           z0mlnd      =max( z0mlnd    ,    5.e-5  )  
 
 1221           z0mbsn =       u2star      *0.536e-3   -  61.8e-6
 
 1222           z0mbsn =   max(z0mbs0      ,z0mbsn)
 
 1311      .                  +z0msea      *(1-lismsk)             
 
 1372           z0hnsv(ikl) =     z0hsea             *(1-lismsk) 
 
 1389  6661     
format(20
x,7f9.6)
 
 1399  6600     
format(/,
' ** SISVAT     *0  ' 
 1400      .            ,
'  za__SV  = ',e12.4,
'  Z0m_SV  = ',e12.4
 
 1401      .            ,
'  sqrCm0  = ',e12.4,
'  Za/Z0m  = ',e12.4
 
 1403      .            ,
'  Z0SaSV  = ',e12.4,
'  Z0h_SV  = ',e12.4
 
 1404      .            ,
'  sqrCh0  = ',e12.4,
'  Za/Z0h  = ',e12.4)
 
 1531  6007   
format(i3,
'  dsn+Buf=',f6.2,6
x,
'q dz *ro =',10f6.2,
 
 1543             z_snsv(ikl)     = z_snsv(ikl) + 
dzsnsv(ikl,isn)
 
 1544             zzsnsv(ikl,isn) = z_snsv(ikl)
 
 1553             z_snsv(ikl)   = max(
zer0,
 
 1590           irupsv(ikl) =      irupsv(ikl)                   
 
 1592           iru_sv(ikl) =     -irupsv(ikl)                   
 
 1649      .         
'        |Net Solar| IR Down | IR Up   | HS/Dwn=+|',
 
 1650      .          
' HL/Dwn=+| Temper. |         |  Snow  |  Rain  |',
 
 1651      .       /,
'        | [W/m2]  | [W/m2]  | [W/m2]  | [W/m2]  |',
 
 1652      .          
' [W/m2]  | [K]     |         | [mm/h] | [mm/h] |',
 
 1653      .       /,
' -------+',7(
'---------+'),2(
'--------+'),
 
 1654      .       /,
' SISVAT |',f8.1,
' |',f8.1,
' |',f8.1,
' |',f8.1,
' |',
 
 1655      .                     f8.1,
' |A',f7.2,
' |', 8
x ,
' |',2(f7.2,
' |'),
 
 1656      .       /,
' Canopy |',f8.1,
' |', 8
x ,
' |',f8.1,
' |',f8.1,
' |',
 
 1657      .                     f8.1,
' |',f8.2,
' |', 8
x ,
' |',2( 7
x ,
' |')
 
 1658      .       /,
' Soil   |',f8.1,
' |', 8
x ,
' |', 8
x ,
' |',f8.1,
' |',
 
 1659      .                     f8.1,
' |',f8.2,
' |', 8
x ,
' |',2( 7
x ,
' |'))
 
 1680      .           
' -----------------+-------------------+', 
 
 1681      .            
'-----------------+-+-----------------+', 
 
 1682      .          
'-------------------+',                     
 
 1683      .         /,
' SOIL/SNOW/VEGET. |                   |', 
 
 1684      .            
' Power,  Forcing |                   |', 
 
 1689      .         /,
'                  |',   11
x ,
'        |', 
 
 1690      .                f9.2,
' [W/m2] |',   11
x ,
'        |', 
 
 1692      .         /,
' -----------------+-------------------+', 
 
 1693      .            
'-----------------+-------------------+', 
 
 1694      .          
'-------------------+',                     
 
 1695      .         /,
' SOIL/SNOW  (TSo) | Energy/dt, Time 0 |', 
 
 1696      .            
' Power,  Forcing |   Sum Tim.0+Forc. |', 
 
 1697      .          
' Energy/dt, Time 1 |',                     
 
 1701      .         /,
'                  |',  f11.2,
' [W/m2] |', 
 
 1702      .                f9.2,
' [W/m2] |',  f11.2,
' [W/m2] |', 
 
 1703      .               f11.2,
' [W/m2] |',                     
 
 1704      .         /,
' -----------------+-------------------+', 
 
 1705      .            
'-----------------+-------------------+', 
 
 1706      .          
'-------------------+',                     
 
 1707      .         /,
'      SNOW  (qSn) | Energy/dt, Time 0 |', 
 
 1708      .            
' Power,  Excess  |   D(Tim.1-0-Forc.)|', 
 
 1709      .          
' Energy/dt, Time 1 |',                     
 
 1713      .         /,
'                  |',  f12.2, 
'[W/m2] |', 
 
 1714      .                f9.2,
' [W/m2] |',  f11.2,
' [W/m2] |', 
 
 1715      .               f12.2, 
'[W/m2] | ',                    
 
 1716      .         /,
' -----------------+-------------------+', 
 
 1717      .            
'-----------------+-------------------+', 
 
 1718      .          
'-------------------+')                     
 
 1732  6001       
format(a18,3i4,
' (EB1'           ,f15.6,
 
 1733      .                
')  - [(EB0           ',f15.6,
')',
 
 1734      .               /,55
x,
'+(ATM->Snow/Soil',f15.6,
')] ',
 
 1735      .                     
'= EBAL'          ,f15.6,
' [W/m2]',
 
 1736      .               /,55
x,
' (ATM->SISVAT'   ,f18.6,
 
 1737      .               /,55
x,
'- Veg. ImBal.',   f18.6,
')  ',
 
 1738      .               /,55
x,
'- ATM->SnoSol',   f18.6,
')  ',
 
 1739      .                     
'= ????'          ,f15.6,
' [W/m2]')
 
 1760  5010     
format(
' SNOW             |   Snow,   Time 0  |',
 
 1761      .            
' Snow,   Forcing |           Sum     |',
 
 1762      .          
'   Snow,   Time 1  |',
 
 1766      .         /,
'                  |',    f13.3,
' [mm] |',
 
 1767      .           
' A',  f9.3,
' [mm] |',    f13.3,
' [mm] |',
 
 1769      .         /,
'                  |',     13
x ,
'      |',
 
 1770      .           
' E',  f9.3,
' [mm] |',     13
x ,
'      |',
 
 1772      .         /,
'                  |',     13
x ,
'      |',
 
 1773      .           
' S',  f9.3,
' [mm] |',     13
x ,
'      |',
 
 1775      .         /,
'                  |',     13
x ,
'      |',
 
 1776      .           
'(M',  f9.3,
' [mm])|  (included in A)  |',
 
 1778      .         /,
'                  |',     13
x ,
'      |',
 
 1779      .           
' R',  f9.3,
' [mm] |',     13
x ,
'      |',
 
 1784      .         /,
' -----------------+-------------------+',
 
 1785      .            
'-----------------+-------------------+',
 
 1786      .          
'-------------------+')
 
 1803  6010       
format(a18,3i4,
' (MB1'        ,f12.6,
 
 1804      .                 
') - [(MB0        ',f12.6,        15
x,
')',
 
 1805      .               /,51
x,
'+(ATM Forcing',f12.6,
' - ',f12.6,
')',
 
 1806      .               /,51
x,
'+(BLS Forcing',f12.6,
' - ',f12.6,
')',
 
 1807      .               /,51
x,
'-(Depo/Sublim',f12.6,        15
x,
')',
 
 1808      .               /,51
x,
' !Melting    ',f12.6,
'  included in A!',
 
 1809      .               /,51
x,
'+(Run  OFF   ',f12.6,        15
x,
')',
 
 1811      .               /,29
x,
'= *BAL'       ,f12.6,      
' [mm w.e.]')
 
 1829  5003     
format(
' SOIL/SNOW  (qSo) |   Water,  Time 0  |',
 
 1830      .            
' Water,  Forcing |           Sum     |',
 
 1831      .          
'   Water,  Time 1  |',
 
 1835      .         /,
'                  |',    f13.3,
' [mm] |',
 
 1836      .                 f11.3,
' [mm] |',    f13.3,
' [mm] |',
 
 1838      .         /,
' -----------------+-------------------+',
 
 1839      .            
'-----------------+-------------------+',
 
 1840      .          
'-------------------+',
 
 1841      .         /,
' SOIL/SNOW/VEGET. |   Water,  Time 0  |',
 
 1842      .            
' Water,  Forcing |           Sum     |',
 
 1843      .          
'   Water,  Time 1  |',
 
 1847      .         /,
'                  |',    f13.3,
' [mm] |',
 
 1848      .                 f11.3,
' [mm] |',    f13.3,
' [mm] |',
 
 1850      .         /,
' -----------------+-------------------+',
 
 1851      .            
'-----------------+-------------------+',
 
 1852      .          
'-------------------+')
 
 1865  6002       
format(30
x,
' NEW Soil Water',3
x,
' Canopy   Water',3
x,
 
 1866      .                 
' OLD SVAT Water',4
x,
' FRC SVAT Water',
 
 1867      .           /,a18,3i4,f15.6,
' + ' ,f15.6,
' - ' ,f15.6,
 
 1868      .                           
' -  ',f15.6,
'    ', 15
x ,
'    ',
 
 1869      .      /,31
x,
'= ',f12.6,
' [mm] (Water Balance)',
 
 1870      .           /,30
x,
' NEW Soil Water',3
x,
'               ',3
x,
 
 1871      .                 
' OLD Soil Water',4
x,
' FRC Soil Water',
 
 1872      .           /,30
x,f15.6,
'   ' , 15
x ,
' - ' ,f15.6,
 
 1873      .                       
' -  ',f15.6,
'    ', 15
x ,
'    ',
 
 1874      .      /,31
x,
'= ',f12.6,
' [mm] (3 terms SUM)')
 
 1884  5004       
format(
' -----+--------+--+-----+--------+----+---+',
 
 1885      .  
'--------+----+---+--------+------+-+--------+--------+',
 
 1886      .           /,
'    n |     z  |     dz |     ro |    eta |',
 
 1887      .  
'     T  |     G1 |     G2 | Extinc |        | HISTORY|',
 
 1888      .           /,
'      |    [m] |    [m] | [kg/m3]| [m3/m3]|',
 
 1889      .  
'    [K] |    [-] |    [-] |    [-] |        |   [-]  |',
 
 1890      .           /,
' -----+--------+--------+--------+--------+',
 
 1891      .  
'--------+--------+--------+--------+--------+--------+')
 
 1893  5005       
format(
'      |        |        |        |W',f6.3,
' |',
 
 1894      .  
'        |        |        |A',f6.3,
' |        |        |')
 
 1902  5015       
format((i5,
' |',2(f7.3,
' |'),               f7.1,
' |',
 
 1903      .           f7.3,
' |' ,  f7.2,
' |', 2(f7.1,
' |'),  f7.3,
' |',
 
 1904      .            7
x ,
' |' ,  i5,
'   |'                          ))
 
 1906  5006       
format(
' -----+--------+--------+--------+--------+',
 
 1907      .  
'--------+--------+--------+--------+--------+--------+')
 
 1911  5007       
format(
' Brgh |',4(8
x,
'|'),  f7.2,
' | [micm] |',4(8
x,
'|'),
 
 1912      .           /,
' VEGE |',4(8
x,
'|'),2(f7.2,
' |'),        2(8
x,
'|'),
 
 1913      .                                   f7.3,
' |',           8
x,
'|' )
 
 1915  5014       
format(
' -----+--------+--------+--------+--------+',
 
 1916      .  
'--------+--------+--------+--------+--------+--------+',
 
 1917      .           /,
'    n |        |     dz |        |    eta |',
 
 1918      .  
'     T  |        |        |        | Root W.| W.Flow |',
 
 1919      .           /,
'      |        |    [m] |        | [m3/m3]|',
 
 1920      .  
'    [K] |        |        |        | [mm/d] | [mm/h] |',
 
 1921      .           /,
' -----+--------+--------+--------+--------+',
 
 1922      .  
'--------+--------+--------+--------+--------+--------+')
 
 1931  5008       
format((i5,
' |',   7
x ,
' |' ,  f7.3,
' |' ,   7
x ,
' |',
 
 1932      .           f7.3,
' |' ,  f7.2,
' |', 2( 7
x ,
' |'),   7
x ,
' |',
 
 1933      .                                     f7.3,
' |' ,  f7.2,
' |'))
 
 1936  5009       
format(
'      |',9(8
x,
'|'),f7.3,
' |')
 
 2000       common/llocal_bsn/blowin
 
 2006       common/rlocal_bsn/facsbs,facubs,       
 
 2007      .                  por_bs,sheabs,rcd10n 
 
 2009       integer  ikl   ,isn   ,isnMAX,is2      
 
 2010       integer  Mobilm,Mobiln                 
 
 2011       integer  Mobile(
klonv)                 
 
 2016       real     SaltM1,SaltM2,SaltMo,SaltMx   
 
 2032       real     dzweqo,dzweqn,bsno_x          
 
 2034       real     PorSno,Salt_f,PorRef,ro_new   
 
 2039       integer  isagr1(
klonv)                 
 
 2040       integer  isagr2(
klonv)                 
 
 2068       data     agblow  / 1.00     /          
 
 2069       data     saltmx  /-5.83e-2  /          
 
 2070       data     facrbs  / 2.868    /          
 
 2071       data     factbs  / 0.085    /          
 
 2072       data     hdrift  / 1.00e+1  /          
 
 2073       data     h_mmwe  / 0.01e00  /          
 
 2074       data     tfv_vk  / 5.10e-1  /          
 
 2082       IF (.NOT.blowin)                                             
THEN 
 2084                facsbs =  1.             / facrbs
 
 2085                facubs =  1.             / factbs
 
 2087                sheabs =                   por_bs/(1.00-por_bs)
 
 2098                write(6,5000)  1./  rcd10n
 
 2099  5000          
format(/,
' Blowing Snow Model  Initialization     ',
 
 2100      .                /,
' Vt / u*t =',f8.2,
' (Neutral Assumption)',
 
 2101      .                /,
'           ', 8
x ,
' (Budd assumes  26.5)',/)
 
 2118       IF (.NOT.blomod)                                     
GO TO 1000
 
 2133         snowok   =  min(1   , max(
isnosv(ikl)      +1 -isn ,0))  
 
 2140         saltok   =  saltok   * snowok
 
 2141         saltm1   = -0.750e-2 * 
g1snsv(ikl,isn) 
 
 2142      .             -0.500e-2 * 
g2snsv(ikl,isn)+ 0.500e00
 
 2146         saltm2   = -0.833d-2 * 
g1snsv(ikl,isn) 
 
 2147      .             -0.583d-2 * 
g2snsv(ikl,isn)+ 0.833d00
 
 2148         saltmo   = (dendok   * saltm1 + (1.-dendok) *     saltm2       )
 
 2154         saltmo   =                                    max(saltmo,min_mo)
 
 2156         saltmo   =  saltok   * saltmo + (1.-saltok) * min(saltmo,saltmx)
 
 2160         saltsu   =     (1.00d0+saltmo)     *facsbs
 
 2172  6010     
format(/,
'SISVAT_BSn',6
x 
 2173      .           ,6
x,i3,2
x,
'G1         =',f6.3,
'   G2         =',f7.3
 
 2174      .           ,      
'   ro [kg/m3] =',f9.3,
'   Age* [Day] =',f9.3
 
 2175      .           ,   /,27
x,
'SaltM1     =',f6.3,
'   SaltM2     =',f7.3
 
 2176      .           ,      
'   Mobility I.=',f9.3,
'   Vt   [m/s] =',f9.3
 
 2177      .           ,   /,27
x,
'            ', 6
x ,
'               ', 7
x 
 2178      .           ,      
'               ', 9
x ,
'   Vn10 [m/s] =',f9.3)
 
 2183         shearx      =                por_bs/max(
eps6,
un_1-por_bs) 
 
 2192         argfac      =  max(
zer0  ,sheabs-shearx)     
 
 2194         fac_mo   =     exp(       argfac       )     
 
 2198         saltsu      =  max(
eps6  ,    saltsu)
 
 2199         saltsu      =  exp(fac_mo*log(saltsu))
 
 2200         argusi      =     -factbs              *
us__sv(ikl)/rcd10n
 
 2201         saltsi(ikl,isn) = (saltsu-exp(argusi)) *facrbs
 
 2207         snowok   =  1 -min(1,iabs(isn-
isnosv(ikl)))
 
 2208         salt_u   =               -log(saltsu)  *facubs
 
 2211         usthsv(ikl) =     snowok *   (salt_u   *rcd10n)
 
 2212      .              + (1.-snowok)*    
usthsv(ikl)
 
 2224  6011     
format(      27
x,
'Fac_Mo     =',f6.3,
'   Por_BS     =',f7.3
 
 2225      .           ,      
'   Drift    I.=',f9.3,
'   ut*_0[m/s] =',f9.3)
 
 2234         mobile(ikl) = 
nsno+1
 
 2238         isnmax      =       max(   1,          
isnosv(ikl)             )
 
 2239         isnmax      =       min( isn,          isnmax                  )
 
 2240         mobiln      = isn * max(
zer0,sign(
un_1,saltsi(ikl     ,isnmax)))
 
 2241         mobilm      =   1 - min(1   ,          mobile(ikl) -1 -mobiln)
 
 2244         mobile(ikl) =                 mobilm * mobiln 
 
 2245      .              +              (1-mobilm)* mobile(ikl)
 
 2256         dbsaux(ikl) = 
dbs_sv(ikl)
 
 2261         zdrift(ikl)     =  zdrift(ikl)
 
 2262      .            + 0.50 * 
dzsnsv(ikl,isn) * (3.25  -saltsi(ikl,isn))
 
 2263         sdrift(ikl,isn) =  saltsi(ikl,isn)
 
 2264      .          *exp(  max(
ea_min, -zdrift(ikl)     *hdrift     ))
 
 2265      .          *min(1,max(0     ,  isn +1          -mobile(ikl)))
 
 2266      .          *min(1,max(0     ,  
isnosv(ikl)     -isn +1     ))
 
 2270         xdrift(ikl)     =           sdrift(ikl,isn) +xdrift(ikl)
 
 2271         zdrift(ikl)     =  zdrift(ikl)
 
 2272      .            + 0.50 * 
dzsnsv(ikl,isn) * (3.25  -saltsi(ikl,isn))
 
 2280         sdrift(ikl,isn) =  sdrift(ikl,isn) /max(
eps6,xdrift(ikl))
 
 2295         zdepos(ikl,isn) =      exp(-zdrift(ikl)   )
 
 2296      .          *min(1,max(0     ,  isn +1          -mobile(ikl)))
 
 2297      .          *min(1,max(0     ,  
isnosv(ikl    ) -isn +1     ))
 
 2299         tdepos(ikl) = tdepos(ikl) + zdepos(ikl,isn)
 
 2300         zdrift(ikl) = zdrift(ikl) + 
dzsnsv(ikl,isn) *
ro__sv(ikl,isn)
 
 2309         zdepos(ikl,isn) = zdepos(ikl,isn) / max(
eps6,tdepos(ikl))
 
 2325      .      18
x,
'MB0',6
x,
'Sno1WE [mm]=',f9.3,19
x,
'0  dbs_SV [mm]=',f9.6)
 
 2336         snowok      = min(1,max(
isnosv(ikl)+1-isn ,0))      
 
 2338         bsno_x      = dbsaux(ikl)     *sdrift(ikl,isn)
 
 2339         dzweqn      = dzweqo          +bsno_x
 
 2340         dzweqn  = max(dzweqn,          h_mmwe *snowok)
 
 2359      .      18
x,
'MB1',6
x,
'Sno1WE [mm]=',f9.3,19
x,
'1  dbs_SV [mm]=',f9.6)
 
 2363      .      18
x,
'MB ',5
x,
'(After  [mm]=',f6.0, 
')-(Erosion[mm]=', f7.3,
 
 2364      .                                         
')-(Before [mm]=', f9.3,
 
 2365      .                                         
')= Budget [mm]=', f9.6)
 
 2382  6003         
format(/,41
x,
'tdepos [-] =',f6.3,40
x,
'Mobil',i3
 
 2383      .              ,/,27
x,
'Salt.Index    sdrift' 
 2384      .              ,      
'    zdepos  ro__snow  ro_bsnow  roN_snow' 
 2385      .              ,                
'  dz__snow  dz_bsnow  dzN_snow' 
 2387      .              ,/,27
x,
'             [kg/m3]   [kg/m3]   [kg/m3]' 
 2394           weagre(ikl) = weagre(ikl) + 
ro__sv(ikl,isn)*
dzsnsv(ikl,isn)
 
 2395           isagr1(ikl) = 
istosv(ikl,isn)
 
 2414           roagr1(ikl) = 
ro__sv(ikl,isn)
 
 2415           roagr2(ikl) = ro_new
 
 2416           hsno_x      = tdepos(ikl)*  zdepos(ikl,isn)
 
 2418           dzagr1(ikl) = 
dzsnsv(ikl,isn)
 
 2419           dzagr2(ikl) = hsno_x     /  ro_new
 
 2426           t_agr1(ikl) =    
tsissv(ikl,isn)
 
 2428           etagr1(ikl) =    
eta_sv(ikl,isn)
 
 2430           g1agr1(ikl) =    
g1snsv(ikl,isn)
 
 2432           g2agr1(ikl) =    
g2snsv(ikl,isn)
 
 2437           agagr1(ikl) =    
agsnsv(ikl,isn)
 
 2446      .                 (isagr1,isagr2,weagre
 
 2447      .                 ,dzagr1,dzagr2,t_agr1,t_agr2
 
 2448      .                 ,roagr1,roagr2,etagr1,etagr2
 
 2449      .                 ,g1agr1,g1agr2,g2agr1,g2agr2
 
 2450      .                 ,agagr1,agagr2,agrege
 
 2465  6004         
format((27
x,i3,f7.2,2f10.6,3f10.3,4f10.6))
 
 2467           istosv(ikl,isn) = isagr1(ikl)
 
 2468           dzsnsv(ikl,isn) = dzagr1(ikl)
 
 2469           tsissv(ikl,isn) = t_agr1(ikl)
 
 2470           ro__sv(ikl,isn) = roagr1(ikl)
 
 2471           eta_sv(ikl,isn) = etagr1(ikl)
 
 2472           g1snsv(ikl,isn) = g1agr1(ikl)
 
 2473           g2snsv(ikl,isn) = g2agr1(ikl)
 
 2474           agsnsv(ikl,isn) = agagr1(ikl)
 
 2485  6008   
format(i3,
'  dsn+Buf=',f6.2,6
x,
'A dz *ro =',10f6.2,
 
 2489           hdrift      =  tdepos(ikl)/
dt__sv 
 2542       common /sisvat_bdu_r/etaust
 
 2544       common /sisvat_bdu_l/logust
 
 2567       data          (f__ust(isot),                       isot=0,
nvgt)
 
 2582       IF (.NOT.logust)                                              
THEN 
 2584                etaust(isot) = 0.0014 * claypc(isot) * claypc(isot) 
 
 2585      .                      + 0.17   * claypc(isot)                
 
 2597         usthdu      = sqrt(
un_1+1.21*exp(0.68*   log(eta_du)    )) 
 
 2602      .                 usthdu     *   max(0,1-
isnosv(ikl))
 
 2670       common/sisvat_sic_r/crodzw,lro__i
 
 2672       common/sisvat_sic_l/sicini
 
 2691       IF (.NOT.sicini)                                              
THEN 
 2695      .                 -(salice/salwat)*(1.-1.e-3*salwat) )   
 
 2707         ocn_ok        =  (1   -
lsmask(ikl)   )                
 
 2708      .             *max(0,1   -
isnosv(ikl)   )                
 
 2837       integer   NLay_s(
klonv)                 
 
 2838       integer   isagr1(
klonv)                 
 
 2839       integer   isagr2(
klonv)                 
 
 2846       integer   isn1  (
klonv)                 
 
 2896       data      dzepsi /   0.0045/            
 
 2897       data      dzxmin /   0.0060/            
 
 2898       data      dz_min /   0.0150/            
 
 2899       data      dz_max /   0.0900/            
 
 2930           isno_n    =             
isnosv(ikl)-isn+1     
 
 2932           iiceok    = min(1,max(0,iice_n         +1))   
 
 2942      .                 *((1-iiceok)*isno_n*isno_n       
 
 2943      .                 +    iiceok *2.   **iice_n)      
 
 2951      .                           dz_dif-dzthin(ikl)))   
 
 2960      .                               *
g1snsv(ikl,max(1,isn-1)))) 
 
 2973           dzthin(ikl) =    (1. - okthin)  * dzthin(ikl) 
 
 2981  4150     
format(/,
'-',a18,i5,
' ',70(
'-'),
 
 2982      .           /,
' Thinest ',i3,
':',f9.3)
 
 3005      .                 *(1 -min(abs(
isnosv(ikl)        
 
 3014      .                               *
g1snsv(ikl,max(1,isn-1)))) 
 
 3028  4151     
format(
' Thinest ',i3,
':',f9.3,
'   Max   =',i3,f12.3)
 
 3032  470    
format(
'Before _zCr1: G1 = ',10f8.1,(/,19
x,10f8.1))
 
 3034  472    
format(
'              G2 = ',10f8.1,(/,19
x,10f8.1))
 
 3050           isagr1(ikl) =    
istosv(ikl,isn)
 
 3052           dzagr1(ikl) =    
dzsnsv(ikl,isn)
 
 3054           t_agr1(ikl) =    
tsissv(ikl,isn)
 
 3056           roagr1(ikl) =    
ro__sv(ikl,isn)
 
 3058           etagr1(ikl) =    
eta_sv(ikl,isn)
 
 3060           g1agr1(ikl) =    
g1snsv(ikl,isn)
 
 3062           g2agr1(ikl) =    
g2snsv(ikl,isn)
 
 3064           agagr1(ikl) =    
agsnsv(ikl,isn)
 
 3066           lstlay      = min(1,max(  0,
isnosv(ikl) -1))  
 
 3068      .     -(1-lstlay)* max(
zer0,                       
 
 3072           agrege(ikl) = max(
zer0,                       
 
 3076      .                * min( max(0   ,
isnosv(ikl)+1     
 
 3085           weagre(ikl) = weagre(ikl) + 
ro__sv(ikl,isn)*
dzsnsv(ikl,isn)
 
 3086      .                                *min(1,max(0,
i_thin(ikl)+1-isn))
 
 3092  410  
format(/,
' Agregation of too THIN Layers')
 
 3097  411  
format(
' dz_ref [cm]:',10f8.2   ,/,(
'             ',10f8.2) )
 
 3098  412  
format(
' dz_dif [cm]:',10f8.2   ,/,(
'             ',10f8.2) )
 
 3099  413  
format(
' dzsnSV [cm]:',10f8.2   ,/,(
'             ',10f8.2) )
 
 3100  414  
format(
'             ',10(i5,3
x),/,(
'             ',10(i5,3
x)))
 
 3107  4111 
format(
' isnoSV     :',  i8  )
 
 3108  4112 
format(
' i_thin     :',  i8  )
 
 3109  4113 
format(
' LIndsv     :',  i8  )
 
 3110  4114 
format(
' Agrege     :',  f8.2)
 
 3111  4115 
format(
' dzagr1     :',  f8.2)
 
 3112  4116 
format(
' dzagr2     :',  f8.2)
 
 3116  471    
format(
'Before _zAg1: G1 = ',10f8.1,(/,19
x,10f8.1))
 
 3125      .                 (isagr1,isagr2,weagre
 
 3126      .                 ,dzagr1,dzagr2,t_agr1,t_agr2
 
 3127      .                 ,roagr1,roagr2,etagr1,etagr2
 
 3128      .                 ,g1agr1,g1agr2,g2agr1,g2agr2
 
 3129      .                 ,agagr1,agagr2,agrege
 
 3141           isn     = min(isn,isn+
lindsv(ikl))
 
 3144      .            -max(0,sign(1,
iicesv(ikl) -isn +icemix))
 
 3146      .            *max(0,sign(1,
iicesv(ikl) -1          ))
 
 3148      .                      +   agrege(ikl) *isagr1(ikl)
 
 3150      .                      +   agrege(ikl) *dzagr1(ikl)
 
 3152      .                      +   agrege(ikl) *t_agr1(ikl)
 
 3154      .                      +   agrege(ikl) *roagr1(ikl)
 
 3156      .                      +   agrege(ikl) *etagr1(ikl)
 
 3158      .                      +   agrege(ikl) *g1agr1(ikl)
 
 3160      .                      +   agrege(ikl) *g2agr1(ikl)
 
 3162      .                      +   agrege(ikl) *agagr1(ikl)
 
 3172             staggr        =  min(1,max(0,i +1 -isn1(ikl)   ))
 
 3174      .            + staggr*((1.-agrege(ikl))*
istosv(ikl,i  )
 
 3175      .                      +   agrege(ikl) *
istosv(ikl,i+1))
 
 3177      .            + staggr*((1.-agrege(ikl))*
dzsnsv(ikl,i  )
 
 3178      .                      +   agrege(ikl) *
dzsnsv(ikl,i+1))
 
 3180      .            + staggr*((1.-agrege(ikl))*
tsissv(ikl,i  )
 
 3181      .                      +   agrege(ikl) *
tsissv(ikl,i+1))
 
 3183      .            + staggr*((1.-agrege(ikl))*
ro__sv(ikl,i  )
 
 3184      .                      +   agrege(ikl) *
ro__sv(ikl,i+1))
 
 3186      .            + staggr*((1.-agrege(ikl))*
eta_sv(ikl,i  )
 
 3187      .                      +   agrege(ikl) *
eta_sv(ikl,i+1))
 
 3189      .            + staggr*((1.-agrege(ikl))*
g1snsv(ikl,i  )
 
 3190      .                      +   agrege(ikl) *
g1snsv(ikl,i+1))
 
 3192      .            + staggr*((1.-agrege(ikl))*
g2snsv(ikl,i  )
 
 3193      .                      +   agrege(ikl) *
g2snsv(ikl,i+1))
 
 3195      .            + staggr*((1.-agrege(ikl))*
agsnsv(ikl,i  )
 
 3196      .                      +   agrege(ikl) *
agsnsv(ikl,i+1))
 
 3217  5991     
format(/,
'First Agregation / Layer',i3,
 
 3218      .           /,
'  i',11
x,
'T',9
x,
'rho',10
x,
'dz',11
x,
'H')
 
 3222  5995     
format(i3,3f12.3,i12)
 
 3238           isno_n    =             
isnosv(ikl)-isn+1     
 
 3240           iiceok    = min(1,max(0,iice_n         +1))   
 
 3241           dz_dif    =(      
dzsnsv(ikl,isn)             
 
 3242      .        - dz_max *((1-iiceok)*isno_n*isno_n       
 
 3243      .                 +    iiceok *2.   **iice_n)  )   
 
 3247      .                           dz_dif-dzthin(ikl)))   
 
 3253           dzthin(ikl) =    (1. - okthin)  * dzthin(ikl) 
 
 3260      .                      sign(
un_1,dzthin(ikl)       
 
 3262      .        * max(0,1-max(0   ,     
isnosv(ikl)       
 
 3264           agrege(ikl) =               thickl            
 
 3265      .        * max(0,1-max(0   ,     
nlaysv(ikl)       
 
 3268           nlay_s(ikl) =               thickl            
 
 3269      .        * max(0,1-max(0   ,     
nlaysv(ikl)       
 
 3273           nlay_s(ikl) = max(0   ,     nlay_s(ikl))      
 
 3278  4152     
format(/,
' Thickest',i3,
':',f9.3,
'   Split =',f4.0)
 
 3286         IF (agrege(ikl).gt.0..AND.
i_thin(ikl).lt.
isnosv(ikl))       
THEN 
 3287           staggr          =  min(1,max(0,isn-
i_thin(ikl)    -1))
 
 3288      .                    *  min(1,max(0,    
isnosv(ikl)-isn+2))
 
 3290      .                    + (1. -  staggr) * 
istosv(ikl ,isn  )
 
 3292      .                    + (1. -  staggr) * 
dzsnsv(ikl ,isn  )
 
 3294      .                    + (1. -  staggr) * 
tsissv(ikl ,isn  )
 
 3296      .                    + (1. -  staggr) * 
ro__sv(ikl ,isn  )
 
 3298      .                    + (1. -  staggr) * 
eta_sv(ikl ,isn  )
 
 3300      .                    + (1. -  staggr) * 
g1snsv(ikl ,isn  )
 
 3302      .                    + (1. -  staggr) * 
g2snsv(ikl ,isn  )
 
 3304      .                    + (1. -  staggr) * 
agsnsv(ikl ,isn  )
 
 3312      .                    + (1.-agrege(ikl))*
dzsnsv(ikl,isn)
 
 3316      .                    + (1.-agrege(ikl))*
istosv(ikl,isn)
 
 3318      .                    + (1.-agrege(ikl))*
dzsnsv(ikl,isn)
 
 3320      .                    + (1.-agrege(ikl))*
tsissv(ikl,isn)
 
 3322      .                    + (1.-agrege(ikl))*
ro__sv(ikl,isn)
 
 3324      .                    + (1.-agrege(ikl))*
eta_sv(ikl,isn)
 
 3326      .                    + (1.-agrege(ikl))*
g1snsv(ikl,isn)
 
 3328      .                    + (1.-agrege(ikl))*
g2snsv(ikl,isn)
 
 3330      .                    + (1.-agrege(ikl))*
agsnsv(ikl,isn)
 
 3333      .                    +     agrege(ikl) *max(0,sign(1,
iicesv(ikl) 
 
 3335      .                                      *max(0,sign(1,
iicesv(ikl) 
 
 3351      .          
'dzsnSV      dz_min      dz_dif      ',
 
 3352      .          
'OKthin      dzthin   i_thin')
 
 3359           isno_n    =             
isnosv(ikl)-isn+1     
 
 3361           iiceok    = min(1,max(0,iice_n         +1))   
 
 3371      .        /max(
eps6,((1-iiceok)*isno_n*isno_n       
 
 3372      .                 +    iiceok *2.   **iice_n))     
 
 3379      .                           dz_dif  - dzthin(ikl)))
 
 3385           dzthin(ikl) =    (1. - okthin) * dzthin(ikl)  
 
 3392  6001     
format(i3,5f12.6,i9)
 
 3402  4153     
format(/,
' Thinest ',i3,
':',f9.3)
 
 3408  473    
format(
'Before _zCr2: G1 = ',10f8.1,(/,19
x,10f8.1))
 
 3425           isagr1(ikl) =    
istosv(ikl,isn)
 
 3427           dzagr1(ikl) =    
dzsnsv(ikl,isn)
 
 3429           t_agr1(ikl) =    
tsissv(ikl,isn)
 
 3431           roagr1(ikl) =    
ro__sv(ikl,isn)
 
 3433           etagr1(ikl) =    
eta_sv(ikl,isn)
 
 3435           g1agr1(ikl) =    
g1snsv(ikl,isn)
 
 3437           g2agr1(ikl) =    
g2snsv(ikl,isn)
 
 3439           agagr1(ikl) =    
agsnsv(ikl,isn)
 
 3441           lstlay      = min(1,max(  0,    
isnosv(ikl)-1   ))
 
 3442           agrege(ikl) = min(1,
 
 3448      .     -(1-lstlay)*max(
zer0,
 
 3458           weagre(ikl) = weagre(ikl) + 
ro__sv(ikl,isn)*
dzsnsv(ikl,isn)
 
 3459      .                                *min(1,max(0,
i_thin(ikl)+1-isn))
 
 3465  4120 
format(
' Agregation of too MUCH Layers')
 
 3477  474    
format(
'Before _zAg2: G1 = ',10f8.1,(/,19
x,10f8.1))
 
 3486      .                 (isagr1,isagr2,weagre
 
 3487      .                 ,dzagr1,dzagr2,t_agr1,t_agr2
 
 3488      .                 ,roagr1,roagr2,etagr1,etagr2
 
 3489      .                 ,g1agr1,g1agr2,g2agr1,g2agr2
 
 3490      .                 ,agagr1,agagr2,agrege
 
 3502           isn     = min(isn,isn+
lindsv(ikl))
 
 3505      .            -max(0,sign(1,
iicesv(ikl) -isn +icemix))
 
 3507      .            *max(0,sign(1,
iicesv(ikl) -1          ))
 
 3509      .                      +   agrege(ikl) *isagr1(ikl)
 
 3511      .                      +   agrege(ikl) *dzagr1(ikl)
 
 3513      .                      +   agrege(ikl) *t_agr1(ikl)
 
 3515      .                      +   agrege(ikl) *roagr1(ikl)
 
 3517      .                      +   agrege(ikl) *etagr1(ikl)
 
 3519      .                      +   agrege(ikl) *g1agr1(ikl)
 
 3521      .                      +   agrege(ikl) *g2agr1(ikl)
 
 3523      .                      +   agrege(ikl) *agagr1(ikl)
 
 3533             staggr        =  min(1,max(0,i +1 -isn1(ikl)   ))
 
 3535      .            + staggr*((1.-agrege(ikl))*
istosv(ikl,i  )
 
 3536      .                      +   agrege(ikl) *
istosv(ikl,i+1))
 
 3538      .            + staggr*((1.-agrege(ikl))*
dzsnsv(ikl,i  )
 
 3539      .                      +   agrege(ikl) *
dzsnsv(ikl,i+1))
 
 3541      .            + staggr*((1.-agrege(ikl))*
tsissv(ikl,i  )
 
 3542      .                      +   agrege(ikl) *
tsissv(ikl,i+1))
 
 3544      .            + staggr*((1.-agrege(ikl))*
ro__sv(ikl,i  )
 
 3545      .                      +   agrege(ikl) *
ro__sv(ikl,i+1))
 
 3547      .            + staggr*((1.-agrege(ikl))*
eta_sv(ikl,i  )
 
 3548      .                      +   agrege(ikl) *
eta_sv(ikl,i+1))
 
 3550      .            + staggr*((1.-agrege(ikl))*
g1snsv(ikl,i  )
 
 3551      .                      +   agrege(ikl) *
g1snsv(ikl,i+1))
 
 3553      .            + staggr*((1.-agrege(ikl))*
g2snsv(ikl,i  )
 
 3554      .                      +   agrege(ikl) *
g2snsv(ikl,i+1))
 
 3556      .            + staggr*((1.-agrege(ikl))*
agsnsv(ikl,i  )
 
 3557      .                      +   agrege(ikl) *
agsnsv(ikl,i+1))
 
 3575  475    
format(
'At End _zSn : G1 = ',10f8.1,(/,19
x,10f8.1))
 
 3655       integer   ikl   ,isn   ,is0   ,is1                                        
 
 3657       real*8    Dtyp_0,Dtyp_1                 
 
 3671       data      dtypmx / 200.0  /             
 
 3674       data      dtypdi /  10.0  /             
 
 3675       data      dtyphi / 100.0  /             
 
 3683           isn         = max(1          ,
i_thin(ikl))                            
 
 3689           is0    = max(1,        
i_thin(ikl)-1 )        
 
 3701      .    *    dendok *((abs(
g1snsv(ikl,isn)            
 
 3704      .                      -
g2snsv(ikl,is0))) *dtypsp  
 
 3706      .                      -
ro__sv(ikl,is0))  *dtypro) 
 
 3708      .    *(1.-dendok)*((abs(
g1snsv(ikl,isn)            
 
 3711      .                      -
g2snsv(ikl,is0))) *dtypdi  
 
 3713      .                      -
ro__sv(ikl,is0))  *dtypro) 
 
 3718      .                      -
istosv(ikl,is0))  *dtyphi) 
 
 3719      .        +             (1 -abs(isn-is0))  * 1.e+6  
 
 3720      .        +  max(0,1-abs(
iicesv(ikl)                
 
 3727           is1    = min(          
i_thin(ikl)+1,         
 
 3740      .    *    dendok *((abs(
g1snsv(ikl,isn)            
 
 3743      .                      -
g2snsv(ikl,is1))) *dtypsp  
 
 3745      .                      -
ro__sv(ikl,is1))  *dtypro) 
 
 3747      .    *(1.-dendok)*((abs(
g1snsv(ikl,isn)            
 
 3750      .                      -
g2snsv(ikl,is1))) *dtypdi  
 
 3752      .                      -
ro__sv(ikl,is1))  *dtypro) 
 
 3757      .                      -
istosv(ikl,is1))  *dtyphi) 
 
 3758      .        +             (1 -abs(isn-is1))  * 1.e+6  
 
 3759      .        +  max(0,1-abs(
iicesv(ikl)                
 
 3768           isno_1      = (1 -min(abs(
isnosv(ikl)        
 
 3770      .                * (1 -min(abs(
isnosv(ikl)        
 
 3782      .                     (isagra,isagrb,weagra                                
 
 3783      .                     ,dzagra,dzagrb,t_agra,t_agrb                         
 
 3784      .                     ,roagra,roagrb,etagra,etagrb                         
 
 3785      .                     ,g1agra,g1agrb,g2agra,g2agrb                         
 
 3786      .                     ,agagra,agagrb,agreg1                                
 
 3847       integer   isagrb(
klonv)                 
 
 3848       real*8    dzagrb(
klonv)                 
 
 3849       real*8    T_agrb(
klonv)                 
 
 3850       real*8    roagrb(
klonv)                 
 
 3851       real*8    etagrb(
klonv)                 
 
 3852       real*8    G1agrb(
klonv)                 
 
 3853       real*8    G2agrb(
klonv)                 
 
 3854       real*8    agagrb(
klonv)                 
 
 3860       integer   isagra(
klonv)                 
 
 3861       real*8    WEagra(
klonv)                 
 
 3862       real*8    Agreg1(
klonv)                 
 
 3863       real*8    dzagra(
klonv)                 
 
 3864       real*8    T_agra(
klonv)                 
 
 3865       real*8    roagra(
klonv)                 
 
 3866       real*8    etagra(
klonv)                 
 
 3867       real*8    G1agra(
klonv)                 
 
 3868       real*8    G2agra(
klonv)                 
 
 3869       real*8    agagra(
klonv)                 
 
 3924           dz      =  dzagra(ikl) + dzagrb(ikl)                                  
 
 3925           dzro_1  =  roagra(ikl) * dzagra(ikl)                                  
 
 3926           dzro_2  =  roagrb(ikl) * dzagrb(ikl)                                  
 
 3927           dzro    =  dzro_1      + dzro_2                                       
 
 3930           wn      = (dzro_1*etagra(ikl) + dzro_2*etagrb(ikl))                   
 
 3932           tn      = (dzro_1*t_agra(ikl) + dzro_2*t_agrb(ikl))                   
 
 3934           ag      = (dzro_1*agagra(ikl) + dzro_2*agagrb(ikl))                   
 
 3939           nh      =  nh__ok *       max(isagra(ikl),isagrb(ikl))                
 
 3940      .          + (1-nh__ok)*       min(isagra(ikl),isagrb(ikl))                
 
 3948  5995     
format(
' WE2,WEa =',2f9.1,
'  nha,b =',2i2,
'  nh__OK,nh =',2i2)        
 
 3958      .                   sign(
unun, g1agra(ikl) *g1agrb(ikl)  - 
eps_21))        
 
 3959           g1same  = (dzro_1*g1agra(ikl) + dzro_2*g1agrb(ikl))                   
 
 3961           g2same  = (dzro_1*g2agra(ikl) + dzro_2*g2agrb(ikl))                   
 
 3967           zronew  =     typ__1  *dzro_1                    
 
 3968      .            + (1.-typ__1) *dzro_2                    
 
 3969           g1_new  =     typ__1  *g1agra(ikl)               
 
 3970      .            + (1.-typ__1) *g1agrb(ikl)               
 
 3971           g2_new  =     typ__1  *g2agra(ikl)               
 
 3972      .            + (1.-typ__1) *g2agrb(ikl)               
 
 3973           zroold  = (1.-typ__1) *dzro_1                    
 
 3975           g1_old  = (1.-typ__1) *g1agra(ikl)               
 
 3976      .            +     typ__1  *g1agrb(ikl)               
 
 3977           g2_old  = (1.-typ__1) *g2agra(ikl)               
 
 3978      .            +     typ__1  *g2agrb(ikl)               
 
 3986           siz_av  = (zronew*siznew+zroold*sizold)      
 
 3988           sph_av  = (zronew*sphnew+zroold*sphold)      
 
 3990           den_av  = (siz_av -(    sph_av *
dscdsv        
 4005           g1diff  =(   -dendok *den_av                                          
 
 4006      .             +(1.-dendok)*sph_av) *
g1_dsv                                  
 4007           g2diff  =     dendok *sph_av  *
g1_dsv                                  
 4008      .             +(1.-dendok)*siz_av                                          
 
 4010      .             +(1.-sameok)*g1diff                                          
 
 4012      .             +(1.-sameok)*g2diff                                          
 
 4018           isagra(ikl)   = agreg1(ikl) *nh +(1.-agreg1(ikl)) *isagra(ikl)        
 
 4019           dzagra(ikl)   = agreg1(ikl) *dz +(1.-agreg1(ikl)) *dzagra(ikl)        
 
 4020           t_agra(ikl)   = agreg1(ikl) *tn +(1.-agreg1(ikl)) *t_agra(ikl)        
 
 4021           roagra(ikl)   = agreg1(ikl) *ro +(1.-agreg1(ikl)) *roagra(ikl)        
 
 4022           etagra(ikl)   = agreg1(ikl) *wn +(1.-agreg1(ikl)) *etagra(ikl)        
 
 4023           g1agra(ikl)   = agreg1(ikl) *g1 +(1.-agreg1(ikl)) *g1agra(ikl)        
 
 4024           g2agra(ikl)   = agreg1(ikl) *g2 +(1.-agreg1(ikl)) *g2agra(ikl)        
 
 4025           agagra(ikl)   = agreg1(ikl) *ag +(1.-agreg1(ikl)) *agagra(ikl)        
 
 4034       subroutine snoptp(jjtime)
 
 4144       integer   isn   ,ikl   ,isn1  ,jjtime
 
 4145       real      sbeta1,sbeta2,sbeta3,sbeta4,sbeta5
 
 4146       real      AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK
 
 4147       real      dalbed,dalbeS,dalbeW
 
 4149       real      RoFrez,DiffRo,SignRo,SnowOK,OpSqrt
 
 4150       real      albSn1,albIc1,a_SnI1,a_SII1
 
 4151       real      albSn2,albIc2,a_SnI2,a_SII2
 
 4152       real      albSn3,albIc3,a_SnI3,a_SII3
 
 4153       real      albSno,albIce,albSnI,albSII,albWIc,albmax
 
 4154       real      doptic,Snow_H,SIce_H,SnownH,SIcenH
 
 4155       real      exarg1,exarg2,exarg3,sign_0,sExt_0
 
 4169       data      sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/
 
 4187       data      doptmx  /2.3e-3/
 
 4205           sph_ok          =  max(
zer0,signg1)
 
 4207           snopsv(ikl,isn) =   1.e-4 *
 
 4227           snopsv(ikl,isn) =  max(
zer0,snopsv(ikl,isn))
 
 4259           snowok =  max(
zer0,signro) 
 
 4261           opsqrt = sqrt(snopsv(ikl,isn))
 
 4263           albsn1 =  0.96-1.580*opsqrt
 
 4264           albsn1 =  max(albsn1,albmin)
 
 4266           albsn1 =  max(albsn1,
zer0)
 
 4267           albsn1 =  min(albsn1,
un_1)
 
 4269           albsn2 =  0.95-15.40*opsqrt
 
 4270           albsn2 =  max(albsn2,
zer0)
 
 4271           albsn2 =  min(albsn2,
un_1)
 
 4273           doptic =  min(snopsv(ikl,isn),doptmx)
 
 4274           albsn3 =  346.3*doptic -32.31*opsqrt +0.88
 
 4275           albsn3 =  max(albsn3,
zer0)
 
 4276           albsn3 =  min(albsn3,
un_1)
 
 4285           albsn1 =  snowok*albsn1+(1.0-snowok)*max(albsno,
ai3dsv)
 
 4286           albsn2 =  snowok*albsn2+(1.0-snowok)*max(albsno,
ai3dsv)
 
 4287           albsn3 =  snowok*albsn3+(1.0-snowok)*max(albsno,
ai3dsv)
 
 4296           snownh =  snow_h  /  hsnosv
 
 4297           snownh =  min(
un_1,  snownh)
 
 4298           sicenh =  sice_h  / (hicesv
 
 4301           sicenh =  min(
un_1,  sicenh)
 
 4309        do isn    =  
isnosv(ikl),1,-1
 
 4310           ro_ave =  ro_ave + 
ro__sv(ikl,isn) * 
dzsnsv(ikl,isn) * snowok
 
 4311           dz_ave =  dz_ave +                   
dzsnsv(ikl,isn) * snowok
 
 4312           snowok =  max(
zer0,sign(
un_1,1.-dz_ave))
 
 4315           ro_ave =  ro_ave / max(dz_ave,
eps6)
 
 4316           snowok =  max(
zer0,sign(
un_1,700.-ro_ave))
 
 4318           snownh =  snowok + snownh * (1. - snowok)   
 
 4329      .           *  (1  -min(1,iabs(isn-
isnosv(ikl)))))    
 
 4333           snowok =  max(
zer0,signro)
 
 4335           albwic = (1. - snowok) * albwic + snowok   
 
 4346           a_sii1      =     albwic      +(albsn1-albwic)     *snownh
 
 4347           a_sii1      = min(a_sii1       ,albsn1)
 
 4349           a_sii2      =     albwic      +(albsn2-albwic)     *snownh
 
 4350           a_sii2      = min(a_sii2       ,albsn2)
 
 4352           a_sii3      =     albwic      +(albsn3-albwic)     *snownh
 
 4353           a_sii3      = min(a_sii3       ,albsn3)
 
 4437      .                  + 
albssv(ikl) *(1.0        - sicenh))
 
 4477           coalbm      =  coalb1(ikl) +coalb2(ikl) +coalb3(ikl)
 
 4478           coalb1(ikl) =  coalb1(ikl)              /coalbm
 
 4479           coalb2(ikl) =  coalb2(ikl)              /coalbm
 
 4480           coalb3(ikl) =  coalb3(ikl)              /coalbm
 
 4487           snowok =  max(
zer0,signro) 
 
 4489           rofrez =  1.e-3      * 
ro__sv(ikl,isn) * (1.0-
eta_sv(ikl,isn))
 
 4491           opsqrt = sqrt(max(
eps6,snopsv(ikl,isn)))
 
 4492           exarg1 =      snowok  *1.e2 *max(sbeta1*rofrez/opsqrt,sbeta2)
 
 4493      .            +(1.0-snowok)           *sbeta5
 
 4494           exarg2 =      snowok  *1.e2 *max(sbeta3*rofrez/opsqrt,sbeta4)
 
 4495      .            +(1.0-snowok)           *sbeta5
 
 4496           exarg3 =      snowok  *1.e2     *sbeta5
 
 4497      .            +(1.0-snowok)           *sbeta5
 
 4520           sext_1(ikl) = sext_1(ikl)
 
 4521      .                          * exp(min(0.0,-exarg1 *
dzsnsv(ikl,isn)))
 
 4522           sign_0      =              sign(
un_1,
epsn   -sext_1(ikl))
 
 4523           sext_0      =               max(
zer0,sign_0)*sext_1(ikl)
 
 4524           sext_1(ikl) = sext_1(ikl)                   -sext_0
 
 4526           sext_2(ikl) = sext_2(ikl)
 
 4527      .                          * exp(min(0.0,-exarg2 *
dzsnsv(ikl,isn)))
 
 4528           sign_0      =              sign(
un_1,
epsn   -sext_2(ikl))
 
 4529           sext_0      =               max(
zer0,sign_0)*sext_2(ikl)
 
 4530           sext_2(ikl) = sext_2(ikl)                   -sext_0
 
 4532           sext_3(ikl) = sext_3(ikl)
 
 4533      .                          * exp(min(0.0,-exarg3 *
dzsnsv(ikl,isn)))
 
 4534           sign_0      =              sign(
un_1,
epsn   -sext_3(ikl))
 
 4535           sext_0      =               max(
zer0,sign_0)*sext_3(ikl)
 
 4536           sext_3(ikl) = sext_3(ikl)                   -sext_0
 
 4538           sex_sv(ikl,isn) = coalb1(ikl) * sext_1(ikl)
 
 4539      .                    + coalb2(ikl) * sext_2(ikl)
 
 4540      .                    + coalb3(ikl) * sext_3(ikl)
 
 4562  460      
format(
'---------------------------------+----+',
 
 4563      .          
'-------+-------+-------+-------+-------+-------+',
 
 4564      .                                  
'-------+-------+-------+',
 
 4565      .         /,
'Snow/Ice Pack ',a18,
' |    |',
 
 4566      .          
' z [m] |0.3/0.8|0.8/1.5|1.5/2.8| Full  |Opt[mm]|',
 
 4567      .                                  
' G1    | G2    | ro    |',
 
 4568      .         /,
'---------------------------------+----+',
 
 4569      .          
'-------+-------+-------+-------+-------+-------+',
 
 4570      .                                  
'-------+-------+-------+')
 
 4575  461      
format(
'Integrated Snow/Ice/Soil  Albedo |',
 
 4576      .            3
x,
' |',  f6.3,
' |' ,4(f6.3,
' |'), 6
x ,
' |',
 
 4580  462      
format(
'Integrated Snow/Ice       Albedo |',
 
 4581      .            i3,
' |',   6
x ,
' |' ,4(f6.3,
' |'), 6
x ,
' |',
 
 4586  463      
format(
'Integrated Water/Bare Ice Albedo |',
 
 4587      .            3
x,
' |',  f6.3,
'w|' ,3( 6
x, 
' |'),  
 
 4588      .                                   f6.3,
' |' ,f6.3,
' |',
 
 4598  465      
format(
'Surficial       Ice Lense        |',
 
 4599      .            i3,
' |', (f6.3,
'i|'),4(f6.3,
' |'),f6.3,
' |',
 
 4609  466      
format(
'Uppermost  Effective Snow Layer  |',
 
 4610      .            i3,
' |', (f6.3,
'*|'),4(f6.3,
' |'),f6.3,
' |',
 
 4685       real      exdRad,k_drad,k___sv(
klonv)
 
 4687       real      zv_fac,zv1fac,deadLF
 
 4688       real      T_Rad0,A_Rad0,A0__sv(
klonv)
 
 4689       real                    r0_Rad,t0_Rad,nu_Rad
 
 4690       real      Tr_Rad,Re_Rad,r__Rad,t__Rad,t1_Rad
 
 4691       real      arggam, gamma,gamasv(
klonv),gammaL
 
 4692       real      denSig,Sig__c,Sigcsv(
klonv)
 
 4693       real      DDifH1,DDifC1,C1__sv(
klonv)
 
 4694       real      DDifH2,DDifC2,C2__sv(
klonv)
 
 4695       real      denS_s,denS_a,den_c1,DDif_L
 
 4696       real      u0_Vis,absg_V,absv_V
 
 4697       real      u0_nIR,absgnI,absvnI
 
 4698       real      argexg,argexk,criLAI(
klonv)
 
 4699       real      residu,d_DDif,dDDifs,dDDifa
 
 4727       DATA (revisl(ivg),renirl(ivg),trvisl(ivg),trnirl(ivg),
 
 4728      .      revisd(ivg),renird(ivg),trvisd(ivg),trnird(ivg),ivg=0,nvgt)
 
 4732      ./0.11,  0.58,  0.07,  0.25,  0.36,  0.58,  0.22,  0.38,
 
 4733      . 0.11,  0.58,  0.07,  0.25,  0.36,  0.58,  0.22,  0.38,
 
 4734      . 0.11,  0.58,  0.07,  0.25,  0.36,  0.58,  0.22,  0.38,
 
 4735      . 0.11,  0.58,  0.07,  0.25,  0.36,  0.58,  0.22,  0.38,
 
 4736      . 0.11,  0.58,  0.07,  0.25,  0.36,  0.58,  0.22,  0.38,
 
 4737      . 0.11,  0.58,  0.07,  0.25,  0.36,  0.58,  0.22,  0.38,
 
 4738      . 0.11,  0.58,  0.07,  0.25,  0.36,  0.58,  0.22,  0.38,
 
 4739      . 0.10,  0.45,  0.05,  0.25,  0.16,  0.39,  0.01,  0.01,
 
 4740      . 0.10,  0.45,  0.05,  0.25,  0.16,  0.39,  0.01,  0.01,
 
 4741      . 0.10,  0.45,  0.05,  0.25,  0.16,  0.39,  0.01,  0.01,
 
 4742      . 0.07,  0.35,  0.05,  0.10,  0.10,  0.39,  0.01,  0.01,
 
 4743      . 0.07,  0.35,  0.05,  0.10,  0.10,  0.39,  0.01,  0.01,
 
 4744      . 0.07,  0.35,  0.05,  0.10,  0.10,  0.39,  0.01,  0.01/
 
 4747      .reviss,renirs,trviss,trnirs
 
 4749      ./0.85,  0.85,  0.00,  0.00/
 
 4764             e_prad = 2.5   *  
coszsv(ikl)       
 
 4767             exdrad =    exp(-k_drad*
lai_sv(ikl))
 
 4768             e1prad = 1.-exp(-e_prad)            
 
 4771             zv_fac =    min( 
sncasv(ikl)/sncamx 
 
 4773             zv1fac = 1.     -       zv_fac      
 
 4774             deadlf = 1.     -       
glf_sv(ikl) 
 
 4780             a_rad0 =      0.25 + 0.697 * e1prad 
 
 4781             t_rad0 = 1. - a_rad0                
 
 4785             re_rad = 
glf_sv(ikl) *revisl(ivg) 
 
 4786      .             + deadlf      *revisd(ivg)
 
 4787             tr_rad = 
glf_sv(ikl) *trvisl(ivg) 
 
 4788      .             + deadlf      *trvisd(ivg)
 
 4792             re_rad = zv1fac      *re_rad      + zv_fac *reviss
 
 4793             tr_rad = zv1fac      *tr_rad      + zv_fac *trviss
 
 4797             r__rad = (2. *re_rad +     tr_rad) / 3.    
 
 4798             t__rad = (    re_rad + 2. *tr_rad) / 3.    
 
 4801             arggam =      t1_rad*t1_rad-r__rad*r__rad  
 
 4802             arggam =  max(arggam,
zer0)                 
 
 4803             gamma  = sqrt(arggam)                      
 
 4804             gammal =  min( gamma*
lai_sv(ikl),40.0)     
 
 4805             ddifh1 =  exp( gammal           )          
 
 4806             ddifh2 =  exp(-gammal           )          
 
 4812             r0_rad = 0.5 *((re_rad+tr_rad) *k_drad     
 
 4813      .                    +(re_rad-tr_rad) /    3.)    
 
 4814             t0_rad = 0.5 *((re_rad+tr_rad) *k_drad     
 
 4815      .                    -(re_rad-tr_rad) /    3.)    
 
 4817             nu_rad = t1_rad-r__rad*
albisv(ikl)         
 
 4818             den_c1 =  gamma*(ddifh1+ddifh2)            
 
 4819      .              +nu_rad*(ddifh1-ddifh2)            
 
 4821             densig =  gamma*gamma - k_drad*k_drad      
 
 4822             dens_s = sign(
un_1,densig)                 
 
 4823             dens_a =  abs(     densig)                 
 
 4824             densig =  max(
eps6,dens_a) * dens_s        
 
 4825             sig__c = (r__rad* r0_rad                   
 
 4826      .               +t0_rad*(k_drad+t1_rad)) / densig 
 
 4828             ddifc1 = ((gamma-nu_rad)*(t_rad0-sig__c*a_rad0)*ddifh2
 
 4829      .             +((k_drad-nu_rad)* sig__c
 
 4830      .               +t0_rad+r__rad * 
albisv(ikl)) *a_rad0 *exdrad)
 
 4832             ddifc2 =  t_rad0        - ddifc1-sig__c*a_rad0
 
 4836             ddif_l =  ddifc1*ddifh1 + ddifc2*ddifh2    
 
 4837      .             +  sig__c*a_rad0 *exdrad            
 
 4838             u0_vis = ((gamma+t1_rad)*ddifc1            
 
 4839      .               -(gamma-t1_rad)*ddifc2            
 
 4840      .             -((k_drad-t1_rad)*sig__c            
 
 4842      .          / max(r__rad,
eps6)                     
 
 4843             u0_vis = min(0.99,max(
eps6,u0_vis))        
 
 4844             absg_v = (1.-
albisv(ikl))*(a_rad0*exdrad   
 
 4846             absv_v = (1.-u0_vis     )- absg_v          
 
 4851             c1__sv(ikl) = ddifc1
 
 4852             c2__sv(ikl) = ddifc2
 
 4853             sigcsv(ikl) = sig__c
 
 4854             k___sv(ikl) = k_drad
 
 4855             a0__sv(ikl) = a_rad0
 
 4861             a_rad0 =      0.80 + 0.185 * e1prad 
 
 4862             t_rad0 = 1. - a_rad0                
 
 4866             re_rad = 
glf_sv(ikl) *renirl(ivg) 
 
 4867      .             + deadlf      *renird(ivg)
 
 4868             tr_rad = 
glf_sv(ikl) *trnirl(ivg) 
 
 4869      .             + deadlf      *trnird(ivg)
 
 4873             re_rad = zv1fac      *re_rad      + zv_fac *renirs
 
 4874             tr_rad = zv1fac      *tr_rad      + zv_fac *trnirs
 
 4878             r__rad = (2. *re_rad +     tr_rad) / 3.    
 
 4879             t__rad = (    re_rad + 2. *tr_rad) / 3.    
 
 4882             arggam =      t1_rad*t1_rad-r__rad*r__rad  
 
 4883             arggam =  max(arggam,
zer0)                 
 
 4884             gamma  = sqrt(arggam)                      
 
 4885             ddifh1 =  exp( gamma*
lai_sv(ikl))          
 
 4886             ddifh2 =  exp(-gamma*
lai_sv(ikl))          
 
 4892             r0_rad = 0.5 *((re_rad+tr_rad) *k_drad     
 
 4893      .                    +(re_rad-tr_rad) /    3.)    
 
 4894             t0_rad = 0.5 *((re_rad+tr_rad) *k_drad     
 
 4895      .                    -(re_rad-tr_rad) /    3.)    
 
 4897             nu_rad = t1_rad-r__rad*
albisv(ikl)         
 
 4898             den_c1 =  gamma*(ddifh1+ddifh2)            
 
 4899      .              +nu_rad*(ddifh1-ddifh2)            
 
 4901             densig =  gamma*gamma - k_drad*k_drad      
 
 4902             dens_s = sign(
un_1,densig)                 
 
 4903             dens_a =  abs(     densig)                 
 
 4904             densig =  max(
eps6,dens_a) * dens_s        
 
 4905             sig__c = (r__rad* r0_rad                   
 
 4906      .               +t0_rad*(k_drad+t1_rad)) / densig 
 
 4908             ddifc1 = ((gamma-nu_rad)*(t_rad0-sig__c*a_rad0)*ddifh2
 
 4909      .             +((k_drad-nu_rad)* sig__c
 
 4910      .               +t0_rad+r__rad * 
albisv(ikl)) *a_rad0 *exdrad)
 
 4912             ddifc2 =  t_rad0        - ddifc1-sig__c*a_rad0
 
 4916             ddif_l =  ddifc1*ddifh1 + ddifc2*ddifh2    
 
 4917      .             +  sig__c*a_rad0 *exdrad            
 
 4918             u0_nir = ((gamma+t1_rad)*ddifc1            
 
 4919      .               -(gamma-t1_rad)*ddifc2            
 
 4920      .             -((k_drad-t1_rad)*sig__c            
 
 4922      .          / max(r__rad,
eps6)                     
 
 4923             u0_nir = min(0.99,max(
eps6,u0_nir))        
 
 4924             absgni = (1.-
albisv(ikl))*(a_rad0*exdrad   
 
 4926             absvni = (1.-u0_nir     )- absgni          
 
 4932             alb_sv(ikl) = (u0_vis+u0_nir)*0.5d0
 
 4933             socasv(ikl) = (absv_v+absvni)*0.5d0
 
 4934             sososv(ikl) = (absg_v+absgni)*0.5d0
 
 4950               argexg      =  min(crilai(ikl)*gamasv(ikl),     
ea_max)
 
 4951               argexk      =  min(crilai(ikl)*k___sv(ikl),     
ea_max)
 
 4952               residu      =      c1__sv(ikl)            *exp( argexg)
 
 4953      .                          +c2__sv(ikl)            *exp(-argexg)
 
 4954      .                          +a0__sv(ikl)*gamasv(ikl)*exp(-argexk)
 
 4957               d_ddif      =      c1__sv(ikl)*gamasv(ikl)*exp( argexg)
 
 4958      .                          -c2__sv(ikl)*gamasv(ikl)*exp(-argexg)
 
 4959      .                          -a0__sv(ikl)*k___sv(ikl)*exp(-argexk)
 
 4960               dddifs      = sign(
un_1,d_ddif)
 
 4961               dddifa      =  abs(     d_ddif)
 
 4962               d_ddif      =  max(
eps6,dddifa) * dddifs
 
 4964               crilai(ikl) =      crilai(ikl)-residu/d_ddif
 
 4965               crilai(ikl) =  max(crilai(ikl),
zer0       )
 
 4966               crilai(ikl) =  min(crilai(ikl),
lai_sv(ikl))
 
 4972               laiesv(ikl) = crilai(ikl) +(exp(-k___sv(ikl)*crilai(ikl))
 
 4973      .                                   -exp(-k___sv(ikl)*
lai_sv(ikl)))
 
 5029       integer   ikl    ,ist    ,ist__s ,ist__w
 
 5031       real      uustar ,thstar ,qqstar ,ssstar
 
 5032       real      thstarv,thstars,thstara
 
 5033       real      zeta   ,zeta_S ,zeta_A
 
 5034       real      fCdCdP ,Cd_min ,cCdUns
 
 5054         ist__s =  min(ist, 1)                        
 
 5062         rapcm0 =      rapcm0     *rapcm0             
 
 5064         cd_m   =  max(cd_min*rapcm0,                 
 
 5065      .                fcdcdp*rapcm0*
vv__sv(ikl)  )   
 
 5066      .          *(1.+max(min(d_tats,
zer0),ccduns)    
 
 5079         us__sv(ikl) = sqrt(uustar)
 
 5092         thstarv     = thstar      + 
tat_sv(ikl) *(0.608*qqstar)
 
 5093         thstars     =     sign(
un_1,thstarv)
 
 5094         thstara     =      abs(     thstarv)
 
 5095         thstarv     =      max(
eps6,thstara)    *thstars
 
 5102         zeta_s      =          sign(
un_1  ,zeta)
 
 5104         zeta        = zeta_s  * max(
eps6  ,zeta_a)
 
 5218       real      VVaSBL(
klonv),VVa_OK          
 
 5224       real      uustar,thstar,qqstar,ssstar,thstarv,thstars,thstara
 
 5225       real      zetam ,zetah ,zeta_S,zeta_A,zeta0m ,zeta0h
 
 5226       real      psim_s,xpsimi,psim_i,psim_z
 
 5227       real      psis_s,psis_z,psis_0
 
 5228       real      psih_s,xpsihi,psih_i,psih_z
 
 5233       real      dustar,u0star,uTstar,usstar
 
 5234       real      sss__F,sss__N,usuth0
 
 5239       real      coef_m,coef_h,stab_s
 
 5242       real      fac_Ri,vuzvun,Kz_vun
 
 5267       data      zetmax/ 0.0428/   
 
 5281         vvasbl(ikl)   = 
vv__sv(ikl)
 
 5283         vvasbl(ikl)   = 
vvmmem(ikl)
 
 5286         dta_ts(ikl)   = 
dtmmem(ikl)
 
 5322           u0star      =      max(
eps6,u0star)      
 
 5323           uustar      = u0star      * u0star       
 
 5324           thstar      = 
uts_sv(ikl) / u0star       
 
 5325           qqstar      = 
uqs_sv(ikl) / u0star       
 
 5326           ssstar      = 
uss_sv(ikl) / u0star       
 
 5334           thstarv     = thstar      + theta0      *(0.608*qqstar)
 
 5336           thstars     =     sign(
un_1,thstarv)
 
 5337           thstara     =      abs(     thstarv)
 
 5338           thstarv     =      max(
eps6,thstara)*thstars
 
 5350           zetam       =           min(zetmax,zetah)
 
 5353           lmomom(ikl) = 
za__sv(ikl) /(max(
eps6,abs(zetam))
 
 5354      .                              *sign(
un_1,    zetam ))
 
 5355           zeta0m      = 
z0m_sv(ikl) / lmomom(ikl)
 
 5360           stab_s      =  max(
zer0,sign(
un_1,zetam))
 
 5363           xpsimi      = sqrt(sqrt(
un_1-coef_m*min(
zer0,zetam)))
 
 5364           psim_i      =   2. *log(
half*(
un_1+xpsimi)) 
 
 5367           psim_z      =    stab_s*psim_s+(1.-stab_s)*psim_i
 
 5373           xpsimi      = sqrt(sqrt(
un_1-coef_m*min(
zer0,zeta0m)))
 
 5374           psim_i      =   2. *log(
half*(
un_1+xpsimi)) 
 
 5377           psim_0      =    stab_s*psim_s+(1.-stab_s)*psim_i
 
 5432  6600     
format(/,
' ** SISVATeSBL *0  ' 
 5433      .            ,
'  Z0m_SV  = ',e12.4,
'  psim_z  = ',e12.4
 
 5434      .            ,
'  LMO_SV  = ',e12.4,
'  uustar  = ',e12.4
 
 5436      .            ,
'  sqrCm0  = ',e12.4,
'  psim_0  = ',e12.4
 
 5437      .            ,
'  LMOmom  = ',e12.4,
'  thstarv = ',e12.4)
 
 5445           vva_ok      =  max(0.000001,       vvasbl(ikl))
 
 5447           sss__f      = (
sqrcm0(ikl) - psim_z + psim_0)
 
 5448           usuth0      =  sss__n /sss__f                
 
 5528  6000     
format(a18,i3,6
x,
'u*   [m/s] =',f6.3,
'   hSalt[mm]='  ,e9.3,
 
 5529      .                  
'   Z0m   [mm] =',f9.3,
'   q   [g/kg] =',f9.3,
 
 5530      .               /,91
x,                    
'   qSa [g/kg] =',f9.3,
 
 5531      .               /,27
x, 
'ut*[m/s]='  ,e9.3,
'   u*-ut*   ='  ,e9.3,
 
 5532      .                  
'   s*  [g/kg] =',f9.3,
'   us* [mm/s] =',f9.3)
 
 5568  6001     
format(18
x,9
x,
'LMO  [m]=',f9.1,
'   zetah[-] =',f9.3)
 
 5577           stab_s      =  max(
zer0,sign(
un_1,zetam))
 
 5580           xpsimi      = sqrt(sqrt(
un_1-coef_m*min(
zer0,zetam)))
 
 5581           psim_i      =   2. *log(
half*(
un_1+xpsimi)) 
 
 5584           psim_z      =    stab_s*psim_s+(1.-stab_s)*psim_i
 
 5587           xpsimi      = sqrt(sqrt(
un_1-coef_m*min(
zer0,zeta0m)))
 
 5588           psim_i      =   2. *log(
half*(
un_1+xpsimi)) 
 
 5591           psim_0      =    stab_s*psim_s+(1.-stab_s)*psim_i
 
 5595           stab_s      =  max(
zer0,sign(
un_1,zetah))
 
 5598           xpsihi      = sqrt(sqrt(
un_1-coef_h*min(
zer0,zetah)))
 
 5599           psih_i      =   2. *log(
half*(
un_1+xpsihi)) 
 
 5600           psih_z      =    stab_s*psih_s+(1.-stab_s)*psih_i
 
 5603           xpsihi      = sqrt(sqrt(
un_1-coef_h*min(
zer0,zeta0h)))
 
 5604           psih_i      =   2. *log(
half*(
un_1+xpsihi)) 
 
 5605           psih_0      =    stab_s*psih_s+(1.-stab_s)*psih_i
 
 5610      .                        /(
sqrch0(ikl)-psih_z+psih_0)
 
 5620           thstar      = 
rcdhsv(ikl) * dta_ts(ikl)
 
 5628           dustar      = max(dustar,abs(
us__sv(ikl)-u0star))
 
 5644  6004     
format(122(
'-'))
 
 5647  6003     
format(
'   V  Ta-Ts  Z0      It' 
 5648      .   ,
' du*     u*    sss__F   CD       Qss       Qs*     ' 
 5649      .   ,
' PseudOL Full-OL zetam   zetah   psim_z  psih_z')
 
 5658  6002     
format(2f6.1,f8.4,i3,f9.6,f6.3,f9.3,3f9.6,2f8.2,2f8.4,2f8.2)
 
 5663  6014     
format(100(
'-'))
 
 5666  6013     
format(
'   V  Ta-Ts  Z0      It' 
 5667      .   ,
' du*     u*    sss__F   W_NUs1   W_NUs2   W_NUs3      ' 
 5668      .   ,
' W_DUs1     W_DUs2 ')
 
 5676  6012     
format(2f6.1,f8.4,i3,f9.6,f6.3,f9.3,3f9.3,2f12.3)
 
 5703           ram_sv(ikl) = 1./(cdm(ikl)*max(vvasbl(ikl),
eps6))
 
 5704           rah_sv(ikl) = 1./(cdh(ikl)*max(vvasbl(ikl),
eps6))
 
 5778       integer   ikl    ,ist    ,ist__s ,ist__w
 
 5779       real      CD_m_0 ,CD_h_0 ,ram0   ,rah0  ,rahMIN
 
 5780       real      d_TaTs ,RiB__D ,RiBulk
 
 5781       real      bmstab ,Am1_FU ,Am2_FU ,Fm_Uns
 
 5782       real      bhstab ,Ah1_FU ,Ah2_FU ,Fh_Uns,dFh_Un
 
 5783       real      Aux_FS ,FStabl ,dFSdRi ,Stabil,Fm_loc
 
 5784       real      uustar ,thstar ,qqstar ,ssstar
 
 5785       real      thstarv,thstars,thstara
 
 5786       real      zeta   ,zeta_S ,zeta_A
 
 5804         ist__s =  min(ist, 1)                        
 
 5811         ram0   =  1.0 /    (cd_m_0     *
vv__sv(ikl)) 
 
 5812         rah0   =  1.0 /    (cd_h_0     *
vv__sv(ikl)) 
 
 5828  6600       
format(/,
'Tem(s,a), Wind  , d_TaTs, RiBulk = ',5e15.6)
 
 5832         bmstab =  ist__s * (13.7 -0.34 /sqrt(cd_m_0))
 
 5834         bmstab =  10.    *  bmstab   * cd_m_0        
 
 5836         am1_fu =  bmstab *    sqrt(abs(ribulk))      
 
 5837         am2_fu =  am1_fu +1.0 +10.*abs(ribulk)       
 
 5838         fm_uns = (am1_fu +1.0)/        am2_fu        
 
 5845  6601       
format(/,
'CD_m_0  , Z0m_SV, bmstab, ist/sw = ',3e15.6,2i15)
 
 5847         bhstab =  ist__s * ( 6.3 -0.18 /sqrt(cd_h_0))
 
 5849         bhstab =  10.    *  bhstab   * cd_h_0        
 
 5851         ah1_fu =  bhstab *    sqrt(abs(ribulk))      
 
 5852         ah2_fu =  ah1_fu +1.0 +10.*abs(ribulk)       
 
 5853         fh_uns = (ah1_fu +1.0)/        ah2_fu        
 
 5854         dfh_un =((ah1_fu +2.0)/(ah2_fu*ah2_fu)) * 5. 
 
 5858         aux_fs =          1.0 + 5.*    ribulk
 
 5859         fstabl =                aux_fs*aux_fs
 
 5860         dfsdri =                aux_fs          *10.
 
 5864         stabil      =          sign(
un_1,d_tats)
 
 5865         fm_loc      =  fstabl * max(
zer0,stabil)
 
 5866      .               - fm_uns * min(
zer0,stabil)
 
 5868      .               - fh_uns * min(
zer0,stabil)
 
 5870      .               - dfh_un * min(
zer0,stabil)
 
 5877  6602       
format(/,
'FStabl  , Stabil, Fm_Uns, Fm_loc = ',4e15.6)
 
 5881         ram_sv(ikl) =  ram0   *          fm_loc
 
 5904  6603       
format(/,
'AeR(m,h), rCD(m,h)               = ',4e15.6)
 
 5913         us__sv(ikl) = sqrt(uustar)
 
 5926         thstarv     = thstar      + 
tat_sv(ikl) *(0.608*qqstar
 
 5928         thstars     =     sign(
un_1,thstarv)
 
 5929         thstara     =      abs(     thstarv)
 
 5930         thstarv     =      max(
eps6,thstara)    *thstars
 
 5937         zeta        =           min(zetmax,zeta)  
 
 5940         zeta_s      =          sign(
un_1  ,zeta)
 
 5942         zeta        = zeta_s  * max(
eps6  ,zeta_a)
 
 5950  6604       
format(/,
'***(m,h), LMO   , zeta           = ',4e15.6)
 
 6081       real      den_qs,arg_qs,qsatvg          
 
 6083       real      FacEvp,FacEvT,Fac_Ev          
 
 6089       real      rrCaOK,snCaOK,dEvpOK          
 
 6114           tveg_0(ikl) = 
tvegsv(ikl)
 
 6120           tau_ca = 1.  - 
tau_sv(ikl)                        
 
 6139           hsv_sv(ikl) = dhsdtv(ikl)                         
 
 6148           den_qs =         
tvegsv(ikl)- 35.8
 
 6149           arg_qs = 17.27 *(
tvegsv(ikl)-273.16)/den_qs
 
 6150           qsatvg = .0038 * exp(arg_qs)
 
 6151           dqs_dt = qsatvg* 4099.2    /(den_qs *den_qs)
 
 6167      .       + (1.-evfrac)*((1-snomsk)*         
rrcasv(ikl)  
 
 6172           facevp = fac_ev      *evfrac     / 
rah_sv(ikl)     
 
 6174           devpdt(ikl) = facevp* dqs_dt                       
 
 6175           facevt = fac_ev * (1.-evfrac)    /(
rah_sv(ikl)     
 
 6178           devtdt(ikl) = facevt* dqs_dt                       
 
 6181           dhldtv(ikl) = 
lhvh2o*(devpdt(ikl)+ devtdt(ikl))    
 
 6182      .                 +
lhfh2o* devpdt(ikl)* snomsk          
 
 6208           dhvdtv    =   dirdtv(ikl) * max(
eps_21,tau_ca)
 
 6216           d_tveg      = hv_tv0       / dhvdtv               
 
 6217           d_tveg      =      sign(
un_1,d_tveg)              
 
 6218      .                       *min( abs(d_tveg)     ,dtvmax) 
 
 6220           hv_max      = max(hv_max,abs(hv_tv0     ))        
 
 6241       IF (                     nit.lt.nitmax)    
GO TO 101                     
 
 6244      .   +dirdtv(ikl) *(
tvegsv(ikl)-tveg_0(ikl))            
 
 6246      .   -dhsdtv(ikl) *(
tvegsv(ikl)-tveg_0(ikl))            
 
 6248      .   +devpdt(ikl) *(
tvegsv(ikl)-tveg_0(ikl))            
 
 6250      .   +devtdt(ikl) *(
tvegsv(ikl)-tveg_0(ikl))            
 
 6252      .   -dhldtv(ikl) *(
tvegsv(ikl)-tveg_0(ikl))            
 
 6267           rrcaok =  max(
rrcasv(ikl), 0.)
 
 6268           sncaok =  max(
sncasv(ikl), 0.)
 
 6269           devpok =     (rrcaok-
rrcasv(ikl)
 
 6274      .    +(1.-snomsk)* 
lhvh2o            * devpok          
 
 6391       integer   ikl   ,isl   ,jsl   ,ist      
 
 6392       integer   ist__s,ist__w                 
 
 6400       real      mu_sno(
klonv),mu_aux          
 
 6410       REAL      TSurf0(
klonv),dTSurf          
 
 6411       real      qsatsg(
klonv),den_qs,arg_qs   
 
 6429       integer   nt_srf,it_srf,itEuBk          
 
 6431       real      agpsrf,xgpsrf,dt_srf,dt_ver   
 
 6454       data      eps__3 /   1.e-3   /          
 
 6455       data      mu_exp /  -0.4343  /          
 
 6456       data      mu_min /   0.172   /          
 
 6457       data      mu_max /   2.000   /          
 
 6458       data      ts_min / 175.      /          
 
 6459       data      ts_max / 300.      /          
 
 6472           mu__dz(ikl,isl) = 0.
 
 6474           dtc_sv(ikl,isl) = 
dtz_sv(isl)                   
 
 6482           ist__s =  min(ist, 1)                           
 
 6488           etamid =  max(etamid,
eps6)
 
 6491           mu_eta =      3.82      *(psimid)**mu_exp       
 
 6492           mu_eta =  min(max(mu_eta, mu_min), mu_max)      
 
 6494           mu_eta =  ist__s *mu_eta +ist__w * 
vk_dsv        
 6496           mu__dz(ikl,isl) = mu_eta/(
dzmisv(isl)           
 
 6499           dtc_sv(ikl,isl) = 
dtz_sv(isl)                   
 
 6515           ist__s =  min(ist, 1)                           
 
 6518           mu_eta =      3.82      *(psimid)**mu_exp       
 
 6519           mu_eta =  min(max(mu_eta, mu_min), mu_max)      
 
 6521           mu_eta =  ist__s *mu_eta +ist__w * 
vk_dsv        
 6527           mu_sno(ikl) =          max(
eps6,mu_sno(ikl))    
 
 6533           mu__dz(ikl,isl) = 2./(
dzsnsv(ikl,isl  )         
 
 6536      .                         *
dz_dsv(    isl-1)/mu_eta) 
 
 6552      .       * max(0,min(
isnosv(ikl)-isl+1,1))            
 
 6564      .      2.                        *mu_aux*mu_sno(ikl) 
 
 6566      .              +
dzsnsv(ikl,isl-1)*mu_aux     )       
 
 6567           mu_sno(ikl)     =            mu_aux             
 
 6571           dtc_sv(ikl,isl) = 
dt__sv/max(eps__3,            
 
 6581           mu__dz(ikl,
isnosv(ikl)+1) = 0.0
 
 6607           elem_a          =  dtc_sv(ikl,isl)         *mu__dz(ikl,isl)
 
 6608           elem_c          =  dtc_sv(ikl,isl)         *mu__dz(ikl,isl+1)
 
 6609           diag_a(ikl,isl) = -elem_a  *
implic 
 6610           diag_c(ikl,isl) = -elem_c  *
implic 
 6611           diag_b(ikl,isl) =  1.0d+0  -diag_a(ikl,isl)-diag_c(ikl,isl)
 
 6613      .                               +elem_c         *
tsissv(ikl,isl+1))
 
 6626           elem_c          =  dtc_sv(ikl,isl)         *mu__dz(ikl,isl+1)
 
 6627           diag_a(ikl,isl) =  0.
 
 6628           diag_c(ikl,isl) = -elem_c  *
implic 
 6629           diag_b(ikl,isl) =  1.0d+0  -diag_a(ikl,isl)-diag_c(ikl,isl)
 
 6641           elem_a          =  dtc_sv(ikl,isl)  *mu__dz(ikl,isl)
 
 6643           diag_a(ikl,isl) = -elem_a  *
implic 
 6644           diag_c(ikl,isl) =  0.
 
 6645           diag_b(ikl,isl) =  1.0d+0  -diag_a(ikl,isl)
 
 6661           irs__d(ikl) = dirsdt(ikl)* 
tsissv(ikl,isl) * 0.75    
 
 6674           f___hl(ikl) = f_hshl(ikl) *    
lx_h2o(ikl)
 
 6685           den_qs      =         
tsissv(ikl,isl)- 35.8          
 
 6686           arg_qs      = 17.27 *(
tsissv(ikl,isl)-273.16)        
 
 6688           qsatsg(ikl) = .0038 *        exp(arg_qs)             
 
 6689           dqs_dt(ikl) = qsatsg(ikl)* 4099.2   /(den_qs *den_qs)
 
 6696               agpsrf       = 
dt__sv*(   1.0-xgpsrf        )    
 
 6697      .                             /(   1.0-xgpsrf**nt_srf)    
 
 6703               etanew(ikl)  =          etabak(ikl)              
 
 6704               eteubk(ikl)  =          etanew(ikl)              
 
 6707               dt_ver       = dt_ver     +dt_srf                
 
 6709               faceta(ikl)  = fac_dt(ikl)*dt_srf                
 
 6723      .           /max(eteubk(ikl),
eps6))                       
 
 6725               psiarg(ikl) = 7.2e-5*psi(ikl)                    
 
 6726               rhusol(ikl) =   exp(-min(
ea_max,psiarg(ikl)))    
 
 6727               shusol(ikl) =     qsatsg(ikl)  *rhusol(ikl)      
 
 6729      .       (etanew(ikl) + faceta(ikl)*(
qat_sv(ikl)           
 
 6733      .      /(1.          + faceta(ikl)* shusol(ikl)           
 
 6737               eteubk(ikl) = eteubk(ikl) -
rootsv(ikl,0)         
 
 6742               etanew(ikl) =  max(eteubk(ikl),
eps6)             
 
 6744               dt_srf      =      dt_srf         * xgpsrf       
 
 6754           d__eta     =  
eta_sv(ikl,isl)-etanew(ikl)            
 
 6756      .                *(etanew(ikl)    -etabak(ikl)) / 
dt__sv   
 6757      .                 +ist__w         *f_hshl(ikl)            
 
 6758      .                *(
qat_sv(ikl)    -qsatsg(ikl))         ) 
 
 6777           tsurf0(ikl)     =  
tsissv(ikl,isl)
 
 6778           elem_a          =  dtc_sv(ikl,isl)*mu__dz(ikl,isl)
 
 6780           diag_a(ikl,isl) = -elem_a *
implic 
 6781           diag_c(ikl,isl) =  0.
 
 6782           diag_b(ikl,isl) =  1.0d+0 -diag_a(ikl,isl)
 
 6783           diag_b(ikl,isl) =  diag_b(ikl,isl)
 
 6784      .  + dtc_sv(ikl,isl) * (dirsdt(ikl)                       
 
 6789           term_d(ikl,isl) =  term_d(ikl,isl)
 
 6807           aux__p(ikl,-
nsol) = diag_b(ikl,-
nsol)
 
 6808           aux__q(ikl,-
nsol) =-diag_c(ikl,-
nsol)/aux__p(ikl,-
nsol)
 
 6813           aux__p(ikl,isl)   = diag_a(ikl,isl)  *aux__q(ikl,isl-1)
 
 6815           aux__q(ikl,isl)   =-diag_c(ikl,isl)  /aux__p(ikl,isl)
 
 6825           tsissv(ikl,isl)   =(term_d(ikl,isl)
 
 6826      .                       -diag_a(ikl,isl)  *
tsissv(ikl,isl-1))
 
 6844           dtsurf            = 
tsissv(ikl,isl) -     tsurf0(ikl)
 
 6845           tsissv(ikl,isl)   = tsurf0(ikl) + sign(1.,dtsurf) 
 
 6846      .              * min(abs(dtsurf),5.e-2*
dt__sv)         
 
 6847           IF (abs(dtsurf) > 5.e-2*
dt__sv) 
THEN 
 6848             write(*,*) 
'abrupt',ikl,
'dTs ',dtsurf,
tsissv(ikl,isl) 
 
 6855           IF (ts_min >  
tsissv(ikl,isl)) 
THEN 
 6856              write(*,*) 
'cold', ikl, 
'couche',isl, 
tsissv(ikl,isl)
 
 6858           IF (ts_max <  
tsissv(ikl,isl)) 
THEN 
 6859              write(*,*) 
'hot ', ikl, 
'couche',isl, 
tsissv(ikl,isl) 
 
 6872           irs_sv(ikl) = irs__d(ikl)                          
 
 6873      .                - dirsdt(ikl) * 
tsissv(ikl,isl)        
 
 6874           hss_sv(ikl) = hs___d(ikl)                          
 
 6876           hls_sv(ikl) = hl___d(ikl)                          
 
 7025       REAL, 
DIMENSION(klonv)   :: zx_mh, zx_nh, zx_oh
 
 7026       REAL, 
DIMENSION(klonv)   :: zx_mq, zx_nq, zx_oq
 
 7027       REAL, 
DIMENSION(klonv)   :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
 
 7028       REAL, 
DIMENSION(klonv)   :: zx_sl, zx_k1
 
 7029       REAL, 
DIMENSION(klonv)   :: d_ts
 
 7030       REAL                     :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
 
 7031       REAL                     :: qsat_new, q1_new
 
 7034       REAL, 
DIMENSION(klonv)   :: IRs__D, dIRsdT
 
 7072         ztherm_i(ig)   = inertie_ice
 
 7099      &          *(1.-d_coef(ig,jk)))
 
 7101      &          (
tsissv(ig,jk)*zdz2(ig,jk)                              &
 
 7102      &          +
dz1_sv(ig,jk)*c_coef(ig,jk)) * z1s
 
 7103           d_coef(ig,jk+1)=
dz1_sv(ig,jk+1)*z1s
 
 7114       mu=1./((2.**1.5-1.)/(2.**(0.5)-1.)-1.)
 
 7127         IF (
isnosv(ig).GT.0) 
THEN 
 7128           IF (
isnosv(ig).GT.1) 
THEN 
 7137         IF (mug(ig)  .LE. 0.05) 
THEN 
 7138          write(*,*)
'Attention mu low', mug(ig)
 
 7140         IF (mug(ig)  .GE. 0.98) 
THEN 
 7141          write(*,*)
'Attention mu high', mug(ig)
 
 7145      &           *min(max(
isnosv(ig),0),1)+                             &
 
 7155      &        (mug(ig)*(1.-d_coef(ig,
isnosv(ig)))+1.)
 
 7161           tsissv(ig,jk-1)=c_coef(ig,jk)+d_coef(ig,jk)                   &
 
 7185           irs__d(ig) = dirsdt(ig)* 
tsissv(ig,isl) * 0.75                  
 
 7202      &          *(1.-d_coef(ig,jk)))
 
 7203           c_coef(ig,jk+1) = (
tsissv(ig,jk)*zdz2(ig,jk)+                 &
 
 7204      &          
dz1_sv(ig,jk)*c_coef(ig,jk)) * z1s
 
 7205           d_coef(ig,jk+1) = 
dz1_sv(ig,jk+1)*z1s
 
 7220         pcapcal(ig)  = ztherm_i(ig)*                                    &
 
 7223         z1s = mug(ig)*(1.-d_coef(ig,
isnosv(ig)))+1.
 
 7224         pcapcal(ig)  = pcapcal(ig)/z1s
 
 7225         pfluxgrd(ig) = ( pfluxgrd(ig)                                   &
 
 7227      &       - mug(ig)* c_coef(ig,
isnosv(ig))                           &
 
 7251         IF (
ps__sv(ig).LT.1.) 
THEN 
 7255         IF (
p1l_sv(ig).LT.1.) 
THEN 
 7259         IF (
tat_sv(ig).LT.180.) 
THEN 
 7263         IF (
qat_sv(ig).LT.1.e-8) 
THEN 
 7267         IF (
tsf_sv(ig).LT.100.) 
THEN 
 7271         IF (
tsf_sv(ig).GT.500.) 
THEN 
 7279          IF (
cdh_sv(ig).LT.1.e-10) 
THEN 
 7300           zcvm5 = r5les*
rlvtt*(1.-zdelta) + r5ies*
rlstt*zdelta
 
 7301           zcvm5 = zcvm5 / 
rcpd / (1.0+rvtmp2*
qat_sv(ig))
 
 7303           zx_qs=min(0.5,zx_qs)
 
 7305           zcor=1./(1.-
retv*zx_qs)
 
 7307           zx_dq_s_dh = foede(
tsf_sv(ig),zdelta,zcvm5,zx_qs,zcor)        &
 
 7308      &          /
rlvtt / zx_pkh(ig)
 
 7310           IF (
tsf_sv(ig).LT.t_coup) 
THEN 
 7320         zx_dq_s_dt(ig) = 
rcpd * zx_pkh(ig) * zx_dq_s_dh
 
 7325         zx_coef(ig) = 
cdh_sv(ig) *                                      &
 
 7339         zx_k1(ig) = zx_coef(ig)
 
 7345         zx_oq(ig) = 1. - (beta(ig) * zx_k1(ig) * 
bcoqsv(ig) * 
dt__sv)
 
 7346         zx_mq(ig) = beta(ig) * zx_k1(ig) *                              &
 
 7347      &       (
acoqsv(ig) - zx_qsat(ig) +                                &
 
 7348      &       zx_dq_s_dt(ig) * 
tsf_sv(ig))                               &
 
 7350         zx_nq(ig) = beta(ig) * zx_k1(ig) * (-1. * zx_dq_s_dt(ig))       &
 
 7355         zx_mh(ig) = zx_k1(ig) * 
acohsv(ig) / zx_oh(ig)
 
 7356         zx_nh(ig) = - (zx_k1(ig) * 
rcpd * zx_pkh(ig))/ zx_oh(ig)
 
 7360      &       (
rsolsv(ig) + zx_mh(ig) + zx_sl(ig) * zx_mq(ig))           & 
 
 7361      &       + dif_grnd(ig) * t_grnd * 
dt__sv)/                         &
 
 7362      &       ( 1. - 
dt__sv * cal(ig)/(
rcpd * zx_pkh(ig)) *              &
 
 7363      &       (zx_nh(ig) + zx_sl(ig) * zx_nq(ig))                        &  
 
 7364      &       + 
dt__sv * dif_grnd(ig))
 
 7380         dldtsv(ig) = zx_sl(ig) * zx_nq(ig)
 
 7386         irs_sv(ig) = irs__d(ig)                                         &
 
 7387      &                - dirsdt(ig) * 
tsfnsv(ig) 
 
 7399         qsat_new=zx_qsat(ig) + zx_dq_s_dt(ig) * d_ts(ig)
 
 7401         qat_sv(ig)=q1_new*(1.-beta(ig)) + beta(ig)*qsat_new 
 
 7512       real      den_qs,arg_qs,qsatvg          
 
 7532       data      bw_min /   4.e-8 /            
 
 7548           psiv_0(ikl) = 
psivsv(ikl)
 
 7567           rootok =      max(
zer0, sign(
un_1,psidif))
 
 7569           plantw(ikl)     = plantw(ikl) + 
rootsv(ikl ,isl)  
 
 7570           dpdpsi(ikl)     = dpdpsi(ikl) + rootok*root_w     
 
 7581           den_qs =         
tvegsv(ikl)- 35.8
 
 7582           arg_qs = 17.27 *(
tvegsv(ikl)-273.16)/den_qs
 
 7583           qsatvg = .0038 * exp(arg_qs)
 
 7590           f___ok =     max(
zer0,sign(
un_1,denomf))
 
 7591           denomf =     max(
eps6,          denomf)           
 
 7593           dfdpsi =             -f_stom  / denomf            
 
 7596           r_stom = numerr *                  f_stom         
 
 7598           drdpsi = r_stom *                  dfdpsi         
 
 7607           facevt = fac_ev * (1.-evfrac)    / denome         
 
 7608           evtran = facevt     *(qsatvg     - 
qat_sv(ikl))   
 
 7609           dedpsi =(evtran     / denome)    * drdpsi         
 
 7615           bwater      =( plantw(ikl)                        
 
 7616      .                 - evtran     )*        f___ok        
 
 7619      .                      sign(
un_1,    abs(bwater)       
 
 7626           dbwdpv      = dpdpsi(ikl)
 
 7628           dbwdpv      = sign(  
un_1,    dbwdpv)             
 
 7629      .                *  max(
eps_21,abs(dbwdpv))            
 
 7635           d_psiv      = bwater       / dbwdpv               
 
 7636           d_psiv      =      sign(
un_1,d_psiv)              
 
 7637      .                       *min( abs(d_psiv)     ,dpvmax) 
 
 7639           bw_max      = max(bw_max,abs(bwater))
 
 7649      .                          /max(
eps_21,plantw(ikl))    
 
 7654       IF (                     nit.lt.nitmax)    
GO TO 101                    
 
 7772       integer  noSnow(
klonv)                 
 
 7806       integer             isnnew,isinew,isnUpD,isnitr
 
 7895           eexdum(ikl) = eexdum(ikl)     -          
eexcsv(ikl)  
 
 7941           dzmelt(ikl) = enmelt      / max(snhlat,    
eps6 )
 
 7942           nosnow(ikl) = nosnow(ikl) 
 
 7943      .      + max(
zer0  ,sign(
un_1,dzmelt(ikl)                
 
 7945      .       *min(1     , max(0 ,1+
isnosv(ikl)-isn))          
 
 7947      .              min(
dzsnsv(ikl, isn),dzmelt(ikl))
 
 7949      .                  
dzsnsv(ikl,isn) -dzmelt(ikl)
 
 7967      .                           *max(0,2-
istosv(ikl,isn)          ) 
 
 7969      .        (1.-okmelt) *               
istosv(ikl,isn)            
 
 7970      .      +     okmelt  *((1-k_face) *  
istdsv(2)                  
 
 7978           layrok      = min(   1, max(0          , 
isnosv(ikl)-isn+1))
 
 7980           wafrez      =   -(     enfrez          * layrok / 
lhfh2o)
 
 7986           rdznew      = wafrez + rdzsno
 
 8001           porvol      = 1.     - rosdry          / 
rhoice           
 8002           porvol      =      max(porvol          , 
zer0  )         
 
 8011      .                         + rosdry          * 
dzsnsv(ikl,isn)
 
 8021      .                +      max(
ispisv(ikl),isn)    * pclose
 
 8043      .             *     max(0     , min(1   ,
isnosv(ikl) +1 -isn ))
 
 8044            isnupd =      max(isnupd,          isnnew)
 
 8046            isinew =      isn*isnupd *max(0, 1-isinew)      
 
 8054            dzsnsv(ikl,isn+isnnew) =(1-isnnew)*
dzsnsv(ikl,isn+isnnew)
 
 8055            ro__sv(ikl,isn+isnnew) =(1-isnnew)*
ro__sv(ikl,isn+isnnew)
 
 8056            eta_sv(ikl,isn+isnnew) =(1-isnnew)*
eta_sv(ikl,isn+isnnew)
 
 8057            g1snsv(ikl,isn+isnnew) =(1-isnnew)*
g1snsv(ikl,isn+isnnew)
 
 8058            g2snsv(ikl,isn+isnnew) =(1-isnnew)*
g2snsv(ikl,isn+isnnew)
 
 8062      .    -isnupd *max(0,min(
ispisv(ikl)-isinew,1))        
 
 8074            nh =    nh + isn* min(
istosv(ikl,isn)-1,1)*max(0,1-nh)
 
 8080      .                     * max(0,min(1,nh+1-isn))
 
 8141  440      
format(
'iSupI    i       dz       ro      eta',
 
 8142      .            
'   PorVol   zSlush     ro_n    eta_n',2
x,a18)
 
 8177  441      
format(2i5,f9.3,f9.1,f9.6,f9.3,f9.6,f9.1,f9.6)  
 
 8190      .  * 
hls_sv(ikl)         * min(isn             , 1   )
 
 8265         z_melt = lastok     *rapdok*thinok
 
 8266         nosnow(ikl)   = nosnow(ikl)   + z_melt
 
 8267         z_melt        =                 z_melt *
dzsnsv(ikl,1) 
 
 8286  6005   
format(i3,
' (noSnow) ')
 
 8290      .     * min(1,iabs(
isnosv(ikl)-nosnow(ikl)))
 
 8330  43     
format(
'SubRoutine SISVAT_qSn: Local Energy and Water Budgets',
 
 8331      .       /,
'=====================================================')
 
 8335  431  
format(
' WARNING (1) in _qSn,',         a18,
 
 8336      .       
': Energy Unbalance in Phase Change = ',e15.6)
 
 8340  432  
format(
' WARNING (2) in _qSn,',         a18,
 
 8341      .       
': Energy Unbalance in Phase Change = ',e15.6)
 
 8348  435    
format(11(
'-'),
'----------+-',3(
'-'),
'----------+-',
 
 8349      .          3(
'-'),
'----------+-',3(
'-'),
'----------+-',
 
 8350      .                 
'----------------+----------------+',
 
 8351      .       /,f8.2,3
x,
'EqSn_0(1) | ',3
x,
'EqSn_d(1) | ',
 
 8352      .              3
x,
'EqSn_1(1) | ',3
x,
'EExcsv(1) | ',
 
 8353      .                 
'E_0+E_d-E_1-EE  |   Water Budget |',
 
 8354      .       /,11(
'-'),
'----------+-',3(
'-'),
'----------+-',
 
 8355      .          3(
'-'),
'----------+-',3(
'-'),
'----------+-',
 
 8356      .                 
'----------------+----------------+')
 
 8365  436    
format(8
x,f12.0,
' + ',f12.0,
' - ',f12.0,
' - ',f12.3,
' = ',f12.3,
 
 8612       real      T3_xOK,T3__OK,T3_nOK          
 
 8614       real      dT1_OK,dT2_OK,dT3xOK,dT3_OK   
 
 8615       real      dT4xOK,dT4_OK,dT4nOK,AngSno   
 
 8616       real      G2_hds,SphrOK,HISupd          
 
 8617       real      H1a_OK,H1b_OK,H1__OK          
 
 8618       real      H23aOK,H23bOK,H23_OK          
 
 8620       real      H45_OK,H4__OK,H5__OK          
 
 8621       real      ViscSn,OK_Liq,OK_Ang,OKxLiq   
 
 8622       real      dSnMas,dzsnew,rosnew,rosmax   
 
 8633       real      vtang1,vtang2,vtang3,vtang4   
 
 8634       real      vtang5,vtang6,vtang7,vtang8   
 
 8635       real      vtang9,vtanga,vtangb,vtangc   
 
 8637       real      vgang1,vgang2,vgang3,vgang4   
 
 8638       real      vgang5,vgang6,vgang7,vgang8   
 
 8639       real      vgang9,vganga,vgangb,vgangc   
 
 8647       real      vvisc1,vvisc2,vvisc3,vvisc4   
 
 8648       real      vvisc5,vvisc6,vvisc7          
 
 8654       real      husi_0,husi_1,husi_2,husi_3   
 
 8680       data       husi_1 / 0.23873 /           
 
 8681       data       husi_2 / 4.18880 /           
 
 8682       data       husi_3 / 0.33333 /           
 
 8683       data       vtail1 / 1.28e-08/           
 
 8684       data       vtail2 / 4.22e-10/           
 
 8686       data       epsi5  / 1.0e-5  /           
 
 8695       data       vsphe2 / 1.0e9   /           
 
 8700       data       vtelv1 / 5.e-1   /           
 
 8709       data       vvisc1 / 0.70    /           
 
 8710       data       vvisc2 / 1.11e5  /           
 
 8711       data       vvisc3 /23.00    /           
 
 8712       data       vvisc4 / 0.10    /           
 
 8713       data       vvisc5 / 1.00    /           
 
 8714       data       vvisc6 / 2.00    /           
 
 8715       data       vvisc7 /10.00    /           
 
 8716       data       rovisc / 0.25    /           
 
 8768           ro_dry(ikl,isn) = 1.e-3 *
ro__sv(ikl,isn)    
 
 8770           etasno(ikl,isn) = 1.e-1 *
dzsnsv(ikl,isn)    
 
 8779           isnp   = min(isn+1,
isnosv(ikl))
 
 8781           dtsndz = abs( (
tsissv(ikl,isnp)-
tsissv(ikl,isn-1)) *2.e-2
 
 8799      .       /max(
eps6,ro_dry(ikl,isn))
 
 8803           exp1wa=   swater**nvdent1
 
 8804           ddendr=max(exp1wa/nvdent2,vdent1*exp(vvap1/
tf_sno))
 
 8814           dendrn= dendrn -ddendr *frac_j  
 
 8815           sphern= sphern +ddendr *frac_j  
 
 8821           g1__wd=ok__de *    (      -dendrn*
g1_dsv)   
 
 8824      .      +(1.-ok__de)*(
adsdsv-min(sphern,vsphe1))  
 
 8837           sphern= sphern +ddendr *frac_j
 
 8845      .            *(husi_2 *( 
g2snsv(ikl,isn)/husi_0)**3
 
 8846      .                      +(vtail1 +vtail2 *exp1wa    )*
dt__sv))
 
 8860           facvap=exp(vvap1/
tsissv(ikl,isn))
 
 8870           dendrn= dendrn-vdent1*facvap*frac_j
 
 8871           sphern= sphern+vsphe2*facvap*frac_j
 
 8877           g1_ldd= ok__de *    (      -dendrn*
g1_dsv)  
 
 8880      .       +(1.-ok__de)*(
adsdsv-min(sphern,vsphe1)) 
 
 8885           diamgx=
g2snsv(ikl,isn)*0.1
 
 8887           istook=min( abs(
istosv(ikl,isn)-
 
 8889           diamok=max(
zer0,  sign(
un_1,vdiam2-diamgx))
 
 8890           no_big=    istook+diamok
 
 8891           no_big=min(no_big,
un_1)
 
 8893           dspher=           vsphe2*facvap*frac_j      
 
 8894           spher0=    sphern+dspher                    
 
 8895           sphbig=    sphern+dspher                    
 
 8897           sphbig=       min(vsphe3,sphbig)            
 
 8898           sphern= no_big *  spher0
 
 8899      .      + (1.-no_big)*  sphbig
 
 8908           okmidt=               ok_mdt  *(1.-oklowt)  
 
 8909           okhigt=          (1. -ok_mdt) *(1.-oklowt)  
 
 8911           facvap=vdent1*exp(vvap1/
tsissv(ikl,isn))
 
 8912      .                 *   (1.e2 *dtsndz)**vvap2
 
 8922           dendrn= dendrn - facvap*frac_j
 
 8923           sphern= sphern - facvap*frac_j
 
 8929           g1_mdd= ok__de *    (      -dendrn*
g1_dsv)  
 
 8937           sphern=         sphern-facvap*frac_j
 
 8942           facvap=vdent1*exp(vvap1/
tsissv(ikl,isn))
 
 8943      .                 *   (1.e2 *dtsndz)**vvap2
 
 8953           dendrn= dendrn - facvap*frac_j              
 
 8954           sphern= sphern - facvap*frac_j              
 
 8960           g1_hdd= ok__de *    (      -dendrn*
g1_dsv)  
 
 8973           sphern= sphern - facvap*frac_j
 
 8982           t3__ok =                    t3_xok  * (1. - t2__ok)
 
 8983           t3_nok =              (1. - t3_xok) * (1. - t2__ok)
 
 8984           ro1_ok = max(
zer0,sign(
un_1,vrang1-ro_dry(ikl,isn)))
 
 8985           ro2_ok = max(
zer0,sign(
un_1,ro_dry(ikl,isn)-vrang2))
 
 8986           dt1_ok = max(
zer0,sign(
un_1,vgang1-dtsndz         ))
 
 8987           dt2_ok = max(
zer0,sign(
un_1,vgang2-dtsndz         ))
 
 8988           dt3xok = max(
zer0,sign(
un_1,vgang3-dtsndz         ))
 
 8989           dt3_ok =                    dt3xok  * (1. - dt2_ok)
 
 8990           dt4xok = max(
zer0,sign(
un_1,vgang4-dtsndz         ))
 
 8991           dt4_ok =                    dt4xok  * (1. - dt3_ok)
 
 8993           dt4nok =              (1. - dt4xok) * (1. - dt3_ok)
 
 9002      .     +t3__ok*(vtang7-vtang8*(
tf_sno-vtang2-
tsissv(ikl,isn))        
 
 9004      .     +t3_nok*(vtanga-vtangb*(
tf_sno-vtang3-
tsissv(ikl,isn))        
 
 9010      .        *(   ro2_ok*(1. - (ro_dry(ikl,isn)-vrang2)                 
 
 9016      .        *(   dt1_ok*(dt2_ok*vgang5*(dtsndz-vgang6)                 
 
 9023      .        *    dt1_ok*(dt3_ok*vgang8*(dtsndz-vgang2)
 
 9025      .                    +dt4_ok*vganga*(dtsndz-vgang3)
 
 9027      .                    +dt4nok*vgangc*(dtsndz-vgang4)
 
 9030           g2_hds = 
g2snsv(ikl,isn) + 1.d2 *angsno*vfi     *frac_j
 
 9039           g1snsv(ikl,isn) = wet_ok * (    ok__wd             *g1__wd    
 
 9040      .                               +(1.-ok__wd)*    ok__ws *g1__ws    
 
 9041      .                               +(1.-ok__wd)*(1.-ok__ws)*g1_bak)   
 
 9043      .                *(    oklowt  *(    ok_ldd             *g1_ldd    
 
 9044      .                               +(1.-ok_ldd)            *g1_lds)   
 
 9045      .                    + okmidt  *(    ok_mdd             *g1_mdd    
 
 9046      .                               +(1.-ok_mdd)            *g1_mds)   
 
 9047      .                    + okhigt  *(    ok_hdd             *g1_hdd    
 
 9048      .                               +(1.-ok_hdd)*    ok_hds *g1_hds    
 
 9049      .                               +(1.-ok_hdd)*(1.-ok_hds)*g1_bak))  
 
 9051           g2snsv(ikl,isn) = wet_ok * (    ok__wd             *g2__wd    
 
 9052      .                               +(1.-ok__wd)*    ok__ws *g2_bak    
 
 9053      .                               +(1.-ok__wd)*(1.-ok__ws)*g2__ws)   
 
 9055      .                *(    oklowt  *(    ok_ldd             *g2_ldd    
 
 9056      .                               +(1.-ok_ldd)            *g2_bak)   
 
 9057      .                    + okmidt  *(    ok_mdd             *g2_mdd    
 
 9058      .                               +(1.-ok_mdd)            *g2_bak)   
 
 9059      .                    + okhigt  *(    ok_hdd             *g2_hdd    
 
 9060      .                               +(1.-ok_hdd)*    ok_hds *g2_bak    
 
 9061      .                               +(1.-ok_hdd)*(1.-ok_hds)*g2_hds))  
 
 9122      .         /,
' isn     =  ',i4,6
x,
'(MAX.:',i4,
')',  
 
 9128      .         /,
' Histor. =  ',i4  ,
 
 9129      .         /,
' Grad(T) =  ',f8.3,
'                   ' ,18i3  ,
 
 9130      ./,         
'                       Current    Case: ',18f3.0,
 
 9131      ./,         
'                       Cases performed: ',18f3.0,
 
 9132      ./,
' ------------------------------------------------------------',
 
 9133      .             
'-----------+------------------+------------------+',
 
 9136      ./,
' ------------------------------------------------------------',
 
 9137      .             
'-----------+------------------+------------------+',
 
 9138      ./,
'    Wet_OK: ',f8.3,
'                  OK__wd: ',f8.3,
'       ',
 
 9139      .             
'           | G1__wd: ',f8.3,
' | G2__wd: ',f8.5,
' |',
 
 9140      ./,
'                                   1.-OK__wd: ',f8.3,
' OK__ws',
 
 9141      .             
': ',f8.3,
' | G1__ws: ',f8.3,
' |                  |',
 
 9143      .             
': ',f8.3,
' |                  | G2__ws: ',f8.5,
' |',
 
 9144      ./,
' 1.-Wet_OK: ',f8.3,
' OKlowT: ',f8.3,
' OK_ldd: ',f8.3,
'       ',
 
 9145      .             
'           | G1_ldd: ',f8.3,
' | G2_ldd: ',f8.5,
' |',
 
 9146      ./,
'                                   1.-OK_ldd: ',f8.3,
'       ',
 
 9147      .             
'           | G1_lds: ',f8.3,
' |                  |',
 
 9148      ./,
'                     OKmidT: ',f8.3,
' OK_mdd: ',f8.3,
'       ',
 
 9149      .             
'           | G1_mdd: ',f8.3,
' | G2_mdd: ',f8.5,
' |',
 
 9150      ./,
'                                   1.-OK_mdd: ',f8.3,
'       ',
 
 9151      .             
'           | G1_mds: ',f8.3,
' |                  |',
 
 9152      ./,
'                     OKhigT: ',f8.3,
' OK_hdd: ',f8.3,
'       ',
 
 9153      .             
'           | G1_hdd: ',f8.3,
' | G2_hdd: ',f8.5,
' |',
 
 9154      ./,
'                                   1.-OK_hdd: ',f8.3,
' OK_hds',
 
 9155      .             
': ',f8.3,
' | G1_hds: ',f8.3,
' |                  |',
 
 9157      .             
': ',f8.3,
' |                  | G2_hds: ',f8.5,
' |',
 
 9158      ./,
' ------------------------------------------------------------',
 
 9159      .             
'-----------+------------------+------------------+',
 
 9161      .             
'           |         ',f8.3,
' |         ',f8.5,
' |',
 
 9162      ./,
' ------------------------------------------------------------',
 
 9163      .             
'-----------+------------------+------------------+')
 
 9178           h1b_ok =     1   - min(1   ,       
istosv(ikl,isn))
 
 9179           h1__ok =                    h1a_ok*h1b_ok
 
 9182           h23bok = max(
zer0,sign(
un_1,etasno(ikl,isn)
 
 9185           h23_ok =                    h23aok*h23bok
 
 9186           h2__ok =     1   - min(1   ,       
istosv(ikl,isn))
 
 9193      .    sphrok*(h1__ok                             *
istdsv(1)
 
 9194      .       +(1.-h1__ok)*    h23_ok         *(h2__ok*
istdsv(2)
 
 9196      .       +(1.-h1__ok)*(1.-h23_ok) *h45_ok*(h4__ok*
istdsv(4)
 
 9198           istosv(ikl,isn) =   hisupd  +
 
 9211           IF  (
g1snsv(ikl,isn).ge.0.)                               
THEN 
 9212             IF(
g1snsv(ikl,isn).lt.vsphe4.and.
istosv(ikl,isn).eq.0)  
THEN 
 9215      .             etasno(ikl,isn)/
dzsnsv(ikl,isn).gt.vtelv1)       
THEN 
 9216               IF  (
istosv(ikl,isn).eq.0)
 
 9240           dsnmas     = 100.*
dzsnsv(ikl,isn)*ro_dry(ikl,isn)
 
 9241           snmass(ikl)=      snmass(ikl)+0.5*dsnmas
 
 9242           viscsn     =      vvisc1         *vvisc2
 
 9243      .               *exp(vvisc3           *ro_dry(ikl,isn)
 
 9245      .                                     *ro_dry(ikl,isn)/rovisc
 
 9264           okxliq     =    max(
zer0,sign(
un_1,vtelv1-etasno(ikl,isn)
 
 9266      .               *    max(0   ,sign(1   ,
istosv(ikl,isn)
 
 9269      .    viscsn*(    ok_liq/(vvisc5+vvisc6*etasno(ikl,isn)
 
 9285      .                     /max(viscsn                      ,
eps6)))
 
 9290           rosnew         =                        min(rosnew ,rosmax)
 
 9294           ro_dry(ikl,isn)= 
ro__sv(ikl,isn)*(1.-
eta_sv(ikl,isn))*1.e-3
 
 9297           snmass(ikl)    = snmass(ikl)+dsnmas*0.5
 
 9307  6600     
format(/,
'WARNING in _GSn: G1,G2 =',2f9.3,
'  (ikl,isn) =',2i4)
 
 9312  6601     
format(/,
'ERROR 1 in _GSn: G1,G2 =',2f9.3,
'  (ikl,isn) =',2i4)
 
 9317  6602     
format(/,
'ERROR 2 in _GSn: G1,G2 =',2f9.3,
'  (ikl,isn) =',2i4)
 
 9323  6603     
format(/,
'ERROR 3 in _GSn: G1,G2 =',2f9.3,
'  (ikl,isn) =',2i4)
 
 9459       integer   isl   ,jsl   ,ist   ,ikl      
 
 9460       integer   ikm   ,ikp   ,ik0   ,ik1      
 
 9461       integer   ist__s,ist__w                 
 
 9479       real      Elem_A,Elem_B,Elem_C          
 
 9548           ist__s = min(ist, 1)                     
 
 9553           dhydif = ist__s    * dhydif              
 
 9556           khydav = ist__s    * 
ks_dsv(ist)         
 
 9628  6001     
format(i3,
'  vRain ,Wg_MAX ',2e12.3,
' mm/hr')
 
 9635  6002     
format(3
x,
'  WgFlow,WExces ',2e12.3,
' mm/hr')
 
 9642  6003     
format(3
x,
'  SoRnOF,drr_SV ',2e12.3,
' mm/hr')
 
 9652  6004     
format(3
x,
'  SatRat,drr_SV ',2e12.3,
' mm/hr')
 
 9663  6005     
format(3
x,
'  eta_SV,       ',e12.3,
' g/kg')
 
 9708           ist__s =      min(ist, 1)                       
 
 9717      .  *(etamid         **(   
bchdsv(ist)+2.))           
 
 9718           dhydtz(ikl,isl) =    dhydif*
dt__sv               
 9721           dhydtz(ikl,isl) =    dhydtz(ikl,isl) * ist__s   
 
 9728           dhydtz(ikl,isl) =    0.0                        
 
 9739           etaaux(ikl,isl) =  
eta_sv(ikl,isl)
 
 9749           elem_a   =         dhydtz(ikl,isl)   
 
 9751           elem_b   =      - (dhydtz(ikl,isl)
 
 9752      .                      +dhydtz(ikl,isl+1)
 
 9755           elem_c   =         dhydtz(ikl,isl+1) 
 
 9782           freedr   =         
iwafsv(ikl)       *  min(ist,1)
 
 9787           elem_b   =      - (dhydtz(ikl,isl+1) 
 
 9790           elem_c   =         dhydtz(ikl,isl+1) 
 
 9792           diag_a(ikl,isl) =  0.
 
 9813           elem_a   =         dhydtz(ikl,isl)   
 
 9815           elem_b   =      - (dhydtz(ikl,isl)
 
 9822           diag_c(ikl,isl) =  0.
 
 9843           aux__p(ikl,-
nsol) = diag_b(ikl,-
nsol)
 
 9844           aux__q(ikl,-
nsol) =-diag_c(ikl,-
nsol)/aux__p(ikl,-
nsol)
 
 9849           aux__p(ikl,isl)   = diag_a(ikl,isl)  *aux__q(ikl,isl-1)
 
 9851           aux__q(ikl,isl)   =-diag_c(ikl,isl)  /aux__p(ikl,isl)
 
 9861           eta_sv(ikl,isl)   =(term_d(ikl,isl)
 
 9862      .                       -diag_a(ikl,isl)  *
eta_sv(ikl,isl-1))
 
 9887           sornof(ikl)     =          sornof(ikl)   
 
 9888      .                    + max(
zer0,satrat)       
 
 9910  6011     
format(2i3,f8.1,i3,3e12.3)
 
 9922           freedr   =          
iwafsv(ikl)       *  min(ist,1)
 
 9924           sornof(ikl) =  sornof(ikl)
 
10060  42     
format(
'SubRoutine SISVAT_qSo: Local Water Budget',
 
10061      .       /,
'=========================================')
 
10067  420    
format(11(
'-'),
'----------+--------------+-',
 
10068      .          3(
'-'),
'----------+--------------+-',
 
10069      .                 
'----------------+----------------+',
 
10070      .       /,f8.2,3
x,
'Wats_0(1) |    Wats_d(1) | ',
 
10071      .              3
x,
'Wats_1(1) | W_0+W_d-W_1  | ',
 
10072      .                 
'   Soil Run OFF |   Soil Evapor. |',
 
10073      .       /,11(
'-'),
'----------+--------------+-',
 
10074      .          3(
'-'),
'----------+--------------+-',
 
10075      .                 
'----------------+----------------+')
 
10080  421    
format(8
x,f12.6,
' + ',f12.6,
' - ',f12.6,
' = ',f12.6,
' | ',f12.6,
 
10118       common/sisvat_weq_l/logweq
 
10131       IF (.NOT.logweq)                                              
THEN 
10133                open(
unit=45,status=
'unknown',file=
'SISVAT_wEq.ve')
 
10155       IF        (
iicesv(1).gt.0)                                    
THEN 
10168       IF (istart.eq.1)                                              
THEN 
10171  45     
format(a18,10(
'-'),
'Pt.',3i4,60(
'-'))
 
10174       write(45,450) labweq,iceweq,
iicesv(ikl),snoweq
 
10175      .                    ,iceweq+snoweq,
isnosv(ikl)
 
10179  450  
format(a6,3
x,
'  I+S =',f11.4,
'(',i2,
') +',f11.4,
' =',
 
10180      .                       f11.4,
'(',i2,
')',
 
real, dimension(:), allocatable, save fh__sv
 
real, dimension(0:nsot, 0:nkhy) akdtsv
 
real, dimension(:), allocatable, save brossv
 
real, dimension(:), allocatable, save sqrcm0
 
real, dimension(:), allocatable, save lx_h2o
 
real, dimension(:), allocatable, save bcohsv
 
real, dimension(:), allocatable, save tau_sv
 
real, dimension(:), allocatable, save vvmmem
 
real, dimension(:), allocatable, save tvegsv
 
real, dimension(:), allocatable, save dsdtsv
 
real, dimension(0:nvgt), parameter pr_dsv
 
real, dimension(:,:), allocatable, save dz1_sv
 
real, dimension(:), allocatable, save irv_sv
 
subroutine sisvat_weq(labWEq, istart)
 
real, dimension(0:nsot) rocssv
 
real, dimension(:), allocatable, save rrmxsv
 
real, dimension(:), allocatable, save socasv
 
real, dimension(:), allocatable, save ps__sv
 
real, dimension(:), allocatable, save p1l_sv
 
real, dimension(:,:), allocatable, save sex_sv
 
real, dimension(-nsol:0) dzavsv
 
real, dimension(-nsol:0) dz_8sv
 
integer, dimension(:), allocatable, save lindsv
 
real, dimension(:), allocatable, save exnrsv
 
real, dimension(:), allocatable, save bufssv
 
integer, dimension(:,:), allocatable, save istosv
 
real, dimension(:), allocatable, save sol_sv
 
real, dimension(:,:), allocatable, save agsnsv
 
real, dimension(:,:), allocatable, save dzsnsv
 
real, dimension(:,:), allocatable, save rootsv
 
real, dimension(:), allocatable, save rsolsv
 
real, dimension(:), allocatable, save alb_sv
 
real, dimension(:), allocatable, save lai0sv
 
real, dimension(:), allocatable, save evp_sv
 
!$Id!common comsoil inertie_sno
 
integer, dimension(:), allocatable, save isnosv
 
integer, dimension(:), allocatable, save isotsv
 
real, dimension(:), allocatable, save dsnbsv
 
real, dimension(:), allocatable, save bcoqsv
 
real, dimension(:), allocatable, save sncasv
 
real, dimension(:), allocatable, save sws_sv
 
real, dimension(:), allocatable, save eexcsv
 
integer, dimension(:), allocatable, save iicesv
 
real, dimension(:,:), allocatable, save g1snsv
 
real, dimension(:), allocatable, save evt_sv
 
real, dimension(:), allocatable, save dfh_sv
 
real, dimension(-nsol:0) dzmisv
 
real, dimension(:,:), allocatable, save tsissv
 
real, dimension(:), allocatable, save uts_sv
 
real, dimension(:), allocatable, save tat_sv
 
real, dimension(:), allocatable, save acoqsv
 
real, dimension(:), allocatable, save tsrfsv
 
real, dimension(:), allocatable, save iru_sv
 
real, dimension(:,:), allocatable, save dz2_sv
 
real, dimension(-nsol:0) dzi_sv
 
integer, dimension(1:5), parameter istdsv
 
real, dimension(:), allocatable, save zwe_sv
 
integer, dimension(nb_wri), save lwrisv
 
real, dimension(:), allocatable, save zwecsv
 
real, dimension(:), allocatable, save psivsv
 
real, dimension(:), allocatable, save cld_sv
 
real, dimension(0:nsot, 0:nkhy) bkdtsv
 
character(len=18), save dahost
 
real, dimension(:), allocatable, save hlv_sv
 
real, dimension(:), allocatable, save lsdzsv
 
real, dimension(:), allocatable, save alb0sv
 
real, dimension(:), allocatable, save rusnsv
 
real, dimension(:), allocatable, save lai_sv
 
real, dimension(:), allocatable, save acohsv
 
real, dimension(:), allocatable, save glf0sv
 
real, dimension(:), allocatable, save bg2ssv
 
real, dimension(:), allocatable, save z0ensv
 
real, dimension(-nsol:0) dziisv
 
subroutine sisvat(SnoMod, BloMod, jjtime)
 
integer, dimension(nb_wri), save j___sv
 
real, dimension(:), allocatable, save usthsv
 
real, dimension(:), allocatable, save z0h_sv
 
integer, dimension(:), allocatable, save iwafsv
 
integer, dimension(:), allocatable, save ivgtsv
 
integer, dimension(:), allocatable, save i_thin
 
real, dimension(0:nvgt,-nsol:0) rf__sv
 
real, dimension(:), allocatable, save rah_sv
 
real, dimension(-nsol:0) dz78sv
 
real, dimension(0:nsot), parameter ks_dsv
 
integer, dimension(nb_wri), save i___sv
 
real, dimension(:), allocatable, save hss_sv
 
real, dimension(:), allocatable, save za__sv
 
real, dimension(:), allocatable, save vv__sv
 
real, dimension(:), allocatable, save sososv
 
real, dimension(:), allocatable, save tsf_sv
 
integer, dimension(:), allocatable, save ispisv
 
real, dimension(:), allocatable, save z0hnsv
 
real, dimension(:), allocatable, save sqrch0
 
real, dimension(:), allocatable, save qsnosv
 
real, dimension(:), allocatable, save rht_sv
 
integer, dimension(nb_wri), save n___sv
 
real, dimension(:), allocatable, save dtmmem
 
real, dimension(:,:), allocatable, save psi_sv
 
real, dimension(:), allocatable, save laiesv
 
real, dimension(0:nvgt), parameter stodsv
 
real, dimension(0:nsot) s2__sv
 
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
 
real, dimension(:), allocatable, save rnofsv
 
subroutine sisvat_bsn(BloMod)
 
real, dimension(:), allocatable, save rcdmsv
 
real, dimension(:), allocatable, save hsv_sv
 
!$Header!c c INCLUDE fxyprim h c c c Fonctions in line c c REAL fyprim REAL rj c c il faut la calculer avant d appeler ces fonctions c c c Fonctions a changer selon x(x) et y(y) choisis.c-----------------------------------------------------------------c c.....ici
 
real, dimension(:), allocatable, save ird_sv
 
real, dimension(0:nsot), parameter etadsv
 
real, dimension(-nsol:0) dz34sv
 
real, dimension(:), allocatable, save dbs_sv
 
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
 
real, dimension(0:nsot), parameter psidsv
 
real, dimension(:), allocatable, save drr_sv
 
real, dimension(:,:), allocatable, save ro__sv
 
real, dimension(:), allocatable, save irs_sv
 
real, dimension(:), allocatable, save lmo_sv
 
subroutine snoptp(jjtime)
 
real, dimension(:), allocatable, save cdh_sv
 
real, dimension(:), allocatable, save z0m_sv
 
real, dimension(:), allocatable, save alb2sv
 
real, dimension(0:nvgt), parameter dh_dsv
 
real, dimension(:), allocatable, save dsn_sv
 
real, dimension(:), allocatable, save glf_sv
 
real, dimension(:), allocatable, save sigmsv
 
real, dimension(-nsol:0) dtz_sv
 
integer, dimension(:), allocatable, save lsmask
 
real, dimension(0:nsot) s1__sv
 
real, dimension(:), allocatable, save albisv
 
real, dimension(:), allocatable, save uqs_sv
 
real, dimension(:), allocatable, save evg_sv
 
real, dimension(:), allocatable, save us__sv
 
real, dimension(:), allocatable, save uss_sv
 
real, dimension(:), allocatable, save slopsv
 
real, dimension(:), allocatable, save rcdhsv
 
real, dimension(:), allocatable, save eso_sv
 
real, dimension(0:nsot), parameter bchdsv
 
real, dimension(-nsol:0) dz_dsv
 
real, dimension(:,:), allocatable, save khydsv
 
real, dimension(:,:), allocatable, save eta_sv
 
real, dimension(:), allocatable, save z0mnsv
 
real, dimension(:), allocatable, save qat_sv
 
real, dimension(:,:), allocatable, save zzsnsv
 
real, dimension(:), allocatable, save tsfnsv
 
real, dimension(:), allocatable, save rrcasv
 
real, dimension(:), allocatable, save alb1sv
 
real, dimension(:), allocatable, save emi_sv
 
subroutine sisvat_zag(isagra, isagrb, WEagra, dzagra, dzagrb, T_agra, T_agrb, roagra, roagrb, etagra, etagrb, G1agra, G1agrb, G2agra, G2agrb, agagra, agagrb, Agreg1)
 
real, dimension(:), allocatable, save alb3sv
 
real, dimension(:), allocatable, save coszsv
 
!$Id!Thermodynamical constants for t0 real clmci real epsi
 
real, dimension(:), allocatable, save albssv
 
!$Header!integer nvarmx s s unit
 
real, dimension(:), allocatable, save hls_sv
 
real, dimension(:), allocatable, save ram_sv
 
real, dimension(:), allocatable, save dldtsv
 
real, dimension(:), allocatable, save swf_sv
 
real, dimension(:,:), allocatable, save g2snsv
 
real, dimension(:), allocatable, save z0e_sv
 
real, dimension(:), allocatable, save esnbsv
 
real, dimension(0:nvgt), parameter z0mdsv
 
real, dimension(:), allocatable, save bg1ssv
 
integer, dimension(:), allocatable, save nlaysv