MODULE module_sf_ssib
REAL, PARAMETER :: CPAIR = 1004.6 &
,STEFAN = 5.669 * 10E-9 &
,GRAV = 9.81 &
,VKC = 0.4 &
,PIE = 3.14159265 &
,TIMCON = PIE/86400. &
,CLAI = 4.2 * 1000. * 0.2 &
,CW = 4.2 * 1000. * 1000. &
,TF = 273.16 &
,GASR = 287.05 &
,HLAT = 2.52E6 &
,SNOMEL = 370518.5 * 1000.
INTEGER, PARAMETER :: ITRUNK = 3
REAL, PARAMETER :: SSISNOW = 0.04 &
,FLMIN = 0.03 &
,FLMAX = 0.10 &
,DZMIN = 0.002 &
,WOMIN = 0.0004 &
,CL = 4212.7 &
,DLM = 3.335d5 &
,RHOWATER = 1000.0 &
,DICE = 920.0 &
,DKSATSNOW= 0.01 &
,SNODEP_CR= 0.07
INTEGER, PARAMETER :: N = 3 &
,N1 = 4 &
,N2 = 4
REAL, DIMENSION (13,2,3,2) :: tran0,ref0
REAL, DIMENSION (13,12,2) :: green0,vcover0,zlt0
REAL, DIMENSION (13,2,3) :: rstpar0
REAL, DIMENSION (13,12) :: z000,d0,z20,z10,rdc0,rbc0
REAL, DIMENSION (13,3) :: depth0,soref0
REAL, DIMENSION (13,2) :: chil0,topt0,tl0,tu0,defac0,ph10,ph20,rootd0
REAL, DIMENSION (13) :: bee0,phsat0,poros0,satco0,slope0
data tran0/ &
0.5000000E-01, 0.5000000E-01, 0.5000000E-01, 0.5000000E-01, &
0.5000000E-01, 0.5000000E-01, 0.7000000E-01, 0.5000000E-01, &
0.5000000E-01, 0.5000000E-01, 0.1000000E-02, 0.5000000E-01, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.7000000E-01, 0.1000000E-02, 0.7000000E-01, &
0.1000000E-02, 0.7000000E-01, 0.1000000E-02, 0.7000000E-01, &
0.1000000E-02, &
0.2500000E+00, 0.2500000E+00, 0.1500000E+00, 0.1000000E+00, &
0.1000000E+00, 0.2500000E+00, 0.2475000E+00, 0.2500000E+00, &
0.2500000E+00, 0.2500000E+00, 0.1000000E-02, 0.2500000E+00, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.2475000E+00, 0.1000000E-02, 0.2475000E+00, &
0.1000000E-02, 0.2475000E+00, 0.1000000E-02, 0.2475000E+00, &
0.1000000E-02, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.2200000E+00, 0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.2200000E+00, 0.1000000E-02, 0.2200000E+00, &
0.1000000E-02, 0.2200000E+00, 0.1000000E-02, 0.2200000E+00, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.3750000E+00, 0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.3750000E+00, 0.1000000E-02, 0.3750000E+00, &
0.1000000E-02, 0.3750000E+00, 0.1000000E-02, 0.3750000E+00, &
0.1000000E-02, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00/
data ref0/ &
0.1000000E+00, 0.1000000E+00, 0.7000000E-01, 0.7000000E-01, &
0.7000000E-01, 0.1000000E+00, 0.1050000E+00, 0.1000000E+00, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-02, 0.1000000E+00, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.1050000E+00, 0.1000000E-02, 0.1050000E+00, &
0.1000000E-02, 0.1050000E+00, 0.1000000E-02, 0.1050000E+00, &
0.1000000E-02, &
0.4500000E+00, 0.4500000E+00, 0.4000000E+00, 0.3500000E+00, &
0.3500000E+00, 0.4500000E+00, 0.5775000E+00, 0.4500000E+00, &
0.4500000E+00, 0.4500000E+00, 0.1000000E-02, 0.4500000E+00, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.5775000E+00, 0.1000000E-02, 0.5775000E+00, &
0.1000000E-02, 0.5775000E+00, 0.1000000E-02, 0.5775000E+00, &
0.1000000E-02, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, &
0.1600000E+00, 0.1600000E+00, 0.1600000E+00, 0.1600000E+00, &
0.1600000E+00, 0.1600000E+00, 0.3600000E+00, 0.1600000E+00, &
0.1600000E+00, 0.1600000E+00, 0.1000000E-02, 0.1600000E+00, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.3600000E+00, 0.1000000E-02, 0.3600000E+00, &
0.1000000E-02, 0.3600000E+00, 0.1000000E-02, 0.3600000E+00, &
0.1000000E-02, &
0.3900000E+00, 0.3900000E+00, 0.3900000E+00, 0.3900000E+00, &
0.3900000E+00, 0.3900000E+00, 0.5775000E+00, 0.3900000E+00, &
0.3900000E+00, 0.3900000E+00, 0.1000000E-02, 0.3900000E+00, &
0.1000000E-02, &
0.1000000E-02, 0.1000000E-02, 0.1000000E-02, 0.1000000E-02, &
0.1000000E-02, 0.5775000E+00, 0.1000000E-02, 0.5775000E+00, &
0.1000000E-02, 0.5775000E+00, 0.1000000E-02, 0.5775000E+00, &
0.1000000E-02, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00/
data green0/ &
0.9050000E+00, 0.2564000E-01, 0.8680600E+00, 0.9132400E+00, &
0.2475200E+00, 0.6319100E+00, 0.5681800E+00, 0.7978700E+00, &
0.8364300E+00, 0.4512600E+00, 0.1000000E-03, 0.2083300E+00, &
0.1000000E-03, &
0.9050000E+00, 0.2564000E-01, 0.8717700E+00, 0.9170300E+00, &
0.2475200E+00, 0.6566600E+00, 0.6218900E+00, 0.5319100E+00, &
0.7172100E+00, 0.4512600E+00, 0.1000000E-03, 0.2083300E+00, &
0.1000000E-03, &
0.9050000E+00, 0.4153800E+00, 0.8847300E+00, 0.9226600E+00, &
0.2475200E+00, 0.5176000E+00, 0.6637200E+00, 0.3623200E+00, &
0.2577300E+00, 0.4512600E+00, 0.1000000E-03, 0.4411800E+00, &
0.1000000E-03, &
0.9050000E+00, 0.7594900E+00, 0.9061000E+00, 0.9247000E+00, &
0.6637200E+00, 0.6527400E+00, 0.6972100E+00, 0.5681800E+00, &
0.7246400E+00, 0.4512600E+00, 0.1000000E-03, 0.7594900E+00, &
0.1000000E-03, &
0.9050000E+00, 0.8875700E+00, 0.9164200E+00, 0.9266400E+00, &
0.8104700E+00, 0.6527400E+00, 0.8104700E+00, 0.5681800E+00, &
0.1736100E+00, 0.4512600E+00, 0.1000000E-03, 0.8875700E+00, &
0.1000000E-03, &
0.9050000E+00, 0.9252000E+00, 0.9259300E+00, 0.9045800E+00, &
0.8680600E+00, 0.7246400E+00, 0.9079900E+00, 0.5681800E+00, &
0.5681800E+00, 0.6218900E+00, 0.1000000E-03, 0.9252000E+00, &
0.1000000E-03, &
0.9050000E+00, 0.8364300E+00, 0.9293700E+00, 0.9021600E+00, &
0.6040900E+00, 0.8712500E+00, 0.8132000E+00, 0.5681800E+00, &
0.5681800E+00, 0.9200800E+00, 0.1000000E-03, 0.8364300E+00, &
0.1000000E-03, &
0.9050000E+00, 0.6967200E+00, 0.8209400E+00, 0.9126500E+00, &
0.5854000E+00, 0.7966000E+00, 0.3943200E+00, 0.8680600E+00, &
0.7246400E+00, 0.6970300E+00, 0.1000000E-03, 0.6967200E+00, &
0.1000000E-03, &
0.9050000E+00, 0.3306900E+00, 0.7123000E+00, 0.8982800E+00, &
0.4990000E+00, 0.7654600E+00, 0.4434600E+00, 0.6505600E+00, &
0.8403400E+00, 0.7567000E-01, 0.1000000E-03, 0.3439200E+00, &
0.1000000E-03, &
0.9050000E+00, 0.1656400E+00, 0.6145700E+00, 0.8548200E+00, &
0.3834400E+00, 0.6146100E+00, 0.5434800E+00, 0.5154600E+00, &
0.8680600E+00, 0.4512600E+00, 0.1000000E-03, 0.1785700E+00, &
0.1000000E-03, &
0.9050000E+00, 0.1538000E-01, 0.8599500E+00, 0.8733600E+00, &
0.2487600E+00, 0.5086500E+00, 0.5531000E+00, 0.6302500E+00, &
0.8875700E+00, 0.4512600E+00, 0.1000000E-03, 0.1470600E+00, &
0.1000000E-03, &
0.9050000E+00, 0.2564000E-01, 0.8599500E+00, 0.9132400E+00, &
0.1984100E+00, 0.7898900E+00, 0.4975100E+00, 0.7978700E+00, &
0.9132400E+00, 0.4512600E+00, 0.1000000E-03, 0.2083300E+00, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03/
data vcover0/ &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.9800000E+00, 0.7500000E+00, 0.7500000E+00, 0.7500000E+00, &
0.5000000E+00, 0.3000000E+00, 0.9000000E+00, 0.1000000E+00, &
0.1000000E+00, 0.3000000E+00, 0.1000000E-04, 0.7500000E-01, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-04, 0.1000000E-03, &
0.1000000E-04/
data chil0/ &
0.1000000E+00, 0.2500000E+00, 0.1300000E+00, 0.1000000E-01, &
0.1000000E-01, 0.1000000E-01, -0.3000000E+00, 0.1000000E-01, &
0.1000000E-01, 0.2000000E+00, 0.1000000E-01, -0.2000000E-01, &
0.1000000E-01, &
0.1000000E+00, 0.2500000E+00, 0.1300000E+00, 0.1000000E-01, &
0.1000000E-01, -0.3000000E+00, -0.3000000E+00, -0.3000000E+00, &
0.1000000E-01, 0.2000000E+00, 0.1000000E-01, -0.2000000E-01, &
0.1000000E-01/
data rstpar0/ &
0.2335900E+04, 0.9802230E+04, 0.6335955E+04, 0.2869680E+04, &
0.2869680E+04, 0.5665121E+05, 0.2582010E+04, 0.9398942E+05, &
0.9398942E+05, 0.9802230E+04, 0.1000000E+04, 0.7459000E+04, &
0.1000000E+04, &
0.2335900E+04, 0.9802230E+04, 0.6335955E+04, 0.2869680E+04, &
0.2869680E+04, 0.2582010E+04, 0.2582010E+04, 0.2582010E+04, &
0.1000000E+01, 0.2582010E+04, 0.1000000E+04, 0.7459000E+04, &
0.1000000E+04, &
0.1450000E-01, 0.1055000E+02, 0.7120000E+01, 0.3690000E+01, &
0.3690000E+01, 0.1083000E+02, 0.1090000E+01, 0.1000000E-01, &
0.1000000E-01, 0.1055000E+02, 0.1000000E+04, 0.5700000E+01, &
0.1000000E+04, &
0.1450000E-01, 0.1055000E+02, 0.7120000E+01, 0.3690000E+01, &
0.3690000E+01, 0.1090000E+01, 0.1090000E+01, 0.1090000E+01, &
0.1000000E+01, 0.1090000E+01, 0.1000000E+04, 0.5700000E+01, &
0.1000000E+04, &
0.1534900E+03, 0.1800000E+03, 0.2065000E+03, 0.2330000E+03, &
0.2330000E+03, 0.1650000E+03, 0.1100000E+03, 0.8550000E+03, &
0.8550000E+03, 0.1800000E+03, 0.1000000E+04, 0.2520000E+02, &
0.1000000E+04, &
0.1534900E+03, 0.1800000E+03, 0.2065000E+03, 0.2330000E+03, &
0.2330000E+03, 0.1100000E+03, 0.1100000E+03, 0.1100000E+03, &
0.1000000E+01, 0.1100000E+03, 0.1000000E+04, 0.2520000E+02, &
0.1000000E+04/
data topt0/ &
0.3030000E+03, 0.3000000E+03, 0.2940000E+03, 0.2880000E+03, &
0.2880000E+03, 0.2970000E+03, 0.3130000E+03, 0.3150000E+03, &
0.3150000E+03, 0.3000000E+03, 0.3100000E+03, 0.3000000E+03, &
0.3100000E+03, &
0.3030000E+03, 0.3000000E+03, 0.2940000E+03, 0.2880000E+03, &
0.2880000E+03, 0.3120000E+03, 0.3130000E+03, 0.3130000E+03, &
0.3150000E+03, 0.2890000E+03, 0.3100000E+03, 0.3000000E+03, &
0.3100000E+03/
data tl0/ &
0.2730000E+03, 0.2730000E+03, 0.2700000E+03, 0.2680000E+03, &
0.2680000E+03, 0.2730000E+03, 0.2830000E+03, 0.2830000E+03, &
0.2830000E+03, 0.2730000E+03, 0.3000000E+03, 0.2730000E+03, &
0.3000000E+03, &
0.2730000E+03, 0.2730000E+03, 0.2700000E+03, 0.2680000E+03, &
0.2680000E+03, 0.2730000E+03, 0.2830000E+03, 0.2830000E+03, &
0.2830000E+03, 0.2730000E+03, 0.3000000E+03, 0.2730000E+03, &
0.3000000E+03/
data tu0/ &
0.3180000E+03, 0.3180000E+03, 0.3150000E+03, 0.3130000E+03, &
0.3130000E+03, 0.3230000E+03, 0.3280000E+03, 0.3230000E+03, &
0.3230000E+03, 0.3230000E+03, 0.3200000E+03, 0.3180000E+03, &
0.3200000E+03, &
0.3180000E+03, 0.3180000E+03, 0.3150000E+03, 0.3130000E+03, &
0.3130000E+03, 0.3230000E+03, 0.3280000E+03, 0.3280000E+03, &
0.3230000E+03, 0.3090000E+03, 0.3200000E+03, 0.3150000E+03, &
0.3200000E+03/
data defac0/ &
0.2730000E-01, 0.3570000E-01, 0.3400000E-01, 0.3100000E-01, &
0.3100000E-01, 0.3570000E-01, 0.2380000E-01, 0.2750000E-01, &
0.2750000E-01, 0.2750000E-01, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, &
0.2730000E-01, 0.3570000E-01, 0.3400000E-01, 0.3100000E-01, &
0.3100000E-01, 0.2380000E-01, 0.2380000E-01, 0.2380000E-01, &
0.2380000E-01, 0.2380000E-01, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00/
data ph10/ &
0.1200000E+01, 0.5350000E+01, 0.1920000E+01, 0.3700000E+01, &
0.7800000E+01, 0.1800000E+01, 0.1730000E+01, 0.1920000E+01, &
0.1390000E+01, 0.9600000E+00, 0.3000000E+01, 0.1800000E+01, &
0.5000000E+01, &
0.1200000E+01, 0.5350000E+01, 0.1920000E+01, 0.3700000E+01, &
0.7800000E+01, 0.1800000E+01, 0.1730000E+01, 0.1920000E+01, &
0.1390000E+01, 0.9600000E+00, 0.3000000E+01, 0.1800000E+01, &
0.5000000E+01/
data ph20/ &
0.6250000E+01, 0.5570000E+01, 0.5730000E+01, 0.5530000E+01, &
0.5660000E+01, 0.5670000E+01, 0.5800000E+01, 0.5610000E+01, &
0.6370000E+01, 0.5370000E+01, 0.6000000E+01, 0.5670000E+01, &
0.6000000E+01, &
0.6250000E+01, 0.5570000E+01, 0.5730000E+01, 0.5530000E+01, &
0.5660000E+01, 0.5670000E+01, 0.5800000E+01, 0.5610000E+01, &
0.6370000E+01, 0.5370000E+01, 0.6000000E+01, 0.5670000E+01, &
0.6000000E+01/
data zlt0/ &
0.5014160E+01, 0.3900000E+00, 0.3456000E+01, 0.6570000E+01, &
0.4040000E+00, 0.1766000E+01, 0.7040000E+00, 0.5780000E+00, &
0.1076000E+01, 0.3776000E+00, 0.1000000E-03, 0.4800000E-01, &
0.1000000E-03, &
0.5014160E+01, 0.3900000E+00, 0.3556000E+01, 0.6870000E+01, &
0.4040000E+00, 0.1546000E+01, 0.8040000E+00, 0.5780000E+00, &
0.9760000E+00, 0.3776000E+00, 0.1000000E-03, 0.4800000E-01, &
0.1000000E-03, &
0.5014160E+01, 0.6500000E+00, 0.3956000E+01, 0.7370000E+01, &
0.4040000E+00, 0.1416000E+01, 0.9040000E+00, 0.4480000E+00, &
0.7760000E+00, 0.3776000E+00, 0.1000000E-03, 0.6800000E-01, &
0.1000000E-03, &
0.5014160E+01, 0.1580000E+01, 0.4856000E+01, 0.7570000E+01, &
0.9040000E+00, 0.1216000E+01, 0.1004000E+01, 0.2880000E+00, &
0.2760000E+00, 0.3776000E+00, 0.1000000E-03, 0.1580000E+00, &
0.1000000E-03, &
0.5014160E+01, 0.3380000E+01, 0.5456000E+01, 0.7770000E+01, &
0.1604000E+01, 0.1186000E+01, 0.1604000E+01, 0.2580000E+00, &
0.5760000E+00, 0.3776000E+00, 0.1000000E-03, 0.3380000E+00, &
0.1000000E-03, &
0.5014160E+01, 0.5080000E+01, 0.6156000E+01, 0.8070000E+01, &
0.2304000E+01, 0.1416000E+01, 0.3304000E+01, 0.2580000E+00, &
0.1760000E+00, 0.5076000E+00, 0.1000000E-03, 0.5080000E+00, &
0.1000000E-03, &
0.5014160E+01, 0.5380000E+01, 0.6456000E+01, 0.7870000E+01, &
0.4304000E+01, 0.2606000E+01, 0.4304000E+01, 0.2580000E+00, &
0.1760000E+00, 0.1737600E+01, 0.1000000E-03, 0.5380000E+00, &
0.1000000E-03, &
0.5014160E+01, 0.4880000E+01, 0.6456000E+01, 0.7670000E+01, &
0.2904000E+01, 0.5206000E+01, 0.3804000E+01, 0.8080000E+00, &
0.2760000E+00, 0.1937600E+01, 0.1000000E-03, 0.4880000E+00, &
0.1000000E-03, &
0.5014160E+01, 0.3780000E+01, 0.5756000E+01, 0.7570000E+01, &
0.2004000E+01, 0.4556000E+01, 0.1804000E+01, 0.1508000E+01, &
0.4760000E+00, 0.1477600E+01, 0.1000000E-03, 0.3780000E+00, &
0.1000000E-03, &
0.5014160E+01, 0.1630000E+01, 0.4556000E+01, 0.7370000E+01, &
0.1304000E+01, 0.3816000E+01, 0.1104000E+01, 0.1148000E+01, &
0.5760000E+00, 0.3776000E+00, 0.1000000E-03, 0.1680000E+00, &
0.1000000E-03, &
0.5014160E+01, 0.6500000E+00, 0.3256000E+01, 0.6870000E+01, &
0.8040000E+00, 0.2806000E+01, 0.9040000E+00, 0.7480000E+00, &
0.6760000E+00, 0.3776000E+00, 0.1000000E-03, 0.6800000E-01, &
0.1000000E-03, &
0.5014160E+01, 0.3900000E+00, 0.3256000E+01, 0.6570000E+01, &
0.5040000E+00, 0.1866000E+01, 0.8040000E+00, 0.5780000E+00, &
0.8760000E+00, 0.3776000E+00, 0.1000000E-03, 0.4800000E-01, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03, 0.1000000E-03, 0.1000000E-03, 0.1000000E-03, &
0.1000000E-03/
data z000/ &
0.2652970E+01, 0.5201000E+00, 0.5706300E+00, 0.1112210E+01, &
0.6414000E+00, 0.8427100E+00, 0.7771000E-01, 0.2446700E+00, &
0.6559000E-01, 0.7524000E-01, 0.1118000E-01, 0.1448500E+00, &
0.1011000E-01, &
0.2652970E+01, 0.5201000E+00, 0.5696600E+00, 0.1102780E+01, &
0.6414000E+00, 0.8087800E+00, 0.7779000E-01, 0.2446700E+00, &
0.6549000E-01, 0.7524000E-01, 0.1118000E-01, 0.1448500E+00, &
0.1011000E-01, &
0.2652970E+01, 0.6664900E+00, 0.5656600E+00, 0.1087660E+01, &
0.6414000E+00, 0.7875000E+00, 0.7785000E-01, 0.2272100E+00, &
0.6521000E-01, 0.7524000E-01, 0.1118000E-01, 0.1752100E+00, &
0.1011000E-01, &
0.2652970E+01, 0.9105300E+00, 0.5654400E+00, 0.1081830E+01, &
0.8633500E+00, 0.7284100E+00, 0.7788000E-01, 0.1998800E+00, &
0.6360000E-01, 0.7524000E-01, 0.1118000E-01, 0.2871900E+00, &
0.1011000E-01, &
0.2652970E+01, 0.1031200E+01, 0.5592300E+00, 0.1076120E+01, &
0.9728300E+00, 0.7284100E+00, 0.7779000E-01, 0.1998800E+00, &
0.6480000E-01, 0.7524000E-01, 0.1118000E-01, 0.4302000E+00, &
0.1011000E-01, &
0.2652970E+01, 0.1043680E+01, 0.5524400E+00, 0.1067790E+01, &
0.1005600E+01, 0.7875000E+00, 0.7712000E-01, 0.1998800E+00, &
0.6331000E-01, 0.7575000E-01, 0.1118000E-01, 0.5087600E+00, &
0.1011000E-01, &
0.2652970E+01, 0.1041940E+01, 0.5497000E+00, 0.1073310E+01, &
0.9967700E+00, 0.9266800E+00, 0.7594000E-01, 0.1998800E+00, &
0.6331000E-01, 0.7767000E-01, 0.1118000E-01, 0.5200300E+00, &
0.1011000E-01, &
0.2652970E+01, 0.1037530E+01, 0.5497000E+00, 0.1078960E+01, &
0.1011190E+01, 0.9715300E+00, 0.7658000E-01, 0.2674000E+00, &
0.6360000E-01, 0.7782000E-01, 0.1118000E-01, 0.5009500E+00, &
0.1011000E-01, &
0.2652970E+01, 0.1036510E+01, 0.5562600E+00, 0.1081830E+01, &
0.9965000E+00, 0.9658800E+00, 0.7776000E-01, 0.2923300E+00, &
0.6446000E-01, 0.7745000E-01, 0.1118000E-01, 0.4503800E+00, &
0.1011000E-01, &
0.2652970E+01, 0.9170700E+00, 0.5686600E+00, 0.1087660E+01, &
0.9386100E+00, 0.9555100E+00, 0.7790000E-01, 0.2803400E+00, &
0.6480000E-01, 0.7524000E-01, 0.1118000E-01, 0.2973700E+00, &
0.1011000E-01, &
0.2652970E+01, 0.6664900E+00, 0.5725100E+00, 0.1102780E+01, &
0.8346400E+00, 0.9204000E+00, 0.7785000E-01, 0.2580600E+00, &
0.6510000E-01, 0.7524000E-01, 0.1118000E-01, 0.1752100E+00, &
0.1011000E-01, &
0.2652970E+01, 0.5201000E+00, 0.5725100E+00, 0.1112210E+01, &
0.7049800E+00, 0.8427100E+00, 0.7779000E-01, 0.2446700E+00, &
0.6537000E-01, 0.7524000E-01, 0.1118000E-01, 0.1448500E+00, &
0.1011000E-01/
data d0/ &
0.2737261E+02, 0.1366377E+02, 0.1813464E+02, 0.1376361E+02, &
0.9193320E+01, 0.1390777E+02, 0.2185200E+00, 0.2812600E+01, &
0.1638000E+00, 0.1062900E+00, 0.6000000E-04, 0.6314240E+01, &
0.4000000E-04, &
0.2737261E+02, 0.1366377E+02, 0.1814677E+02, 0.1380041E+02, &
0.9193320E+01, 0.1376090E+02, 0.2265800E+00, 0.2812600E+01, &
0.1548100E+00, 0.1062900E+00, 0.6000000E-04, 0.6314240E+01, &
0.4000000E-04, &
0.2737261E+02, 0.1461883E+02, 0.1819051E+02, 0.1385740E+02, &
0.9193320E+01, 0.1367074E+02, 0.2332800E+00, 0.2662290E+01, &
0.1343400E+00, 0.1062900E+00, 0.6000000E-04, 0.7639520E+01, &
0.4000000E-04, &
0.2737261E+02, 0.1569677E+02, 0.1825890E+02, 0.1387880E+02, &
0.9903400E+01, 0.1344527E+02, 0.2389500E+00, 0.2390910E+01, &
0.6191000E-01, 0.1062900E+00, 0.6000000E-04, 0.1070958E+02, &
0.4000000E-04, &
0.2737261E+02, 0.1632865E+02, 0.1829956E+02, 0.1389946E+02, &
0.1030010E+02, 0.1344527E+02, 0.2605400E+00, 0.2390910E+01, &
0.1096800E+00, 0.1062900E+00, 0.6000000E-04, 0.1278272E+02, &
0.4000000E-04, &
0.2737261E+02, 0.1662263E+02, 0.1833903E+02, 0.1392915E+02, &
0.1053455E+02, 0.1367074E+02, 0.2988000E+00, 0.2390910E+01, &
0.5103000E-01, 0.1229900E+00, 0.6000000E-04, 0.1356813E+02, &
0.4000000E-04, &
0.2737261E+02, 0.1666297E+02, 0.1835387E+02, 0.1390953E+02, &
0.1091967E+02, 0.1425275E+02, 0.3251800E+00, 0.2390910E+01, &
0.5103000E-01, 0.2152100E+00, 0.6000000E-04, 0.1366182E+02, &
0.4000000E-04, &
0.2737261E+02, 0.1660123E+02, 0.1835387E+02, 0.1388922E+02, &
0.1068047E+02, 0.1459719E+02, 0.3130700E+00, 0.2974600E+01, &
0.6191000E-01, 0.2289700E+00, 0.6000000E-04, 0.1349985E+02, &
0.4000000E-04, &
0.2737261E+02, 0.1641343E+02, 0.1831739E+02, 0.1387880E+02, &
0.1044517E+02, 0.1452246E+02, 0.2649800E+00, 0.3137710E+01, &
0.9547000E-01, 0.1996100E+00, 0.6000000E-04, 0.1301951E+02, &
0.4000000E-04, &
0.2737261E+02, 0.1572679E+02, 0.1823553E+02, 0.1385740E+02, &
0.1016423E+02, 0.1443002E+02, 0.2438100E+00, 0.3062460E+01, &
0.1096800E+00, 0.1062900E+00, 0.6000000E-04, 0.1090759E+02, &
0.4000000E-04, &
0.2737261E+02, 0.1461883E+02, 0.1810866E+02, 0.1380041E+02, &
0.9814290E+01, 0.1422050E+02, 0.2332800E+00, 0.2907360E+01, &
0.1225000E+00, 0.1062900E+00, 0.6000000E-04, 0.7639520E+01, &
0.4000000E-04, &
0.2737261E+02, 0.1366377E+02, 0.1810866E+02, 0.1376361E+02, &
0.9417390E+01, 0.1390777E+02, 0.2265800E+00, 0.2812600E+01, &
0.1450200E+00, 0.1062900E+00, 0.6000000E-04, 0.6314240E+01, &
0.4000000E-04/
data z10/ &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03, &
0.1000000E+01, 0.1150000E+02, 0.1600000E+02, 0.8500000E+01, &
0.7000000E+01, 0.1000000E+02, 0.1000000E+00, 0.2000000E+01, &
0.1000000E+00, 0.1000000E+00, 0.1000000E-01, 0.1150000E+02, &
0.1000000E-03/
data z20/ &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00, &
0.3500000E+02, 0.2000000E+02, 0.2000000E+02, 0.1700000E+02, &
0.1400000E+02, 0.1800000E+02, 0.6000000E+00, 0.5000000E+01, &
0.5000000E+00, 0.6000000E+00, 0.1000000E+00, 0.2000000E+02, &
0.1000000E+00/
data rdc0/ &
0.2858700E+03, 0.2113200E+03, 0.2985200E+03, 0.5654100E+03, &
0.1852000E+03, 0.2301300E+03, 0.2443000E+02, 0.1036000E+03, &
0.2311000E+02, 0.2286000E+02, 0.2376000E+02, 0.1949000E+03, &
0.2850000E+02, &
0.2858700E+03, 0.2113200E+03, 0.3013500E+03, 0.5870500E+03, &
0.1852000E+03, 0.2244200E+03, 0.2463000E+02, 0.1036000E+03, &
0.2294000E+02, 0.2286000E+02, 0.2376000E+02, 0.1949000E+03, &
0.2850000E+02, &
0.2858700E+03, 0.2187800E+03, 0.3124600E+03, 0.6234600E+03, &
0.1852000E+03, 0.2215700E+03, 0.2480000E+02, 0.1023500E+03, &
0.2262000E+02, 0.2286000E+02, 0.2376000E+02, 0.1964400E+03, &
0.2850000E+02, &
0.2858700E+03, 0.2434000E+03, 0.3312300E+03, 0.6381300E+03, &
0.2048700E+03, 0.2164100E+03, 0.2496000E+02, 0.1007200E+03, &
0.2189000E+02, 0.2286000E+02, 0.2376000E+02, 0.2014400E+03, &
0.2850000E+02, &
0.2858700E+03, 0.2948700E+03, 0.3458300E+03, 0.6528600E+03, &
0.2330100E+03, 0.2164100E+03, 0.2572000E+02, 0.1007200E+03, &
0.2230000E+02, 0.2286000E+02, 0.2376000E+02, 0.2071300E+03, &
0.2850000E+02, &
0.2858700E+03, 0.3459000E+03, 0.3619400E+03, 0.6750500E+03, &
0.2620800E+03, 0.2215700E+03, 0.2774000E+02, 0.1007200E+03, &
0.2182000E+02, 0.2301000E+02, 0.2376000E+02, 0.2107900E+03, &
0.2850000E+02, &
0.2858700E+03, 0.3551800E+03, 0.3685400E+03, 0.6602400E+03, &
0.3443100E+03, 0.2500700E+03, 0.3006000E+02, 0.1007200E+03, &
0.2182000E+02, 0.2436000E+02, 0.2376000E+02, 0.2113100E+03, &
0.2850000E+02, &
0.2858700E+03, 0.3418400E+03, 0.3685400E+03, 0.6454900E+03, &
0.2870900E+03, 0.2885700E+03, 0.2886000E+02, 0.1053000E+03, &
0.2189000E+02, 0.2469000E+02, 0.2376000E+02, 0.2104200E+03, &
0.2850000E+02, &
0.2858700E+03, 0.3072200E+03, 0.3528500E+03, 0.6381300E+03, &
0.2495800E+03, 0.2780300E+03, 0.2590000E+02, 0.1079400E+03, &
0.2216000E+02, 0.2404000E+02, 0.2376000E+02, 0.2081500E+03, &
0.2850000E+02, &
0.2858700E+03, 0.2448400E+03, 0.3236500E+03, 0.6231300E+03, &
0.2211200E+03, 0.2668400E+03, 0.2511000E+02, 0.1065900E+03, &
0.2230000E+02, 0.2286000E+02, 0.2376000E+02, 0.2018800E+03, &
0.2850000E+02, &
0.2858700E+03, 0.2187800E+03, 0.2927900E+03, 0.5870500E+03, &
0.2008900E+03, 0.2475700E+03, 0.2480000E+02, 0.1044900E+03, &
0.2244000E+02, 0.2286000E+02, 0.2376000E+02, 0.1964400E+03, &
0.2850000E+02, &
0.2858700E+03, 0.2113200E+03, 0.2927900E+03, 0.5654100E+03, &
0.1892600E+03, 0.2301300E+03, 0.2464000E+02, 0.1036000E+03, &
0.2277000E+02, 0.2286000E+02, 0.2376000E+02, 0.1949000E+03, &
0.2850000E+02/
data rbc0/ &
0.5430000E+01, 0.6936000E+02, 0.8590000E+01, 0.8800000E+00, &
0.7850000E+01, 0.2661000E+02, 0.2207000E+02, 0.2188000E+02, &
0.1761000E+02, 0.4351000E+02, 0.3592951E+05, 0.5600000E+03, &
0.3546177E+05, &
0.5430000E+01, 0.6936000E+02, 0.8450000E+01, 0.8600000E+00, &
0.7850000E+01, 0.3044000E+02, 0.2053000E+02, 0.2188000E+02, &
0.1942000E+02, 0.4351000E+02, 0.3592951E+05, 0.5600000E+03, &
0.3546177E+05, &
0.5430000E+01, 0.4257000E+02, 0.7980000E+01, 0.8400000E+00, &
0.7850000E+01, 0.3295000E+02, 0.1934000E+02, 0.2673000E+02, &
0.2446000E+02, 0.4351000E+02, 0.3592951E+05, 0.4019700E+03, &
0.3546177E+05, &
0.5430000E+01, 0.1897000E+02, 0.7180000E+01, 0.8300000E+00, &
0.3810000E+01, 0.4003000E+02, 0.1838000E+02, 0.3712000E+02, &
0.6928000E+02, 0.4351000E+02, 0.3592951E+05, 0.1855200E+03, &
0.3546177E+05, &
0.5430000E+01, 0.1035000E+02, 0.6810000E+01, 0.8200000E+00, &
0.2400000E+01, 0.4003000E+02, 0.1516000E+02, 0.3712000E+02, &
0.3303000E+02, 0.4351000E+02, 0.3592951E+05, 0.9801000E+02, &
0.3546177E+05, &
0.5430000E+01, 0.7880000E+01, 0.6480000E+01, 0.8100000E+00, &
0.1860000E+01, 0.3295000E+02, 0.1068000E+02, 0.3712000E+02, &
0.8702000E+02, 0.3568000E+02, 0.3592951E+05, 0.7224000E+02, &
0.3546177E+05, &
0.5430000E+01, 0.7610000E+01, 0.6360000E+01, 0.8200000E+00, &
0.1290000E+01, 0.1870000E+02, 0.8300000E+01, 0.3712000E+02, &
0.8702000E+02, 0.1449000E+02, 0.3592951E+05, 0.6938000E+02, &
0.3546177E+05, &
0.5430000E+01, 0.8090000E+01, 0.6360000E+01, 0.8300000E+00, &
0.1600000E+01, 0.1318000E+02, 0.9330000E+01, 0.1722000E+02, &
0.6928000E+02, 0.1281000E+02, 0.3592951E+05, 0.7434000E+02, &
0.3546177E+05, &
0.5430000E+01, 0.9570000E+01, 0.6660000E+01, 0.8300000E+00, &
0.2040000E+01, 0.1420000E+02, 0.1457000E+02, 0.1317000E+02, &
0.4003000E+02, 0.1669000E+02, 0.3592951E+05, 0.8988000E+02, &
0.3546177E+05, &
0.5430000E+01, 0.1847000E+02, 0.7400000E+01, 0.8400000E+00, &
0.2820000E+01, 0.1559000E+02, 0.1760000E+02, 0.1497000E+02, &
0.3303000E+02, 0.4351000E+02, 0.3592951E+05, 0.1757600E+03, &
0.3546177E+05, &
0.5430000E+01, 0.4257000E+02, 0.8880000E+01, 0.8600000E+00, &
0.4210000E+01, 0.1933000E+02, 0.1934000E+02, 0.1906000E+02, &
0.2810000E+02, 0.4351000E+02, 0.3592951E+05, 0.4019700E+03, &
0.3546177E+05, &
0.5430000E+01, 0.6936000E+02, 0.8880000E+01, 0.8800000E+00, &
0.6400000E+01, 0.2661000E+02, 0.2053000E+02, 0.2188000E+02, &
0.2165000E+02, 0.4351000E+02, 0.3592951E+05, 0.5600000E+03, &
0.3546177E+05/
data rootd0/ &
0.1000000E+01, 0.1000000E+01, 0.1000000E+01, 0.5000000E+00, &
0.5000000E+00, 0.5000000E+00, 0.5000000E+00, 0.5000000E+00, &
0.5000000E+00, 0.2000000E+00, 0.1000000E+00, 0.1000000E+01, &
0.1000000E+01, &
0.1000000E+01, 0.1000000E+01, 0.1000000E+01, 0.5000000E+00, &
0.5000000E+00, 0.5000000E+00, 0.5000000E+00, 0.5000000E+00, &
0.5000000E+00, 0.2000000E+00, 0.1000000E+00, 0.1000000E+01, &
0.1000000E+01/
data soref0/ &
0.1100000E+00, 0.1100000E+00, 0.1100000E+00, 0.1100000E+00, &
0.1100000E+00, 0.1100000E+00, 0.1000000E+00, 0.1000000E+00, &
0.3000000E+00, 0.1000000E+00, 0.3000000E+00, 0.1000000E+00, &
0.1000000E+00, &
0.2250000E+00, 0.2250000E+00, 0.2250000E+00, 0.2250000E+00, &
0.2250000E+00, 0.2250000E+00, 0.2000000E+00, 0.2000000E+00, &
0.3500000E+00, 0.2000000E+00, 0.3500000E+00, 0.1500000E+00, &
0.1500000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
0.0000000E+00/
data bee0/ &
0.7120000E+01, 0.7120000E+01, 0.7120000E+01, 0.7120000E+01, &
0.7120000E+01, 0.7120000E+01, 0.7120000E+01, 0.4050000E+01, &
0.4050000E+01, 0.7120000E+01, 0.4050000E+01, 0.7797000E+01, &
0.4804000E+01/
data phsat0/ &
-0.8600000E-01, -0.8600000E-01, -0.8600000E-01, -0.8600000E-01, &
-0.8600000E-01, -0.8600000E-01, -0.8600000E-01, -0.3500000E-01, &
-0.3500000E-01, -0.8600000E-01, -0.3500000E-01, -0.1980000E+00, &
-0.1670000E+00/
data poros0/ &
0.4200000E+00, 0.4200000E+00, 0.4200000E+00, 0.4200000E+00, &
0.4200000E+00, 0.4200000E+00, 0.4200000E+00, 0.4352000E+00, &
0.4352000E+00, 0.4200000E+00, 0.4352000E+00, 0.4577000E+00, &
0.4352000E+00/
data satco0/ &
0.2000000E-04, 0.2000000E-04, 0.2000000E-04, 0.2000000E-04, &
0.2000000E-04, 0.2000000E-04, 0.2000000E-04, 0.1760000E-03, &
0.1760000E-03, 0.2000000E-04, 0.1760000E-03, 0.3500000E-05, &
0.7620000E-04/
data slope0/ &
0.1736000E+00, 0.1736000E+00, 0.1736000E+00, 0.1736000E+00, &
0.1736000E+00, 0.1736000E+00, 0.1736000E+00, 0.8720000E-01, &
0.8720000E-01, 0.1736000E+00, 0.8720000E-01, 0.3420000E+00, &
0.8720000E-01/
data depth0/ &
0.2000000E-01, 0.2000000E-01, 0.2000000E-01, 0.2000000E-01, &
0.2000000E-01, 0.2000000E-01, 0.2000000E-01, 0.2000000E-01, &
0.2000000E-01, 0.2000000E-01, 0.2000000E-01, 0.2000000E-01, &
0.1000000E+01, &
0.1480000E+01, 0.1480000E+01, 0.1480000E+01, 0.1480000E+01, &
0.1480000E+01, 0.1480000E+01, 0.4700000E+00, 0.4700000E+00, &
0.4700000E+00, 0.1700000E+00, 0.1700000E+00, 0.1480000E+01, &
0.1000000E+01, &
0.2000000E+01, 0.2000000E+01, 0.2000000E+01, 0.2000000E+01, &
0.2000000E+01, 0.2000000E+01, 0.1000000E+01, 0.1000000E+01, &
0.1000000E+01, 0.1000000E+01, 0.3000000E+00, 0.2000000E+01, &
0.1000000E+01/
CONTAINS
SUBROUTINE SSIB( II, JJ, DDTT, ITIME, ZLAT, SUNANGLE, &
PPL, PPC, RLWDOWN, ZWIND2, &
WWW1, WWW2, WWW3, &
TC, TGS, TD, &
SNOA, ROFF, &
UMM, VMM, QM, TM, &
PM, PSUR, ivgtyp, &
SWDOWN1, SNOB, &
SALB11, SALB12, SALB21, SALB22, &
RADFRAC11, RADFRAC12, RADFRAC21, RADFRAC22, &
XHSFLX, ELATEN, GHTFLX, XHLFLX, TGEFF, &
USTAR, RIB, FM, FH, CM, &
XLHF, XSHF, XGHF, XEGS, XECI, XECT, &
XEGI, XEGT, XSDN, XSUP, XLDN, XLUP, &
XWAT, XHCX, XHGX, XZLT, XVCF, XXZ0, &
XVEG, XDD, &
ISNOW,SWE,SNOWDEN,SNOWDEPTH,TKAIR, &
DZO1,WO1,TSSN1,TSSNO1,BWO1,BTO1,CTO1,FIO1,FLO1,BIO1,BLO1,HO1, &
DZO2,WO2,TSSN2,TSSNO2,BWO2,BTO2,CTO2,FIO2,FLO2,BIO2,BLO2,HO2, &
DZO3,WO3,TSSN3,TSSNO3,BWO3,BTO3,CTO3,FIO3,FLO3,BIO3,BLO3,HO3, &
DZO4,WO4,TSSN4,TSSNO4,BWO4,BTO4,CTO4,FIO4,FLO4,BIO4,BLO4,HO4, &
DAY, CLOUD, Q2M, TA, BEDO, &
sw_physics, MMINLU &
)
INTEGER, DIMENSION (12) :: IDAYS
REAL, DIMENSION (2) :: CAPAC, SATCAP, GREEN, VCOVER, ZLT, CHIL, TOPT, TL, &
TU, DEFAC, PH1, PH2, RST, ROOTD, RADT, PAR, PD
REAL, DIMENSION (3) :: WWW, SOREF, ZDEPTH, ROOTP, PHSOIL, YMATT, YMATQ
REAL, DIMENSION (2,2) :: RADFRAC, SALB
REAL, DIMENSION (2,3) :: RSTPAR
REAL, DIMENSION (2,4) :: RSTFAC
REAL, DIMENSION (3,2) :: RADN
REAL, DIMENSION (2,3,2) :: TRAN, REF, ALBEDO, EXTK
REAL, DIMENSION (2,2,2) :: RADFAC
INTEGER, DIMENSION (24) :: IVUSGS
REAL, DIMENSION (13) :: TD_DEPTH
INTEGER :: sw_physics
CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
REAL, DIMENSION (N2) :: SS,SSO,POROSITY,H,HO,BI,BIO,DZ,DZO,BW,BWO,BL
REAL, DIMENSION (N2) :: BLO,TSSN,TSSNO,W,WO,WF,FI,FIO, FL,FLO,DMLT
REAL, DIMENSION (N2) :: DMLTO,BT,BTO,S,SO,CT,CTO,DLIQVOL,DICEVOL
REAL, DIMENSION (N2) :: QK,PDZDTC,DMASS,DSOL,DHP,THK
DATA IDAYS/31,59,90,120,151,181,212,243,273,304,334,366/
DATA TD_DEPTH/1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.0, 1.0, 1.0 &
& , 0.5, 0.5, 1.5, 1.5/
DATA IVUSGS / 7, 12, 12, 12, 12, 12, 7, 9, &
8, 6, 2, 5, 1, 4, 3, 0, &
10, 3, 11, 10, 10, 10, 10, 13/
IF(MMINLU.EQ.'SSIB') THEN
ITYPE=IVGTYP
ELSEIF(MMINLU.EQ.'USGS') THEN
ITYPE=IVUSGS(IVGTYP)
ELSE
CALL wrf_error_fatal3("<stdin>",1021,&
'SSIB LSM only works with SSIB or USGS vegetation (landuse) map' )
ENDIF
if(itype.le.0.or.itype.gt.13) then
print *,"veg type: ",itype
CALL wrf_error_fatal3("<stdin>",1029,&
'module_sf_ssib: ERROR in vegetation/landuse map' )
endif
INTG=1
XADJ=0.
CTLPA=1.
NROOT=1
WFSOIL=0.
ZWIND=ZWIND2*0.5
IMONTH=1
IDAY=INT(DAY)
DO I=1,12
IF(IDAY.LE.IDAYS(I)) THEN
IMONTH=I
EXIT
ENDIF
ENDDO
IF(ZLAT.LT.0.0) THEN
MON_COR=IMONTH+6
IF(MON_COR.GT.12) MON_COR=MON_COR-12
ELSE
MON_COR=IMONTH
ENDIF
IF (ITIME.EQ.1) TA=TC
PSURF=PSUR*0.01
DTT =DDTT*FLOAT(INTG)
CALL VEGOUT(TRAN,REF,GREEN,VCOVER,CHIL, &
RSTPAR,TOPT,TL,TU,DEFAC,PH1,PH2, &
ZLT,Z0,XDD,Z2,Z1,RDC,RBC,ROOTD,SOREF, &
BEE, PHSAT, POROS, SATCO,SLOPE, &
ZDEPTH,MON_COR,ITYPE)
IF (ITYPE.EQ.12) CALL CROPS(ZLAT,DAY,CHIL, &
ZLT,GREEN,VCOVER,RSTPAR,TOPT,TL,TU,DEFAC,PH2,PH1)
IF (ITIME.EQ.1) THEN
STLEV1=0.05
STLEV2=1.05
DEPTH = TD_DEPTH(ITYPE)
IF (DEPTH.GT.STLEV1.AND.DEPTH.LE.STLEV2)THEN
TD = ( (STLEV2-DEPTH)*TGS + (DEPTH-STLEV1)*TD ) &
& /(STLEV2-STLEV1)
ELSE IF(DEPTH.GT.STLEV2)THEN
TD = ( (DEPTH-STLEV1)*TD - (DEPTH-STLEV2)*TGS) &
& /(STLEV2-STLEV1)
ENDIF
ENDIF
WWW(1) = WWW1 / POROS
WWW(2) = WWW2 / POROS
WWW(3) = WWW3 / POROS
SNOA = SNOA/1000.
SNOB = SNOB/1000.
CALL CONVDIM(0, &
DZO1,WO1,TSSN1,TSSNO1,BWO1,BTO1,CTO1,FIO1,FLO1,BIO1,BLO1,HO1, &
DZO2,WO2,TSSN2,TSSNO2,BWO2,BTO2,CTO2,FIO2,FLO2,BIO2,BLO2,HO2, &
DZO3,WO3,TSSN3,TSSNO3,BWO3,BTO3,CTO3,FIO3,FLO3,BIO3,BLO3,HO3, &
DZO4,WO4,TSSN4,TSSNO4,BWO4,BTO4,CTO4,FIO4,FLO4,BIO4,BLO4,HO4, &
DZO, WO, TSSN, TSSNO, BWO, BTO, CTO, FIO, FLO, BIO, BLO, HO )
IF (ITIME.EQ.1) THEN
ISNOW = 1
SNOWDEN = 3.75
SWE = SNOA
SNOWDEPTH = SWE * SNOWDEN
TGG=AMIN1(273.15,TGS)
ENDIF
CAPAC(1)=SNOB
CAPAC(2)=SNOA
IF (ITIME.EQ.1) THEN
IF (SNOA.GT.0.) THEN
CAPAC(1) = ZLT(1) * 0.0001
TC = AMIN1(TC ,TF-0.01)
TGS = AMIN1(TGS,TF-0.01)
ENDIF
ENDIF
UM=SQRT(UMM**2+VMM**2)
RHOAIR=100./GASR*(PSURF+0.01*PM)/(TM+TA)
AKAPPA = GASR/CPAIR
BPS =1.0 / EXP ( AKAPPA * ALOG (0.01*PM/PSURF) )
IF (ISNOW.EQ.0) THEN
TSOIL=TGS
TGS=TSSNO(N)
CAPAC(2)=SWE
IPTYPE=2
IF(TM.ge.TF) IPTYPE=1
END IF
EM=(PSURF*QM)/0.6220
IF (ITIME.EQ.1) EA=EM
SUNANG=AMAX1(SUNANGLE,0.01746)
IF (sw_physics.eq.3 .or. sw_physics.eq.4) THEN
radfrac11 = amax1(radfrac11,0.025)
radfrac12 = amax1(radfrac12,0.025)
radfrac21 = amax1(radfrac21,0.025)
radfrac22 = amax1(radfrac22,0.025)
swdown = radfrac11 + radfrac12 + radfrac21 + radfrac22
RADFRAC(1,1) = radfrac11/swdown
RADFRAC(1,2) = radfrac12/swdown
RADFRAC(2,1) = radfrac21/swdown
RADFRAC(2,2) = radfrac22/swdown
ELSE
swdown = amax1(swdown1,0.1)
CLOUD = AMAX1(CLOUD,0.0)
CLOUD = AMIN1(CLOUD,1.0)
CLOUD = AMAX1(0.58,CLOUD)
DIFRAT = 0.0604 / ( SUNANG-0.0223 ) + 0.0683
IF ( DIFRAT .LT. 0.0 ) DIFRAT = 0.0
IF ( DIFRAT .GT. 1.0 ) DIFRAT = 1.0
DIFRAT = DIFRAT + ( 1.0 - DIFRAT ) * CLOUD
VNRAT = ( 580.0 - CLOUD*464.0 ) / ( ( 580.0-CLOUD*499.0) &
& + ( 580.0 - CLOUD*464.0 ) )
RADFRAC(1,1) = (1.0-DIFRAT)*VNRAT
RADFRAC(1,2) = DIFRAT*VNRAT
RADFRAC(2,1) = (1.0-DIFRAT)*(1.0-VNRAT)
RADFRAC(2,2) = DIFRAT*(1.0-VNRAT)
ENDIF
RADN(1,1) = RADFRAC(1,1) * SWDOWN
RADN(1,2) = RADFRAC(1,2) * SWDOWN
RADN(2,1) = RADFRAC(2,1) * SWDOWN
RADN(2,2) = RADFRAC(2,2) * SWDOWN
RADN(3,1) = 0.
RADN(3,2) = RLWDOWN
CALL RADAB (TRAN,REF,GREEN,VCOVER,CHIL,ZLT,Z2,Z1,SOREF,TC, &
TGS,SATCAP,EXTK,RADFAC,THERMK,RADT,PAR,PD,ALBEDO,SALB, &
TGEFF,SUNANG,XADJ,CAPAC,RADN,ZLWUP,RADFRAC, &
ISNOW,SNOWDEN,SNOWDEPTH,SWDOWN,BEDO,SNOCV,0, &
fsdown,fldown,fsup,flup)
CALL ROOT1(PHSAT,BEE,WWW,PHSOIL)
CALL STOMA1(GREEN,VCOVER,CHIL,ZLT,PAR,PD,EXTK,SUNANG,RST, &
RSTPAR, CTLPA)
RSTUN = RST(1)
CALL INTERCS (DTT,VCOVER,ZLT,TM,TC,TGS,CAPAC,WWW,PPC,PPL, &
ROFF,ZDEPTH,POROS,CCX,CG,SATCO,SATCAP,SPWET, &
EXTK,ISNOW,P0,CSOIL,dzsoil,CHISL,SMELT)
CALL SET0(TSSNO,BWO,BLO,BIO,HO,FLO,FIO,WO,DZO, &
SSO,CTO,BTO,DMLTO,WF,DHP)
IF (ISNOW.EQ.0) THEN
PRCP=P0
TKAIR=TM
CALL GETMET(IPTYPE,PRCP,TKAIR, &
PRCPS,PRCPW,FIFALL,FLFALL,BIFALL,BLFALL)
SOLAR=0.
DO 1100 IVEG = 2, 2
DO 1100 IWAVE = 1, 2
DO 1100 IRAD = 1, 2
SOLAR=SOLAR+RADFAC(IVEG,IWAVE,IRAD)*RADN(IWAVE,IRAD)
1100 CONTINUE
CALL SNOW_1ST (DTT,TM,SOLAR,PRCPW,PRCPS,BIO,BLO,DICEVOL, &
DLIQVOL,TSSNO,PDZDTC,POROSITY,SO,SSO,WF,DHP,DZO,WO, &
BWO,BTO,CTO,DMASS,DSOL,SNROFF,HROFF,SNOWDEPTH,SOLSOIL, &
FLO,FIO,DMLTO,HO,BIFALL,BLFALL,FLFALL)
CALL TEMRS2(DTT,TC,TGS,TD,TA,TM,QM,EM,PSUR,WWW,CAPAC,SATCAP, &
DTC,DTG,RA,RST,ZDEPTH,BEE,PHSAT,POROS,XDD,Z0,RDC,RBC,VCOVER, &
Z2,ZLT,DEFAC,TU,TL,TOPT,RSTFAC,NROOT,ROOTD,PHSOIL,ROOTP, &
PH1,PH2,ECT,ECI,EGT,EGI,EGS,HC,HG,EC,EG,EA,RADT,CHF,SHF, &
ALBEDO,ZLWUP,THERMK,RHOAIR,ZWIND,UM,USTAR,DRAG,CCX,CG, &
ISNOW,CHISL,TSOIL,SOLSOIL,CSOIL,WFSOIL,POROSITY, &
DZ,DZO,W,WO,WF,TSSN,TSSNO,BW,BWO,CT,CTO,FI,FIO,FL,FLO, &
BL,BLO,BI,BIO,BT,BTO,DMLT,DMLTO,H,HO,S,SO,SS,SSO,DSOL,DHP, &
DICEVOL,DLIQVOL,THK,QK,SWE,SNOWDEN,SNOWDEPTH,TKAIR,SNROFF, &
DZSOIL,BPS,rib,CU,XCT,flup,ii,jj)
CALL OLD(TSSN,BW,BL,BI,H,FL,FI,W,DZ,SS,CT,BT,DMLT, &
TSSNO,BWO,BLO,BIO,HO,FLO,FIO,WO,DZO,SSO,CTO,BTO,DMLTO)
ELSE
CALL TEMRS1(DTT,TC,TGS,TD,TA,TM,QM,EM,PSUR,WWW,CAPAC,SATCAP, &
DTC,DTG,RA,RST,ZDEPTH,BEE,PHSAT,POROS,XDD,Z0,RDC,RBC,VCOVER,Z2, &
ZLT,DEFAC,TU,TL,TOPT,RSTFAC,NROOT,ROOTD,PHSOIL,ROOTP,PH1,PH2, &
ECT,ECI,EGT,EGI,EGS,HC,HG,EC,EG,EA,RADT,CHF,SHF,ALBEDO,ZLWUP, &
THERMK,RHOAIR,ZWIND,UM,USTAR,DRAG,CCX,CG, ISNOW,SNOWDEN, &
BPS,rib,CU,XCT,flup,ii,jj)
SWE=CAPAC(2)
SNOWDEPTH=SWE*SNOWDEN
SNROFF=0.
END IF
CALL UPDAT1(DTT,TC,TGS,TD,CAPAC,DTC,DTG,ECT,ECI,EGT,EGI, &
EGS,EG,HC,HG,HFLUX,ETMASS,ROFF, &
1,ROOTD,ROOTP,POROS,BEE,SATCO,SLOPE, &
PHSAT,ZDEPTH,WWW,CCX,CG,CHF,SHF, ISNOW,WFSOIL,SWE,SNROFF,SMELT)
IF (ISNOW.EQ.0) THEN
CAPAC(2)=SWE
IF (SNOWDEPTH.LT.SNODEP_CR) THEN
ISNOW=1
CALL LAYER1 (CSOIL,TGS,DZSOIL,H,W,SNOWDEPTH,SWE,STEMP,N2)
ELSE
ISNOW=0
CALL MODNODE(SNOWDEPTH,DZO,WO,HO,TSSNO,BWO,BIO, &
BLO,BTO,FIO,FLO,CTO,DLIQVOL,DICEVOL)
END IF
ELSE IF(ISNOW.GT.0) THEN
IF (CAPAC(2)*SNOWDEN.GT.SNODEP_CR) THEN
SWE=CAPAC(2)
SNOWDEPTH=CAPAC(2)*SNOWDEN
ISNOW=0
CALL LAYERN (TGS,SWE,SNOWDEPTH, DZO,BWO,WO,BTO,CTO,FLO,FIO, &
HO,BLO,BIO,DLIQVOL,DICEVOL,TSSNO,DMLTO)
ELSE
ISNOW=1
END IF
END IF
ROFF=ROFF+SNROFF
SOILM=(WWW(1)*ZDEPTH(1)+WWW(2)*ZDEPTH(2)+WWW(3)*ZDEPTH(3))*POROS
UMOM=RHOAIR*CU*USTAR*UMM
VMOM=RHOAIR*CU*USTAR*VMM
HLFLX= ETMASS/RHOAIR/DTT
HSFLX= HFLUX/CPAIR/RHOAIR/DTT
ULWSF1=TGEFF*TGEFF*TGEFF*TGEFF*STEFAN
Q2M=0.622*EA/(PSURF-EA)
EVAP=ETMASS*HLAT
CM=(USTAR*USTAR)/(UM*UM)
CH=1/(UM*RA)
FM=VKC/CU
FH=VKC/XCT
EVAPSOIL=EGS /DTT
EVAPWC=ECI /DTT
EVAPDC=ECT /DTT
EVAPSN=EGI /DTT
EVAPGX=EGT /DTT
ELATEN=EVAPSOIL+EVAPWC+EVAPDC+EVAPSN+EVAPGX
XHLFLX=ELATEN/HLAT
GHTFLX=CHF+SHF
xhsflx=(hc+hg)/dtt
CALL CONVDIM(1, &
DZO1,WO1,TSSN1,TSSNO1,BWO1,BTO1,CTO1,FIO1,FLO1,BIO1,BLO1,HO1, &
DZO2,WO2,TSSN2,TSSNO2,BWO2,BTO2,CTO2,FIO2,FLO2,BIO2,BLO2,HO2, &
DZO3,WO3,TSSN3,TSSNO3,BWO3,BTO3,CTO3,FIO3,FLO3,BIO3,BLO3,HO3, &
DZO4,WO4,TSSN4,TSSNO4,BWO4,BTO4,CTO4,FIO4,FLO4,BIO4,BLO4,HO4, &
DZO, WO, TSSN, TSSNO, BWO, BTO, CTO, FIO, FLO, BIO, BLO, HO )
WWW1=WWW(1)*POROS
WWW2=WWW(2)*POROS
WWW3=WWW(3)*POROS
SNOA = CAPAC(2)
SNOB = CAPAC(1)
SNOA = SNOA*1000.
SNOB = SNOB*1000.
SALB11=SALB(1,1)
SALB12=SALB(1,2)
SALB21=SALB(2,1)
SALB22=SALB(2,2)
xlhf = elaten
xshf = xhsflx
xghf = ghtflx
xegs = evapsoil
xeci = evapwc
xect = evapdc
xegi = evapsn
xegt = evapgx
xsdn = fsdown
xsup = fsup
xldn = fldown
xlup = flup
xwat = soilm
xhcx = hc/dtt
xhgx = hg/dtt
xzlt = zlt(1)
xvcf = vcover(1)
xxz0 = z0
xveg = float(itype)
END SUBROUTINE SSIB
SUBROUTINE SSIB_SEAICE &
( II, JJ, DDTT, ITIME, ZLAT, SUNANGLE, &
PPL, PPC, RLWDOWN, ZWIND2, &
WWW1, WWW2, WWW3, &
TC, TGS, TD, &
SNOA, ROFF, YICE, &
UMM, VMM, QM, TM, &
PM, PSUR, &
SWDOWN1, SNOB, &
SALB11, SALB12, SALB21, SALB22, &
RADFRAC11, RADFRAC12, RADFRAC21, RADFRAC22, &
XHSFLX, ELATEN, GHTFLX, XHLFLX, TGEFF, &
USTAR, RIB, FM, FH, CM, &
XLHF, XSHF, XGHF, &
XSDN, XSUP, XLDN, XLUP, &
XWAT, XXZ0, &
XVEG, &
DAY, CLOUD, Q2M, TA, BEDO, &
sw_physics,ice_threshold &
)
INTEGER, DIMENSION (12) :: IDAYS
REAL, DIMENSION (2) :: CAPAC, SATCAP, GREEN, VCOVER, ZLT, CHIL, TOPT, TL, &
TU, DEFAC, PH1, PH2, RST, ROOTD, RADT, PAR, PD
REAL, DIMENSION (3) :: WWW, SOREF, ZDEPTH, ROOTP, PHSOIL, YMATT, YMATQ
REAL, DIMENSION (2,2) :: RADFRAC, SALB
REAL, DIMENSION (2,3) :: RSTPAR
REAL, DIMENSION (2,4) :: RSTFAC
REAL, DIMENSION (3,2) :: RADN
REAL, DIMENSION (2,3,2) :: TRAN, REF, ALBEDO, EXTK
REAL, DIMENSION (13) :: TD_DEPTH
REAL :: ice_threshold
INTEGER :: sw_physics
DATA IDAYS/31,59,90,120,151,181,212,243,273,304,334,366/
DATA TD_DEPTH/1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.0, 1.0, 1.0 &
& , 0.5, 0.5, 1.5, 1.5/
XADJ=0.
CTLPA=1.
NROOT=1
INTG=1
ITYPE=13
ZWIND=ZWIND2*0.5
IMONTH=1
IDAY=INT(DAY)
DO I=1,12
IF(IDAY.LE.IDAYS(I)) THEN
IMONTH=I
EXIT
ENDIF
ENDDO
IF(ZLAT.LT.0.0) THEN
MON_COR=IMONTH+6
IF(MON_COR.GT.12) MON_COR=MON_COR-12
ELSE
MON_COR=IMONTH
ENDIF
IF (ITIME.EQ.1) TA=TC
PSURF=PSUR*0.01
DTT =DDTT*FLOAT(INTG)
CALL VEGOUT(TRAN,REF,GREEN,VCOVER,CHIL, &
RSTPAR,TOPT,TL,TU,DEFAC,PH1,PH2, &
ZLT,Z0,XDD,Z2,Z1,RDC,RBC,ROOTD,SOREF, &
BEE, PHSAT, POROS, SATCO,SLOPE, &
ZDEPTH,MON_COR,ITYPE)
IF (ITIME.EQ.1) THEN
STLEV1=0.05
STLEV2=1.05
DEPTH = TD_DEPTH(ITYPE)
IF (DEPTH.GT.STLEV1.AND.DEPTH.LE.STLEV2)THEN
TD = ( (STLEV2-DEPTH)*TGS + (DEPTH-STLEV1)*TD ) &
& /(STLEV2-STLEV1)
ELSE IF(DEPTH.GT.STLEV2)THEN
TD = ( (DEPTH-STLEV1)*TD - (DEPTH-STLEV2)*TGS) &
& /(STLEV2-STLEV1)
ENDIF
ENDIF
WWW(1) = 1.
WWW(2) = 1.
WWW(3) = 1.
SNOA = SNOA/1000.
SNOB = SNOB/1000.
CAPAC(1)=SNOB
CAPAC(2)=SNOA
SNOWDEN = 3.75
ISNOW = 1
IF (ITIME.EQ.1) THEN
TA=TGS
CAPAC(1)=0.
CAPAC(2)=0.
IF (SNOA.GT.0.) CAPAC(1) = ZLT(1) * 0.0001
TC = AMIN1(TC ,273.15)
TGS= AMIN1(TGS,273.15)
TD = AMIN1(TD ,272.50)
ELSE
IF( YICE .LT. ice_threshold ) THEN
CAPAC(1)= 0.
CAPAC(2)= 0.
XADIA = EXP(GASR/CPAIR*LOG(PSUR/PM))
XX = MIN(TM*XADIA,273.15)
TC = MIN(TM*XADIA,273.15)
TGS= MIN(TM*XADIA,273.15)
IF(TD.EQ.0.) TD=272.5
TD = MIN(TD,272.5)
ENDIF
ENDIF
UM=SQRT(UMM**2+VMM**2)
RHOAIR=100./GASR*(PSURF+0.01*PM)/(TM+TA)
AKAPPA = GASR/CPAIR
BPS =1.0 / EXP ( AKAPPA * ALOG (0.01*PM/PSURF) )
EM=(PSURF*QM)/0.6220
IF (ITIME.EQ.1) EA=EM
SUNANG=AMAX1(SUNANGLE,0.01746)
IF (sw_physics.eq.3 .or. sw_physics.eq.4) THEN
radfrac11 = amax1(radfrac11,0.025)
radfrac12 = amax1(radfrac12,0.025)
radfrac21 = amax1(radfrac21,0.025)
radfrac22 = amax1(radfrac22,0.025)
swdown = radfrac11 + radfrac12 + radfrac21 + radfrac22
RADFRAC(1,1) = radfrac11/swdown
RADFRAC(1,2) = radfrac12/swdown
RADFRAC(2,1) = radfrac21/swdown
RADFRAC(2,2) = radfrac22/swdown
ELSE
swdown = amax1(swdown1,0.1)
CLOUD = AMAX1(CLOUD,0.0)
CLOUD = AMIN1(CLOUD,1.0)
CLOUD = AMAX1(0.58,CLOUD)
DIFRAT = 0.0604 / ( SUNANG-0.0223 ) + 0.0683
IF ( DIFRAT .LT. 0.0 ) DIFRAT = 0.0
IF ( DIFRAT .GT. 1.0 ) DIFRAT = 1.0
DIFRAT = DIFRAT + ( 1.0 - DIFRAT ) * CLOUD
VNRAT = ( 580.0 - CLOUD*464.0 ) / ( ( 580.0-CLOUD*499.0) &
& + ( 580.0 - CLOUD*464.0 ) )
RADFRAC(1,1) = (1.0-DIFRAT)*VNRAT
RADFRAC(1,2) = DIFRAT*VNRAT
RADFRAC(2,1) = (1.0-DIFRAT)*(1.0-VNRAT)
RADFRAC(2,2) = DIFRAT*(1.0-VNRAT)
ENDIF
RADN(1,1) = RADFRAC(1,1) * SWDOWN
RADN(1,2) = RADFRAC(1,2) * SWDOWN
RADN(2,1) = RADFRAC(2,1) * SWDOWN
RADN(2,2) = RADFRAC(2,2) * SWDOWN
RADN(3,1) = 0.
RADN(3,2) = RLWDOWN
CALL RADAB_ICE(TRAN,REF,GREEN,VCOVER,CHIL,ZLT,Z2,Z1,SOREF, &
TC,TGS,SATCAP,EXTK,CLOSS,GLOSS,THERMK,P1F,P2F, &
RADT,PAR,PD,SALB,ALBEDO,TGEFF,SUNANG,XADJ,CAPAC, &
RADN,BEDO,ZLWUP,RADFRAC,SWDOWN,SNOCV,1, &
fsdown,fldown,fsup,flup)
CALL ROOT1(PHSAT,BEE,WWW,PHSOIL)
CALL STOMA1(GREEN,VCOVER,CHIL,ZLT,PAR,PD,EXTK,SUNANG,RST, &
RSTPAR, CTLPA)
POROSAVE=POROS
POROS=0.95
CALL INTERC(DTT ,VCOVER,ZLT,TM,TC,TGS,CAPAC,WWW,PPC,PPL,ROFF, &
ZDEPTH,POROS,CCX,CG,SATCO,SATCAP,SPWET,EXTK,RNOFFS,FILTR, &
SMELT)
CALL TEMRS1(DTT,TC,TGS,TD,TA,TM,QM,EM,PSUR,WWW,CAPAC,SATCAP, &
DTC,DTG,RA,RST,ZDEPTH,BEE,PHSAT,POROS,XDD,Z0,RDC,RBC,VCOVER,Z2, &
ZLT,DEFAC,TU,TL,TOPT,RSTFAC,NROOT,ROOTD,PHSOIL,ROOTP,PH1,PH2, &
ECT,ECI,EGT,EGI,EGS,HC,HG,EC,EG,EA,RADT,CHF,SHF,ALBEDO,ZLWUP, &
THERMK,RHOAIR,ZWIND,UM,USTAR,DRAG,CCX,CG, ISNOW,SNOWDEN, &
BPS,rib,CU,XCT,flup,ii,jj)
CALL UPDAT1_ICE(DTT ,TC,TGS,TD,CAPAC,DTC,DTG,ECT,ECI,EGT,EGI, &
EGS,EG,HC,HG,HFLUX,ETMASS,FILTR,SOILDIF,SOILDRA,ROFF, &
RNOFFB,RNOFFS,NROOT,ROOTD,ROOTP,POROS,BEE,SATCO,SLOPE, &
PHSAT,ZDEPTH,WWW,CCX,CG,CHF,SHF,SMELT)
POROS=POROSAVE
TD = AMIN1(TD ,273.15)
TC = AMIN1(TC ,273.15)
TGS = AMIN1(TGS,273.15)
SOILM=(WWW(1)*ZDEPTH(1)+WWW(2)*ZDEPTH(2)+WWW(3)*ZDEPTH(3))*POROS
UMOM=RHOAIR*CU*USTAR*UMM
VMOM=RHOAIR*CU*USTAR*VMM
HLFLX= ETMASS/RHOAIR/DTT
HSFLX= HFLUX/CPAIR/RHOAIR/DTT
ULWSF1=TGEFF*TGEFF*TGEFF*TGEFF*STEFAN
Q2M=0.622*EA/(PSURF-EA)
EVAP=ETMASS*HLAT
CM=(USTAR*USTAR)/(UM*UM)
CH=1/(UM*RA)
FM=VKC/CU
FH=VKC/XCT
ELATEN=EVAP/DTT
XHLFLX=ELATEN/HLAT
GHTFLX=CHF+SHF
xhsflx=(hc+hg)/dtt
WWW1=WWW(1)*POROS
WWW2=WWW(2)*POROS
WWW3=WWW(3)*POROS
SNOA = CAPAC(2)
SNOB = CAPAC(1)
SALB11=SALB(1,1)
SALB12=SALB(1,2)
SALB21=SALB(2,1)
SALB22=SALB(2,2)
xlhf = elaten
xshf = xhsflx
xghf = ghtflx
xsdn = fsdown
xsup = fsup
xldn = fldown
xlup = flup
xwat = soilm
xxz0 = z0
xveg = float(itype)
END SUBROUTINE SSIB_SEAICE
SUBROUTINE CROPS(XLAT,DAY,CHIL,ZLT,GREEN,XCOVER &
,RSTPAR,TOPT,TL,TU,DEFAC,PH2,PH1)
REAL, DIMENSION (2) :: GREEN, XCOVER, CHIL, ZLT, TOPT, TL, TU, DEFAC, PH1, PH2
REAL, DIMENSION (2,3) :: RSTPAR
REAL, DIMENSION (9) :: PHENST, WLAI, WGRN
DATA WLAI/1.0, 2.0, 6.0, 4.0, 3.0, 1.0, 0.01, 0.01, 1.0/
DATA WGRN/0.6, 0.9, 0.8, 0.5, 0.2, 0.1, 0.01, 0.01, 0.6/
DATA IHEAD,IEND,DEND,IWHEAT/3,9,244.,12/,SYR/365.25E0/
IF (XLAT.LT.0.) THEN
RDAY= DAY+184
IF (RDAY.GT.365) RDAY=RDAY-365
ELSE
RDAY= DAY
END IF
JULDAY=INT(RDAY+0.2)
PHI=XLAT
APHI = ABS(PHI)
IF (APHI.GT.55.) PHI=SIGN(55.,PHI)
IF (APHI.LT.20.) PHI=SIGN(20.,PHI)
FLIP = 0.0
PHENST(2) = 4.50 *ABS(PHI) - 64.0 + FLIP
PHENST(3) = 4.74 *ABS(PHI) - 46.2 + FLIP
PHENST(4) = 4.86 *ABS(PHI) - 30.8 + FLIP
PHENST(5) = 4.55 *ABS(PHI) - 3.0 + FLIP
PHENST(6) = 4.35 *ABS(PHI) + 11.5 + FLIP
PHENST(7) = PHENST(6) + 3.0
DEMG = ABS( 5.21 *ABS(PHI) - 0.3 )
PHENST(1) = PHENST(2) - DEMG
PHENST(9) = PHENST(1)
PHENST(8) = PHENST(9) - 5.0
DO 10 NS = 1,9
IF(PHENST(NS) .LT. 0.0E0)PHENST(NS) = PHENST(NS) + 365.
IF(PHENST(NS) .GT. 365. )PHENST(NS) = PHENST(NS) - 365.
10 CONTINUE
ROOTGC = 1.0
CHILW =-0.02
TLAI = 0.5
GRLF = 0.6
DO 50 NS = 1,8
TOP = PHENST(NS+1)
BOT = PHENST(NS)
DIFF1 = TOP-BOT
DIFF2 = RDAY-BOT
IF(RDAY.GE. BOT .AND. RDAY .LE. TOP ) GO TO 40
IF(BOT .LT. TOP ) GO TO 50
ICOND = 0
IF(RDAY .GE. BOT .AND. RDAY .LE. 365.) ICOND = 1
IF(RDAY .GE. 0.0 .AND. RDAY .LE. TOP ) ICOND = 2
IF(ICOND .EQ. 0)GO TO 50
IF(ICOND .EQ. 2)GO TO 35
DIFF1 = 365. - BOT + TOP
DIFF2 = RDAY - BOT
GO TO 40
35 CONTINUE
DIFF1 = 365. - BOT + TOP
DIFF2 = 365. - BOT + RDAY
40 CONTINUE
IF ((RDAY.GT.PHENST(IHEAD)).AND.(RDAY.LE.DEND)) THEN
TLAI=WLAI(IHEAD)
GRLF=WGRN(IHEAD)
GO TO 77
END IF
IF ((RDAY.GT.DEND).AND.(RDAY.LE.PHENST(IEND))) THEN
DIFF1=PHENST(IEND)-DEND
DIFF2=RDAY-DEND
PERC = DIFF2/DIFF1
TLAI = PERC*(WLAI(IEND)-WLAI(IHEAD)) + WLAI(IHEAD)
GRLF = PERC*(WGRN(IEND)-WGRN(IHEAD)) + WGRN(IHEAD)
GO TO 77
END IF
PERC = DIFF2/DIFF1
TLAI = PERC*(WLAI(NS+1)-WLAI(NS)) + WLAI(NS)
GRLF = PERC*(WGRN(NS+1)-WGRN(NS)) + WGRN(NS)
77 CONTINUE
GO TO 95
50 CONTINUE
95 CONTINUE
XCOVER(1)=0.90*(1.0 - EXP(-TLAI))
ZLTGMX = WLAI(IHEAD)
ROOTGC = 2910.0 * (0.5 +0.5 *TLAI/ZLTGMX * GRLF)
IF (NS.NE.1.AND.NS.NE.2) CHILW=-0.2
ZLT (1) = TLAI
GREEN (1) = GRLF
CHIL (1) = CHILW
END SUBROUTINE CROPS
SUBROUTINE ROOT1(PHSAT,BEE,WWW,PHSOIL)
REAL, DIMENSION (3) :: WWW, PHSOIL
DO 1000 IL = 1, 3
PHSOIL(IL) = PHSAT * AMAX1( 0.05, WWW(IL) ) ** ( - BEE )
1000 CONTINUE
END SUBROUTINE ROOT1
SUBROUTINE STOMA1(GREEN,VCOVER,CHIL,ZLT,PAR,PD,EXTK,SUNANG,RST, &
RSTPAR,CTLPA)
REAL, DIMENSION (2) :: GREEN, VCOVER, ZLT, CHIL, PAR, PD, RST
REAL, DIMENSION (2,3) :: RSTPAR
REAL, DIMENSION (2,3,2) :: EXTK
DO 1000 IVEG = 1, 2
AT = ZLT(IVEG) / VCOVER(IVEG)
IF (SUNANG .LE. 0.02) THEN
XABC = RSTPAR(IVEG,1) / RSTPAR(IVEG,2) + RSTPAR(IVEG,3)
RST(IVEG) = 0.5 / XABC * AT
IF (RST(IVEG) .LT. 0.) RST(IVEG) = 0.00001
GO TO 1010
END IF
GAMMA = ( RSTPAR(IVEG,1) + RSTPAR(IVEG,2) * RSTPAR(IVEG,3) ) / &
RSTPAR(IVEG,3)
POWER1 = AMIN1( 50., AT * EXTK(IVEG,1,1) )
POWER2 = AMIN1( 50., AT * EXTK(IVEG,1,2) )
AA = 0.5 - 0.633 * CHIL(IVEG)- 0.33 * CHIL(IVEG)* CHIL(IVEG)
BB = 0.877 * ( 1. - 2. * AA )
ZAT = ALOG( ( EXP(-POWER1) + 1. )/2. ) * PD(IVEG) &
/ ( POWER1/AT )
ZAT = ZAT + ALOG( ( EXP(-POWER2) + 1. )/2. ) &
* ( 1. - PD(IVEG) ) / ( POWER2/AT )
POW1 = AMIN1( 50., (POWER1*ZAT/AT) )
POW2 = AMIN1( 50., (POWER2*ZAT/AT) )
ZK = 1. / ZAT * ALOG( PD(IVEG) * EXP ( POW1 ) &
+ ( 1. - PD(IVEG) ) * EXP ( POW2 ) )
POW = AMIN1( 50., ZK*AT )
EKAT = EXP ( POW )
AVFLUX = PAR(IVEG) * ( PD(IVEG) / SUNANG * ( AA + BB * SUNANG ) &
+ ( 1. - PD(IVEG) )*( BB / 3. + AA * 1.5 &
+ BB / 4. * PIE ))
RHO4 = GAMMA / AVFLUX
RST(IVEG) = RSTPAR(IVEG,2)/GAMMA * ALOG(( RHO4 * EKAT + 1. ) / &
( RHO4 + 1. ) )
RST(IVEG) = RST(IVEG) - ALOG (( RHO4 + 1. / EKAT ) / &
( RHO4 + 1. ) )
RST(IVEG) = RST(IVEG) / ( ZK * RSTPAR(IVEG,3) )
1010 RST(IVEG) = 1. / ( RST(IVEG) * GREEN(IVEG) + 0.0000001)
1000 CONTINUE
RST(1) = RST(1) * CTLPA
END SUBROUTINE STOMA1
SUBROUTINE VEGOUT(XTRAN,XREF,XGREEN,XVCOVER,XCHIL, &
XRSTPAR,XTOPT,XTL,XTU,XDEFAC,XPH1,XPH2, &
XZLT,XZ0,XDD,XZ2,XZ1,XRDC,XRBC,XROOTD,XSOREF, &
XBEE, XPHSAT, XPOROS, XSATCO,XSLOPE, &
XDEPTH,MONTH,ITYPE)
REAL, DIMENSION (2) :: XGREEN, XVCOVER, XZLT, XCHIL, XTOPT, XTL, &
XTU, XDEFAC, XPH1, XPH2, XROOTD
REAL, DIMENSION (3) :: XSOREF, XDEPTH
REAL, DIMENSION (2,3) :: XRSTPAR
REAL, DIMENSION (2,3,2) :: XTRAN, XREF
DO IW=1,3
XTRAN(1,IW,1)=TRAN0(ITYPE,1,IW,1)
XTRAN(1,IW,2)=TRAN0(ITYPE,1,IW,2)
XTRAN(2,IW,1)=TRAN0(ITYPE,2,IW,1)
XTRAN(2,IW,2)=TRAN0(ITYPE,2,IW,2)
XREF (1,IW,1)= REF0(ITYPE,1,IW,1)
XREF (1,IW,2)= REF0(ITYPE,1,IW,2)
XREF (2,IW,1)= REF0(ITYPE,2,IW,1)
XREF (2,IW,2)= REF0(ITYPE,2,IW,2)
XRSTPAR(1,IW)=RSTPAR0(ITYPE,1,IW)
XRSTPAR(2,IW)=RSTPAR0(ITYPE,2,IW)
XSOREF (IW) =SOREF0(ITYPE,IW)
END DO
DO IV=1,2
XCHIL(IV)=CHIL0(ITYPE,IV)
XTOPT(IV)=TOPT0(ITYPE,IV)
XTL(IV)=TL0(ITYPE,IV)
XTU(IV)=TU0(ITYPE,IV)
XDEFAC(IV)=DEFAC0(ITYPE,IV)
XPH1(IV)=PH10(ITYPE,IV)
XPH2(IV)=PH20(ITYPE,IV)
XROOTD(IV)=ROOTD0(ITYPE,IV)
XZLT(IV)=ZLT0(ITYPE,MONTH,IV)
XGREEN(IV)=GREEN0(ITYPE,MONTH,IV)
XVCOVER(IV)=VCOVER0(ITYPE,MONTH,IV)
END DO
DO IDEP=1,3
XDEPTH(IDEP)=DEPTH0(ITYPE,IDEP)
END DO
XBEE=BEE0(ITYPE)
XPHSAT=PHSAT0(ITYPE)
XSATCO=SATCO0(ITYPE)
XPOROS=POROS0(ITYPE)
XSLOPE=SLOPE0(ITYPE)
XZ2=Z20(ITYPE,MONTH)
XZ1=Z10(ITYPE,MONTH)
XZ0= Z000(ITYPE,MONTH)
XDD= D0(ITYPE,MONTH)
XRBC=RBC0 (ITYPE,MONTH)
XRDC=RDC0 (ITYPE,MONTH)
END SUBROUTINE VEGOUT
SUBROUTINE COMBO (DDZ2,DZP,DZM,WP,WM,HP,HM,TP,TM,BWP,BWM,BIP, &
BIM,BLP,BLM,BTP,BTM,FIP,FIM,FLP,FLM,CTP,CTM, &
DLIQVOLP,DLIQVOLM,DICEVOLP,DICEVOLM)
RATIO= DDZ2/dzm
dzp=dzp + RATIO*dzm
wp = wp + RATIO*wm
hp = hp + RATIO*hm
bwp= wp*rhowater/dzp
btp= bwp
ctp= (1.9e6)*(bwp/920.0)
dmlt=wp*rhowater*dlm
if(hp.ge.(-1.0)*dmlt)then
tp=273.16
fip=(-1.0)*hp/dmlt
flp=1.0-fip
blp=bwp*flp
bip=bwp*fip
dliqvolp = blp/rhowater
dicevolp = bip/dice
else
flp=0.0
fip=1.0
tp=(hp+dmlt)/(ctp*dzp)+273.16
bip=bwp
blp=0.0
dliqvolp = 0.0
dicevolp = bip/dice
endif
dzm=dzm - RATIO*dzm
wm = wm - RATIO*wm
hm = hm - RATIO*hm
END SUBROUTINE COMBO
SUBROUTINE COMPACT(BI,T,BL,OVERBURDEN,PDZDT,SS,DICE)
data c2,c3,c4,c5/23d-3,2.777d-6,0.04,2.0/
data dm/150/
data eta0/0.9d6/
if(bi .ge. dice .or. ss .ge. 1.) return
ddz1=-c3*exp(-c4*(273.15-t))
if(bi .gt. dm) ddz1=ddz1*exp(-46.0d-3*(bi-dm))
if(bl .gt. 0.01) ddz1=ddz1*c5
ddz2=-overburden*exp(-0.08*(273.15-t)-c2*bi)/eta0
ddz3=0d0
pdzdt=ddz1+ddz2+ddz3
END SUBROUTINE COMPACT
SUBROUTINE GETMET(IPTYPE,PRCP_TOTAL,TAIR, &
PRCP_S,PRCP_W,FI_FALL,FL_FALL,BI_FALL,BL_FALL)
IF (PRCP_TOTAL.gt.0.) THEN
IF(IPTYPE.EQ.2)THEN
PRCP_S=PRCP_TOTAL
PRCP_W=0.0
ELSE IF(IPTYPE.EQ.1)THEN
PRCP_W=PRCP_TOTAL
PRCP_S=0.0
FL_FALL=1.0
FI_FALL=0.
BL_FALL=1000.0
BI_FALL=0.0
ENDIF
ELSE
PRCP_W=0.0
PRCP_S=0.0
IPTYPE = 0
RETURN
END IF
IF (IPTYPE.NE.1) THEN
IF (TAIR .GT. 275.15) THEN
BI_FALL =189
ELSE IF (TAIR.GT.258.16)THEN
BI_FALL=50+1.7*(TAIR-258.16)**1.5d0
ELSE
BI_FALL=50
ENDIF
FL_FALL = 0
FI_FALL=1.0
BL_FALL=0.0
ENDIF
END SUBROUTINE GETMET
SUBROUTINE INTERCS (DTT,VCOVER,ZLAI,TM,TC,TGS,CAPAC,WWW,PPC,PPL, &
ROFF,ZDEPTH,POROS,CCX,CG,SATCO,SATCAP,SPWET, &
EXTK,ISNOW,P0,CSOIL,DZSOIL, &
CHISL,SMELT)
DIMENSION CAPACP(2), SNOWP(2), PCOEFS(2,2)
DATA PCOEFS(1,1)/ 20. /, PCOEFS(1,2)/ .206E-8 /, &
PCOEFS(2,1)/ 0.0001 /, PCOEFS(2,2)/ 0.9999 /, BP /20. /
DIMENSION VCOVER(2),ZLAI(2),WWW(3),CAPAC(2),SATCAP(2),EXTK(2,3,2)
DIMENSION ZDEPTH(3),SNOWW(2)
AP = PCOEFS(2,1)
CP = PCOEFS(2,2)
TOTALP = PPC + PPL
IF(TOTALP.LT.1.E-8)GO TO 6000
AP = PPC/TOTALP * PCOEFS(1,1) + PPL/TOTALP * PCOEFS(2,1)
CP = PPC/TOTALP * PCOEFS(1,2) + PPL/TOTALP * PCOEFS(2,2)
6000 CONTINUE
ROFF = 0.
THRU = 0.
FPI = 0.
THETA=WWW(1)*POROS
CHISL=( 9.8E-4+1.2E-3*THETA )/( 1.1-0.4*THETA )
CHISL=CHISL*4.186E2
DIFSL=5.E-7
ROCS =CHISL/DIFSL
D1 =SQRT(DIFSL*86400.0)
CSOIL=ROCS*D1/SQRT(PIE)/2.0
dzsoil=D1/SQRT(PIE)/2.0
THALAS=0.
OCEANS=0.
POLAR=0.
CSOIL=CSOIL*(1.0-THALAS)+10.E10*OCEANS+POLAR*3.6*4.2E4
P0 = TOTALP * 0.001
DO 1000 IVEG = 1, 2
SPWET1 = AMIN1 ( 0.05, CAPAC(IVEG))*CW
TS = TC
SPECHT = ZLAI(1) * CLAI
IF ( IVEG .EQ. 1 ) GO TO 1100
TS = TGS
SPECHT = CSOIL
1100 CONTINUE
XSC = AMAX1(0., CAPAC(IVEG) - SATCAP(IVEG) )
IF(IVEG.EQ.2 .AND. TS.LE.TF )GO TO 1170
CAPAC(IVEG) = CAPAC(IVEG) - XSC
ROFF = ROFF + XSC
1170 CONTINUE
CAPACP(IVEG) = 0.
SNOWP(IVEG) = 0.
IF( TS .GT. TF ) CAPACP(IVEG) = CAPAC(IVEG)
IF( TS .LE. TF ) SNOWP(IVEG) = CAPAC(IVEG)
CAPAC(IVEG) = CAPACP(IVEG)
SNOWW(IVEG) = SNOWP(IVEG)
ZLOAD = CAPAC(IVEG) + SNOWW(IVEG)
FPI = ( 1.-EXP( - EXTK(IVEG,3,1) * ZLAI(IVEG)/VCOVER(IVEG) ) ) &
* VCOVER(IVEG)
TTI = P0 * ( 1.-FPI )
XS = 1.
IF ( P0 .LT. 1.E-9 ) GO TO 1150
ARG = ( SATCAP(IVEG)-ZLOAD )/( P0*FPI*AP ) -CP/AP
IF ( ARG .LT. 1.E-9 ) GO TO 1150
XS = -1./BP * ALOG( ARG )
XS = AMIN1( XS, 1. )
XS = AMAX1( XS, 0. )
1150 TEX = P0*FPI * ( AP/BP*( 1.- EXP( -BP*XS )) + CP*XS ) - &
( SATCAP(IVEG) - ZLOAD ) * XS
TEX = AMAX1( TEX, 0. )
THRU = TTI + TEX
IF(IVEG.EQ.2.AND.TGS.LE.TF)THRU = 0.
PINF = P0 - THRU
IF( TM .GT. TF ) CAPAC(IVEG) = CAPAC(IVEG) + PINF
IF( TM .LE. TF ) SNOWW(IVEG) = SNOWW(IVEG) + PINF
IF( IVEG .EQ. 1 ) GO TO 1300
IF( TM .GT. TF ) GO TO 1200
SNOWW(IVEG) = SNOWP(IVEG) + P0
THRU = 0.
GO TO 1300
1200 EQUDEP = SATCO * DTT
XS = 1.
IF ( THRU .LT. 1.E-9 ) GO TO 1250
ARG = EQUDEP / ( THRU * AP ) -CP/AP
IF ( ARG .LT. 1.E-9 ) GO TO 1250
XS = -1./BP * ALOG( ARG )
XS = AMIN1( XS, 1. )
XS = AMAX1( XS, 0. )
1250 ROFFO = THRU * ( AP/BP * ( 1.-EXP( -BP*XS )) + CP*XS ) &
-EQUDEP*XS
ROFFO = AMAX1 ( ROFFO, 0. )
ROFF = ROFF + ROFFO
WWW(1) = WWW(1) + (THRU - ROFFO) / ( POROS*ZDEPTH(1) )
1300 CONTINUE
DIFF = ( CAPAC(IVEG)+SNOWW(IVEG) - CAPACP(IVEG)-SNOWP(IVEG) )*CW
CCP = SPECHT + SPWET1
CCT = SPECHT + SPWET1 + DIFF
TSD = ( TS * CCP + TM * DIFF ) / CCT
FREEZE = 0.
IF ( TS .GT. TF .AND. TM .GT. TF ) GO TO 2000
IF ( TS .LE. TF .AND. TM .LE. TF ) GO TO 2000
TTA = TS
TTB = TM
CCA = CCP
CCB = DIFF
IF ( TSD .GT. TF ) GO TO 2100
CCC = CAPACP(IVEG) * SNOMEL
IF ( TS .LT. TM ) CCC = DIFF * SNOMEL / CW
TSD = ( TTA * CCA + TTB * CCB + CCC ) / CCT
FREEZE = ( TF * CCT - ( TTA * CCA + TTB * CCB ) )
FREEZE = (AMIN1 ( CCC, FREEZE )) / SNOMEL
IF(TSD .GT. TF)TSD = TF - 0.1
GO TO 2000
2100 CONTINUE
CCC = - SNOWW(IVEG) * SNOMEL
IF ( TS .GT. TM ) CCC = - DIFF * SNOMEL / CW
TSD = ( TTA * CCA + TTB * CCB + CCC ) / CCT
FREEZE = ( TF * CCT - ( TTA * CCA + TTB * CCB ) )
FREEZE = (AMAX1( CCC, FREEZE )) / SNOMEL
IF(TSD .LE. TF)TSD = TF - 0.1
2000 CONTINUE
SMELT = FREEZE
SNOWW(IVEG) = SNOWW(IVEG) + FREEZE
CAPAC(IVEG) = CAPAC(IVEG) - FREEZE
IF( IVEG .EQ. 1 ) TC = TSD
IF( IVEG .EQ. 2 ) TGS = TSD
IF( SNOWW(IVEG) .LT. 0.0000001 ) GO TO 3000
ZMELT = CAPAC(IVEG)
CAPAC(IVEG) = 0.
WWW(1) = WWW(1) + ZMELT / ( POROS * ZDEPTH(1) )
3000 CONTINUE
CAPAC(IVEG) = CAPAC(IVEG) + SNOWW(IVEG)
SNOWW(IVEG) = 0.
P0 = THRU
IF (ISNOW.eq.0) go to 1001
1000 CONTINUE
1001 CCX = ZLAI(1) * CLAI + CAPAC(1) * CW
SPWET = AMIN1 ( 0.05, CAPAC(2)) * CW
CG = (CSOIL + SPWET)
END SUBROUTINE INTERCS
SUBROUTINE INTERC(DTT,VCOVER,ZLT,TM,TC,TGS,CAPAC,WWW,PPC,PPL,ROFF, &
ZDEPTH,POROS,CCX,CG,SATCO,SATCAP,SPWET,EXTK,RNOFFS,FILTR,SMELT)
REAL, DIMENSION (2) :: VCOVER, ZLT, CAPAC, SATCAP, SNOWW, CAPACP, SNOWP
REAL, DIMENSION (3) :: WWW, ZDEPTH
REAL, DIMENSION (2,2) :: PCOEFS
REAL, DIMENSION (2,3,2) :: EXTK
DATA PCOEFS(1,1)/ 20. /, PCOEFS(1,2)/ .206E-8 /, &
PCOEFS(2,1)/ 0.0001 /, PCOEFS(2,2)/ 0.9999 /, BP /20. /
AP = PCOEFS(2,1)
CP = PCOEFS(2,2)
TOTALP = PPC + PPL
IF(TOTALP.LT.1.E-8)GO TO 6000
AP = PPC/TOTALP * PCOEFS(1,1) + PPL/TOTALP * PCOEFS(2,1)
CP = PPC/TOTALP * PCOEFS(1,2) + PPL/TOTALP * PCOEFS(2,2)
6000 CONTINUE
ROFF = 0.
THRU = 0.
FPI = 0.
THETA=WWW(1)*POROS
CHISL=( 9.8E-4+1.2E-3*THETA )/( 1.1-0.4*THETA )
CHISL=CHISL*4.186E2
DIFSL=5.E-7
ROCS =CHISL/DIFSL
D1 =SQRT(DIFSL*86400.0)
CSOIL=ROCS*D1/SQRT(PIE)/2.0
THALAS=0.
OCEANS=0.
POLAR=0.
CSOIL=CSOIL*(1.0-THALAS)+10.E10*OCEANS+POLAR*3.6*4.2E4
P0 = TOTALP * 0.001
DO 1000 IVEG = 1, 2
SPWET1 = AMIN1 ( 0.05, CAPAC(IVEG))*CW
TS = TC
SPECHT = ZLT(1) * CLAI
IF ( IVEG .EQ. 1 ) GO TO 1100
TS = TGS
SPECHT = CSOIL
1100 CONTINUE
XSC = AMAX1(0., CAPAC(IVEG) - SATCAP(IVEG) )
IF(IVEG.EQ.2 .AND. TS.LE.TF )GO TO 1170
CAPAC(IVEG) = CAPAC(IVEG) - XSC
ROFF = ROFF + XSC
RNOFFS = XSC*1000. + RNOFFS
1170 CONTINUE
CAPACP(IVEG) = 0.
SNOWP(IVEG) = 0.
IF( TS .GT. TF ) CAPACP(IVEG) = CAPAC(IVEG)
IF( TS .LE. TF ) SNOWP(IVEG) = CAPAC(IVEG)
CAPAC(IVEG) = CAPACP(IVEG)
SNOWW(IVEG) = SNOWP(IVEG)
ZLOAD = CAPAC(IVEG) + SNOWW(IVEG)
FPI = ( 1.-EXP( - EXTK(IVEG,3,1) * ZLT(IVEG)/VCOVER(IVEG) ) ) &
* VCOVER(IVEG)
TTI = P0 * ( 1.-FPI )
XS = 1.
IF ( P0 .LT. 1.E-9 ) GO TO 1150
ARG = ( SATCAP(IVEG)-ZLOAD )/( P0*FPI*AP ) -CP/AP
IF ( ARG .LT. 1.E-9 ) GO TO 1150
XS = -1./BP * ALOG( ARG )
XS = AMIN1( XS, 1. )
XS = AMAX1( XS, 0. )
1150 TEX = P0*FPI * ( AP/BP*( 1.- EXP( -BP*XS )) + CP*XS ) - &
( SATCAP(IVEG) - ZLOAD ) * XS
TEX = AMAX1( TEX, 0. )
THRU = TTI + TEX
IF(IVEG.EQ.2.AND.TGS.LE.TF)THRU = 0.
PINF = P0 - THRU
IF( TM .GT. TF ) CAPAC(IVEG) = CAPAC(IVEG) + PINF
IF( TM .LE. TF ) SNOWW(IVEG) = SNOWW(IVEG) + PINF
IF( IVEG .EQ. 1 ) GO TO 1300
IF( TM .GT. TF ) GO TO 1200
SNOWW(IVEG) = SNOWP(IVEG) + P0
THRU = 0.
GO TO 1300
1200 EQUDEP = SATCO * DTT
XS = 1.
IF ( THRU .LT. 1.E-9 ) GO TO 1250
ARG = EQUDEP / ( THRU * AP ) -CP/AP
IF ( ARG .LT. 1.E-9 ) GO TO 1250
XS = -1./BP * ALOG( ARG )
XS = AMIN1( XS, 1. )
XS = AMAX1( XS, 0. )
1250 ROFFO = THRU * ( AP/BP * ( 1.-EXP( -BP*XS )) + CP*XS ) &
-EQUDEP*XS
ROFFO = AMAX1 ( ROFFO, 0. )
ROFF = ROFF + ROFFO
RNOFFS = RNOFFS + ROFFO*1000.
FILTR = FILTR + (THRU - ROFFO)
WWW(1) = WWW(1) + (THRU - ROFFO) / ( POROS*ZDEPTH(1) )
1300 CONTINUE
DIFF = ( CAPAC(IVEG)+SNOWW(IVEG) - CAPACP(IVEG)-SNOWP(IVEG) )*CW
CCP = SPECHT + SPWET1
CCT = SPECHT + SPWET1 + DIFF
TSD = ( TS * CCP + TM * DIFF ) / CCT
FREEZE = 0.
IF ( TS .GT. TF .AND. TM .GT. TF ) GO TO 2000
IF ( TS .LE. TF .AND. TM .LE. TF ) GO TO 2000
TTA = TS
TTB = TM
CCA = CCP
CCB = DIFF
IF ( TSD .GT. TF ) GO TO 2100
CCC = CAPACP(IVEG) * SNOMEL
IF ( TS .LT. TM ) CCC = DIFF * SNOMEL / CW
TSD = ( TTA * CCA + TTB * CCB + CCC ) / CCT
FREEZE = ( TF * CCT - ( TTA * CCA + TTB * CCB ) )
FREEZE = (AMIN1 ( CCC, FREEZE )) / SNOMEL
IF(TSD .GT. TF)TSD = TF - 0.1
GO TO 2000
2100 CONTINUE
CCC = - SNOWW(IVEG) * SNOMEL
IF ( TS .GT. TM ) CCC = - DIFF * SNOMEL / CW
TSD = ( TTA * CCA + TTB * CCB + CCC ) / CCT
FREEZE = ( TF * CCT - ( TTA * CCA + TTB * CCB ) )
FREEZE = (AMAX1( CCC, FREEZE )) / SNOMEL
IF(TSD .LE. TF)TSD = TF - 0.1
2000 CONTINUE
SMELT = FREEZE
SNOWW(IVEG) = SNOWW(IVEG) + FREEZE
CAPAC(IVEG) = CAPAC(IVEG) - FREEZE
IF( IVEG .EQ. 1 ) TC = TSD
IF( IVEG .EQ. 2 ) TGS = TSD
IF( SNOWW(IVEG) .LT. 0.0000001 ) GO TO 3000
ZMELT = 0.
ZMELT = CAPAC(IVEG)
CAPAC(IVEG) = 0.
WWW(1) = WWW(1) + ZMELT / ( POROS * ZDEPTH(1) )
FILTR = FILTR + ZMELT
3000 CONTINUE
CAPAC(IVEG) = CAPAC(IVEG) + SNOWW(IVEG)
SNOWW(IVEG) = 0.
freeze=0.0
P0 = THRU
1000 CONTINUE
CCX = ZLT(1) * CLAI + CAPAC(1) * CW
SPWET = AMIN1 ( 0.05, CAPAC(2))*CW
CG = (CSOIL + SPWET)
END SUBROUTINE INTERC
SUBROUTINE LAYER1 (CSOIL,TGS,DZSOIL,H,W,SNOWDEPTH,SWE,STEMP,ND)
parameter (dice=920.0, rhowater=1000.0,dlm=3.335d5)
dimension h(nd),w(nd)
swe=w(1)+w(2)+w(3)
snh=h(1)+h(2)+h(3)+csoil*(tgs-273.16)
dmlto=swe*dlm*rhowater
scv=1.9e+6*(swe/snowdepth)/dice
if (snh.gt.0.0) then
stemp=snh/(swe*4.18*10**6.+csoil)+273.16
else if (snh.gt.-dmlto) then
stemp=273.16
else
stemp=(snh+dmlto)/(scv*snowdepth+csoil)+273.16
end if
END SUBROUTINE LAYER1
SUBROUTINE LAYERN (TG,SNOW_WE,SNOW_DEPTH, DZ0,BW0,W0,BT0,CT0, &
FL0,FI0,H0,BL0,BI0,DLIQV0,DICEV0,TSSN0,DMLT0)
DIMENSION DZ0(4),W0(4),BW0(4),BT0(4),CT0(4),FL0(4),FI0(4),H0(4), &
BL0(4),BI0(4),DLIQV0(4),DICEV0(4),TSSN0(4),DMLT0(4)
IF(SNOW_DEPTH.GT.0.05.AND.SNOW_DEPTH.LE.0.06) THEN
DZ0(1)=0.02
DZ0(2)=0.02
DZ0(3)=SNOW_DEPTH- DZ0(1)- DZ0(2)
ELSE IF ( SNOW_DEPTH.GT.0.06.AND.SNOW_DEPTH.LE.0.08) THEN
DZ0(3)=0.02
DZ0(2)=0.02
DZ0(1)=SNOW_DEPTH- DZ0(3)- DZ0(2)
ELSE IF ( SNOW_DEPTH.GT.0.08.AND.SNOW_DEPTH.LE.0.62) THEN
DZ0(3)=0.02
DZ0(2)=(SNOW_DEPTH- DZ0(3))*0.33333333
DZ0(1)=(SNOW_DEPTH- DZ0(3))*0.66666667
ELSE IF ( SNOW_DEPTH.GT.0.62) THEN
DZ0(3)=0.02
DZ0(2)=0.20
DZ0(1)=SNOW_DEPTH- DZ0(3)- DZ0(2)
End IF
do 777 i=1,N
TSSN0(I)=TG
BW0(I)=SNOW_WE*RHOWATER/SNOW_DEPTH
777 continue
do 666 i=1,N
W0(I)=(BW0(I)*DZ0(I))/RHOWATER
BT0(I)=BW0(I)
CT0(I)=(BW0(I)/920.0)*1.9e+6
IF (TSSN0(I).EQ.273.16)THEN
FL0(I)= FLMIN
FI0(I)=1.0- FLMIN
H0(I)=(-1.0)*W0(I)*FI0(I)*DLM*RHOWATER
BL0(I)=BW0(I)*FL0(I)
BI0(I)=BW0(I)*FI0(I)
DLIQV0(I) = BL0(I)/RHOWATER
DICEV0(I) = BI0(I)/DICE
ELSE IF(TSSN0(I).LT.273.16) THEN
FL0(I)=0.0
FI0(I)=1.0
DMLT0(I)=W0(I)*DLM*RHOWATER
H0(I)=(TSSN0(I)-273.16)*CT0(I)*DZ0(I)-DMLT0(I)
BL0(I)=0.0
BI0(I)=BW0(I)
DLIQV0(I)=0.0
DICEV0(I) = BI0(I)/DICE
ENDIF
666 continue
END SUBROUTINE LAYERN
SUBROUTINE MODNODE(SNOWDEPTH,DZO,WO,HO,TSSNO,BWO,BIO, &
BLO,BTO,FIO,FLO,CTO,DLIQVOL,DICEVOL)
DIMENSION DZO(N1),WO(N1),HO(N1),TSSNO(N1),BWO(N1),BIO(N1),BLO(N1), &
BTO(N1),FIO(N1),FLO(N1),CTO(N1),DLIQVOL(N1),DICEVOL(N1)
IF (SNOWDEPTH.le.0.06) then
DZ1=0.02
DZ2=0.02
DZ3=SNOWDEPTH-( DZ2+DZ1)
ELSE IF (SNOWDEPTH.gt.0.06) then
DZ3=0.02
ENDIF
DDZ3=DZ3-dzo(3)
IF (DDZ3.GT.0.0) THEN
DDZ3=MIN(DDZ3,dzo(2))
CALL COMBO (DDZ3,dzo(3),dzo(2),wo(3),wo(2),ho(3),ho(2), &
tssno(3),tssno(2),bwo(3),bwo(2),bio(3),bio(2),blo(3),blo(2), &
bto(3),bto(2),fio(3),fio(2),flo(3),flo(2),cto(3),cto(2), &
dliqvol(3),dliqvol(2),dicevol(3),dicevol(2))
ELSE
DDZ3=-DDZ3
CALL COMBO (DDZ3,dzo(2),dzo(3),wo(2),wo(3),ho(2),ho(3), &
tssno(2),tssno(3),bwo(2),bwo(3),bio(2),bio(3),blo(2),blo(3), &
bto(2),bto(3),fio(2),fio(3),flo(2),flo(3),cto(2),cto(3), &
dliqvol(2),dliqvol(3),dicevol(2),dicevol(3))
END IF
SUM12=dzo(1)+dzo(2)
IF (SNOWDEPTH.le.0.06) THEN
DZ2=0.5*SUM12
ELSE IF (SNOWDEPTH.gt.0.06.and.SNOWDEPTH.le.0.08) THEN
DZ2=0.02
ELSE IF (SNOWDEPTH.gt.0.08.and.SNOWDEPTH.le.0.62) THEN
DZ2=0.33333333*SUM12
ELSE IF (SNOWDEPTH.gt.0.62) THEN
DZ2=0.20
ENDIF
DDZ2=DZ2-dzo(2)
IF (DDZ2.GT.0.0) THEN
CALL COMBO (DDZ2,dzo(2),dzo(1),wo(2),wo(1),ho(2),ho(1), &
tssno(2),tssno(1),bwo(2),bwo(1),bio(2),bio(1),blo(2),blo(1), &
bto(2),bto(1),fio(2),fio(1),flo(2),flo(1),cto(2),cto(1), &
dliqvol(2),dliqvol(1),dicevol(2),dicevol(1))
ELSE
DDZ2=-DDZ2
CALL COMBO (DDZ2,dzo(1),dzo(2),wo(1),wo(2),ho(1),ho(2), &
tssno(1),tssno(2),bwo(1),bwo(2),bio(1),bio(2),blo(1),blo(2), &
bto(1),bto(2),fio(1),fio(2),flo(1),flo(2),cto(1),cto(2), &
dliqvol(1),dliqvol(2),dicevol(1),dicevol(2))
END IF
SNOWDEPTH=dzo(1)+dzo(2)+dzo(3)
END SUBROUTINE MODNODE
SUBROUTINE NEWSNOW(PRCP,BIFALL,BLFALL,FLFALL,TKAIR, &
DZO,WO,BWO,CTO,HO,DMLTO,FIO,FLO,BIO,BLO,DLIQVOL,DICEVOL,TSSNO,WF)
dzfall=prcp*rhowater/bifall
dzo=dzo+dzfall
wo=wo+prcp
bwo=(wo*rhowater)/dzo
cto=1.9e+6*(bwo/920.0)
dum=(tkair-273.16)*cto*dzfall &
-(1.0-flfall)*(blfall+bifall)*dlm*dzfall
ho=ho+dum
dmlto=wo*rhowater*dlm
if (ho.ge.-dmlto) then
tssno=273.16
fio=-ho/dmlto
flo=1.0-fio
blo=bwo*flo
bio=bwo*fio
dliqvol=blo/rhowater
dicevol=bio/dice
else
fio=1.0
flo=0.0
bio=bwo
blo=0.0
dliqvol=0.0
dicevol=bio/dice
wf=0.0
tssno=(ho+dmlto)/(cto*dzo)+273.16
end if
END SUBROUTINE NEWSNOW
SUBROUTINE NEWTON(A1,Y,FINC,NOX,NONPOS,IWOLK,L,ZINC,A2,Y1,ITER)
REAL, DIMENSION (3) :: IWALK, NEX, ITER
REAL, DIMENSION (3) :: ZINC, A2, Y1
DATA CONS/1.0/
ERTOL = 0.05 * FINC
IWALK(L) = IWOLK
NEX(L)=NOX
IF ( ITER(L) .GE. 490 ) GO TO 160
IF (ERTOL .LT. 0.00000001) ERTOL=0.000001
IF (ABS(Y) .LE. ERTOL) GO TO 150
IF((ABS(Y-Y1(L))).LE.0.01*ERTOL .AND. IWALK(L).EQ.0 ) GO TO 8
IF(ABS(Y1(L)).GT.ERTOL) GO TO 1
A2(L)=A1
A1=A1-Y
NEX(L)=0
Y1(L)=Y
ITER(L)=1
IF (IWALK(L) .EQ. 3) GO TO 101
IWALK(L)=0
GO TO 101
1 ITER(L)=ITER(L)+1
IF(ITER(L) .EQ. 10) IWALK(L)=1
IF(IWALK(L) .NE. 0) GO TO 2
IF(ABS(Y) .GT. ERTOL) GO TO 3
NEX(L)=1
GO TO 150
3 A=A1-Y*(A1-A2(L))/(Y-Y1(L))
IF(ABS(A-A1).GT.(10.0*FINC)) &
A=A1+10.0*FINC*SIGN(CONS,(A-A1))
A2(L)=A1
A1=A
Y1(L)=Y
GO TO 101
2 IF(IWALK(L).EQ.2)GO TO 4
IF(IWALK(L).EQ.3) GO TO 6
IF(SIGN(CONS,Y).EQ.SIGN(CONS,Y1(L))) GO TO 3
ZINC(L)=(A1-A2(L))/4.0
A1=A2(L)+ZINC(L)
IWALK(L)=2
NEX(L)=0
GO TO 101
4 IF(SIGN(CONS,Y) .EQ.SIGN(CONS,Y1(L))) GO TO 5
ZINC(L)=-ZINC(L)/4.0
A2(L)=A1
A1=A1+ZINC(L)
NEX(L)=0
Y1(L)=Y
GO TO 101
5 A2(L)=A1
A1=A1+ZINC(L)
Y1(L)=Y
NEX(L)=0
GO TO 101
6 IF(SIGN(CONS,Y).EQ.SIGN(CONS,Y1(L))) GO TO 7
IWALK(L)=1
GO TO 2
7 A2(L) = A1
A1 = A1+FINC
Y1(L)=Y
NEX(L) = 0
GO TO 101
8 A1 = A1 + FINC*2.0
NEX(L)=0
GO TO 101
160 CONTINUE
900 FORMAT ( 3X,' FAILURE TO CONVERGE AFTER 490 ITERATIONS', &
/, 3X,' Y = ',2G12.5,2X,I14)
150 NEX(L) = 1
ZINC(L)=0.0
ITER(L) = 0
IWALK(L)=0
Y1(L)=0.0
Y=0.0
A2(L)=0.0
101 CONTINUE
IF(NONPOS.EQ.1.AND.A1.LT.0.0) A1=A2(L)/2.0
NOX = NEX(L)
IWOLK = IWALK(L)
END SUBROUTINE NEWTON
SUBROUTINE OLD(TSSN,BW,BL,BI,H,FL,FI,W,DZ,SS,CT,BT,DMLT, &
TSSNO,BWO,BLO,BIO,HO,FLO,FIO,WO,DZO,SSO,CTO,BTO,DMLTO)
DIMENSION TSSN(N1),BW(N1),BL(N1),BI(N1),H(N1),FL(N1),FI(N1), &
W(N1),DZ(N1),SS(N1),CT(N1),BT(N1),DMLT(N1), TSSNO(N1), &
BWO(N1),BLO(N1),BIO(N1),HO(N1),FLO(N1),FIO(N1), &
WO(N1),DZO(N1),SSO(N1),CTO(N1),BTO(N1),DMLTO(N1)
DO 20 I=1,N
TSSNO(I)=TSSN(I)
BWO(I)=BW(I)
BLO(I)=BL(I)
BIO(I)=BI(I)
HO(I)=H(I)
FLO(I)=FL(I)
FIO(I)=FI(I)
WO(I)=W(I)
DZO(I)=DZ(I)
SSO(I)=SS(I)
CTO(I)=CT(I)
BTO(I)=BT(I)
DMLTO(I)=DMLT(I)
20 CONTINUE
END SUBROUTINE OLD
SUBROUTINE RADAB (TRAN,REF,GREEN,VCOVER,CHIL,ZLAI,Z2,Z1,SOREF,TC, &
TGS,SATCAP,EXTK,RADFAC,THERMK,RADT,PAR,PD,ALBEDO,SALB, &
TGEFF,SUNANG,XADJ,CAPAC,RADN,ZLWUP,FRAC, &
ISNOW,SNOWDEN,SNOWDEPTH,SWDOWN,XALBEDO,SCOV2,ISICE, &
fsdown,fldown,fsup,flup)
DIMENSION TRANC1(2), TRANC2(2), TRANC3(2)
DIMENSION CAPAC(2), SATCAP(2), TRAN(2,3,2), REF(2,3,2), SOREF(3)
DIMENSION GREEN(2), VCOVER(2), ZLAI(2), CHIL(2), RADN(3,2),RADT(2)
DIMENSION RADFAC(2,2,2), RADSAV(12), PAR(2), PD(2), ALBEDO(2,3,2)
DIMENSION SALB(2,2), EXTK(2,3,2), FRAC(2,2)
DIMENSION sr(2)
data sr/0.85,0.65/
f=max(sunang,0.01746)
xref1=1.05
xref2=0.20
FMELT = 1.
IF ( ABS(TF-TGS) .LT. 0.5 ) FMELT = 0.6
SATCAP(1) = ZLAI(1) * 0.0001
SATCAP(2) = ZLAI(2) * 0.0001
IF (ISNOW.eq.0) THEN
DEPCOV = AMAX1( 0., (SNOWDEPTH-Z1))
ELSE
DEPCOV = AMAX1( 0., (CAPAC(2)*SNOWDEN-Z1))
END IF
DEPCOV = AMIN1( DEPCOV, (Z2-Z1)*0.95 )
SATCAP(1) = SATCAP(1) * ( 1. - DEPCOV / ( Z2 - Z1 ) )
do 202 iveg = 1, 2
do 202 iwave = 1, 3
do 202 irad = 1, 2
albedo(iveg,iwave,irad)=0.
202 continue
DO 1000 IWAVE = 1, 2
DO 2000 IVEG = 2, 1,-1
SCOV = 0.
IF( IVEG .EQ. 2 ) GO TO 100
IF( TC .LE. TF ) SCOV = AMIN1( 0.5, CAPAC(1) / SATCAP(1) )
100 CONTINUE
REFF1 = ( 1. - SCOV ) * REF(IVEG,IWAVE,1) + SCOV * ( xref1 - &
IWAVE * xref2 ) * FMELT
REFF2 = ( 1. - SCOV ) * REF(IVEG,IWAVE,2) + SCOV * ( xref1 - &
IWAVE * xref2 ) * FMELT
TRAN1 = TRAN(IVEG,IWAVE,1) * ( 1. - SCOV ) &
+ SCOV * ( 1.- ( xref1 - IWAVE * xref2 ) * FMELT ) &
* TRAN(IVEG,IWAVE,1)
TRAN2 = TRAN(IVEG,IWAVE,2) * ( 1. - SCOV ) &
+ SCOV * ( 1.- ( xref1 - IWAVE * xref2 ) * FMELT ) * 0.9 &
* TRAN(IVEG,IWAVE,2)
SCAT = GREEN(IVEG)*( TRAN1 + REFF1 ) +( 1. - GREEN(IVEG) ) * &
( TRAN2 + REFF2)
CHIV = CHIL(IVEG)
IF ( ABS(CHIV) .LE. 0.01 ) CHIV = 0.01
AA = 0.5 - 0.633 * CHIV - 0.33 * CHIV * CHIV
BB = 0.877 * ( 1. - 2. * AA )
PROJ = AA + BB * F
EXTKB = ( AA + BB * F ) / F
ZMEW = 1. / BB * ( 1. - AA / BB * ALOG ( ( AA + BB ) / AA ) )
ACSS = SCAT / 2. * PROJ / ( PROJ + F * BB )
ACSS = ACSS * ( 1. - F * AA / ( PROJ + F * BB ) * ALOG ( ( PROJ &
+ F * BB + F * AA ) / ( F * AA ) ) )
EXTK( IVEG, IWAVE, 1 ) = PROJ / F * SQRT( 1.-SCAT )
EXTK( IVEG, IWAVE, 2 ) = 1. / ZMEW * SQRT( 1.-SCAT )
EXTK( IVEG, 3, 1 ) = AA + BB
EXTK( IVEG, 3, 2 ) = 1./ZMEW
UPSCAT = GREEN(IVEG) * TRAN1 + ( 1. - GREEN(IVEG) ) * TRAN2
UPSCAT = 0.5 * ( SCAT + ( SCAT - 2. * UPSCAT ) * &
(( 1. - CHIV ) / 2. ) ** 2 )
BETAO = ( 1. + ZMEW * EXTKB ) / ( SCAT * ZMEW * EXTKB ) * ACSS
BE = 1. - SCAT + UPSCAT
CE = UPSCAT
BOT = ( ZMEW * EXTKB ) ** 2 + ( CE**2 - BE**2 )
IF ( ABS(BOT) .GT. 1.E-10) GO TO 200
SCAT = SCAT* 0.98
BE = 1. - SCAT + UPSCAT
BOT = ( ZMEW * EXTKB ) ** 2 + ( CE**2 - BE**2 )
200 CONTINUE
DE = SCAT * ZMEW * EXTKB * BETAO
FE = SCAT * ZMEW * EXTKB * ( 1. - BETAO )
CCE = DE * BE - ZMEW * DE * EXTKB + CE * FE
FFE = BE * FE + ZMEW * FE * EXTKB + CE * DE
TORE = -CCE / BOT
SIGE = -FFE / BOT
PSI = SQRT(BE**2 - CE**2)/ZMEW
IF (ISNOW.eq.0) THEN
SDEP=SNOWDEPTH
ELSE
SDEP = CAPAC(2) *SNOWDEN
END IF
FAC = ( SDEP - Z1 ) / ( Z2 - Z1 )
FAC = AMAX1( 0., FAC )
FAC = AMIN1( 0.99, FAC )
ZAT = ZLAI(IVEG) / VCOVER(IVEG)
IF ( IVEG .EQ. 1 ) ZAT = ZAT * (1.-FAC)
POWER1 = AMIN1( PSI*ZAT, 50. )
POWER2 = AMIN1( EXTKB*ZAT, 50. )
EPSI = EXP( - POWER1 )
EK = EXP ( - POWER2 )
ROSB = SOREF(IWAVE)
ROSD = SOREF(IWAVE)
IF ( IVEG .EQ. 2 ) GO TO 300
ROSB = ALBEDO(2,IWAVE,1)
ROSD = ALBEDO(2,IWAVE,2)
300 CONTINUE
GE = ROSB / ROSD
F1 = BE - CE / ROSD
ZP = ZMEW * PSI
DEN = ( BE + ZP ) * ( F1 - ZP ) / EPSI - &
( BE - ZP ) * ( F1 + ZP ) * EPSI
ALPHA = CE * ( F1 - ZP ) / EPSI / DEN
BETA = -CE * ( F1 + ZP ) * EPSI / DEN
F1 = BE - CE * ROSD
DEN = ( F1 + ZP ) / EPSI - ( F1 - ZP ) * EPSI
GAMMA = ( F1 + ZP ) / EPSI / DEN
DELTA = - ( F1 - ZP ) * EPSI / DEN
ALBEDO(IVEG,IWAVE,2) = ALPHA + BETA
IF ( IVEG .EQ. 1 ) GO TO 400
SCOV2 = 0.
IF ( TGS .LE. TF ) SCOV2 = AMIN1( 1., CAPAC(2) / 0.004 )
IF (ISICE.EQ.1) SCOV2=1.
ALBEDO(2,IWAVE,2) = &
ROSD * ( 1. - VCOVER(2) ) + ALBEDO(2,IWAVE,2) * VCOVER(2)
ALBEDO(2,IWAVE,2) = &
( 1. - SCOV2 ) * ALBEDO(2,IWAVE,2) + SCOV2 * &
( xref1-IWAVE*xref2 ) * &
FMELT
400 CONTINUE
TRANC2(IWAVE) = GAMMA * EPSI + DELTA / EPSI
F1 = BE - CE / ROSD
ZMK = ZMEW * EXTKB
DEN = ( BE + ZP ) * ( F1 - ZP ) / EPSI - &
( BE - ZP ) * ( F1 + ZP ) * EPSI
ALPHA = ( DE - TORE * ( BE + ZMK ) ) * ( F1 - ZP ) / EPSI - &
( BE - ZP ) * ( DE - CE*GE - TORE * ( F1 + ZMK ) ) * EK
ALPHA = ALPHA / DEN
BETA = ( BE + ZP ) * (DE - CE*GE - TORE * ( F1 + ZMK ))* EK - &
( DE - TORE * ( BE + ZMK ) ) * ( F1 + ZP ) * EPSI
BETA = BETA / DEN
F1 = BE - CE * ROSD
DEN = ( F1 + ZP ) / EPSI - ( F1 - ZP ) * EPSI
GAMMA = - SIGE * ( F1 + ZP ) / EPSI - &
( FE + CE * GE * ROSD + SIGE * ( ZMK - F1 ) ) * EK
GAMMA = GAMMA / DEN
DELTA = ( CE * GE * ROSD + FE + SIGE * ( ZMK - F1 ) ) * EK &
+ SIGE * ( F1 - ZP ) * EPSI
DELTA = DELTA / DEN
ALBEDO(IVEG,IWAVE,1) = TORE + ALPHA + BETA
IF( IVEG .EQ. 1 ) GO TO 500
ALBEDO(2,IWAVE,1) = ROSB * ( 1. - VCOVER(2) ) &
+ ALBEDO(2,IWAVE,1) * VCOVER(2)
ALBEDO(2,IWAVE,1) = ( 1. - SCOV2 ) * ALBEDO(2,IWAVE,1) + &
SCOV2 * ( xref1-IWAVE*xref2 ) * FMELT
500 CONTINUE
TRANC1(IWAVE) = EK
TRANC3(IWAVE) = SIGE * EK + GAMMA * EPSI + DELTA / EPSI
2000 CONTINUE
RADFAC(2,IWAVE,1) = ( 1.-VCOVER(1) ) * ( 1.-ALBEDO(2,IWAVE,1) ) &
+ VCOVER(1) * ( TRANC1(IWAVE) * ( 1.-ALBEDO(2,IWAVE,1) ) &
+ TRANC3(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) ) )
RADFAC(2,IWAVE,2) = ( 1.-VCOVER(1) ) * ( 1.-ALBEDO(2,IWAVE,2) ) &
+ VCOVER(1) * TRANC2(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) )
RADFAC(1,IWAVE,1) = VCOVER(1) * ( ( 1.-ALBEDO(1,IWAVE,1) ) &
- TRANC1(IWAVE) * ( 1.-ALBEDO(2,IWAVE,1) ) &
- TRANC3(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) ) )
RADFAC(1,IWAVE,2) = VCOVER(1) * ( ( 1.-ALBEDO(1,IWAVE,2) ) &
- TRANC2(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) ) )
DO 3000 IRAD = 1, 2
SALB(IWAVE,IRAD) = ( 1.-VCOVER(1) ) * ALBEDO(2,IWAVE,IRAD) + &
VCOVER(1) * ALBEDO(1,IWAVE,IRAD)
3000 CONTINUE
IF ( IWAVE .EQ. 2 ) GO TO 600
RADSAV(1) = 1. - VCOVER(1) &
+ VCOVER(1) * ( TRANC1(IWAVE) + TRANC3(IWAVE) )
RADSAV(2) = 1. - VCOVER(1) + VCOVER(1) * TRANC2(IWAVE)
600 CONTINUE
1000 CONTINUE
if (xadj.eq.0.) go to 730
xx = radfac(1,1,2) + radsav(2)
xy = radfac(1,1,1) + radsav(1)
ssum = salb(1,1)*frac(1,1) + salb(1,2)*frac(1,2)+ &
salb(2,1)*frac(2,1) + salb(2,2)*frac(2,2)
do 650 iwave = 1, 2
salb(iwave,2) = salb(iwave,2) + xadj * salb(iwave,2) / ssum
x0 = 1. - salb(iwave,2)
x1 = radfac(1,iwave,2) + radfac(2,iwave,2)
x2 = radfac(1,iwave,2) / x1
x3 = radfac(2,iwave,2) / x1
radfac(1,iwave,2) = x0 * x2
radfac(2,iwave,2) = x0 * x3
650 continue
640 format(1x,'unrealistic value, dif',2i12,4e11.4)
do 750 iwave = 1, 2
salb(iwave,1) = salb(iwave,1) + xadj * salb(iwave,1) / ssum
x0 = 1. - salb(iwave,1)
x1 = radfac(1,iwave,1) + radfac(2,iwave,1)
x2 = radfac(1,iwave,1) / x1
x3 = radfac(2,iwave,1) / x1
radfac(1,iwave,1) = x0 * x2
radfac(2,iwave,1) = x0 * x3
radsav(1) = xy - radfac(1,1,1)
radsav(2) = xx - radfac(1,1,2)
750 continue
740 format(1x,'unrealistic value',2i12,4e11.4)
730 continue
swup = radn(1,1)*salb(1,1) + radn(1,2)*salb(1,2) &
+ radn(2,1)*salb(2,1) + radn(2,2)*salb(2,2)
if ((swdown.gt.0.01).and.(swup.gt.0.01)) then
xalbedo = swup / swdown
if (xalbedo.gt.1.) then
swup = 0.
xalbedo = 999.
write (6, *) 'albebo incorrect',xalbedo
endif
else
swup = 0.0
xalbedo = .1
endif
TC4 = TC * TC * TC * TC
TG4 = TGS * TGS * TGS * TGS
ZKAT = EXTK(1,3,2) * ZLAI(1) / VCOVER(1)
ZKAT = AMIN1( 50. , ZKAT )
ZKAT = AMAX1( 1.E-5, ZKAT )
THERMK = EXP(-ZKAT)
FAC1 = VCOVER(1) * ( 1.-THERMK )
FAC2 = 1.
CLOSS = 2. * FAC1 * STEFAN * TC4
CLOSS = CLOSS - FAC2 * FAC1 * STEFAN * TG4
GLOSS = FAC2 * STEFAN * TG4
GLOSS = GLOSS - FAC1 * FAC2 * STEFAN * TC4
ZLWUP = FAC1 * STEFAN * TC4 + (1. - FAC1 ) * FAC2 * STEFAN * TG4
TGEFF = SQRT( SQRT ( ( ZLWUP / STEFAN ) ) )
RADSAV(3) = EXTK(1,1,1)
RADSAV(4) = EXTK(1,1,2)
RADSAV(5) = EXTK(2,1,1)
RADSAV(6) = EXTK(2,1,2)
RADSAV(7) = THERMK
RADSAV(8) = EXTK(1,3,1)
RADSAV(9) = EXTK(2,3,1)
RADSAV(10)= CLOSS
RADSAV(11)= GLOSS
RADSAV(12)= TGEFF
P1F = RADSAV(1)
P2F = RADSAV(2)
RADT(1) = 0.
RADT(2) = 0.
DO 7000 IVEG = 1, 2
DO 7000 IWAVE = 1, 2
DO 7000 IRAD = 1, 2
RADT(IVEG) = RADT(IVEG)+RADFAC(IVEG,IWAVE,IRAD)*RADN(IWAVE,IRAD)
7000 CONTINUE
fsdown = radn(1,1)+radn(1,2)+radn(2,1)+radn(2,2)
fsup = fsdown-radt(1)-radt(2)
SWCAN=RADT(1)
SWGND=RADT(2)
RADT(1) = RADT(1) + RADN(3,2)*VCOVER(1)*(1.- THERMK) &
- CLOSS
RADT(2) = RADT(2) + RADN(3,2)*( 1.-VCOVER(1)*(1-THERMK) ) &
- GLOSS
fldown = radn(3,2)
flup = closs+gloss
PAR(1) = RADN(1,1) + RADN(1,2) + 0.001
PD(1) = ( RADN(1,1) + 0.001 ) / PAR(1)
P1 = P1F * RADN(1,1) + 0.001
P2 = P2F * RADN(1,2)
PAR(2) = P1 + P2
PD(2) = P1 / PAR(2)
END SUBROUTINE RADAB
SUBROUTINE RADAB_ICE(TRAN,REF,GREEN,VCOVER,CHIL,ZLT,Z2,Z1,SOREF, &
TC,TGS,SATCAP,EXTK,CLOSS,GLOSS,THERMK,P1F,P2F, &
RADT,PAR,PD,SALB,ALBEDO,TGEFF,SUNANG,XADJ,CAPAC, &
RADN,bedo,ZLWUP,RADFRAC,SWDOWN,SCOV2,ISICE, &
fsdown,fldown,fsup,flup)
REAL, DIMENSION (2) :: TRANC1, TRANC2, TRANC3, CAPAC, SATCAP, &
GREEN, VCOVER, ZLT, CHIL, RADT, PAR, PD
REAL, DIMENSION (3) :: SOREF
REAL, DIMENSION (2,2) :: RADFRAC, SALB
REAL, DIMENSION (3,2) :: RADN
REAL, DIMENSION (2,2,2) :: RADFAC
REAL, DIMENSION (2,3,2) :: TRAN, REF, ALBEDO, EXTK
REAL, DIMENSION (12) :: RADSAV
f=max(sunang,0.01746)
FMELT = 1.
IF ( ABS(TF-TGS) .LT. 0.5 ) FMELT = 0.6
SATCAP(1) = ZLT(1) * 0.0001
SATCAP(2) = ZLT(2) * 0.0001
DEPCOV = AMAX1( 0., (CAPAC(2)*5.-Z1) )
DEPCOV = AMIN1( DEPCOV, (Z2-Z1)*0.95 )
SATCAP(1) = SATCAP(1) * ( 1. - DEPCOV / ( Z2 - Z1 ) )
do 202 iveg = 1, 2
do 202 iwave = 1, 3
do 202 irad = 1, 2
albedo(iveg,iwave,irad)=0.
202 continue
DO 1000 IWAVE = 1,2
DO 2000 IVDUM = 1,2
IF ( IVDUM .EQ. 1 ) IVEG = 2
IF ( IVDUM .EQ. 2 ) IVEG = 1
SCOV = 0.
IF( IVEG .EQ. 2 ) GO TO 100
IF( TC .LE. TF ) SCOV = AMIN1( 0.5, CAPAC(1) / SATCAP(1) )
100 CONTINUE
REFF1 = ( 1. - SCOV ) * REF(IVEG,IWAVE,1) + SCOV * ( 1.2 - &
IWAVE * 0.4 ) * FMELT
REFF2 = ( 1. - SCOV ) * REF(IVEG,IWAVE,2) + SCOV * ( 1.2 - &
IWAVE * 0.4 ) * FMELT
TRAN1 = TRAN(IVEG,IWAVE,1) * ( 1. - SCOV ) &
+ SCOV * ( 1.- ( 1.2 - IWAVE * 0.4 ) * FMELT ) &
* TRAN(IVEG,IWAVE,1)
TRAN2 = TRAN(IVEG,IWAVE,2) * ( 1. - SCOV ) &
+ SCOV * ( 1.- ( 1.2 - IWAVE * 0.4 ) * FMELT ) * 0.9 &
* TRAN(IVEG,IWAVE,2)
SCAT = GREEN(IVEG)*( TRAN1 + REFF1 ) +( 1. - GREEN(IVEG) ) * &
( TRAN2 + REFF2)
CHIV = CHIL(IVEG)
IF ( ABS(CHIV) .LE. 0.01 ) CHIV = 0.01
AA = 0.5 - 0.633 * CHIV - 0.33 * CHIV * CHIV
BB = 0.877 * ( 1. - 2. * AA )
PROJ = AA + BB * F
EXTKB = ( AA + BB * F ) / F
ZMEW = 1. / BB * ( 1. - AA / BB * ALOG ( ( AA + BB ) / AA ) )
ACSS = SCAT / 2. * PROJ / ( PROJ + F * BB )
ACSS = ACSS * ( 1. - F * AA / ( PROJ + F * BB ) * ALOG ( ( PROJ &
+ F * BB + F * AA ) / ( F * AA ) ) )
EXTK( IVEG, IWAVE, 1 ) = PROJ / F * SQRT( 1.-SCAT )
EXTK( IVEG, IWAVE, 2 ) = 1. / ZMEW * SQRT( 1.-SCAT )
EXTK( IVEG, 3, 1 ) = AA + BB
EXTK( IVEG, 3, 2 ) = 1./ZMEW
UPSCAT = GREEN(IVEG) * TRAN1 + ( 1. - GREEN(IVEG) ) * TRAN2
UPSCAT = 0.5 * ( SCAT + ( SCAT - 2. * UPSCAT ) * &
(( 1. - CHIV ) / 2. ) ** 2 )
BETAO = ( 1. + ZMEW * EXTKB ) / ( SCAT * ZMEW * EXTKB ) * ACSS
BE = 1. - SCAT + UPSCAT
CE = UPSCAT
BOT = ( ZMEW * EXTKB ) ** 2 + ( CE**2 - BE**2 )
IF ( ABS(BOT) .GT. 1.E-10) GO TO 200
SCAT = SCAT* 0.98
BE = 1. - SCAT + UPSCAT
BOT = ( ZMEW * EXTKB ) ** 2 + ( CE**2 - BE**2 )
200 CONTINUE
DE = SCAT * ZMEW * EXTKB * BETAO
FE = SCAT * ZMEW * EXTKB * ( 1. - BETAO )
CCE = DE * BE - ZMEW * DE * EXTKB + CE * FE
FFE = BE * FE + ZMEW * FE * EXTKB + CE * DE
TORE = -CCE / BOT
SIGE = -FFE / BOT
PSI = SQRT(BE**2 - CE**2)/ZMEW
SDEP = CAPAC(2) * 5.
FAC = ( SDEP - Z1 ) / ( Z2 - Z1 )
FAC = AMAX1( 0., FAC )
FAC = AMIN1( 0.99, FAC )
ZAT = ZLT(IVEG) / VCOVER(IVEG)
IF ( IVEG .EQ. 1 ) ZAT = ZAT * (1.-FAC)
POWER1 = AMIN1( PSI*ZAT, 50. )
POWER2 = AMIN1( EXTKB*ZAT, 50. )
EPSI = EXP( - POWER1 )
EK = EXP ( - POWER2 )
ROSB = SOREF(IWAVE)
ROSD = SOREF(IWAVE)
IF ( IVEG .EQ. 2 ) GO TO 300
ROSB = ALBEDO(2,IWAVE,1)
ROSD = ALBEDO(2,IWAVE,2)
300 CONTINUE
GE = ROSB / ROSD
F1 = BE - CE / ROSD
ZP = ZMEW * PSI
DEN = ( BE + ZP ) * ( F1 - ZP ) / EPSI - &
( BE - ZP ) * ( F1 + ZP ) * EPSI
ALPHA = CE * ( F1 - ZP ) / EPSI / DEN
BETA = -CE * ( F1 + ZP ) * EPSI / DEN
F1 = BE - CE * ROSD
DEN = ( F1 + ZP ) / EPSI - ( F1 - ZP ) * EPSI
GAMMA = ( F1 + ZP ) / EPSI / DEN
DELTA = - ( F1 - ZP ) * EPSI / DEN
ALBEDO(IVEG,IWAVE,2) = ALPHA + BETA
IF ( IVEG .EQ. 1 ) GO TO 400
SCOV2 = 0.
IF (ISICE.EQ.1) SCOV2=1.
IF ( TGS .LE. TF ) SCOV2 = AMIN1( 1., CAPAC(2) / 0.004 )
ALBEDO(2,IWAVE,2)= ROSD * ( 1. - VCOVER(2) ) + ALBEDO(2,IWAVE,2) * VCOVER(2)
ALBEDO(2,IWAVE,2) = &
( 1. - SCOV2 ) * ALBEDO(2,IWAVE,2) + SCOV2 * &
( 1.2-IWAVE*0.4 ) * FMELT
400 CONTINUE
TRANC2(IWAVE) = GAMMA * EPSI + DELTA / EPSI
F1 = BE - CE / ROSD
ZMK = ZMEW * EXTKB
DEN = ( BE + ZP ) * ( F1 - ZP ) / EPSI - &
( BE - ZP ) * ( F1 + ZP ) * EPSI
ALPHA = ( DE - TORE * ( BE + ZMK ) ) * ( F1 - ZP ) / EPSI - &
( BE - ZP ) * ( DE - CE*GE - TORE * ( F1 + ZMK ) ) * EK
ALPHA = ALPHA / DEN
BETA = ( BE + ZP ) * (DE - CE*GE - TORE * ( F1 + ZMK ))* EK - &
( DE - TORE * ( BE + ZMK ) ) * ( F1 + ZP ) * EPSI
BETA = BETA / DEN
F1 = BE - CE * ROSD
DEN = ( F1 + ZP ) / EPSI - ( F1 - ZP ) * EPSI
GAMMA = - SIGE * ( F1 + ZP ) / EPSI - &
( FE + CE * GE * ROSD + SIGE * ( ZMK - F1 ) ) * EK
GAMMA = GAMMA / DEN
DELTA = ( CE * GE * ROSD + FE + SIGE * ( ZMK - F1 ) ) * EK &
+ SIGE * ( F1 - ZP ) * EPSI
DELTA = DELTA / DEN
ALBEDO(IVEG,IWAVE,1) = TORE + ALPHA + BETA
IF( IVEG .EQ. 1 ) GO TO 500
ALBEDO(2,IWAVE,1) = ROSB * ( 1. - VCOVER(2) ) &
+ ALBEDO(2,IWAVE,1) * VCOVER(2)
ALBEDO(2,IWAVE,1) = ( 1. - SCOV2 ) * ALBEDO(2,IWAVE,1) + &
SCOV2 * ( 1.2-IWAVE*0.4 ) * FMELT
500 CONTINUE
TRANC1(IWAVE) = EK
TRANC3(IWAVE) = SIGE * EK + GAMMA * EPSI + DELTA / EPSI
2000 CONTINUE
RADFAC(2,IWAVE,1) = ( 1.-VCOVER(1) ) * ( 1.-ALBEDO(2,IWAVE,1) ) &
+ VCOVER(1) * ( TRANC1(IWAVE) * ( 1.-ALBEDO(2,IWAVE,1) ) &
+ TRANC3(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) ) )
RADFAC(2,IWAVE,2) = ( 1.-VCOVER(1) ) * ( 1.-ALBEDO(2,IWAVE,2) ) &
+ VCOVER(1) * TRANC2(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) )
RADFAC(1,IWAVE,1) = VCOVER(1) * ( ( 1.-ALBEDO(1,IWAVE,1) ) &
- TRANC1(IWAVE) * ( 1.-ALBEDO(2,IWAVE,1) ) &
- TRANC3(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) ) )
RADFAC(1,IWAVE,2) = VCOVER(1) * ( ( 1.-ALBEDO(1,IWAVE,2) ) &
- TRANC2(IWAVE) * ( 1.-ALBEDO(2,IWAVE,2) ) )
DO 3000 IRAD = 1, 2
SALB(IWAVE,IRAD) = ( 1.-VCOVER(1) ) * ALBEDO(2,IWAVE,IRAD) + &
VCOVER(1) * ALBEDO(1,IWAVE,IRAD)
3000 CONTINUE
IF ( IWAVE .EQ. 2 ) GO TO 600
RADSAV(1) = 1. - VCOVER(1) &
+ VCOVER(1) * ( TRANC1(IWAVE) + TRANC3(IWAVE) )
RADSAV(2) = 1. - VCOVER(1) + VCOVER(1) * TRANC2(IWAVE)
600 CONTINUE
1000 CONTINUE
if (xadj.eq.0.) go to 730
xx = radfac(1,1,2) + radsav(2)
xy = radfac(1,1,1) + radsav(1)
ssum = salb(1,1)*radfrac(1,1) + salb(1,2)*radfrac(1,2)+ &
salb(2,1)*radfrac(2,1) + salb(2,2)*radfrac(2,2)
do 650 iwave = 1, 2
salb(iwave,2) = salb(iwave,2) + xadj * salb(iwave,2) / ssum
x0 = 1. - salb(iwave,2)
x1 = radfac(1,iwave,2) + radfac(2,iwave,2)
x2 = radfac(1,iwave,2) / x1
x3 = radfac(2,iwave,2) / x1
radfac(1,iwave,2) = x0 * x2
radfac(2,iwave,2) = x0 * x3
if (salb(iwave,2).gt.1..or.radfac(1,iwave,2).gt.1..or. &
radfac(2,iwave,2).gt.1..or.salb(iwave,2).lt.0..or. &
radfac(1,iwave,2).lt.0..or.radfac(2,iwave,2).lt.0.) then
stop 999
end if
650 continue
640 format(1x,'unrealistic value, dif',2i12,4e11.4)
do 750 iwave = 1, 2
salb(iwave,1) = salb(iwave,1) + xadj * salb(iwave,1) / ssum
x0 = 1. - salb(iwave,1)
x1 = radfac(1,iwave,1) + radfac(2,iwave,1)
x2 = radfac(1,iwave,1) / x1
x3 = radfac(2,iwave,1) / x1
radfac(1,iwave,1) = x0 * x2
radfac(2,iwave,1) = x0 * x3
radsav(1) = xy - radfac(1,1,1)
radsav(2) = xx - radfac(1,1,2)
if (salb(iwave,1).gt.1..or.radfac(1,iwave,1).gt.1..or. &
radfac(2,iwave,1).gt.1..or.salb(iwave,1).lt.0..or. &
radfac(1,iwave,1).lt.0..or.radfac(2,iwave,1).lt.0.) then
write(7,740) nymdh,iwave,salb(iwave,1),radfac(1,iwave,1), &
radfac(2,iwave,1)
stop 999
end if
750 continue
740 format(1x,'unrealistic value',2i12,4e11.4)
730 continue
sibsu = radn(1,1)*salb(1,1) + radn(1,2)*salb(1,2) &
+ radn(2,1)*salb(2,1) + radn(2,2)*salb(2,2)
if ((swdown.gt.0.01).and.(sibsu.gt.0.01)) then
bedo = sibsu / swdown
if (bedo.gt.1.) then
sibsu = 0.
bedo = .1
print*,'albebo incorrect',ii,jj,bedo,sibsu,swdown, &
radn(1,1),radn(1,2),radn(2,1),radn(2,2)
endif
else
sibsu = 0.0
bedo = .1
endif
TC4 = TC * TC * TC * TC
TG4 = TGS * TGS * TGS * TGS
ZKAT = EXTK(1,3,2) * ZLT(1) / VCOVER(1)
ZKAT = AMIN1( 50. , ZKAT )
ZKAT = AMAX1( 1.E-5, ZKAT )
THERMK = EXP(-ZKAT)
FAC1 = VCOVER(1) * ( 1.-THERMK )
FAC2 = 1.
CLOSS = 2. * FAC1 * STEFAN * TC4
CLOSS = CLOSS - FAC2 * FAC1 * STEFAN * TG4
GLOSS = FAC2 * STEFAN * TG4
GLOSS = GLOSS - FAC1 * FAC2 * STEFAN * TC4
ZLWUP = FAC1 * STEFAN * TC4 + (1. - FAC1 ) * FAC2 * STEFAN * TG4
TGEFF = SQRT( SQRT ( ( ZLWUP / STEFAN ) ) )
RADSAV(3) = EXTK(1,1,1)
RADSAV(4) = EXTK(1,1,2)
RADSAV(5) = EXTK(2,1,1)
RADSAV(6) = EXTK(2,1,2)
RADSAV(7) = THERMK
RADSAV(8) = EXTK(1,3,1)
RADSAV(9) = EXTK(2,3,1)
RADSAV(10)= CLOSS
RADSAV(11)= GLOSS
RADSAV(12)= TGEFF
P1F = RADSAV(1)
P2F = RADSAV(2)
EXTK(1,1,1) = RADSAV(3)
EXTK(1,1,2) = RADSAV(4)
EXTK(2,1,1) = RADSAV(5)
EXTK(2,1,2) = RADSAV(6)
THERMK = RADSAV(7)
EXTK(1,3,1) = RADSAV(8)
EXTK(2,3,1) = RADSAV(9)
CLOSS = RADSAV(10)
GLOSS = RADSAV(11)
TGEFF = RADSAV(12)
RADT(1) = 0.
RADT(2) = 0.
DO 7000 IVEG = 1, 2
DO 7000 IWAVE = 1, 2
DO 7000 IRAD = 1, 2
RADT(IVEG) = RADT(IVEG)+RADFAC(IVEG,IWAVE,IRAD)*RADN(IWAVE,IRAD)
7000 CONTINUE
fsdown = radn(1,1)+radn(1,2)+radn(2,1)+radn(2,2)
fsup = fsdown-radt(1)-radt(2)
SWCAN=RADT(1)
SWGND=RADT(2)
RADT(1) = RADT(1) + RADN(3,2)*VCOVER(1)*(1.- THERMK) &
- CLOSS
RADT(2) = RADT(2) + RADN(3,2)*( 1.-VCOVER(1)*(1-THERMK) ) &
- GLOSS
fldown = radn(3,2)
flup = closs+gloss
PAR(1) = RADN(1,1) + RADN(1,2) + 0.001
PD(1) = ( RADN(1,1) + 0.001 ) / PAR(1)
P1 = P1F * RADN(1,1) + 0.001
P2 = P2F * RADN(1,2)
PAR(2) = P1 + P2
PD(2) = P1 / PAR(2)
END SUBROUTINE RADAB_ICE
SUBROUTINE RASIT5(TRIB,CTNI,CUNI,RA,Z2,Z0,D,ZZWIND,UMM1, &
RHOA,TMM,U2,USTAR,DRAG,TA,bps,rib,CU,CT,iii,jjj)
FS(X) = 66.85 * X
FT(X) = 0.904 * X
FV(X) = 0.315 * X
G2= 0.75
G3= 0.75
Z22 = Z2
ZL = Z2 + 11.785 * Z0
ZWIND = ZZWIND
TM = TMM
UMM = UMM1
if(zwind.le.d.or.zl.le.d) d=min(zwind,zl)-0.1
Z2 = D + Z0
CUNI = ALOG((ZWIND-D)/Z0)/VKC
IF (ZL.LT.ZWIND) THEN
XCT1 = ALOG((ZWIND-D)/(ZL-D))
XCT2 = ALOG((ZL-D)/(Z2-D))
XCTU2 = ALOG((ZL-D)/(Z22-D))
CTNI = (XCT1 + G3 * XCT2) / VKC
ELSE
XCT2 = ALOG((ZWIND-D)/(Z2-D))
XCTU2 = ALOG((ZWIND-D)/(Z22-D))
CTNI = G3 * XCT2 /VKC
END IF
UM=AMAX1(UMM,2.)
USTARN=UM/CUNI
VENTN =RHOA /CTNI*USTARN
IF (ZL.LT.ZWIND) THEN
U2 = UM - 1. / VKC * USTARN * (XCT1 + G2 * XCTU2)
ELSE
U2 = UM - 1. / VKC * USTARN * G2 * XCTU2
END IF
if(u2.lt.0.01) u2=0.01
THM=TM*bps
THVGM=TRIB-THM
IF (TA.EQ.0.) THVGM = 0.
RIB = -THVGM*GRAV*(ZWIND-D) / (THM*(UM-U2)**2)
RIB = MAX(-10.E0,RIB)
RIB = MIN(0.1643E0,RIB)
IF(RIB.LT.0.0)THEN
GRIB = +RIB
GRZL = +RIB*(ZL-D)/(ZWIND-D)
GRZ2 = +RIB*(Z2-D)/(ZWIND-D)
FVV = FV(GRIB)
IF (ZL.LT.ZWIND) THEN
FTT = FT(GRIB) + (G3-1.) * FT(GRZL) - G3 * FT(GRZ2)
ELSE
FTT = G3*(FT(GRIB) - FT(GRZ2))
END IF
CUI = CUNI + FVV
CTI = CTNI + FTT
ELSE
RZL = RIB/(ZWIND-D)*(ZL-D)
RZ2 = RIB/(ZWIND-D)*(Z2-D)
FVV = FS(RIB)
IF (ZL.LT.ZWIND) THEN
FTT = FS(RIB) + (G3-1) * FS(RZL) - G3 * FS(RZ2)
ELSE
FTT = G3 * (FS(RIB) - FS(RZ2))
END IF
312 CUI = CUNI + FVV
CTI = CTNI + FTT
ENDIF
310 CONTINUE
CU=1./CUI
CT=1./CTI
USTAR =UM*CU
RAF = CTI / USTAR
IF (RAF.LT.0.80) RAF = 0.80
RA = RAF
UEST = USTAR
DRAG = RHOA * UEST*UEST
Z2 = Z22
END SUBROUTINE RASIT5
SUBROUTINE SDSOL(DSOL,DMASS,N,SOLAR,SOLSOIL)
parameter(nd = 4)
integer n
real dsol(nd),dmass(nd),fext(nd)
gsize = 5.d-4
bext = 400.0
cv = 3.795d-3
depth = 30
do i=1,n
fext(i) = 0.0
enddo
tmass = 0.0
do 10 i=1,n
j=n+1-i
tmass=tmass+dmass(j)
if(tmass.gt.depth) goto 30
fext(j)=exp(-cv*dmass(j)/sqrt(gsize))
if(j .eq. n) fext(n)=exp(-bext*2d-3)*fext(n)
10 continue
30 tsolt = solar
do 20 i=1,n
j=n+1-i
if(tsolt .le. 0d0)then
dsol(j)=0d0
tsolb=0.0
else
tsolb=tsolt*fext(j)
dsol(j)=tsolt-tsolb
tsolt=tsolb
end if
20 continue
solsoil = tsolb
END SUBROUTINE SDSOL
SUBROUTINE SET0(TSSNO,BWO,BLO,BIO,HO,FLO,FIO,WO,DZO, &
SSO,CTO,BTO,DMLTO,WF,DHP)
DIMENSION WF(N1),DHP(N1),TSSNO(N1),BWO(N1),BLO(N1),BIO(N1),HO(N1), &
FLO(N1),FIO(N1),WO(N1),DZO(N1),SSO(N1),CTO(N1),BTO(N1),DMLTO(N1)
DO 100 I=N+1,N1
TSSNO(I)=0.0
BWO(I)=0.0
BLO(I)=0.0
BIO(I)=0.0
HO(I)=0.0
FLO(I)=0.0
FIO(I)=0.0
WO(I)=0.0
DZO(I)=0.0
SSO(I)=0.0
CTO(I)=0.0
BTO(I)=0.0
DMLTO(I)=0.0
100 CONTINUE
DO 200 I=1,N1
WF(I)=0.0
DHP(I)=0.0
200 CONTINUE
END SUBROUTINE SET0
SUBROUTINE SNOW_1ST (DTT,TM,SOLAR,PRCPW,PRCPS,BIO,BLO,DICEVOL, &
DLIQVOL,TSSNO,PDZDTC,POROSITY,SO,SSO,WF,DHP,DZO,WO, &
BWO,BTO,CTO,DMASS,DSOL,SNROFF,HROFF,SNOWDEPTH,SOLSOIL, &
FLO,FIO,DMLTO,HO,BIFALL,BLFALL,FLFALL)
DIMENSION BIO(N1),BLO(N1),DICEVOL(N1),DLIQVOL(N1),TSSNO(N1), &
PDZDTC(N1),POROSITY(N1),SO(N1),SSO(N1),WF(N1),DHP(N1), &
DZO(N1),WO(N1),BWO(N1),BTO(N1),CTO(N1),DMASS(N1), &
DSOL(N1),FLO(N1),FIO(N1),DMLTO(N1),HO(N1)
tkair = TM
prcp = prcpw+prcps
snroff = 0.0
hroff = 0.0
if(prcpw.gt.0.0)then
wf(n+1)=amin1(prcpw, dksatsnow*dtt)
dhp(n+1)=(wf(n+1)/dtt)*cl*rhowater*(tkair-273.16)
snroff =snroff+(prcpw-wf(n+1))
hroff=hroff+(prcpw-wf(n+1))*cl*rhowater*(tkair-273.16)
else if(prcps.gt.0.0)then
wf(n+1)=0.0
dhp(n+1)=0.0
CALL NEWSNOW(PRCP,BIFALL,BLFALL,FLFALL,TKAIR, &
DZO(N),WO(N),BWO(N),CTO(N),HO(N),DMLTO(N),FIO(N),FLO(N), &
BIO(N),BLO(N),DLIQVOL(N),DICEVOL(N),TSSNO(N),WF(N))
endif
do 277 i=n,1,-1
dicevol(i) = bio(i)/dice
dliqvol(i) = blo(i)/rhowater
porosity(i)=1.0-dicevol(i)
porosity(i)=amin1(porosity(i),1.0)
porosity(i)=amax1(porosity(i),0.0)
so(i)=ssisnow
if(porosity(i).ne.0.0) so(i)=dliqvol(i)/porosity(i)
sso(i)=(so(i)-ssisnow)/(1.0-ssisnow)
277 continue
overburden=0.0
do 377 i=n,1,-1
overburden=overburden+ wo(i)*rhowater
call COMPACT(BIO(I),TSSNO(I),BLO(I),OVERBURDEN,PDZDTC(I), &
SSO(I),DICE)
377 continue
do 390 i = 1,n
if((sso(i).lt.1.0.and.porosity(i).gt.0.0))then
dzot=dzo(i)*(1d0+pdzdtc(i)*dtt)
dzo(i)=amax1(dzot,dzmin)
if(wo(i).gt.womin)then
bwo(i)=(wo(i)*rhowater)/dzo(i)
if (bwo(i).gt.920.0) then
bwo(i)=920.0
dzo(i)=(wo(i)*rhowater)/bwo(i)
end if
endif
blo(i)=bwo(i)*flo(i)
bio(i)=bwo(i)*fio(i)
bto(i)=bwo(i)
end if
dicevol(i) = bio(i)/dice
dliqvol(i) = blo(i)/rhowater
dummy = dliqvol(i) + dicevol(i)
if(dummy.gt.1.0)then
dliqvol(i) = 1.0 - dicevol(i)
blo(i) = dliqvol(i)*rhowater
bwo(i) = blo(i) + bio(i)
dzo(i)=(wo(i)*rhowater)/bwo(i)
endif
cto(i)=(bwo(i)/920.0)*1.9e+6
porosity(i)=1.0-dicevol(i)
if(porosity(i) .gt. 1.0)porosity(i)=1.0
if(porosity(i) .lt. 0.0)porosity(i)=0.0
if(porosity(i).gt.0.0)then
so(i)=dliqvol(i)/porosity(i)
else
so(i)=ssisnow
endif
if(so(i).gt.ssisnow)then
sso(i)=(so(i)-ssisnow)/(1.0-ssisnow)
else
sso(i)=0.0
endif
dmass(i)=bto(i)*dzo(i)
390 continue
SNOWDEPTH=dzo(1)+dzo(2)+dzo(3)
IF (solar .gt. 0d0 ) THEN
call sdsol(dsol,dmass,n,solar,solsoil)
ELSE
do 112 i=1,n
dsol(i)=0d0
112 continue
solsoil=0.0
END IF
END SUBROUTINE SNOW_1ST
SUBROUTINE SNRESULT(DTT,I,ICASE,WFSOIL,TSOIL,B1,B2,FFF,DELTH,WWW, &
ZDEPTH,POROS,BI,BIO,DZ,DZO,W,BW,BWO,H,HO,QK,BL,BT,CT, &
FI,FL,WF,TSSN,DLIQVOL,DICEVOL,SNROFF,HROFF,QSOIL)
DIMENSION BIO(N1),DZO(N1),W(N1),BWO(N1),HO(N1), &
DZ(N1),BI(N1),BW(N1),BL(N1),BT(N1),CT(N1),FI(N1),FL(N1), &
WF(N1),H(N1),TSSN(N1),DLIQVOL(N1),DICEVOL(N1),QK(N1), &
WWW(3),ZDEPTH(3)
DIMENSION DELTH(20)
DATA BWE/200.0/
hx=0.0
IF (ICASE.EQ.1) THEN
fi(i)=1.0
fl(i)=0.0
dz(i)=dzo(i)
bw(i)=(w(i)*rhowater)/dz(i)
if((w(i)/dz(i)).lt.0.05.or.(w(i)/dz(i)) &
.gt.(dice/1000.0))then
bw(i) = bwo(i)
dz(i) = (w(i)*rhowater)/bw(i)
endif
bi(i)=bw(i)
bl(i)=0.0
bt(i)=bw(i)
wf(i)=0.0
if (i.eq.1) wfsoil=0.0
dliqvol(i)=0.0
dicevol(i)=bi(i)/dice
ct(i)=(bw(i)/920.0)*1.9e+6
if (i.eq.n) then
h(i)=ct(i)*dz(i)*(tssn(i)-273.16)-rhowater*dlm*w(n)*fi(n)
else
tssn(i) = ( ho(i) + ct(i)*dz(i)*273.16 + b1*dtt &
+ rhowater*dlm*w(i) ) &
/ ( ct(i)*dz(i) - b2*dtt )
h(i) = ho(i) + (b1+b2*tssn(i))*dtt
end if
if(tssn(i).gt.273.16) then
WRITE( message,* ) 'Warning: Snow Temp above freezing',i,tssn(i)
tssn(i)=273.16
CALL wrf_message ( message )
endif
ELSE IF (ICASE.EQ.2) THEN
fl(i)=1.0-fi(i)
tssn(i)=273.16
wf(i)=0.0
If(bwo(i).ge.bwe) Then
if(fl(i).gt.flmin)then
wf(i) = w(i)-(fi(i)/(1.0-flmin))*w(i)
w(i) = (fi(i)/(1.0-flmin))*w(i)
dum = wf(i)
fl(i)=flmin
fi(i)=1.0-fl(i)
endif
Else
flm = flmin+(flmax-flmin)*((bwe-bwo(i))/bwe)
if(fl(i).gt.flm)then
wf(i) = w(i)-(fi(i)/(1.0-flm))*w(i)
w(i) = (fi(i)/(1.0-flm))*w(i)
dum = wf(i)
fl(i)=flm
fi(i)=1.0-fl(i)
endif
Endif
If( wf(i).gt.0.0) Then
if(i.ne.1)then
wf(i)=amin1(dum, dksatsnow*dtt)
snroff = snroff + (dum - wf(i))
hroff=hroff+(dum-wf(i))*cl*rhowater*(tssn(i)-273.16)
else
if(www(1).ge.1.0) then
snroff = snroff + wf(i)
wfsoil=0.0
else
slwet=www(1)*poros*zdepth(1)
www(1)=(slwet+wf(i))/(poros*zdepth(1))
if(www(1).gt.1.0) then
snroff = snroff + (www(1)-1.0)*poros*zdepth(1)
www(1)=1.0
endif
wfsoil=0.0
endif
hroff=hroff + wf(i)*cl*rhowater*(tssn(i)-273.16)
endif
Endif
xnodalmelt=bio(i)*dzo(i)-w(i)*rhowater*fi(i)
If(xnodalmelt.gt.0.0.and.bio(i)*dzo(i).gt.0.0 &
.and.(bio(i).lt.250.0.or.(i.eq.n.and. &
bio(i).lt.400.0))) Then
ddz3=-xnodalmelt/(bio(i)*dzo(i))
dz(i)=dzo(i)*(1.0+ddz3)
Else
dz(i)=dzo(i)
Endif
bw(i)=(w(i)*rhowater)/dz(i)
If((w(i)/dz(i)).lt.0.05.or.(w(i)/dz(i)) &
.gt.(dice/1000.0)) Then
bw(i) = bwo(i)
dz(i) = (w(i)*rhowater)/bw(i)
Endif
bi(i)=bw(i)*fi(i)
bl(i)=bw(i)-bi(i)
bt(i)=bw(i)
ct(i)=(bw(i)/920.0)*1.9e+6
dliqvol(i)=bl(i)/rhowater
dicevol(i)=bi(i)/dice
h(i)=(-1.0)*w(i)*fi(i)*dlm*rhowater
ELSE IF (ICASE.EQ.3) THEN
fl(i) = 1.0
fi(i) = 0.0
wf(i) = w(i)
dum= wf(i)
dz(i) = 10e-15
w(i) = 10e-15
bw(i) =rhowater
bl(i)=bw(i)
bi(i)=0.0
dliqvol(i) = 1.0
dicevol(i) = 0.0
ct(i)=(bw(i)/920.0)*1.9e+6
tssn(i) = 273.16
h(i) = 0.0
If (i.eq.n) Then
if (i.eq.1) then
hx=(-1.0)*w(i)*fff*dlm*rhowater/dtt
snroff=wf(1)+snroff
wfsoil=0.0
else
wf(i)=amin1(dum, dksatsnow*dtt)
snroff = snroff + (dum - wf(i))
hroff=hroff+(dum-wf(i))*cl*rhowater*(tssn(i)-273.16)
delth(n-1) = (-1.0)*w(i)*fff*dlm*rhowater/dtt
end if
Else
if(i.eq.1)then
hx = ho(i)/dtt + b1+b2*tssn(i)
if(www(1).ge.1.0) then
snroff = snroff + wf(i)
wfsoil=0.0
else
slwet=www(1)*poros*zdepth(1)
www(1)=(slwet+wf(i))/(poros*zdepth(1))
if(www(1).gt.1.0) then
snroff = snroff + (www(1)-1.0)*poros*zdepth(1)
www(1)=1.0
endif
wfsoil=0.0
endif
else
wf(i)=amin1(dum, dksatsnow*dtt)
snroff = snroff + (dum - wf(i))
hroff=hroff+(dum-wf(i))*cl*rhowater*(tssn(i)-273.16)
delth(i-1) = ho(i)/dtt + b1+b2*tssn(i)
endif
End if
END IF
if (i.eq.1) qsoil = qk(1)*(tssn(1) - tsoil) + hx
END SUBROUTINE SNRESULT
SUBROUTINE STRES1 (IFIRST,RSTM,ROOTP, &
RSTFAC,RST,TC,ETC,RB,TGS,ETGS,RD,TU,TL,TOPT,EA, &
DEFAC,PH1,PH2,NROOT,ZDEPTH,PHSOIL,ROOTD,VCOVER,DROP)
DIMENSION TOPT(2), TL(2), TU(2), DEFAC(2), VCOVER(2)
DIMENSION PH1(2), PH2(2), RST(2), RSTFAC(2,4),XDRR(3)
DIMENSION ROOTD(2),ROOTP(3),ZDEPTH(3),PHSOIL(3), RSTM(2), DEP(3)
DO 1000 IVEG = 1, 2
TV = TC
ETV = ETC
RAIR = RB * 2.
IF ( IVEG .EQ. 1 ) GO TO 100
TV = TGS
ETV = ETGS
RAIR = RD
100 CONTINUE
TV = AMIN1 ( ( TU(IVEG) - 0.1 ), TV )
TV = AMAX1 ( ( TL(IVEG) + 0.1 ), TV )
IF( IFIRST .EQ. 0 ) GO TO 200
RSTM(IVEG) = RST(IVEG)
D2 = ( TU(IVEG) - TOPT(IVEG) ) / ( TOPT(IVEG) - TL(IVEG) )
D1 = 1. /(( TOPT(IVEG) - TL(IVEG) )* &
EXP( ALOG( TU(IVEG) - TOPT(IVEG))*D2))
RSTFAC(IVEG,3) = D1*( TV-TL(IVEG)) * EXP(ALOG(TU(IVEG)-TV)*D2)
IF (RSTFAC(IVEG,3).LT.0.) RSTFAC(IVEG,3) = 0.
IF (RSTFAC(IVEG,3).GT.1.) RSTFAC(IVEG,3) = 1.
XDRR(1)=-PHSOIL(1)
XDRR(2)=-PHSOIL(2)
XDRR(3)=-PHSOIL(3)
IF(XDRR(1).le.0.001) XDRR(1)=0.001
IF(XDRR(2).le.0.001) XDRR(2)=0.001
IF(XDRR(3).le.0.001) XDRR(3)=0.001
XDRR(1)=ALOG(XDRR(1))
XDRR(2)=ALOG(XDRR(2))
XDRR(3)=ALOG(XDRR(3))
IF (NROOT.EQ.1) THEN
XROT = ROOTD(1)
DO 7400 I = 1, 3
7400 DEP(I) = 0.
DO 7500 I = 1, 3
DEP(I) = AMIN1(ZDEPTH(I), XROT)
XROT = XROT - ZDEPTH(I)
IF (XROT.LE.0.) GO TO 7410
7500 CONTINUE
7410 CONTINUE
XDR=(XDRR(1)*DEP(1)+XDRR(2)*DEP(2)+XDRR(3)*DEP(3))/ROOTD(1)
ELSE
XDR=XDRR(1)*ROOTP(1)+XDRR(2)*ROOTP(2)+XDRR(3)*ROOTP(3)
END IF
RSTFAC(IVEG,2) = 1. - EXP(- PH1(IVEG) * (PH2(IVEG) - XDR))
IF (RSTFAC(IVEG,2).GT.1.) RSTFAC(IVEG,2) = 1.
IF (RSTFAC(IVEG,2).LT.0.) RSTFAC(IVEG,2) = 0.
200 RST(IVEG) = RSTM(IVEG)
EPOT = ETV - EA
EPOT = AMAX1(0.0001,(ETV-EA))
RSTFAC(IVEG,1) = 1./ ( 1 + DEFAC(IVEG)*DROP )
IF (RSTFAC(IVEG,1).LT.0.) RSTFAC(IVEG,1) = 0.
IF (RSTFAC(IVEG,1).GT.1.) RSTFAC(IVEG,1) = 1.
300 FTPD = RSTFAC(IVEG,1) * RSTFAC(IVEG,2) * RSTFAC(IVEG,3)
RSTFAC(IVEG,4) = AMAX1( FTPD, 0.00001 )
RST(IVEG) = RST(IVEG) / RSTFAC(IVEG,4) / VCOVER(IVEG)
RST(IVEG) = AMIN1( RST(IVEG), 100000. )
1000 CONTINUE
END SUBROUTINE STRES1
SUBROUTINE TEMRS1 &
(DTT,TC,TGS,TD,TA,TM,QM,EM,PSURF,WWW,CAPAC,SATCAP, &
DTC,DTG,RA,RST,ZDEPTH,BEE,PHSAT,POROS,D,Z0,RDC,RBC,VCOVER,Z2, &
ZLAI,DEFAC,TU,TL,TOPT,RSTFAC,NROOT,ROOTD,PHSOIL,ROOTP,PH1,PH2, &
ECT,ECI,EGT,EGI,EGS,HC,HG,EC,EG,EA,RADT,CHF,SHF,ALBEDO,ZLWUP, &
THERMK,RHOAIR,ZWIND,UM,USTAR,DRAG,CCX,CG, ISNOW,SNOWDEN, &
BPS,rib,CU,XCT,flup,iii,jjj)
REAL ZINC(3), A2(3), Y1(3), ITEX(3),RSTM(2)
DIMENSION WWW(3), CAPAC(2), SATCAP(2), ZDEPTH(3)
DIMENSION VCOVER(2), ZLAI(2), RADT(2),ALBEDO(2,3,2)
DIMENSION TOPT(2), TL(2), TU(2), DEFAC(2)
DIMENSION PH1(2), PH2(2), RST(2), RSTFAC(2,4)
DIMENSION ROOTD(2), ROOTP(3), PHSOIL(3)
E(X) = EXP( 21.18123 - 5418. / X ) / .622
GE(X) = EXP( 21.18123 - 5418. / X ) * 5418. &
/ (X*X) / .622
ETC = E(TC)
ETGS = E(TGS)
GETC = GE(TC)
GETGS = GE(TGS)
PSY = CPAIR / HLAT * PSURF/100. / .622
RCP = RHOAIR * CPAIR
WC = AMIN1( 1., CAPAC(1)/SATCAP(1) )
WG = AMIN1( 1., CAPAC(2)/SATCAP(2) )
FAC = AMIN1( www(1), 0.99 )
FAC = AMAX1( FAC, 0.02 )
RSOIL = 101840. * (1. - FAC ** 0.0027)
PSIT = PHSAT * FAC ** (- BEE )
ARGG = AMAX1(-10.,(PSIT*GRAV/461.5/TGS))
HR = EXP(ARGG)
PILPHR = HR
RESD = D
RESZ0 = Z0
RESRDC = RDC
RESRBC = RBC
RESV2 = VCOVER(2)
IF ( TGS .GT. TF ) GO TO 100
SDEP = CAPAC(2) *SNOWDEN
SDEP = AMIN1( SDEP, (Z2*0.95) )
D = Z2 - ( Z2-D ) / Z2 * ( Z2 - SDEP )
Z0 = Z0 / ( Z2-RESD ) * ( Z2-D )
RDC = RDC * ( Z2-SDEP ) / Z2
RBC = RBC * Z2 / ( Z2-SDEP )
VCOVER(2) = 1.
WG = AMIN1( 1., CAPAC(2) / 0.004 )
RST(2) = RSOIL
100 CONTINUE
IFIRST = 1
ICOUNT = 0
TGEN = TGS
TCEN = TC
FC = 1.
FG = 1.
TRIB = TA
EA = EM
HT = 0.
IONCE = 0
1000 CONTINUE
ICOUNT = ICOUNT + 1
CALL RASIT5(TRIB,CTNI,CUNI,RA,Z2,Z0,D,ZWIND,UM, &
RHOAIR,TM,U2,USTAR,DRAG,TA,bps,rib,CU,XCT,iii,jjj)
IF ( IFIRST .EQ. 1 ) THEN
RB = 1.0/(SQRT(U2)/RBC+ZLAI(1)*.004)
TGTA = TGS- TA
TEMDIF = ( TGTA + SQRT(TGTA*TGTA) ) / 2. + 0.1
FIH = SQRT( 1. + 9. * GRAV *TEMDIF * Z2 / TGS / ( U2*U2) )
RD = RDC / U2 / FIH
ENDIF
D1 = 1./RA + 1./RB + 1./RD
TA = ( TGS/RD + TC/RB + TM/RA *bps) / D1
HT = ( TA - TM ) * RCP / RA
RCC = RST(1)*FC + 2. * RB
COC = (1.-WC)/RCC + WC/(2.*RB)
RG = RST(2)*FG
RSURF = RSOIL*FG
COG1 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD)*HR &
+ VCOVER(2)/(RSURF+RD+44.)*HR
COG2 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD) &
+ VCOVER(2)/(RSURF+RD+44.)
COG1 = COG1 + WG/RD * VCOVER(2)
COG2 = COG2 + WG/RD * VCOVER(2)
D2 = 1./RA + COC + COG2
TOP = COC * ETC + COG1 * ETGS + EM / RA
EA = TOP / D2
DROP = AMAX1( 0., (E(TA)-EA) )
CALL STRES1 (IFIRST, RSTM,ROOTP, &
RSTFAC,RST,TC,ETC,RB,TGS,ETGS,RD,TU,TL,TOPT,EA, &
DEFAC,PH1,PH2,NROOT,ZDEPTH,PHSOIL,ROOTD,VCOVER,DROP)
IFIRST = 0
ERIB = EA
TRIB = TA
IF ( ICOUNT .LE. 4 ) GO TO 1000
TC3 = TC * TC * TC
TG3 = TGS * TGS * TGS
FAC1 = ( 1. - ALBEDO(1,3,2) ) * ( 1.-THERMK ) * VCOVER(1)
FAC2 = 1. - ALBEDO(2,3,2)
RNCDTC = - 2. * 4. * FAC1 * STEFAN * TC3
RNCDTG = 4. * FAC1 * FAC2 * STEFAN * TG3
RNGDTG = - 4. * FAC2 * STEFAN * TG3
RNGDTC = 4. * FAC1 * FAC2 * STEFAN * TC3
IF ( EA .GT. ETC ) FC = 0.
IF ( EA .GT. ETGS) FG = 0.
I = 0
NOX = 0
NONPOS = 1
IWALK = 0
LX = 2
FINC = 1.
ITEX(LX) = 0.
ZINC(LX) = 0.
A2(LX) = 0.
Y1(LX) = 0.
2000 CONTINUE
CALL RASIT5(TRIB,CTNI,CUNI,RA,Z2,Z0,D,ZWIND,UM, &
RHOAIR,TM,U2,USTAR,DRAG,TA,bps,rib,CU,XCT,iii,jjj)
RCP = RHOAIR * CPAIR
D1 = 1./RA + 1./RB + 1./RD
TA = ( TGS/RD + TC/RB + TM/RA *bps) / D1
HC = RCP * ( TC - TA ) / RB * DTT
HG = RCP * ( TGS - TA ) / RD * DTT
HCDTC = RCP / RB * ( 1./RA + 1./RD ) / D1
HCDTG = - RCP / ( RB * RD ) / D1
HCDTM = - RCP / ( RB * RA ) / D1 * BPS
HGDTG = RCP / RD * ( 1./RA + 1./RB ) / D1
HGDTC = - RCP / ( RD * RB ) / D1
HGDTM = - RCP / ( RD * RA ) / D1 *BPS
HRR = HR
IF ( FG .LT. .5 ) HRR = 1.
RCC = RST(1)*FC + 2. * RB
COC = (1.-WC)/RCC + WC/(2.*RB)
RG = RST(2)*FG
IF (ISNOW.eq.0) THEN
RSURF=RSOIL
ELSE
RSURF = RSOIL*FG
END IF
COG1 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD)*HRR &
+ VCOVER(2)/(RSURF+RD+44.)*HRR
COG2 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD) &
+ VCOVER(2)/(RSURF+RD+44.)
COG1 = COG1 + WG/RD * VCOVER(2)
COG2 = COG2 + WG/RD * VCOVER(2)
D2 = 1./RA + COC + COG2
TOP = COC * ETC + COG1 * ETGS + EM/RA
EA = TOP / D2
EC = ( ETC - EA ) * COC * RCP/PSY * DTT
EG = ( ETGS*COG1 - EA*COG2 ) * RCP/PSY * DTT
DEADTC = GETC * COC / D2
DEADTG = GETGS * COG1 / D2
ECDTC = ( GETC - DEADTC ) * COC * RCP / PSY
ECDTG = - DEADTG * COC * RCP / PSY
EGDTG = ( GETGS*COG1 - DEADTG*COG2 ) * RCP / PSY
EGDTC = - DEADTC * COG2 * RCP / PSY
DEADQM = 0.622 * PSURF /( (0.622+QM)**2 * RA * D2 )
ECDQM = -DEADQM * COC * RCP / PSY
EGDQM = -DEADQM * COG2 * RCP / PSY
AK = 1/ RCP / BPS
AH = 1/ (HLAT*RHOAIR)
CCODTC = CCX / DTT - RNCDTC + HCDTC + ECDTC
CCODTG = - RNCDTG + HCDTG + ECDTG
CCORHS = RADT(1) - ( HC + EC ) / DTT
GCODTG = CG / DTT + TIMCON*CG*2. - RNGDTG + HGDTG + EGDTG
GCODTC = - RNGDTC + HGDTC + EGDTC
GCORHS = RADT(2) - TIMCON*CG*2. * ( TGS -TD ) - ( HG + EG ) / DTT
DENOM = CCODTC * GCODTG - CCODTG * GCODTC
DTC = ( CCORHS * GCODTG - CCODTG * GCORHS ) / DENOM
DTG = ( CCODTC * GCORHS - CCORHS * GCODTC ) / DENOM
ECPOT = ( (ETC - EA) + (GETC - DEADTC)*DTC - DEADTG*DTG )
ECI = ECPOT * WC /(2.*RB) * RCP/PSY * DTT
ECIDIF=AMAX1(0.0,(ECI-CAPAC(1)*1.E3*HLAT))
ECI =AMIN1(ECI,( CAPAC(1)*1.E3*HLAT))
EGPOT = ( (ETGS - EA) + (GETGS - DEADTG)*DTG - DEADTC*DTC )
EGI = EGPOT * VCOVER(2) * WG/RD * RCP/PSY * DTT
EGIDIF=AMAX1(0.0,(EGI-CAPAC(2)*1.E3*HLAT))
EGI =AMIN1(EGI,( CAPAC(2)*1.E3*HLAT))
TGEN = TGS + DTG
TCEN = TC + DTC
D1 = 1./RA + 1./RB + 1./RD
TAEN = ( TGEN / RD + TCEN / RB + TM / RA *bps) / D1
HEND = ( TAEN - TM ) * RCP / RA + (ECIDIF + EGIDIF)/DTT
Y= TRIB - TAEN
I = I + 1
HT = HEND
IF ( I .GT. 20 ) GO TO 200
CALL NEWTON(TRIB,Y,FINC,NOX,NONPOS,IWALK,LX,ZINC,A2,Y1,ITEX)
IF(NOX.NE.1)GO TO 2000
200 CONTINUE
HRR = HR
IF ( FG .LT. .5 ) HRR = 1.
RSURF = RSOIL*FG
COCT = (1.-WC)/RCC
COGT = VCOVER(2) * (1.-WG)/( RG + RD )
COGS1 = (1.-VCOVER(2)) / ( RD + RSURF ) * HRR &
+ VCOVER(2) / ( RD + RSURF + 44.) * HRR
COGS2 = COGS1 / HRR
ECT = ECPOT * COCT * RCP/PSY * DTT
EGT = EGPOT * COGT * RCP/PSY * DTT
EGS = (ETGS + GETGS*DTG ) * COGS1 &
- ( EA + DEADTG*DTG + DEADTC*DTC ) * COGS2
EGS = EGS * RCP/PSY * DTT
EGSMAX = WWW(1) / 2. * ZDEPTH(1) * POROS * HLAT * 1000.
EGIADD = AMAX1( 0., EGS - EGSMAX )
EGS = AMIN1 ( EGS, EGSMAX )
EGIDIF = EGIDIF + EGIADD
HC = HC + (HCDTC*DTC + HCDTG*DTG)*DTT + ECIDIF
HG = HG + (HGDTC*DTC + HGDTG*DTG)*DTT + EGIDIF
ECF = SIGN( 1., ECPOT )
EGF = SIGN( 1., EGPOT )
DEWC = FC * 2. - 1.
DEWG = FG * 2. - 1.
IF(DEWC*ECF.GT.0.0) GO TO 300
HC = HC + ECI + ECT
ECI = 0.
ECT = 0.
300 IF(DEWG*EGF.GT.0.0) GO TO 400
HG = HG + EGS + EGI + EGT
EGS = 0.
EGI = 0.
EGT = 0.
400 CONTINUE
EC = ECI + ECT
EG = EGT + EGS + EGI
TC = TCEN
TGS = TGEN
TA = TAEN
EA = EA + DEADTC*DTC + DEADTG*DTG
RADT(1) = RADT(1) + RNCDTC*DTC + RNCDTG*DTG
RADT(2) = RADT(2) + RNGDTC*DTC + RNGDTG*DTG
FLUP = FLUP - (RNCDTC+RNGDTC)*DTC - (RNCDTG+RNGDTG)*DTG
CHF = CCX / DTT * DTC
SHF = CG / DTT * DTG + TIMCON*CG*2. * ( TGS - TD )
ZLWUP = ZLWUP - RNCDTC * DTC / 2. &
- RNGDTG * DTG * (1.-VCOVER(1)*(1.-THERMK) )
IF ( TGS .GT. TF ) GO TO 500
EGS = EG - EGI
EGT = 0.
500 CONTINUE
VCOVER(2) = RESV2
D = RESD
Z0 = RESZ0
RDC = RESRDC
RBC = RESRBC
END SUBROUTINE TEMRS1
SUBROUTINE TEMRS2 &
(DTT,TC,TGS,TD,TA,TM,QM,EM,PSURF,WWW,CAPAC,SATCAP, &
DTC,DTG,RA,RST,ZDEPTH,BEE,PHSAT,POROS,D,Z0,RDC,RBC,VCOVER, &
Z2,ZLAI,DEFAC,TU,TL,TOPT,RSTFAC,NROOT,ROOTD,PHSOIL,ROOTP, &
PH1,PH2,ECT,ECI,EGT,EGI,EGS,HC,HG,EC,EG,EA,RADT,CHF,SHF, &
ALBEDO,ZLWUP,THERMK,RHOAIR,ZWIND,UM,USTAR,DRAG,CCX,CG, &
ISNOW,CHISL,TSOIL,SOLSOIL,CSOIL,WFSOIL,POROSITY, &
DZ,DZO,W,WO,WF,TSSN,TSSNO,BW,BWO,CT,CTO,FI,FIO,FL,FLO, &
BL,BLO,BI,BIO,BT,BTO,DMLT,DMLTO,H,HO,S,SO,SS,SSO,DSOL,DHP, &
DICEVOL,DLIQVOL,THK,QK,SWE,SNOWDEN,SNOWDEPTH,TKAIR,SNROFF, &
DZSOIL,BPS,rib,CU,XCT,flup,iii,jjj)
REAL WORK(N1),WORK1(N1),DELTH(20)
data DELTH/20*0.0/
REAL ZINC(3), A2(3), Y1(3), ITEX(3),RSTM(2)
DIMENSION SSO(N1),POROSITY(N1),H(N1),HO(N1),DZ(N1),DZO(N1),CT(N1), &
BI(N1),BIO(N1),BW(N1),BWO(N1),BL(N1),BLO(N1),CTO(N1), &
TSSN(N1),TSSNO(N1),DLIQVOL(N1),DICEVOL(N1),DSOL(N1), &
W(N1),WO(N1),WF(N1),FI(N1),FIO(N1),FL(N1),FLO(N1), &
DMLT(N1),DMLTO(N1),BT(N1),BTO(N1),S(N1),SO(N1),SS(N1), &
PDZDTC(N1),DMASS(N1),THK(N1),DHP(N1),QK(N1)
DIMENSION WWW(3),CAPAC(2),SATCAP(2),ZDEPTH(3),VCOVER(2),ZLAI(2), &
RADT(2),ALBEDO(2,3,2),TOPT(2),TL(2),TU(2),DEFAC(2), &
PH1(2),PH2(2),RST(2),RSTFAC(2,4), &
ROOTD(2),ROOTP(3),PHSOIL(3)
E(X) = EXP( 21.18123 - 5418. / X ) / .622
GE(X) = EXP( 21.18123 - 5418. / X ) * 5418. &
/ (X*X) / .622
ETC = E(TC)
ETGS = E(TGS)
GETC = GE(TC)
GETGS = GE(TGS)
PSY = CPAIR / HLAT * PSURF/ 100. / .622
RCP = RHOAIR * CPAIR
WC = AMIN1( 1., CAPAC(1)/SATCAP(1) )
IF (ISNOW.eq.0) THEN
WG=1.0
ELSE
WG = AMIN1( 1., CAPAC(2)/SATCAP(2) )
END IF
FAC = AMIN1( www(1), 0.99 )
FAC = AMAX1( FAC, 0.02 )
IF (ISNOW.eq.0) THEN
RSOIL=10000000000.
ELSE
RSOIL = 101840. * (1. - FAC ** 0.0027)
END IF
PSIT = PHSAT * FAC ** (- BEE )
ARGG = AMAX1(-10.,(PSIT*GRAV/461.5/TGS))
HR = EXP(ARGG)
PILPHR = HR
RESD = D
RESZ0 = Z0
RESRDC = RDC
RESRBC = RBC
RESV2 = VCOVER(2)
IF ( TGS .GT. TF ) GO TO 100
IF (ISNOW.eq.0) THEN
SDEP = SNOWDEPTH
ELSE
SDEP = CAPAC(2) * SNOWDEN
END IF
SDEP = AMIN1( SDEP, (Z2*0.95) )
D = Z2 - ( Z2-D ) / Z2 * ( Z2 - SDEP )
Z0 = Z0 / ( Z2-RESD ) * ( Z2-D )
RDC = RDC * ( Z2-SDEP ) / Z2
RBC = RBC * Z2 / ( Z2-SDEP )
VCOVER(2) = 1.
IF (ISNOW.eq.0) THEN
WG=1.0
ELSE
WG = AMIN1( 1., CAPAC(2) / 0.004 )
END IF
RST(2) = RSOIL
100 CONTINUE
IFIRST = 1
ICOUNT = 0
TGEN = TGS
TCEN = TC
FC = 1.
FG = 1.
TRIB = TA
EA = EM
HT = 0.
IONCE = 0
1000 CONTINUE
ICOUNT = ICOUNT + 1
CALL RASIT5(TRIB,CTNI,CUNI,RA,Z2,Z0,D,ZWIND,UM, &
RHOAIR,TM,U2,USTAR,DRAG,TA,bps,rib,CU,XCT,iii,jjj)
IF ( IFIRST .EQ. 1 ) THEN
RB = 1.0/(SQRT(U2)/RBC+ZLAI(1)*.004)
TGTA = TGS- TA
TEMDIF = ( TGTA + SQRT(TGTA*TGTA) ) / 2. + 0.1
FIH = SQRT( 1. + 9. * GRAV *TEMDIF * Z2 / TGS / ( U2*U2) )
RD = RDC / U2 / FIH
ENDIF
D1 = 1./RA + 1./RB + 1./RD
TA = ( TGS/RD + TC/RB + TM/RA *bps) / D1
HT = ( TA - TM ) * RCP / RA
RCC = RST(1)*FC + 2. * RB
COC = (1.-WC)/RCC + WC/(2.*RB)
RG = RST(2)*FG
IF (ISNOW.eq.0) THEN
RSURF = RSOIL
ELSE
RSURF = RSOIL*FG
END IF
COG1 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD)*HR &
+ VCOVER(2)/(RSURF+RD+44.)*HR
COG2 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD) &
+ VCOVER(2)/(RSURF+RD+44.)
COG1 = COG1 + WG/RD * VCOVER(2)
COG2 = COG2 + WG/RD * VCOVER(2)
D2 = 1./RA + COC + COG2
TOP = COC * ETC + COG1 * ETGS + EM / RA
EA = TOP / D2
DROP = AMAX1( 0., (E(TA)-EA) )
CALL STRES1 (IFIRST, RSTM,ROOTP, &
RSTFAC,RST,TC,ETC,RB,TGS,ETGS,RD,TU,TL,TOPT,EA, &
DEFAC,PH1,PH2,NROOT,ZDEPTH,PHSOIL,ROOTD,VCOVER,DROP)
IFIRST = 0
ERIB = EA
TRIB = TA
IF ( ICOUNT .LE. 4 ) GO TO 1000
TC3 = TC * TC * TC
TG3 = TGS * TGS * TGS
FAC1 = ( 1. - ALBEDO(1,3,2) ) * ( 1.-THERMK ) * VCOVER(1)
FAC2 = 1. - ALBEDO(2,3,2)
RNCDTC = - 2. * 4. * FAC1 * STEFAN * TC3
RNCDTG = 4. * FAC1 * FAC2 * STEFAN * TG3
RNGDTG = - 4. * FAC2 * STEFAN * TG3
RNGDTC = 4. * FAC1 * FAC2 * STEFAN * TC3
IF ( EA .GT. ETC ) FC = 0.
IF ( EA .GT. ETGS) FG = 0.
II = 0
NOX = 0
NONPOS = 1
IWALK = 0
LX = 2
FINC = 1.
ITEX(LX) = 0.
ZINC(LX) = 0.
A2(LX) = 0.
Y1(LX) = 0.
IF (ISNOW.eq.0) THEN
CALL TPROPTY(CHISL,BWO,DZO,TKAIR,DZSOIL, THK,QK)
tssn(n+1) = tkair
icount = 0
do i=1,n
work(i) = tssno(i)
work1(i) = dliqvol(i)
end do
hx = 0.0
NK=n
ELSE
NK=1
END IF
RADDWN=solsoil
RADDWN=RADDWN+dsol(1)+dsol(2)
RNG = RADT(2) - RADDWN
RADT(2)=RNG
do 57 ik = NK , 1 , -1
IF (ISNOW.ne.0) go to 2000
If((sso(ik).lt.1d0.and.porosity(ik).gt.0d0))then
udum0 = dzo(ik)*(porosity(ik) -work1(ik))
if(udum0.lt.0.0) then
print*,' udum0 is WRONG in thermal.f'
STOP
endif
if(wf(ik+1).gt.udum0)then
uuu=udum0
snroff = snroff + (wf(ik+1)-udum0)
hroff=hroff+(wf(ik+1)-udum0)*cl*rhowater*(tssn(ik+1)-273.16)
else
uuu=wf(ik+1)
endif
dhp(ik+1)=(uuu*cl*rhowater*(tssn(ik+1)-273.16))/dtt
w(ik)=wo(ik)+ uuu
bwo(ik)=rhowater*w(ik)/dzo(ik)
cto(ik)=(bwo(ik)/920.0)*1.9e+6
dmlto(ik)=w(ik)*dlm*rhowater
if (ho(ik).lt.-dmlto(ik)) then
fio(ik)=1.0
flo(ik)=0.0
tssno(ik)=( ho(ik)+dmlto(ik))/(cto(ik)*dzo(ik))+273.16
else
tssno(ik)=273.16
fio(ik)=-ho(ik)/dmlto(ik)
flo(ik)=1.0-fio(ik)
end if
blo(ik)=bwo(ik)*flo(ik)
bio(ik)=bwo(ik)*fio(ik)
dliqvol(ik)=blo(ik)/rhowater
dicevol(ik)=bio(ik)/dice
Else
w(ik)=wo(ik)
snroff = snroff +wf(ik+1)
hroff=hroff+wf(ik+1)*cl*rhowater*(tssn(ik+1)-273.16)
dhp(ik+1) = 0.0
End if
TGS=tssno(NK)
If (ik.lt.Nk) Then
if(ik.ne.1) then
b1 = dsol(ik) + delth(ik) &
+ qk(ik+1)*(tssn(ik+1)-tssno(ik)) + qk(ik)*work(ik-1)
else
b1 = dsol(ik) + delth(ik) &
+ qk(ik+1)*(tssn(ik+1)-tssno(ik)) + qk(ik)*tsoil
endif
b2 = - qk(ik)
delth(ik) = 0.0
End if
dmlt(ik)=w(ik)*dlm*rhowater
If (ik.lt.NK.and.ik.ge.1) Then
fff = -( ho(ik) + (b1+b2*273.16)*dtt ) &
/ ( rhowater*dlm*w(ik) )
if(fff.gt.0.0.and.fff.le.1.0) then
ICASE=2
fi(ik)=fff
else if (fff.gt.1.0) then
ICASE=1
else if (fff.le.0.0) then
ICASE=3
end if
End if
If (ik.lt.NK) go to 3000
2000 CONTINUE
CALL RASIT5(TRIB,CTNI,CUNI,RA,Z2,Z0,D,ZWIND,UM, &
RHOAIR,TM,U2,USTAR,DRAG,TA,bps,rib,CU,XCT,iii,jjj)
RCP = RHOAIR * CPAIR
D1 = 1./RA + 1./RB + 1./RD
TA = ( TGS/RD + TC/RB + TM/RA *bps) / D1
HC = RCP * ( TC - TA ) / RB * DTT
HG = RCP * ( TGS - TA ) / RD * DTT
HCDTC = RCP / RB * ( 1./RA + 1./RD ) / D1
HCDTG = - RCP / ( RB * RD ) / D1
HCDTM = - RCP / ( RB * RA ) / D1 * BPS
HGDTG = RCP / RD * ( 1./RA + 1./RB ) / D1
HGDTC = - RCP / ( RD * RB ) / D1
HGDTM = - RCP / ( RD * RA ) / D1 *BPS
HRR = HR
IF ( FG .LT. .5 ) HRR = 1.
RCC = RST(1)*FC + 2. * RB
COC = (1.-WC)/RCC + WC/(2.*RB)
RG = RST(2)*FG
IF (ISNOW.eq.0) THEN
RSURF=RSOIL
ELSE
RSURF = RSOIL*FG
END IF
COG1 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD)*HRR &
+ VCOVER(2)/(RSURF+RD+44.)*HRR
COG2 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD) &
+ VCOVER(2)/(RSURF+RD+44.)
COG1 = COG1 + WG/RD * VCOVER(2)
COG2 = COG2 + WG/RD * VCOVER(2)
D2 = 1./RA + COC + COG2
TOP = COC * ETC + COG1 * ETGS + EM/RA
EA = TOP / D2
EC = ( ETC - EA ) * COC * RCP/PSY * DTT
EG = ( ETGS*COG1 - EA*COG2 ) * RCP/PSY * DTT
DEADTC = GETC * COC / D2
DEADTG = GETGS * COG1 / D2
ECDTC = ( GETC - DEADTC ) * COC * RCP / PSY
ECDTG = - DEADTG * COC * RCP / PSY
EGDTG = ( GETGS*COG1 - DEADTG*COG2 ) * RCP / PSY
EGDTC = - DEADTC * COG2 * RCP / PSY
DEADQM = 0.622 * PSURF /( (0.622+QM)**2 * RA * D2 )
ECDQM = -DEADQM * COC * RCP / PSY
EGDQM = -DEADQM * COG2 * RCP / PSY
AK = 1/ RCP / BPS
AH = 1/ (HLAT*RHOAIR)
CCODTC = CCX / DTT - RNCDTC + HCDTC + ECDTC
CCODTG = - RNCDTG + HCDTG + ECDTG
CCORHS = RADT(1) - ( HC + EC ) / DTT
IF (ISNOW.eq.0) THEN
GCODTG= cto(n)*dzo(n)/DTT - RNGDTG + HGDTG + EGDTG + qk(n)
ELSE
GCODTG = CG / DTT + TIMCON*CG*2. - RNGDTG + HGDTG + EGDTG
END IF
GCODTC = - RNGDTC + HGDTC + EGDTC
IF (ISNOW.ne.0) THEN
GCORHS = RADT(2)-TIMCON*CG*2.*( TGS -TD )-( HG + EG )/ DTT
ELSE
fi(n)=1.0
GCORHS1 = ho(n)/DTT+RNG - ( HG + EG ) / DTT +dhp(n+1) &
- qk(n)*(TGS -tssno(n-1))-cto(n)*dzo(n)*(tssno(n)-273.16)/DTT
GCORHS = GCORHS1+ rhowater*dlm*w(n)*fi(n)/DTT
END IF
DENOM = CCODTC * GCODTG - CCODTG * GCODTC
DTC = ( CCORHS * GCODTG - CCODTG * GCORHS ) / DENOM
DTG = ( CCODTC * GCORHS - CCORHS * GCODTC ) / DENOM
IF (ISNOW.eq.0) THEN
If ((TGS+DTG).le.273.16) Then
TGSNEW=(TGS+DTG)
ICASE=1
h(NK)=( TGSNEW-273.16)*cto(n)*dzo(n)-fi(NK)*w(NK)*dlm*rhowater
Dh_DTT_DTG=(h(NK)-ho(NK))/DTT
tonm1=tssno(NK-1)
qkn=qk(n)
Else
DTG=273.16-TGS
DTC= (CCORHS - CCODTG*DTG)/CCODTC
fi(NK)=(GCODTC*DTC+GCODTG*DTG-GCORHS1)/(rhowater*dlm*w(n))*DTT
if (fi(NK).ge.0.0.and.fi(NK).le.1.0) then
h(NK)=-fi(n)*w(n)*dlm*rhowater
Dh_DTT_DTG=(h(NK)-ho(NK))/DTT
ICASE=2
tonm1=tssno(NK-1)
qkn=qk(n)
else if (fi(NK).lt.0.)then
h(NK)= -fi(NK)*w(NK)*dlm*rhowater
Dh_DTT_DTG=(h(NK)-ho(NK))/DTT
ICASE=3
tonm1=tssno(NK-1)
qkn=qk(n)
fff=fi(NK)
fi(NK)=0.0
end if
End if
END IF
ECPOT = ( (ETC - EA) + (GETC - DEADTC)*DTC - DEADTG*DTG )
ECI = ECPOT * WC /(2.*RB) * RCP/PSY * DTT
ECIDIF=AMAX1(0.0,(ECI-CAPAC(1)*1.E3*HLAT))
ECI =AMIN1(ECI,( CAPAC(1)*1.E3*HLAT))
EGPOT = ( (ETGS - EA) + (GETGS - DEADTG)*DTG - DEADTC*DTC )
EGI = EGPOT * VCOVER(2) * WG/RD * RCP/PSY * DTT
EGIDIF=AMAX1(0.0,(EGI-CAPAC(2)*1.E3*HLAT))
EGI =AMIN1(EGI,( CAPAC(2)*1.E3*HLAT))
TGEN = TGS + DTG
TCEN = TC + DTC
D1 = 1./RA + 1./RB + 1./RD
TAEN = ( TGEN / RD + TCEN / RB + TM / RA *bps) / D1
HEND = ( TAEN - TM ) * RCP / RA + (ECIDIF + EGIDIF)/DTT
Y= TRIB - TAEN
II = II + 1
HT = HEND
IF ( II .GT. 20 ) GO TO 200
CALL NEWTON(TRIB,Y,FINC,NOX,NONPOS,IWALK,LX,ZINC,A2,Y1,ITEX)
IF(NOX.NE.1)GO TO 2000
200 CONTINUE
3000 IF (ISNOW.eq.0) THEN
IF (ICASE.eq.1.and.ik.eq.NK) THEN
tssn(NK)=TGS+DTG
END IF
If (ik.eq.NK) then
SNOFAC = HLAT / (HLAT + SNOMEL /1000.)
egidw = EGI*SNOFAC /HLAT/1000.
w(n)=w(n)-egidw
swe=swe-egidw
dzold=dzo(n)
dzo(n)=dzo(n)-egidw*rhowater/bwo(n)
ho(n)=ho(n)*dzo(n)/dzold
capac(2)=swe
End if
CALL SNRESULT(DTT,IK,ICASE,WFSOIL,TSOIL,B1,B2,FFF,DELTH,WWW, &
ZDEPTH,POROS,BI,BIO,DZ,DZO,W,BW,BWO,H,HO,QK,BL,BT,CT, &
FI,FL,WF,TSSN,DLIQVOL,DICEVOL,SNROFF,HROFF,QSOIL)
END IF
57 CONTINUE
SNOWDEPTH=DZO(1)+DZO(2)+DZO(3)
IF (ISNOW.eq.0) THEN
SWE=W(1)+W(2)+W(3)
CAPAC(2)=SWE
IF((DZ(1)+DZ(2)+DZ(3)).NE.0.0) THEN
SNOWDEN=(BW(1)*DZ(1)+BW(2)*DZ(2)+BW(3)*DZ(3)) &
/(DZ(1)+DZ(2)+DZ(3))
SNOWDEN=1000./SNOWDEN
ENDIF
ENDIF
HRR = HR
IF ( FG .LT. .5 ) HRR = 1.
IF (ISNOW.eq.0) THEN
RSURF = RSOIL
ELSE
RSURF = RSOIL*FG
END IF
COCT = (1.-WC)/RCC
COGT = VCOVER(2) * (1.-WG)/( RG + RD )
COGS1 = (1.-VCOVER(2)) / ( RD + RSURF ) * HRR &
+ VCOVER(2) / ( RD + RSURF + 44.) * HRR
COGS2 = COGS1 / HRR
ECT = ECPOT * COCT * RCP/PSY * DTT
EGT = EGPOT * COGT * RCP/PSY * DTT
EGS = (ETGS + GETGS*DTG ) * COGS1 &
- ( EA + DEADTG*DTG + DEADTC*DTC ) * COGS2
EGS = EGS * RCP/PSY * DTT
IF (ISNOW.eq.0) EGS=0.0
EGSMAX = WWW(1) / 2. * ZDEPTH(1) * POROS * HLAT * 1000.
EGIADD = AMAX1( 0., EGS - EGSMAX )
EGS = AMIN1 ( EGS, EGSMAX )
EGIDIF = EGIDIF + EGIADD
HC = HC + (HCDTC*DTC + HCDTG*DTG)*DTT + ECIDIF
HG = HG + (HGDTC*DTC + HGDTG*DTG)*DTT + EGIDIF
IF (ISNOW.eq.0) go to 401
ECF = SIGN( 1., ECPOT )
EGF = SIGN( 1., EGPOT )
DEWC = FC * 2. - 1.
DEWG = FG * 2. - 1.
IF(DEWC*ECF.GT.0.0) GO TO 300
HC = HC + ECI + ECT
ECI = 0.
ECT = 0.
300 IF(DEWG*EGF.GT.0.0) GO TO 400
HG = HG + EGS + EGI + EGT
EGS = 0.
EGI = 0.
EGT = 0.
400 CONTINUE
401 CONTINUE
EC = ECI + ECT
EG = EGT + EGS + EGI
TGSOLD=TGS
TC = TCEN
TGS = TGEN
IF (ISNOW.eq.0) tssn(n)=TGS
TA = TAEN
EA = EA + DEADTC*DTC + DEADTG*DTG
RADT(1) = RADT(1) + RNCDTC*DTC + RNCDTG*DTG
RADT(2) = RADT(2) + RNGDTC*DTC + RNGDTG*DTG
FLUP = FLUP - (RNCDTC+RNGDTC)*DTC - (RNCDTG+RNGDTG)*DTG
CHF = CCX / DTT * DTC
IF (ISNOW.eq.0) THEN
SHF= Dh_DTT_DTG - dhp(n+1)+ qkn*(TGSOLD-tonm1) &
+qkn*DTG
ELSE
SHF = CG / DTT * DTG + TIMCON*CG*2. * ( TGS - TD )
END IF
ZLWUP = ZLWUP - RNCDTC * DTC / 2. &
- RNGDTG * DTG * (1.-VCOVER(1)*(1.-THERMK) )
IF ( TGS .GT. TF ) GO TO 500
EGS = EG - EGI
EGT = 0.
500 CONTINUE
VCOVER(2) = RESV2
D = RESD
Z0 = RESZ0
RDC = RESRDC
RBC = RESRBC
IF (ISNOW.eq.0) THEN
TGS=TSOIL
ATMP= (QSOIL+SOLSOIL)/CSOIL
BTMP=2.*3.1416/86400.
CTMP=CSOIL*BTMP/CSOIL/(365.*3.1416)**0.5
TGS=(TSOIL+ATMP*DTT+BTMP*DTT*TD/(1.+CTMP*DTT))/ &
(1.+BTMP*DTT*(1.-CTMP*DTT/(1.+CTMP*DTT)))
TD=(CTMP*DTT*TGS+TD)/(1.+CTMP*DTT)
END IF
END SUBROUTINE TEMRS2
SUBROUTINE TPROPTY(THKSOIL,BWO,DZO,TKAIR,DZSOIL, THK,QK)
DIMENSION BWO(N1),DZO(N1),THK(N1),QK(N1)
do 37 i=1,n
thkice=2.290d0
thkair=2.30d-2
thk(i) = thkair+(7.75d-5 *bwo(i)+ 1.105d-6* &
bwo(i)*bwo(i))*(thkice -thkair)+0.1
37 continue
do 47 i=2,n
qk(i)=2.0*thk(i)*thk(i-1)/(thk(i)*dzo(i-1)+thk(i-1)*dzo(i))
47 continue
qk(1)= 2.0*thk(1)*thksoil/(thk(1)*dzsoil+thksoil*dzo(1))
END SUBROUTINE TPROPTY
SUBROUTINE UPDAT1(DTT,TC,TGS,TD,CAPAC,DTC,DTG,ECT,ECI,EGT,EGI, &
EGS,EG,HC,HG,HFLUX,ETMASS,ROFF, &
NROOT,ROOTD,ROOTP,POROS,BEE,SATCO,SLOPE, &
PHSAT,ZDEPTH,WWW,CCX,CG,CHF,SHF,ISNOW,WFSOIL,SWE,SNROFF,smelt)
DIMENSION EF(3)
DIMENSION WWW(3), CAPAC(2),SNOWW(2),ROOTD(2), ZDEPTH(3), ROOTP(3)
DIMENSION TEMW(3),TEMWP(3),TEMWPP(3),AAA(2),BBB(2),CCC(2),QQQ(2)
SNOFAC = HLAT / ( HLAT + SNOMEL /1000. )
FACKS = 1.
IF ( (TC-DTC) .LE. TF ) FACKS = SNOFAC
IF ( (ECT+ECI) .GT. 0.) GO TO 100
ECI = ECT + ECI
ECT = 0.
FACKS = 1. / FACKS
100 CAPAC(1)=CAPAC(1) - ECI*FACKS/HLAT/1000.
ECMASS = ( ECT + ECI * FACKS ) / HLAT
IF (ISNOW.eq.0) FACKS = SNOFAC
IF (ISNOW.EQ.0) go to 201
FACKS = 1.
IF ( (TGS-DTG) .LE. TF ) FACKS = SNOFAC
IF ( (EGT+EGI) .GT. 0. ) GO TO 200
EGI = EGT + EGI
EGT = 0.
FACKS = 1. / FACKS
200 CAPAC(2)=CAPAC(2) - EGI*FACKS/HLAT/1000.
201 EGMASS = ( EGT + EGS + EGI * FACKS ) / HLAT
ETMASS = ECMASS + EGMASS
HFLUX = ( HC + HG ) / DTT
DO 1000 IVEG = 1, 2
IF ( CAPAC(IVEG) .GT. 0.000001 ) GO TO 1000
WWW(1) = WWW(1) + CAPAC(IVEG) / ( POROS*ZDEPTH(1) )
CAPAC(IVEG) = 0.
1000 CONTINUE
DO 6000 IVEG = 1, 2
IF (ISNOW.EQ.0.and.IVEG.EQ.2) THEN
ZMELT= WFSOIL
WWW(1) = WWW(1) + ZMELT / ( POROS * ZDEPTH(1) )
CAPAC(2)= SWE
GO TO 6000
END IF
CCT = CCX
TS = TC
DTS = DTC
FLUX = CHF
IF ( IVEG .EQ. 1 ) GO TO 110
CCT = CG
TS = TGS
DTS = DTG
FLUX = CCT * DTG / DTT
110 CONTINUE
TTA = TS - DTS
TTB = TS
SNOWW(IVEG) = 0.
IF ( TTA .LE. TF ) SNOWW(IVEG) = CAPAC(IVEG)
CAPAC(IVEG) = CAPAC(IVEG) - SNOWW(IVEG)
IF ( TTA .GT. TF .AND. TTB .GT. TF ) GO TO 120
IF ( TTA .LE. TF .AND. TTB .LE. TF ) GO TO 120
DTF = TF - TTA
DTIME1 = CCT * DTF / FLUX
HF = FLUX*(DTT-DTIME1)
FCAP = - CAPAC(IVEG) * SNOMEL
SPWET = AMIN1( 5. , SNOWW(IVEG) )
IF ( DTS .GT. 0. ) FCAP = SPWET * SNOMEL
DTIME2 = FCAP / FLUX
DTF2 = FLUX * (DTT-DTIME1-DTIME2)/CCT
TN = TF + DTF2
TS = TF - 0.1
IF (ABS(HF) .GE.ABS(FCAP) ) TS = TN
CHANGE = HF
IF (ABS(CHANGE) .GE.ABS(FCAP) ) CHANGE = FCAP
CHANGE = CHANGE / SNOMEL
IF (CHANGE.GT.0.0) SMELT=CHANGE+SMELT
SNOWW(IVEG) = SNOWW(IVEG) - CHANGE
CAPAC(IVEG) = CAPAC(IVEG) + CHANGE
IF ( IVEG .EQ. 1 ) TC = TS
IF ( IVEG .EQ. 2 ) TGS = TS
IF ( SNOWW(IVEG) .LT. 0.00001 ) GO TO 120
ZMELT = CAPAC(IVEG)
WWW(1) = WWW(1) + ZMELT / ( POROS * ZDEPTH(1) )
CAPAC(IVEG) = 0.
120 CONTINUE
CAPAC(IVEG) = CAPAC(IVEG) + SNOWW(IVEG)
6000 CONTINUE
IF (ISNOW.NE.0) THEN
FLUXEF = SHF - CCT*DTG/DTT
TD = TD + FLUXEF / ( CG * 2. * SQRT ( PIE*365. ) ) * DTT
END IF
change=0.0
WWW(1) = WWW(1) - EGS / HLAT / 1000. / ( POROS * ZDEPTH(1) )
DO 2000 IVEG = 1, 2
IF ( IVEG .EQ. 1 ) ABSOIL = ECT / HLAT / 1000.
IF ( IVEG .EQ. 2 ) ABSOIL = EGT / HLAT / 1000.
IF (NROOT.EQ.1) THEN
EF(2) = 0.
EF(3) = 0.
TOTDEP = ZDEPTH(1)
DO 3000 IL = 2, 3
TOTDEP = TOTDEP + ZDEPTH(IL)
IF ( ROOTD(IVEG) .LT. TOTDEP ) GO TO 400
EF(IL) = ZDEPTH(IL) / ROOTD(IVEG)
GO TO 500
400 CONTINUE
EF(IL) = ROOTD(IVEG) - TOTDEP + ZDEPTH(IL)
EF(IL) = EF(IL) / ROOTD(IVEG)
GO TO 600
500 CONTINUE
3000 CONTINUE
600 EFT = EF(2) + EF(3)
EFT = MAX(EFT,0.1E-5)
EF(2) = EF(2) / EFT
EF(3) = EF(3) / EFT
DO 4000 IL = 2, 3
WWW(IL) = WWW(IL) - ABSOIL * EF(IL) / ( POROS * ZDEPTH(IL) )
4000 CONTINUE
ELSE
EF(1) = ROOTP(1)
EF(2) = ROOTP(2)
EF(3) = ROOTP(3)
DO 4004 IL = 1, 3
WWW(IL) = WWW(IL) - ABSOIL * EF(IL) / ( POROS * ZDEPTH(IL) )
4004 CONTINUE
END IF
2000 CONTINUE
DO 5000 IL = 1, 2
IF ( WWW(IL) .GT. 0. ) GO TO 5000
WWW(IL+1) = WWW(IL+1) + WWW(IL) * ZDEPTH(IL)/ZDEPTH(IL+1)
WWW(IL) = 0.
5000 CONTINUE
do 8000 i = 1, 3
TEMW(I) = AMAX1( 0.03, WWW(I) )
TEMWP(I) = TEMW(I) ** ( -BEE )
TEMWPP(I) = AMIN1( 1., TEMW(I)) ** ( 2.*BEE+ 3. )
8000 CONTINUE
POWS = 2.*BEE+2.
Q3G = TEMW(3)**(-POWS) + SATCO/ZDEPTH(3)/POROS*SLOPE*POWS*DTT
Q3G = Q3G ** ( 1. / POWS )
Q3G = - ( 1. / Q3G - WWW(3) ) * POROS * ZDEPTH(3) / DTT
Q3G = AMAX1( 0., Q3G )
Q3G = AMIN1( Q3G, WWW(3)*POROS*ZDEPTH(3)/DTT )
Q3G = Q3G + 0.002*POROS*ZDEPTH(3)*0.5 / 86400. * WWW(3)
WMAX = AMAX1( WWW(1), WWW(2), WWW(3), 0.05 )
WMAX = AMIN1( WMAX, 1. )
PMAX = WMAX**(-BEE)
WMIN = (PMAX-2./( PHSAT*(ZDEPTH(1)+2.*ZDEPTH(2)+ZDEPTH(3)))) &
**(-1./BEE)
WMIN = AMIN1( WWW(1), WWW(2), WWW(3), WMIN )
WMIN = AMAX1( WMIN, 0.02 )
PMIN = WMIN**(-BEE)
DPDW = PHSAT*( PMAX-PMIN )/( WMAX-WMIN )
do 8200 i = 1, 2
RSAME = 0.
AVK = TEMWP(I)*TEMWPP(I) - TEMWP(I+1)*TEMWPP(I+1)
DIV = TEMWP(I+1) - TEMWP(I)
IF ( ABS(DIV) .LT. 1.E-6 ) RSAME = 1.
AVK = SATCO*AVK / ( ( 1. + 3./BEE ) * DIV + RSAME )
AVKMIN = SATCO * AMIN1( TEMWPP(I), TEMWPP(I+1) )
AVKMAX = SATCO * AMAX1( TEMWPP(I), TEMWPP(I+1) )*1.01
AVK = AMAX1( AVK, AVKMIN )
AVK = AMIN1( AVK, AVKMAX )
TSNOW = AMIN1 ( TF-0.01, TGS )
AREAS = AMIN1 (0.999,13.2*SNOWW(2))
TGG = TSNOW*AREAS + TGS*(1.-AREAS)
TS = TGG*(2-I) + TD*(I-1)
PROPS = ( TS-(TF-10.) ) / 10.
PROPS = AMAX1( 0.05, AMIN1( 1.0, PROPS ) )
AVK = AVK * PROPS
Q3G = Q3G * PROPS
DPDWDZ = DPDW * 2./( ZDEPTH(I) + ZDEPTH(I+1) )
AAA(I) = 1. + AVK*DPDWDZ*( 1./ZDEPTH(I)+1./ZDEPTH(I+1) ) &
*DTT/POROS
BBB(I) =-AVK * DPDWDZ * 1./ZDEPTH(2)*DTT/POROS
CCC(I) = AVK * ( DPDWDZ * ( WWW(I)-WWW(I+1) ) + 1. + &
(I-1)*DPDWDZ*Q3G*1./ZDEPTH(3)*DTT/POROS )
8200 CONTINUE
DENOM = ( AAA(1)*AAA(2) - BBB(1)*BBB(2) )
RDENOM = 0.
IF ( ABS(DENOM) .LT. 1.E-6 ) RDENOM = 1.
RDENOM = ( 1.-RDENOM)/( DENOM + RDENOM )
QQQ(1) = ( AAA(2)*CCC(1) - BBB(1)*CCC(2) ) * RDENOM
QQQ(2) = ( AAA(1)*CCC(2) - BBB(2)*CCC(1) ) * RDENOM
WWW(3) = WWW(3) - Q3G*DTT/(POROS*ZDEPTH(3))
ROFF = ROFF + Q3G * DTT
do 8300 i = 1, 2
QMAX = WWW(I) * (POROS*ZDEPTH(I) /DTT)
QMIN = -WWW(I+1) * (POROS*ZDEPTH(I+1)/DTT)
QQQ(I) = AMIN1( QQQ(I),QMAX)
QQQ(I) = AMAX1( QQQ(I),QMIN)
WWW(I) = WWW(I) - QQQ(I)/(POROS*ZDEPTH(I) /DTT)
WWW(I+1) = WWW(I+1) + QQQ(I)/(POROS*ZDEPTH(I+1)/DTT)
8300 continue
do 8400 i = 1, 3
EXCESS = AMAX1(0.,(WWW(I) - 1.))
WWW(I) = WWW(I) - EXCESS
ROFF = ROFF + EXCESS * POROS*ZDEPTH(I)
8400 continue
do 8402 i = 1,2
DEFICIT = AMAX1 (0.,(1.E-12 - WWW(I)))
WWW (I) = WWW(I) + DEFICIT
WWW (I+1) = WWW(I+1) - DEFICIT * ZDEPTH(I) / ZDEPTH (I+1)
8402 CONTINUE
WWW(3) = AMAX1 (WWW(3),1.E-12)
800 CONTINUE
IF (WWW(1) .GT.1.) THEN
WWW(2) = WWW(2) + (WWW(1)-1.) * ZDEPTH(1) / ZDEPTH(2)
WWW(1) = 1.
END IF
If (WWW(2) .GT.1.) THEN
WWW(3) = WWW(3) + (WWW(2)-1.) * ZDEPTH(2) / ZDEPTH(3)
WWW(2) = 1.
END IF
IF (WWW(3) .GT.1.) THEN
ROFF = ROFF + (WWW(3)-1.)* ZDEPTH(3) * POROS
WWW(3) = 1.
END IF
END SUBROUTINE UPDAT1
SUBROUTINE UPDAT1_ICE(DTT,TC,TGS,TD,CAPAC,DTC,DTG,ECT,ECI,EGT,EGI, &
EGS,EG,HC,HG,HFLUX,ETMASS,FILTR,SOILDIF,SOILDRA,ROFF, &
RNOFFB,RNOFFS,NROOT,ROOTD,ROOTP,POROS,BEE,SATCO,SLOPE, &
PHSAT,ZDEPTH,WWW,CCX,CG,CHF,SHF,SMELT)
REAL, DIMENSION (2) :: CAPAC, SNOWW, ROOTD, aaa, bbb, ccc, qqq
REAL, DIMENSION (3) :: WWW, EF, ZDEPTH, ROOTP, temw, temwp, temwpp
SNOFAC = HLAT / ( HLAT + SNOMEL /1000. )
FACKS = 1.
IF ( (TC-DTC) .LE. TF ) FACKS = SNOFAC
IF ( (ECT+ECI) .GT. 0.) GO TO 100
ECI = ECT + ECI
ECT = 0.
FACKS = 1. / FACKS
100 CAPAC(1)=CAPAC(1) - ECI*FACKS/HLAT/1000.
ECMASS = ( ECT + ECI * FACKS ) / HLAT
FACKS = 1.
IF ( (TGS-DTG) .LE. TF ) FACKS = SNOFAC
IF ( (EGT+EGI) .GT. 0. ) GO TO 200
EGI = EGT + EGI
EGT = 0.
FACKS = 1. / FACKS
200 CAPAC(2)=CAPAC(2) - EGI*FACKS/HLAT/1000.
EGMASS = ( EGT + EGS + EGI * FACKS ) / HLAT
ETMASS = ECMASS + EGMASS
HFLUX = ( HC + HG )
DO 1000 IVEG = 1, 2
IF ( CAPAC(IVEG) .GT. 0.000001 ) GO TO 300
FILTR = FILTR + CAPAC(IVEG)
WWW(1) = WWW(1) + CAPAC(IVEG) / ( POROS*ZDEPTH(1) )
CAPAC(IVEG) = 0.
300 CONTINUE
1000 CONTINUE
DO 7000 IVEG = 1, 2
CCT = CCX
TS = TC
DTS = DTC
FLUX = CHF
IF ( IVEG .EQ. 1 ) GO TO 7100
CCT = CG
TS = TGS
DTS = DTG
FLUX = CCT * DTG / DTT
7100 CONTINUE
TTA = TS - DTS
TTB = TS
SNOWW(IVEG) = 0.
IF ( TTA .LE. TF ) SNOWW(IVEG) = CAPAC(IVEG)
CAPAC(IVEG) = CAPAC(IVEG) - SNOWW(IVEG)
IF ( TTA .GT. TF .AND. TTB .GT. TF ) GO TO 7200
IF ( TTA .LE. TF .AND. TTB .LE. TF ) GO TO 7200
DTF = TF - TTA
DTIME1 = CCT * DTF / FLUX
HF = FLUX*(DTT-DTIME1)
FCAP = - CAPAC(IVEG) * SNOMEL
SPWET = AMIN1( 5. , SNOWW(IVEG) )
IF ( DTS .GT. 0. ) FCAP = SPWET * SNOMEL
DTIME2 = FCAP / FLUX
DTF2 = FLUX * (DTT-DTIME1-DTIME2)/CCT
TN = TF + DTF2
TS = TF - 0.1
IF (ABS(HF) .GE.ABS(FCAP) ) TS = TN
CHANGE = HF
IF (ABS(CHANGE) .GE.ABS(FCAP) ) CHANGE = FCAP
CHANGE = CHANGE / SNOMEL
IF (CHANGE.GT.0.0) SMELT=CHANGE+SMELT
SNOWW(IVEG) = SNOWW(IVEG) - CHANGE
CAPAC(IVEG) = CAPAC(IVEG) + CHANGE
IF ( IVEG .EQ. 1 ) TC = TS
IF ( IVEG .EQ. 2 ) TGS = TS
IF ( SNOWW(IVEG) .LT. 0.00001 ) GO TO 7200
ZMELT = 0.
ZMELT = CAPAC(IVEG)
FILTR = FILTR+ ZMELT
WWW(1) = WWW(1) + ZMELT / ( POROS * ZDEPTH(1) )
CAPAC(IVEG) = 0.
7200 CONTINUE
CAPAC(IVEG) = CAPAC(IVEG) + SNOWW(IVEG)
7000 CONTINUE
FLUXEF = SHF - CCT*DTG/DTT
TD = TD + FLUXEF / ( CG * 2. * SQRT ( PIE*365. ) ) * DTT
change=0.0
FILTR = FILTR - EGS / HLAT / 1000.
WWW(1) = WWW(1) - EGS / HLAT / 1000. / ( POROS * ZDEPTH(1) )
DO 2000 IVEG = 1, 2
IF ( IVEG .EQ. 1 ) ABSOIL = ECT / HLAT / 1000.
IF ( IVEG .EQ. 2 ) ABSOIL = EGT / HLAT / 1000.
IF (NROOT.EQ.1) THEN
EF(2) = 0.
EF(3) = 0.
TOTDEP = ZDEPTH(1)
DO 3000 IL = 2, 3
TOTDEP = TOTDEP + ZDEPTH(IL)
IF ( ROOTD(IVEG) .LT. TOTDEP ) GO TO 400
EF(IL) = ZDEPTH(IL) / ROOTD(IVEG)
GO TO 500
400 CONTINUE
EF(IL) = ROOTD(IVEG) - TOTDEP + ZDEPTH(IL)
EF(IL) = EF(IL) / ROOTD(IVEG)
GO TO 600
500 CONTINUE
3000 CONTINUE
600 EFT = EF(2) + EF(3)
EFT = MAX(EFT,0.1E-5)
EF(2) = EF(2) / EFT
EF(3) = EF(3) / EFT
DO 4000 IL = 2, 3
WWW(IL) = WWW(IL) - ABSOIL * EF(IL) / ( POROS * ZDEPTH(IL) )
4000 CONTINUE
ELSE
EF(1) = ROOTP(1)
EF(2) = ROOTP(2)
EF(3) = ROOTP(3)
DO 4004 IL = 1, 3
WWW(IL) = WWW(IL) - ABSOIL * EF(IL) / ( POROS * ZDEPTH(IL) )
4004 CONTINUE
END IF
2000 CONTINUE
DO 5000 IL = 1, 2
IF ( WWW(IL) .GT. 0. ) GO TO 700
WWW(IL+1) = WWW(IL+1) + WWW(IL) * ZDEPTH(IL)/ZDEPTH(IL+1)
WWW(IL) = 0.
700 CONTINUE
5000 CONTINUE
do 8000 i = 1, 3
TEMW(I) = AMAX1( 0.03, WWW(I) )
TEMWP(I) = TEMW(I) ** ( -BEE )
TEMWPP(I) = AMIN1( 1., TEMW(I)) ** ( 2.*BEE+ 3. )
8000 CONTINUE
POWS = 2.*BEE+2.
Q3G = TEMW(3)**(-POWS) + SATCO/ZDEPTH(3)/POROS*SLOPE*POWS*DTT
Q3G = Q3G ** ( 1. / POWS )
Q3G = - ( 1. / Q3G - WWW(3) ) * POROS * ZDEPTH(3) / DTT
Q3G = AMAX1( 0., Q3G )
Q3G = AMIN1( Q3G, WWW(3)*POROS*ZDEPTH(3)/DTT )
Q3G = Q3G + 0.002*POROS*ZDEPTH(3)*0.5 / 86400. * WWW(3)
WMAX = AMAX1( WWW(1), WWW(2), WWW(3), 0.05 )
WMAX = AMIN1( WMAX, 1. )
PMAX = WMAX**(-BEE)
WMIN = (PMAX-2./( PHSAT*(ZDEPTH(1)+2.*ZDEPTH(2)+ZDEPTH(3)))) &
**(-1./BEE)
WMIN = AMIN1( WWW(1), WWW(2), WWW(3), WMIN )
WMIN = AMAX1( WMIN, 0.02 )
PMIN = WMIN**(-BEE)
DPDW = PHSAT*( PMAX-PMIN )/( WMAX-WMIN )
DO 8200 I = 1, 2
RSAME = 0.
AVK = TEMWP(I)*TEMWPP(I) - TEMWP(I+1)*TEMWPP(I+1)
DIV = TEMWP(I+1) - TEMWP(I)
IF ( ABS(DIV) .LT. 1.E-6 ) RSAME = 1.
AVK = SATCO*AVK / ( ( 1. + 3./BEE ) * DIV + RSAME )
AVKMIN = SATCO * AMIN1( TEMWPP(I), TEMWPP(I+1) )
AVKMAX = SATCO * AMAX1( TEMWPP(I), TEMWPP(I+1) )*1.01
AVK = AMAX1( AVK, AVKMIN )
AVK = AMIN1( AVK, AVKMAX )
TSNOW = AMIN1 ( TF-0.01, TGS )
AREAS = AMIN1 (0.999,13.2*SNOWW(2))
TGG = TSNOW*AREAS + TGS*(1.-AREAS)
TS = TGG*(2-I) + TD*(I-1)
PROPS = ( TS-(TF-10.) ) / 10.
PROPS = AMAX1( 0.05, AMIN1( 1.0, PROPS ) )
AVK = AVK * PROPS
Q3G = Q3G * PROPS
DPDWDZ = DPDW * 2./( ZDEPTH(I) + ZDEPTH(I+1) )
AAA(I) = 1. + AVK*DPDWDZ*( 1./ZDEPTH(I)+1./ZDEPTH(I+1) ) &
*DTT/POROS
BBB(I) =-AVK * DPDWDZ * 1./ZDEPTH(2)*DTT/POROS
CCC(I) = AVK * ( DPDWDZ * ( WWW(I)-WWW(I+1) ) + 1. + &
(I-1)*DPDWDZ*Q3G*1./ZDEPTH(3)*DTT/POROS )
8200 CONTINUE
DENOM = ( AAA(1)*AAA(2) - BBB(1)*BBB(2) )
RDENOM = 0.
IF ( ABS(DENOM) .LT. 1.E-6 ) RDENOM = 1.
RDENOM = ( 1.-RDENOM)/( DENOM + RDENOM )
QQQ(1) = ( AAA(2)*CCC(1) - BBB(1)*CCC(2) ) * RDENOM
QQQ(2) = ( AAA(1)*CCC(2) - BBB(2)*CCC(1) ) * RDENOM
WWW(3) = WWW(3) - Q3G*DTT/(POROS*ZDEPTH(3))
ROFF = ROFF + Q3G * DTT
DO 8300 I = 1, 2
QMAX = WWW(I) * (POROS*ZDEPTH(I) /DTT)
QMIN = -WWW(I+1) * (POROS*ZDEPTH(I+1)/DTT)
QQQ(I) = AMIN1( QQQ(I),QMAX)
QQQ(I) = AMAX1( QQQ(I),QMIN)
WWW(I) = WWW(I) - QQQ(I)/(POROS*ZDEPTH(I) /DTT)
WWW(I+1) = WWW(I+1) + QQQ(I)/(POROS*ZDEPTH(I+1)/DTT)
8300 CONTINUE
SOILDIF=SOILDIF+QQQ(1)*DTT*1000.
SOILDRA=SOILDRA+Q3G*DTT*1000.
DO 8400 I = 1, 3
EXCESS = AMAX1(0.,(WWW(I) - 1.))
WWW(I) = WWW(I) - EXCESS
ROFF = ROFF + EXCESS * POROS*ZDEPTH(I)
IF (I.LT.2) THEN
RNOFFS= RNOFFS+ 1000.*EXCESS*POROS*ZDEPTH(I)
ELSE
RNOFFB= RNOFFB+ 1000.*EXCESS*POROS*ZDEPTH(I)
ENDIF
8400 CONTINUE
DO 8402 I = 1,2
DEFICIT = AMAX1 (0.,(1.E-12 - WWW(I)))
IF (I.EQ.1) SOILDIF=SOILDIF-DEFICIT* &
ZDEPTH(1)*POROS
WWW (I) = WWW(I) + DEFICIT
WWW (I+1) = WWW(I+1) - DEFICIT * ZDEPTH(I) / ZDEPTH (I+1)
8402 CONTINUE
WWW(3) = AMAX1 (WWW(3),1.E-12)
800 CONTINUE
IF (WWW(1) .GT.1.) THEN
WWW(2) = WWW(2) + (WWW(1)-1.) * ZDEPTH(1)/ &
ZDEPTH(2)
SOILDIF=SOILDIF+(WWW(1)-1.)*ZDEPTH(1) &
*POROS*1000.
WWW(1) = 1.
END IF
If (WWW(2) .GT.1.) WWW(3) = WWW(3) + (WWW(2)-1.) * &
ZDEPTH(2) / ZDEPTH(3)
IF (WWW(2) .GT.1.) WWW(2) = 1.
IF (WWW(3) .GT.1.) THEN
ROFF = ROFF + (WWW(3)-1.)*POROS*ZDEPTH(3)
RNOFFB=RNOFFB+((WWW(3)-1.)*ZDEPTH(3)* &
POROS*1000.)
WWW(3) = 1.
END IF
END SUBROUTINE UPDAT1_ICE
SUBROUTINE CONVDIM(IOFLAG, &
DZO1,WO1,TSSN1,TSSNO1,BWO1,BTO1,CTO1,FIO1,FLO1,BIO1,BLO1,HO1, &
DZO2,WO2,TSSN2,TSSNO2,BWO2,BTO2,CTO2,FIO2,FLO2,BIO2,BLO2,HO2, &
DZO3,WO3,TSSN3,TSSNO3,BWO3,BTO3,CTO3,FIO3,FLO3,BIO3,BLO3,HO3, &
DZO4,WO4,TSSN4,TSSNO4,BWO4,BTO4,CTO4,FIO4,FLO4,BIO4,BLO4,HO4, &
DZO, WO, TSSN, TSSNO, BWO, BTO, CTO, FIO, FLO, BIO, BLO, HO )
REAL, DIMENSION (4) :: DZO,WO,TSSN,TSSNO,BWO,BTO,CTO,FIO,FLO,BIO,BLO,HO
IF (IOFLAG.EQ.0) THEN
DZO (1) = DZO1
WO (1) = WO1
TSSN (1) = TSSN1
TSSNO (1) = TSSNO1
BWO (1) = BWO1
BTO (1) = BTO1
CTO (1) = CTO1
FIO (1) = FIO1
FLO (1) = FLO1
BIO (1) = BIO1
BLO (1) = BLO1
HO (1) = HO1
DZO (2) = DZO2
WO (2) = WO2
TSSN (2) = TSSN2
TSSNO (2) = TSSNO2
BWO (2) = BWO2
BTO (2) = BTO2
CTO (2) = CTO2
FIO (2) = FIO2
FLO (2) = FLO2
BIO (2) = BIO2
BLO (2) = BLO2
HO (2) = HO2
DZO (3) = DZO3
WO (3) = WO3
TSSN (3) = TSSN3
TSSNO (3) = TSSNO3
BWO (3) = BWO3
BTO (3) = BTO3
CTO (3) = CTO3
FIO (3) = FIO3
FLO (3) = FLO3
BIO (3) = BIO3
BLO (3) = BLO3
HO (3) = HO3
DZO (4) = DZO4
WO (4) = WO4
TSSN (4) = TSSN4
TSSNO (4) = TSSNO4
BWO (4) = BWO4
BTO (4) = BTO4
CTO (4) = CTO4
FIO (4) = FIO4
FLO (4) = FLO4
BIO (4) = BIO4
BLO (4) = BLO4
HO (4) = HO4
ELSEIF (IOFLAG.EQ.1) THEN
DZO1 = DZO(1)
WO1 = WO(1)
TSSN1 = TSSN(1)
TSSNO1 = TSSNO(1)
BWO1 = BWO(1)
BTO1 = BTO(1)
CTO1 = CTO(1)
FIO1 = FIO(1)
FLO1 = FLO(1)
BIO1 = BIO(1)
BLO1 = BLO(1)
HO1 = HO(1)
DZO2 = DZO(2)
WO2 = WO(2)
TSSN2 = TSSN(2)
TSSNO2 = TSSNO(2)
BWO2 = BWO(2)
BTO2 = BTO(2)
CTO2 = CTO(2)
FIO2 = FIO(2)
FLO2 = FLO(2)
BIO2 = BIO(2)
BLO2 = BLO(2)
HO2 = HO(2)
DZO3 = DZO(3)
WO3 = WO(3)
TSSN3 = TSSN(3)
TSSNO3 = TSSNO(3)
BWO3 = BWO(3)
BTO3 = BTO(3)
CTO3 = CTO(3)
FIO3 = FIO(3)
FLO3 = FLO(3)
BIO3 = BIO(3)
BLO3 = BLO(3)
HO3 = HO(3)
DZO4 = DZO(4)
WO4 = WO(4)
TSSN4 = TSSN(4)
TSSNO4 = TSSNO(4)
BWO4 = BWO(4)
BTO4 = BTO(4)
CTO4 = CTO(4)
FIO4 = FIO(4)
FLO4 = FLO(4)
BIO4 = BIO(4)
BLO4 = BLO(4)
HO4 = HO(4)
ELSE
print*,'something wrong in CONVDIM',IOFLAG
STOP
ENDIF
END SUBROUTINE CONVDIM
END MODULE module_sf_ssib