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