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( IX, JX, 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, UV10, & 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 INTEGER, DIMENSION (20) :: IVMODIS REAL, DIMENSION (13) :: TD_DEPTH INTEGER :: sw_physics CHARACTER(LEN=*), INTENT(IN ) :: MMINLU CHARACTER*256 :: message 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/ DATA IVMODIS / 4, 1, 5, 2, 3, 8, 9, 6, & 6, 7, 7, 12, 11, 12, 13, 11, & 0, 10, 10, 10/ IF(MMINLU.EQ.'SSIB') THEN ITYPE=IVGTYP ELSEIF(MMINLU.EQ.'USGS') THEN ITYPE=IVUSGS(IVGTYP) ELSEIF(MMINLU.EQ.'MODIS') THEN ITYPE=IVMODIS(IVGTYP) ELSE IF (MMINLU .EQ. 'MODIFIED_IGBP_MODIS_NOAH') THEN ITYPE=IVMODIS(IVGTYP) ELSE CALL wrf_error_fatal3("",1013,& '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("",1020,& '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) TOTWB = WWW(1) * POROS * ZDEPTH(1) & + WWW(2) * POROS * ZDEPTH(2) & + WWW(3) * POROS * ZDEPTH(3) & + CAPAC(1) + CAPAC(2) 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,UV10) 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,UV10) 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 ENDWB = WWW(1) * POROS * ZDEPTH(1) & + WWW(2) * POROS * ZDEPTH(2) & + WWW(3) * POROS * ZDEPTH(3) & + CAPAC(1) + CAPAC(2) - (PPC+PPL)/1000. + ETMASS/1000. + ROFF ERRW = TOTWB - ENDWB IF(ABS(ERRW) .GT. 0.0001) THEN WRITE(message,*) 'SSIB WATER BALANCE WARNING: ',ERRW CALL wrf_message ( message ) ENDIF ZLHS = RADT(1) + RADT(2) - CHF - SHF ZRHS = HFLUX + (ECT + ECI + EGT + EGI + EGS)/DTT ERRH = ZLHS - ZRHS IF(ABS(ERRH) .GT. 1.) THEN WRITE(message,*) 'SSIB ENERGY BALANCE WARNING: ',ERRH CALL wrf_message ( message ) ENDIF 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 & ( IX, JX, 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, UV10, & 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,UV10) 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 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,UV10) 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 USTAR =UM*CU RAF = CTI / USTAR IF (RAF.LT.0.80) THEN RAF = 0.80 CTI = RAF*USTAR ENDIF CT = 1./CTI RA = RAF UEST = USTAR DRAG = RHOA * UEST*UEST Z2 = Z22 CUNI10 = ALOG((z2+10-D)/Z0)/VKC UV10 = USTAR * (CUNI10+FVV) 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,UV10) 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,UV10) 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,UV10) 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,UV10) 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,UV10) 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,UV10) 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