11 use machine ,
only : kind_phys, kind_dbl_prec
13 use physcons,
only : rhowater, con_t0c, con_hfus, con_hvap, &
14 con_pi, con_rv, con_g, con_csol, con_tice
25 real (kind_phys),
parameter :: tfrz = con_t0c
26 real (kind_phys),
parameter :: xls = con_hvap + con_hfus
27 real (kind_phys),
parameter :: piconst = con_pi
28 real (kind_phys),
parameter :: r_v = con_rv
29 real (kind_phys),
parameter :: grav = con_g
30 real (kind_phys),
parameter :: sheatice = con_csol
32 real (kind_phys),
parameter :: rhoice = 917._kind_phys
33 real (kind_phys),
parameter :: sheatsn = 2090._kind_phys
34 real (kind_phys),
parameter :: P1000mb = 100000._kind_phys
36 real (kind_phys),
parameter :: zero = 0._kind_dbl_prec
37 real (kind_phys),
parameter :: one = 1._kind_dbl_prec
59 integer,
PARAMETER :: NLUS=50
66 INTEGER,
PARAMETER :: NSLTYPE=30
73 INTEGER,
PARAMETER :: NSLOPE=30
74 real (kind_phys) :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
75 REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
89 DT,init,lsm_cold_start,KTAU,iter,NSL, &
90 graupelncv,snowncv,rainncv,raincv, &
91 ZS,RAINBL,SNOW,SNOWH,SNOWC,FRZFRAC,frpcpn, &
92 rhosnf,precipfr,exticeden, hgt,stdev, &
93 Z3D,P8W,T3D,QV3D,QC3D,RHO3D,EMISBCK, &
94 GLW,GSWdn,GSW,EMISS,CHKLOWQ, CHS, &
95 FLQC,FLHC,rhonewsn_ex,mosaic_lu, &
96 mosaic_soil,isncond_opt,isncovr_opt, &
97 MAVAIL,CANWAT,VEGFRA, &
98 ALB,ZNT,Z0,SNOALB,ALBBCK,LAI, &
99 landusef, nlcat, soilctop, nscat, &
101 QSFC,QSG,QVG,QCG,DEW,SOILT1,TSNAV, &
102 TBOT,IVGTYP,ISLTYP,XLAND, &
103 ISWATER,ISICE,XICE,XICE_THRESHOLD, &
104 CP,RV,RD,G0,PI,LV,STBOLT, &
105 SOILMOIS,SH2O,SMAVAIL,SMMAX, &
106 TSO,SOILT,EDIR,EC,ETT,SUBLIM,SNOH, &
107 HFX,QFX,LH,INFILTR, &
108 RUNOFF1,RUNOFF2,ACRUNOFF,SFCEXC, &
109 SFCEVP,GRDFLX,SNOWFALLAC,ACSNOW,SNOM, &
110 SMFR3D,KEEPFR3DFLAG, &
111 add_fire_heat_flux,fire_heat_flux, &
112 myj,shdmin,shdmax,rdlai2d, &
113 ims,ime, jms,jme, kms,kme, &
114 its,ite, jts,jte, kts,kte, &
199 real (kind_phys),
INTENT(IN ) :: xlat,xlon
200 real (kind_phys),
INTENT(IN ) :: dt
201 LOGICAL,
INTENT(IN ) :: myj,frpcpn,init,lsm_cold_start,exticeden
202 INTEGER,
INTENT(IN ) :: nlcat, nscat
203 INTEGER,
INTENT(IN ) :: mosaic_lu,mosaic_soil
204 INTEGER,
INTENT(IN ) :: isncond_opt,isncovr_opt
205 INTEGER,
INTENT(IN ) :: ktau, iter, nsl, isice, iswater, &
206 ims,ime, jms,jme, kms,kme, &
207 its,ite, jts,jte, kts,kte
211 real (kind_phys),
DIMENSION( ims:ime, kms:kme, jms:jme ) , &
212 INTENT(IN ) :: qv3d, &
219 real (kind_phys),
DIMENSION( ims:ime , jms:jme ), &
220 INTENT(IN ) :: rainbl, &
233 real (kind_phys),
DIMENSION( ims:ime , jms:jme ), &
234 INTENT(IN ) :: graupelncv, &
238 real (kind_phys),
DIMENSION( ims:ime),
INTENT(IN ) :: rhonewsn_ex
240 real (kind_phys),
DIMENSION( ims:ime , jms:jme ),
INTENT(IN ):: shdmax
241 real (kind_phys),
DIMENSION( ims:ime , jms:jme ),
INTENT(IN ):: shdmin
242 real (kind_phys),
DIMENSION( ims:ime , jms:jme ),
INTENT(IN ):: hgt
243 real (kind_phys),
DIMENSION( ims:ime , jms:jme ),
INTENT(IN ):: stdev
244 LOGICAL,
intent(in) :: add_fire_heat_flux
245 real (kind_phys),
DIMENSION( ims:ime , jms:jme ),
INTENT(IN ):: fire_heat_flux
246 LOGICAL,
intent(in) :: rdlai2d
248 real (kind_phys),
DIMENSION( 1:nsl),
INTENT(IN ) :: zs
250 real (kind_phys),
DIMENSION( ims:ime , jms:jme ), &
268 real (kind_phys),
DIMENSION( ims:ime , jms:jme ), &
272 INTEGER,
DIMENSION( ims:ime , jms:jme ), &
273 INTENT(IN ) :: ivgtyp, &
275 real (kind_phys),
DIMENSION( ims:ime , 1:nlcat, jms:jme ),
INTENT(IN):: landusef
276 real (kind_phys),
DIMENSION( ims:ime , 1:nscat, jms:jme ),
INTENT(IN):: soilctop
278 real (kind_phys),
INTENT(IN ) :: cp,g0,lv,stbolt,rv,rd,pi, &
281 real (kind_phys),
DIMENSION( ims:ime , 1:nsl, jms:jme ) , &
282 INTENT(INOUT) :: soilmois,sh2o,tso
284 real (kind_phys),
DIMENSION( ims:ime, jms:jme ) , &
285 INTENT(INOUT) :: soilt, &
310 real (kind_phys),
DIMENSION( ims:ime, jms:jme ) , &
311 INTENT(INOUT) :: smavail, &
314 real (kind_phys),
DIMENSION( its:ite, jts:jte ) :: &
332 real (kind_phys),
DIMENSION( its:ite, jts:jte ) :: &
342 real (kind_phys),
DIMENSION( ims:ime, 1:nsl, jms:jme) &
346 real (kind_phys),
DIMENSION( ims:ime, jms:jme ),
INTENT(OUT) :: &
369 real (kind_phys) :: cn, &
378 real (kind_phys),
DIMENSION(1:NSL) :: zsmain, &
382 real (kind_phys),
DIMENSION(1:2*(nsl-2)) :: dtdzs
384 real (kind_phys),
DIMENSION(1:5001) :: tbq
387 real (kind_phys),
DIMENSION( 1:nsl ) :: soilm1d, &
393 real (kind_phys),
DIMENSION( 1:nsl ) :: keepfr
395 real (kind_phys),
DIMENSION( 1:nlcat ) :: lufrac
396 real (kind_phys),
DIMENSION( 1:nscat ) :: soilfrac
398 real (kind_phys) :: rsm, &
402 real (kind_phys) :: prcpms, &
423 real (kind_phys) :: cq,r61,r273,arp,brp,x,evs,eis
424 real (kind_phys) :: cropfr, cropsm, newsm, factor
426 real (kind_phys) :: meltfactor, ac,as, wb,rovcp
428 INTEGER :: iland,isoil,iforest
430 INTEGER :: i,j,k,nzs,nzs1,nddzs
432 logical :: debug_print
435 real (kind_phys) :: testptlat, testptlon
437 character(len=*),
intent(out) :: errmsg
438 integer,
intent(out) :: errflg
446 debug_print = .false.
460 cq=173.15_kind_dbl_prec-.05_kind_dbl_prec
461 r273=1._kind_dbl_prec/tfrz
462 r61=6.1153_kind_dbl_prec*0.62198_kind_dbl_prec
463 arp=77455._kind_dbl_prec*41.9_kind_dbl_prec/461.525_kind_dbl_prec
464 brp=64._kind_dbl_prec*41.9_kind_dbl_prec/461.525_kind_dbl_prec
467 cq=cq+.05_kind_dbl_prec
468 evs=exp(17.67_kind_dbl_prec*(cq-tfrz)/(cq-29.65_kind_dbl_prec))
469 eis=exp(22.514_kind_dbl_prec-6.15e3_kind_dbl_prec/cq)
481 if(init .and. iter == 1)
then
483 if( lsm_cold_start )
then
489 IF((soilt1(i,j) .LT. 170._kind_phys) .or. (soilt1(i,j) .GT.400._kind_phys))
THEN
490 IF(snowc(i,j).gt.zero)
THEN
491 soilt1(i,j)=min(tfrz,0.5_kind_phys*(soilt(i,j)+tso(i,1,j)) )
492 IF (debug_print )
THEN
494 'Temperature inside snow is initialized in RUCLSM ', soilt1(i,j),i,xlat,xlon
497 soilt1(i,j) = tso(i,1,j)
500 tsnav(i,j) =min(zero,0.5_kind_phys*(soilt(i,j)+tso(i,1,j))-tfrz)
503 if(hgt(i,j) > 2500._kind_phys)
then
504 snoalb(i,j) = min(0.65_kind_phys,snoalb(i,j))
507 patmb=p8w(i,kms,j)*1.e-2_kind_phys
508 qsg(i,j) = qsn(soilt(i,j),tbq)/patmb
510 if((qcg(i,j) < zero) .or. (qcg(i,j) > 0.1_kind_phys))
then
511 qcg(i,j) = qc3d(i,1,j)
512 if (debug_print )
then
513 print *,
'QCG is initialized in RUCLSM ', qcg(i,j),qc3d(i,1,j),i,xlat,xlon
517 if((qvg(i,j) .LE. zero) .or. (qvg(i,j) .GT.0.1_kind_phys))
then
518 qvg(i,j) = qv3d(i,1,j)
519 if (debug_print )
then
520 print *,
'QVG is initialized in RUCLSM ', qvg(i,j),mavail(i,j),qsg(i,j),i,xlat,xlon
523 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
528 snowfallac(i,j) = zero
530 rhosnf(i,j) = -1.e3_kind_phys
537 sfcrunoff(i,j) = zero
545 waterbudget(i,j) = zero
546 acwaterbudget(i,j) = zero
588 IF (debug_print )
THEN
589 if (abs(xlat-testptlat).lt.0.2 .and. &
590 abs(xlon-testptlon).lt.0.2)
then
591 print 100,
'(RUC start) i=',i,
' lat,lon=',xlat,xlon, &
592 'mavail ', mavail(i,j),
' soilt',soilt(i,j),
'qvg ',qvg(i,j),&
593 'p8w',p8w(i,1,j),
'sflay qfx',qfx(i,j),
'sflay hfx',hfx(i,j),&
594 'gsw ',gsw(i,j),
'glw ',glw(i,j),
'soilt ',soilt(i,j), &
595 'chs ',chs(i,j),
'flqc ',flhc(i,j),
'alb ',alb(i,j), &
596 'rainbl ',rainbl(i,j),
'dt ',dt
597 print *,
'nzs',nzs,
'ivgtyp ',ivgtyp(i,j),
'isltyp ',isltyp(i,j)
604 qvatm = qv3d(i,kms,j)
605 qcatm = qc3d(i,kms,j)
606 patm = p8w(i,kms,j)*1.e-5_kind_phys
610 conflx = z3d(i,kms,j)*0.5_kind_phys
618 prcpncliq = rainncv(i,j)*(1.-frzfrac(i,j))
619 prcpncfr = rainncv(i,j)*frzfrac(i,j)
623 if(frzfrac(i,j) > zero .and. tabs < tfrz)
then
624 prcpculiq = max(zero,raincv(i,j)*(one-frzfrac(i,j)))
625 prcpcufr = max(zero,raincv(i,j)*frzfrac(i,j))
628 prcpcufr = max(zero,raincv(i,j))
632 prcpculiq = max(zero,raincv(i,j))
636 prcpms = (prcpncliq + prcpculiq)/dt*1.e-3_kind_phys
637 newsnms = (prcpncfr + prcpcufr)/dt*1.e-3_kind_phys
639 if((prcpncfr + prcpcufr) > zero)
then
641 snowrat=min(one,max(zero,snowncv(i,j)/(prcpncfr + prcpcufr)))
642 grauprat=min(one,max(zero,graupelncv(i,j)/(prcpncfr + prcpcufr)))
643 icerat=min(one,max(zero,(prcpncfr-snowncv(i,j)-graupelncv(i,j)) &
644 /(prcpncfr + prcpcufr)))
645 curat=min(one,max(zero,(prcpcufr/(prcpncfr + prcpcufr))))
649 if (tabs.le.tfrz)
then
651 newsnms = rainbl(i,j)/dt*1.e-3_kind_phys
656 prcpms = rainbl(i,j)/dt*1.e-3_kind_phys
663 precipfr(i,j) = newsnms * dt *1.e3_kind_phys
670 qkms=flqc(i,j)/rho/mavail(i,j)
672 tkms=flhc(i,j)/rho/(cp*(one+0.84_kind_phys*qvatm))
675 snwe=snow(i,j)*1.e-3_kind_phys
677 canwatr=canwat(i,j)*1.e-3_kind_phys
680 rhosnfall=rhosnf(i,j)
688 zshalf(k)=0.5_kind_phys*(zsmain(k-1) + zsmain(k))
692 lufrac(k) = landusef(i,k,j)
695 soilfrac(k) = soilctop(i,k,j)
703 IF (debug_print )
THEN
704 print *,
' DT,NZS1, ZSMAIN, ZSHALF --->', dt,nzs1,zsmain,zshalf
710 x=dt/2./(zshalf(k+1)-zshalf(k))
711 dtdzs(k1)=x/(zsmain(k)-zsmain(k-1))
713 dtdzs(k2)=x/(zsmain(k+1)-zsmain(k))
716 cw =4.183e6_kind_dbl_prec
722 kqwrtz=7.7_kind_dbl_prec
723 kice=2.2_kind_dbl_prec
724 kwt=0.57_kind_dbl_prec
729 c1sn=0.026_kind_dbl_prec
730 c2sn=21._kind_dbl_prec
737 rhonewsn = 200._kind_phys
738 if(snow(i,j).gt.zero .and. snowh(i,j).gt.0.02_kind_phys)
then
739 rhosn = snow(i,j)/snowh(i,j)
741 rhosn = 300._kind_phys
744 IF (debug_print )
THEN
746 if (abs(xlat-testptlat).lt.0.2 .and. &
747 abs(xlon-testptlon).lt.0.2)
then
748 print*,
' lat,lon=',xlat,xlon
749 print *,
'before SOILVEGIN - z0,znt',i,z0(i,j),znt(i,j)
750 print *,
'ILAND, ISOIL =',i,iland,isoil
757 CALL soilvegin ( debug_print, mosaic_lu, mosaic_soil, &
758 soilfrac,nscat,shdmin(i,j),shdmax(i,j), &
759 nlcat,iland,isoil,iswater,myj,iforest,lufrac,vegfra(i,j), &
760 emissl(i,j),pc(i,j),msnf(i,j),facsnf(i,j), &
761 znt(i,j),lai(i,j),rdlai2d, &
762 qwrtz,rhocs,bclh,dqm,ksat,psis,qmin,ref,wilt,i,j,errmsg, errflg)
765 emisbck(i,j) = emissl(i,j)
769 IF (debug_print )
THEN
771 if (abs(xlat-testptlat).lt.0.2 .and. &
772 abs(xlon-testptlon).lt.0.2)
then
773 print*,
' lat,lon=',xlat,xlon
774 print *,
'after SOILVEGIN - z0,znt,lai',i,z0(i,j),znt(i,j),lai(i,j)
775 print *,
'NLCAT,iland,EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J)', &
776 nlcat,iland,emissl(i,j),pc(i,j),znt(i,j),lai(i,j),i,j
777 print *,
'NSCAT,QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT',&
778 nscat,qwrtz,rhocs,bclh,dqm,ksat,psis,qmin,ref,wilt,i,j
784 sat = 5.e-4_kind_phys
787 IF(iforest.gt.2)
THEN
793 meltfactor = 2.0_kind_phys
796 if(zsmain(k).ge.0.4_kind_phys)
then
808 meltfactor = 0.85_kind_phys
811 if(zsmain(k).ge.1.1_kind_phys)
then
819 IF (debug_print )
THEN
820 if (abs(xlat-testptlat).lt.0.2 .and. &
821 abs(xlon-testptlon).lt.0.2)
then
822 print*,
' lat,lon=',xlat,xlon
823 print *,
' ZNT, LAI, VEGFRA, SAT, EMIS, PC --->', &
824 znt(i,j),lai(i,j),vegfra(i,j),sat,emissl(i,j),pc(i,j)
825 print *,
' ZS, ZSMAIN, ZSHALF, CONFLX, CN, SAT, --->', zs,zsmain,zshalf,conflx,cn,sat
826 print *,
'NROOT, meltfactor, iforest, ivgtyp, i,j ', nroot,meltfactor,iforest,ivgtyp(i,j),i,j
830 IF((xland(i,j)-1.5).GE.0._kind_phys)
THEN
839 acsnow(i,j)=acsnow(i,j)+precipfr(i,j)
844 patmb=p8w(i,1,j)*1.e-2_kind_phys
845 qvg(i,j) = qsn(soilt(i,j),tbq)/patmb
846 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
848 q2sat=qsn(tabs,tbq)/patmb
853 tso(i,k,j)= soilt(i,j)
856 IF (debug_print )
THEN
857 if (abs(xlat-testptlat).lt.0.2 .and. &
858 abs(xlon-testptlon).lt.0.2)
then
859 print*,
' water point'
860 print*,
' lat,lon=',xlat,xlon,
'SOILT=', soilt(i,j)
867 if(xice(i,j).ge.xice_threshold)
then
873 IF(seaice(i,j).GT.0.5_kind_phys)
THEN
875 IF (debug_print )
THEN
876 if (abs(xlat-testptlat).lt.0.2 .and. &
877 abs(xlon-testptlon).lt.0.2)
then
878 print*,
' sea-ice at water point'
879 print*,
' lat,lon=',xlat,xlon
888 znt(i,j) = 0.011_kind_phys
890 emissl(i,j) = emisbck(i,j)
896 patmb=p8w(i,1,j)*1.e-2_kind_phys
897 qvg(i,j) = qsn(soilt(i,j),tbq)/patmb
899 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
902 soilmois(i,k,j) = one
905 keepfr3dflag(i,k,j) = zero
906 tso(i,k,j) = min(con_tice,tso(i,k,j))
915 soilm1d(k) = min(max(zero,soilmois(i,k,j)-qmin),dqm)
916 tso1d(k) = tso(i,k,j)
917 soiliqw(k) = min(max(zero,sh2o(i,k,j)-qmin),soilm1d(k))
918 soilice(k) =(soilm1d(k) - soiliqw(k))/0.9_kind_phys
922 smfrkeep(k) = smfr3d(i,k,j)
923 keepfr(k) = keepfr3dflag(i,k,j)
926 lmavail(i,j)=max(0.00001_kind_phys,min(one,soilm1d(1)/(ref-qmin)))
928 IF (debug_print )
THEN
929 if (abs(xlat-testptlat).lt.0.2 .and. &
930 abs(xlon-testptlon).lt.0.2)
then
931 print*,
' lat,lon=',xlat,xlon
932 print *,
'LAND, i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO', &
933 i,j,tso1d,soilm1d,patm,tabs,qvatm,qcatm,rho
934 print *,
'CONFLX =',conflx
935 print *,
'SMFRKEEP,KEEPFR ',smfrkeep,keepfr
942 smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(k))* &
943 (zshalf(k+1)-zshalf(k))
946 if (debug_print .and. abs(xlat-testptlat).lt.0.2 &
947 .and. abs(xlon-testptlon).lt.0.2)
then
948 print *,
'Old soilm1d ',i,soilm1d
951 canwatold(i,j) = canwatr
953 CALL sfctmp (debug_print, dt,ktau,conflx,i,j, &
954 xlat, xlon, testptlat, testptlon, &
956 nzs,nddzs,nroot,meltfactor, &
957 isncond_opt,isncovr_opt, &
958 iland,isoil,ivgtyp(i,j),isltyp(i,j), &
959 prcpms, newsnms,snwe,snhei,snowfrac, &
960 exticeden,rhosn,rhonewsn_ex(i),rhonewsn, &
961 rhosnfall,snowrat,grauprat,icerat,curat, &
962 patm,tabs,qvatm,qcatm,rho, &
963 glw(i,j),gswdn(i,j),gsw(i,j), &
964 emissl(i,j),emisbck(i,j), &
965 msnf(i,j), facsnf(i,j), &
966 qkms,tkms,pc(i,j),lmavail(i,j), &
967 canwatr,vegfra(i,j),alb(i,j),znt(i,j), &
968 snoalb(i,j),albbck(i,j),lai(i,j), &
969 hgt(i,j),stdev(i,j), &
970 myj,seaice(i,j),isice, &
971 add_fire_heat_flux,fire_heat_flux(i,j), &
974 rhocs,dqm,qmin,ref, &
975 wilt,psis,bclh,ksat, &
976 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
978 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
981 snweprint,snheiprint,rsm, &
982 soilm1d,tso1d,smfrkeep,keepfr, &
983 soilt(i,j),soilt1(i,j),tsnav(i,j),dew(i,j), &
984 qvg(i,j),qsg(i,j),qcg(i,j),smelt(i,j), &
985 snoh(i,j),snflx(i,j),snom(i,j),snowfallac(i,j), &
986 acsnow(i,j),edir(i,j),ec(i,j),ett(i,j),qfx(i,j), &
987 lh(i,j),hfx(i,j),
sflx(i,j),sublim(i,j), &
988 evapl(i,j),prcpl(i,j),budget(i,j),runoff1(i,j), &
989 runoff2(i,j),soilice,soiliqw,infiltrp,smf(i,j))
999 if(mosaic_lu == 1)
then
1001 factor = max(zero,min(one,(vegfra(i,j)-shdmin(i,j))/max(one,(shdmax(i,j)-shdmin(i,j)))))
1002 if (debug_print )
then
1003 if (abs(xlat-testptlat).lt.0.1 .and. &
1004 abs(xlon-testptlon).lt.0.1)
then
1005 print *,
' lat,lon=',xlat,xlon,
' factor=',factor
1009 if((ivgtyp(i,j) == natural .or. ivgtyp(i,j) == crop) .and. factor > 0.75)
then
1014 cropsm=1.05_kind_phys*wilt - qmin
1015 cropfr = min(one,lufrac(crop) + 0.4*lufrac(natural))
1016 newsm = cropsm*cropfr + (1.-cropfr)*soilm1d(k)
1017 if(soilm1d(k) < newsm)
then
1018 IF (debug_print )
THEN
1019 if (abs(xlat-testptlat).lt.0.1 .and. &
1020 abs(xlon-testptlon).lt.0.1)
then
1021 print * ,
'Soil moisture is below wilting in cropland areas at time step',ktau
1022 print * ,
' lat,lon=',xlat,xlon
1023 print * ,
' lufrac=',lufrac,
'factor=',factor &
1024 ,
'lai,ivgtyp,lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm,', &
1025 lai(i,j),ivgtyp(i,j),lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm
1029 IF (debug_print )
THEN
1030 if (abs(xlat-testptlat).lt.0.1 .and. &
1031 abs(xlon-testptlon).lt.0.1)
then
1032 print*,
' lat,lon=',xlat,xlon
1033 print * ,
'Added soil water to cropland areas, k,soilm1d(k)',k,soilm1d(k)
1051 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* &
1052 (zshalf(k+1)-zshalf(k))
1053 smmax(i,j) =smmax(i,j)+(qmin+dqm)* &
1054 (zshalf(k+1)-zshalf(k))
1057 if (debug_print)
then
1058 if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
1059 print 100,
'(RUC runoff) i=',i,
' lat,lon=',xlat,xlon, &
1060 'RUNOFF1', runoff1(i,j),
'RUNOFF2 ',runoff2(i,j), &
1061 'edir ',edir(i,j),
'ec ',ec(i,j),
'ett ',ett(i,j)
1067 acrunoff(i,j) = (runoff1(i,j)+runoff2(i,j))*dt*rhowater
1068 smavail(i,j) = smavail(i,j) * rhowater
1069 smmax(i,j) = smmax(i,j) * rhowater
1070 smtotold(i,j) = smtotold(i,j) * rhowater
1074 soilmois(i,k,j) = soilm1d(k) + qmin
1075 sh2o(i,k,j) = min(soiliqw(k) + qmin,soilmois(i,k,j))
1076 tso(i,k,j) = tso1d(k)
1080 smfr3d(i,k,j) = smfrkeep(k)
1081 keepfr3dflag(i,k,j) = keepfr(k)
1089 patmb=p8w(i,1,j)*1.e-2_kind_phys
1090 q2sat=qsn(tabs,tbq)/patmb
1091 qsfc(i,j) = qvg(i,j)/(one+qvg(i,j))
1095 IF((qvatm.GE.q2sat*0.95_kind_phys).AND.qvatm.LT.qvg(i,j))
THEN
1101 if(snow(i,j)==zero) emissl(i,j) = emisbck(i,j)
1102 emiss(i,j) = emissl(i,j)
1106 snow(i,j) = min(3._kind_phys,snwe)*1000._kind_phys
1107 snowh(i,j) = min(7.5_kind_phys,snhei)
1108 canwat(i,j) = canwatr*1000._kind_phys
1110 if (debug_print)
then
1111 if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
1112 print *,
'snow(i,j),soilt(i,j),xice(i,j),tso(i,:,j)',snow(i,j),soilt(i,j),xice(i,j),tso(i,:,j)
1115 infiltr(i,j) = infiltrp
1117 mavail(i,j) = lmavail(i,j)
1118 IF (debug_print )
THEN
1119 if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
1120 print *,
' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
1123 sfcevp(i,j) = sfcevp(i,j) + qfx(i,j) * dt
1124 grdflx(i,j) = -one *
sflx(i,j)
1134 rhosnf(i,j)=rhosnfall
1137 sfcevp(i,j) = sfcevp(i,j) + qfx(i,j) * dt
1145 if (debug_print )
then
1146 if (abs(xlat-testptlat).lt.0.2 .and. &
1147 abs(xlon-testptlon).lt.0.2)
then
1153 ac=canwat(i,j)-canwatold(i,j)*rhowater
1154 as=snwe-snowold(i,j)
1155 wb = smavail(i,j)-smtotold(i,j)
1156 waterbudget(i,j)=rainbl(i,j)+smelt(i,j)*dt*rhowater &
1158 -runoff1(i,j)*dt*rhowater-runoff2(i,j)*dt*rhowater &
1161 print *,
'soilm1d ',i,soilm1d
1162 print 100,
'(RUC budgets) i=',i,
' lat,lon=',xlat,xlon, &
1163 'budget ',budget(i,j),
'waterbudget',waterbudget(i,j), &
1164 'rainbl ',rainbl(i,j),
'runoff1 ',runoff1(i,j), &
1165 'smelt ',smelt(i,j)*dt*1.e3_kind_phys,
'smc change ',wb, &
1166 'snwe change ',as,
'canw change ',ac,
'runoff2 ',runoff2(i,j), &
1167 'qfx*dt ',qfx(i,j)*dt,
'smavail ',smavail(i,j),
'smcold',smtotold(i,j)
1169 print *,
'Smf=',smf(i,j),i,j
1170 print *,
'SNOW,SNOWold',i,j,snwe,snowold(i,j)
1171 print *,
'SNOW-SNOWold',i,j,max(zero,snwe-snowold(i,j))
1172 print *,
'CANWATold, canwat ',i,j,canwatold(i,j),canwat(i,j)
1173 print *,
'canwat(i,j)-canwatold(i,j)',max(zero,canwat(i,j)-canwatold(i,j))
1177 100
format (
";;; ",a,i4,a,2f14.7/(4(a10,
'='es14.7)))
1180 IF (debug_print )
THEN
1181 if (abs(xlat-testptlat).lt.0.2 .and. &
1182 abs(xlon-testptlon).lt.0.2)
then
1183 print *,
'LAND, i,tso1d,soilm1d,soilt - end of time step', &
1184 i,tso1d,soilm1d,soilt(i,j)
1185 print *,
'LAND, QFX, HFX after SFCTMP', i,lh(i,j),hfx(i,j)
1210 SUBROUTINE sfctmp (debug_print, delt,ktau,conflx,i,j, & !--- input variables
1211 xlat,xlon,testptlat,testptlon, &
1212 nzs,nddzs,nroot,meltfactor, &
1213 isncond_opt,isncovr_opt, &
1214 ILAND,ISOIL,IVGTYP,ISLTYP,PRCPMS, &
1215 NEWSNMS,SNWE,SNHEI,SNOWFRAC, &
1216 exticeden,RHOSN,RHONEWSN_ex,RHONEWSN,RHOSNFALL, &
1217 snowrat,grauprat,icerat,curat, &
1218 PATM,TABS,QVATM,QCATM,rho, &
1219 GLW,GSWdn,GSW,EMISS,EMISBCK,msnf,facsnf, &
1220 QKMS,TKMS,PC,MAVAIL,CST,VEGFRA,ALB,ZNT, &
1221 ALB_SNOW,ALB_SNOW_FREE,lai,hgt,stdev, &
1223 add_fire_heat_flux,fire_heat_flux, &
1224 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & !--- soil fixed fields
1225 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1226 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, & !--- constants
1228 snweprint,snheiprint,rsm, & !---output variables
1229 soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
1230 tsnav,dew,qvg,qsg,qcg, &
1231 SMELT,SNOH,SNFLX,SNOM,SNOWFALLAC,ACSNOW, &
1232 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
1233 evapl,prcpl,fltot,runoff1,runoff2,soilice, &
1234 soiliqw,infiltr,smf)
1241 INTEGER,
INTENT(IN ) :: isice,i,j,nroot,ktau,nzs , &
1243 integer,
intent(in ) :: isncond_opt,isncovr_opt
1245 real (kind_phys),
INTENT(IN ) :: delt,conflx,meltfactor,xlat,xlon
1246 real (kind_phys),
INTENT(IN ) :: testptlat,testptlon
1247 real (kind_phys),
INTENT(IN ) :: c1sn,c2sn,rhonewsn_ex
1248 LOGICAL,
INTENT(IN ) :: myj, debug_print, exticeden
1250 real (kind_phys) , &
1251 INTENT(IN ) :: patm, &
1255 real (kind_phys) , &
1256 INTENT(IN ) :: glw, &
1270 LOGICAL,
INTENT(IN ) :: add_fire_heat_flux
1272 INTEGER,
INTENT(IN ) :: ivgtyp, isltyp
1274 real (kind_phys) , &
1275 INTENT(INOUT) :: emiss, &
1284 real (kind_phys) :: &
1296 real (kind_phys),
INTENT(IN ) :: cn, &
1307 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
1312 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
1314 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: tbq
1319 real (kind_phys),
DIMENSION( 1:nzs ) , &
1320 INTENT(INOUT) :: ts1d, &
1323 real (kind_phys),
DIMENSION( 1:nzs ) , &
1324 INTENT(INOUT) :: keepfr
1326 real (kind_phys),
DIMENSION(1:NZS),
INTENT(INOUT) :: soilice, &
1330 INTEGER,
INTENT(INOUT) :: iland,isoil
1334 real (kind_phys) , &
1335 INTENT(INOUT) :: dew, &
1374 real (kind_phys),
DIMENSION(1:NZS) :: &
1386 real (kind_phys) :: &
1411 real (kind_phys),
INTENT(INOUT) :: rsm, &
1418 real (kind_phys) :: bsn, xsn , &
1419 rainf, snth, newsn, prcpms, newsnms , &
1420 t3, upflux, xinet, snowfrac2, m
1421 real (kind_phys) :: snhei_crit, snhei_crit_newsn, keep_snow_albedo, snowfracnewsn
1422 real (kind_phys) :: newsnowratio, dd1
1424 real (kind_phys) :: rhonewgr,rhonewice
1426 real (kind_phys) :: rnet,gswnew,gswin,emissn,zntsn,emiss_snowfree
1427 real (kind_phys) :: vegfrac, snow_mosaic, snfr, vgfr
1428 real :: cice, albice, albsn, drip, dripsn, dripliq
1429 real :: interw, intersn, infwater, intwratio
1432 integer,
parameter :: ilsnow=99
1434 IF (debug_print )
THEN
1435 print *,
' in SFCTMP',i,j,nzs,nddzs,nroot, &
1436 snwe,rhosn,snom,smelt,ts1d
1447 snhei_crit=0.01601_kind_phys*rhowater/rhosn
1448 snhei_crit_newsn=0.0005_kind_phys*rhowater/rhosn
1450 zntsn = z0tbl(isice)
1457 rhonewsn = 100._kind_phys
1458 if(snhei == zero) snowfrac=zero
1471 vegfrac=0.01_kind_phys*vegfra
1491 albice=alb_snow_free
1493 emissn = 0.99_kind_phys
1494 emiss_snowfree = emisbck
1500 if(seaice.ge.0.5_kind_phys)
then
1502 tice(k) = ts1d(k) - tfrz
1503 rhosice(k) = 917.6_kind_phys/(one-0.000165_kind_phys*tice(k))
1504 cice = 2115.85_kind_phys +7.7948_kind_phys*tice(k)
1505 capice(k) = cice*rhosice(k)
1506 thdifice(k) = 2.260872_kind_phys/capice(k)
1512 albice = min(alb_snow_free,max(alb_snow_free - 0.05_kind_phys, &
1513 alb_snow_free - 0.1_kind_phys*(tice(1)+10._kind_phys)/10._kind_phys ))
1516 IF (debug_print )
THEN
1517 print *,
'alb_snow_free',alb_snow_free
1518 print *,
'GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE',&
1519 gsw,gswnew,glw,soilt,emiss,alb,albice,snwe
1522 if(snhei.gt.0.0081_kind_phys*rhowater/rhosn)
then
1524 bsn=delt/3600._kind_phys*c1sn*exp(0.08_kind_phys*min(zero,tsnav)-c2sn*rhosn*1.e-3_kind_phys)
1525 if(bsn*snwe*100._kind_phys.ge.1.e-4_kind_phys)
then
1526 xsn=rhosn*(exp(bsn*snwe*100._kind_phys)-one)/(bsn*snwe*100._kind_phys)
1527 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
1532 if(snowfrac < 0.75_kind_phys) snow_mosaic = one
1536 acsnow=acsnow+newsn*rhowater
1538 IF(newsn.GT.zero)
THEN
1540 IF (debug_print )
THEN
1541 print *,
'THERE IS NEW SNOW, newsn', newsn
1544 newsnowratio = min(one,newsn/(snwe+newsn))
1557 rhonewsn = rhonewsn_ex
1559 rhonewsn=min(125._kind_phys,rhowater/max(8._kind_phys,(17._kind_phys*tanh((276.65_kind_phys-tabs)*0.15_kind_phys))))
1560 rhonewgr=min(500._kind_phys,rhowater/max(2._kind_phys,(3.5_kind_phys*tanh((274.15_kind_phys-tabs)*0.3333_kind_phys))))
1566 rhosnfall = min(500._kind_phys,max(58.8_kind_phys,(rhonewsn*snowrat + &
1567 rhonewgr*grauprat + rhonewice*icerat + rhonewgr*curat)))
1569 if (debug_print)
then
1571 print *,
' xlat, xlon', xlat, xlon
1572 print *,
'snow_mosaic = ',snow_mosaic
1573 print *,
'new snow,newsnowratio,rhosnfall =',newsn,newsnowratio,rhosnfall
1574 print *,
'snowrat,grauprat,icerat,curat,rhonewgr,rhonewsn,rhonewice',snowrat,grauprat,icerat,curat,rhonewgr,rhonewsn,rhonewice
1583 xsn=(rhosn*snwe+rhonewsn*newsn)/ &
1585 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
1588 IF(prcpms > zero)
THEN
1601 if(vegfrac > 0.01_kind_phys)
then
1604 interw=0.25_kind_phys*delt*prcpms*(one-exp(-0.5_kind_phys*lai))*vegfrac
1605 intersn=0.25_kind_phys*newsn*(one-exp(-0.5_kind_phys*lai))*vegfrac
1606 infwater=prcpms - interw/delt
1607 if((interw+intersn) > zero)
then
1608 intwratio=interw/(interw+intersn)
1612 dd1=cst + interw + intersn
1627 IF(newsn.GT.zero)
THEN
1629 snwe=max(zero,snwe+newsn-intersn)
1631 if(drip > zero)
then
1632 if (snow_mosaic==one)
then
1633 dripliq=drip*intwratio
1634 dripsn = drip - dripliq
1636 infwater=infwater+dripliq
1643 snhei=snwe*rhowater/rhosn
1644 newsn=newsn*rhowater/rhonewsn
1647 IF(snhei.GT.zero)
THEN
1664 if(isncovr_opt == 1)
then
1665 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
1666 elseif(isncovr_opt == 2)
then
1667 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
1668 if(ivgtyp == glacier .or. ivgtyp == bare)
then
1670 snowfrac2 = tanh( snhei/(2.5_kind_phys * 0.2_kind_phys *(rhosn/rhonewsn)**m))
1675 snowfrac2 = tanh( snhei/(2.5_kind_phys *min(0.2_kind_phys,znt) *(rhosn/rhonewsn)**m))
1678 snowfrac = 0.5_kind_phys*(snowfrac+snowfrac2)
1685 snowfrac = tanh( snhei/(10._kind_phys * facsnf *(rhosn/rhonewsn)**m))
1688 if(newsn > zero )
then
1689 snowfracnewsn=min(one,snowfallac*1.e-3_kind_phys/snhei_crit_newsn)
1694 if(hgt > 2500._kind_phys .and. ivgtyp == glacier) snowfrac=min(0.85_kind_phys,snowfrac)
1697 if(ivgtyp == urban) snowfrac=min(0.75_kind_phys,snowfrac)
1699 if(snowfrac < 0.75_kind_phys) snow_mosaic = one
1701 keep_snow_albedo = zero
1702 IF (snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys)
THEN
1704 keep_snow_albedo = one
1709 IF (debug_print )
THEN
1710 print *,
'SNHEI_CRIT,SNOWFRAC,SNHEI_CRIT_newsn,SNOWFRACnewsn', &
1711 snhei_crit,snowfrac,snhei_crit_newsn,snowfracnewsn
1716 IF(newsn.eq.zero .and. znt.le.0.2_kind_phys .and. ivgtyp.ne.isice)
then
1717 if( snhei .le. 2._kind_phys*znt)
then
1719 znt=0.55_kind_phys*znt+0.45_kind_phys*z0tbl(iland)
1720 elseif( snhei .gt. 2._kind_phys*znt .and. snhei .le. 4._kind_phys*znt)
then
1721 znt=0.2_kind_phys*znt+0.8_kind_phys*z0tbl(iland)
1722 elseif(snhei > 4._kind_phys*znt)
then
1729 IF(seaice .LT. 0.5_kind_phys)
THEN
1735 if( snow_mosaic == one)
then
1737 if(keep_snow_albedo > 0.9_kind_phys .and. albsn < 0.4_kind_phys)
then
1742 albsn = 0.7_kind_phys
1747 albsn = max(keep_snow_albedo*alb_snow, &
1748 min((alb_snow_free + &
1749 (alb_snow - alb_snow_free) * snowfrac), alb_snow))
1750 if(newsn > zero .and. keep_snow_albedo > 0.9_kind_phys .and. albsn < 0.4_kind_phys)
then
1755 albsn = 0.7_kind_phys
1759 emiss = max(keep_snow_albedo*emissn, &
1760 min((emiss_snowfree + &
1761 (emissn - emiss_snowfree) * snowfrac), emissn))
1764 IF (debug_print )
THEN
1766 print *,
'Snow on soil ALBsn,emiss,snow_mosaic',i,j,albsn,emiss,snow_mosaic
1777 if(albsn.lt.0.4_kind_phys .or. keep_snow_albedo==1)
then
1781 alb = min(albsn,max(albsn - 0.1_kind_phys*(soilt - 263.15_kind_phys)/ &
1782 (tfrz-263.15_kind_phys)*albsn, albsn - 0.05_kind_phys))
1786 if( snow_mosaic == one)
then
1790 albsn = max(keep_snow_albedo*alb_snow, &
1791 min((albice + (alb_snow - albice) * snowfrac), alb_snow))
1792 emiss = max(keep_snow_albedo*emissn, &
1794 min((emiss_snowfree + &
1795 (emissn - emiss_snowfree) * snowfrac), emissn))
1798 IF (debug_print )
THEN
1799 print *,
'Snow on ice snow_mosaic,ALBsn,emiss',i,j,albsn,emiss,snow_mosaic
1804 if(albsn.lt.alb_snow .or. keep_snow_albedo .eq.one)
then
1808 alb = min(albsn,max(albsn - 0.15_kind_phys*albsn*(soilt - 263.15_kind_phys)/ &
1809 (tfrz-263.15_kind_phys), albsn - 0.1_kind_phys))
1814 if (snow_mosaic==one)
then
1817 if(seaice .LT. 0.5_kind_phys)
then
1822 gswnew=gswin*(one-alb_snow_free)
1824 t3 = stbolt*soilt*soilt*soilt
1826 xinet = emiss_snowfree*(glw-upflux)
1827 rnet = gswnew + xinet
1828 IF ( add_fire_heat_flux .and. fire_heat_flux >0 )
then
1829 IF (debug_print )
THEN
1830 print *,
'RNET snow-free, fire_heat_flux, xlat/xlon',rnet, fire_heat_flux,xlat,xlon
1832 rnet = rnet + fire_heat_flux
1835 IF (debug_print )
THEN
1836 print *,
'Fractional snow - snowfrac=',snowfrac
1837 print *,
'Snowfrac<1 GSWin,GSWnew -',gswin,gswnew,
'SOILT, RNET',soilt,rnet
1840 soilm1ds(k) = soilm1d(k)
1842 smfrkeeps(k) = smfrkeep(k)
1843 keepfrs(k) = keepfr(k)
1844 soilices(k) = soilice(k)
1845 soiliqws(k) = soiliqw(k)
1859 CALL soil(debug_print,xlat, xlon, testptlat, testptlon,&
1861 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1862 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew,gswin, &
1863 emiss_snowfree,rnet,qkms,tkms,pc,csts,dripliq, &
1864 infwater,rho,vegfrac,lai,myj, &
1866 qwrtz,rhocs,dqm,qmin,ref,wilt, &
1867 psis,bclh,ksat,sat,cn, &
1868 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1870 lv,cp,rovcp,g0,cw,stbolt,tabs, &
1873 soilm1ds,ts1ds,smfrkeeps,keepfrs, &
1874 dews,soilts,qvgs,qsgs,qcgs,edir1s,ec1s, &
1875 ett1s,eetas,qfxs,hfxs,ss,evapls,prcpls,fltots,runoff1s, &
1876 runoff2s,mavails,soilices,soiliqws, &
1883 gswnew=gswin*(one-albice)
1885 t3 = stbolt*soilt*soilt*soilt
1887 xinet = emiss_snowfree*(glw-upflux)
1888 rnet = gswnew + xinet
1889 IF (debug_print )
THEN
1890 print *,
'Fractional snow - snowfrac=',snowfrac
1891 print *,
'Snowfrac<1 GSWin,GSWnew -',gswin,gswnew,
'SOILT, RNET',soilt,rnet
1904 CALL sice(debug_print,xlat,xlon, &
1906 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1907 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew, &
1908 0.98_kind_phys,rnet,qkms,tkms,rho,myj, &
1910 tice,rhosice,capice,thdifice, &
1911 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1913 lv,cp,rovcp,cw,stbolt,tabs, &
1915 ts1ds,dews,soilts,qvgs,qsgs,qcgs, &
1916 eetas,qfxs,hfxs,ss,evapls,prcpls,fltots &
1918 edir1 = eeta*1.e-3_kind_phys
1939 gswnew=gswin*(one-alb)
1941 t3 = stbolt*soilt*soilt*soilt
1943 xinet = emiss*(glw-upflux)
1944 rnet = gswnew + xinet
1945 IF (debug_print )
THEN
1947 print *,
'RNET=',rnet
1948 print *,
'SNOW - I,J,newsn,snwe,snhei,GSW,GSWnew,GLW,UPFLUX,ALB',&
1949 i,j,newsn,snwe,snhei,gsw,gswnew,glw,upflux,alb
1950 print *,
'GSWnew',gswnew,
'alb=',alb
1953 if (seaice .LT. 0.5_kind_phys)
then
1955 IF ( add_fire_heat_flux .and. fire_heat_flux>0 )
then
1956 IF (debug_print )
THEN
1957 print *,
'RNET snow, fire_heat_flux, xlat/xlon',rnet, fire_heat_flux,xlat,xlon
1959 rnet = rnet + fire_heat_flux
1961 if(snow_mosaic==one)
then
1966 CALL snowsoil (debug_print,xlat,xlon,testptlat,testptlon, &
1967 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1968 isncond_opt,isncovr_opt, &
1969 meltfactor,rhonewsn,snhei_crit, &
1970 iland,prcpms,rainf,newsn,snhei,snwe,snfr, &
1971 rhosn,patm,qvatm,qcatm, &
1972 glw,gswnew,gswin,emiss,rnet,ivgtyp, &
1973 qkms,tkms,pc,cst,dripsn,infwater, &
1974 rho,vegfrac,alb,znt,lai, &
1977 qwrtz,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1978 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1980 lv,cp,rovcp,g0,cw,stbolt,tabs, &
1983 ilnb,snweprint,snheiprint,rsm, &
1984 soilm1d,ts1d,smfrkeep,keepfr, &
1985 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1986 smelt,snoh,snflx,snom,edir1,ec1,ett1,eeta, &
1987 qfx,hfx,s,sublim,prcpl,fltot,runoff1,runoff2, &
1988 mavail,soilice,soiliqw,infiltr )
1991 if(snow_mosaic==one)
then
1997 CALL snowseaice (debug_print,xlat,xlon, &
1998 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
1999 isncond_opt,isncovr_opt, &
2000 meltfactor,rhonewsn,snhei_crit, &
2001 iland,prcpms,rainf,newsn,snhei,snwe,snfr, &
2002 rhosn,patm,qvatm,qcatm, &
2003 glw,gswnew,emiss,rnet, &
2004 qkms,tkms,rho,myj, &
2007 tice,rhosice,capice,thdifice, &
2008 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2010 lv,cp,rovcp,cw,stbolt,tabs, &
2012 ilnb,snweprint,snheiprint,rsm,ts1d, &
2013 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
2014 smelt,snoh,snflx,snom,eeta, &
2015 qfx,hfx,s,sublim,prcpl,fltot &
2017 edir1 = eeta*1.e-3_kind_phys
2035 if (snow_mosaic==one)
then
2038 if(seaice .LT. 0.5_kind_phys)
then
2040 IF (debug_print )
THEN
2042 print *,
' xlat, xlon', xlat, xlon
2043 print *,
' snowfrac = ',snowfrac
2044 print *,
' SOILT snow on land', ktau, i,j,soilt
2045 print *,
' SOILT on snow-free land', i,j,soilts
2046 print *,
' ts1d,ts1ds',i,j,ts1d,ts1ds
2047 print *,
' SNOW flux',i,j, snflx
2048 print *,
' Ground flux on snow-covered land',i,j, s
2049 print *,
' Ground flux on snow-free land', i,j,ss
2050 print *,
' CSTS, CST', i,j,csts,cst
2054 soilm1d(k) = soilm1ds(k)*(1.-snowfrac) + soilm1d(k)*snowfrac
2055 ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac
2056 smfrkeep(k) = smfrkeeps(k)*(1.-snowfrac) + smfrkeep(k)*snowfrac
2057 if(snowfrac > 0.5_kind_phys)
then
2058 keepfr(k) = keepfr(k)
2060 keepfr(k) = keepfrs(k)
2062 soilice(k) = soilices(k)*(1.-snowfrac) + soilice(k)*snowfrac
2063 soiliqw(k) = soiliqws(k)*(1.-snowfrac) + soiliqw(k)*snowfrac
2065 dew = dews*(one-snowfrac) + dew*snowfrac
2066 soilt = soilts*(one-snowfrac) + soilt*snowfrac
2067 qvg = qvgs*(one-snowfrac) + qvg*snowfrac
2068 qsg = qsgs*(one-snowfrac) + qsg*snowfrac
2069 qcg = qcgs*(one-snowfrac) + qcg*snowfrac
2070 edir1 = edir1s*(one-snowfrac) + edir1*snowfrac
2071 ec1 = ec1s*(one-snowfrac) + ec1*snowfrac
2072 cst = csts*(one-snowfrac) + cst*snowfrac
2073 ett1 = ett1s*(one-snowfrac) + ett1*snowfrac
2074 eeta = eetas*(one-snowfrac) + eeta*snowfrac
2075 qfx = qfxs*(one-snowfrac) + qfx*snowfrac
2076 hfx = hfxs*(one-snowfrac) + hfx*snowfrac
2077 s = ss*(one-snowfrac) + s*snowfrac
2078 evapl = evapls*(one-snowfrac)
2079 prcpl = prcpls*(one-snowfrac) + prcpl*snowfrac
2080 fltot = fltots*(one-snowfrac) + fltot*snowfrac
2081 alb = max(keep_snow_albedo*alb, &
2082 min((alb_snow_free + (alb - alb_snow_free) * snowfrac), alb))
2084 emiss = max(keep_snow_albedo*emissn, &
2085 min((emiss_snowfree + &
2086 (emissn - emiss_snowfree) * snowfrac), emissn))
2088 runoff1 = runoff1s*(one-snowfrac) + runoff1*snowfrac
2089 runoff2 = runoff2s*(one-snowfrac) + runoff2*snowfrac
2090 mavail = mavails*(one-snowfrac) + one*snowfrac
2091 infiltr = infiltrs*(one-snowfrac) + infiltr*snowfrac
2093 IF (debug_print )
THEN
2095 print *,
' Ground flux combined', xlat,xlon, s
2096 print *,
' SOILT combined on land', soilt
2097 print *,
' TS combined on land', ts1d
2102 IF (debug_print )
THEN
2103 print *,
'SOILT snow on ice', soilt
2106 ts1d(k) = ts1ds(k)*(one-snowfrac) + ts1d(k)*snowfrac
2108 dew = dews*(one-snowfrac) + dew*snowfrac
2109 soilt = soilts*(one-snowfrac) + soilt*snowfrac
2110 qvg = qvgs*(one-snowfrac) + qvg*snowfrac
2111 qsg = qsgs*(one-snowfrac) + qsg*snowfrac
2112 qcg = qcgs*(one-snowfrac) + qcg*snowfrac
2114 eeta = eetas*(one-snowfrac) + eeta*snowfrac
2115 qfx = qfxs*(one-snowfrac) + qfx*snowfrac
2116 hfx = hfxs*(one-snowfrac) + hfx*snowfrac
2117 s = ss*(one-snowfrac) + s*snowfrac
2118 prcpl = prcpls*(one-snowfrac) + prcpl*snowfrac
2119 fltot = fltots*(one-snowfrac) + fltot*snowfrac
2120 alb = max(keep_snow_albedo*alb, &
2121 min((albice + (alb - alb_snow_free) * snowfrac), alb))
2122 emiss = max(keep_snow_albedo*emissn, &
2123 min((emiss_snowfree + &
2124 (emissn - emiss_snowfree) * snowfrac), emissn))
2125 runoff1 = runoff1s*(one-snowfrac) + runoff1*snowfrac
2126 runoff2 = runoff2s*(one-snowfrac) + runoff2*snowfrac
2127 IF (debug_print )
THEN
2128 print *,
'SOILT combined on ice', soilt
2164 IF(snhei == zero)
then
2169 emiss = emiss_snowfree
2173 if(isncovr_opt == 1)
then
2174 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
2175 elseif(isncovr_opt == 2)
then
2177 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
2178 if(ivgtyp == glacier .or. ivgtyp == bare)
then
2180 snowfrac2 = tanh( snhei/(2.5_kind_phys * 0.2_kind_phys *(rhosn/rhonewsn)**m))
2185 snowfrac2 = tanh( snhei/(2.5 *min(0.2,znt) *(rhosn/rhonewsn)**m))
2188 snowfrac = 0.5_kind_phys*(snowfrac+snowfrac2)
2195 snowfrac = tanh( snhei/(10._kind_phys * facsnf *(rhosn/rhonewsn)**m))
2200 if(hgt > 2500._kind_phys .and. ivgtyp == glacier) snowfrac=min(0.85_kind_phys,snowfrac)
2202 if(ivgtyp == urban) snowfrac=min(0.75_kind_phys,snowfrac)
2206 IF (debug_print )
then
2208 print *,
'Snowfallac xlat, xlon',xlat,xlon
2209 print *,
'newsn [m],rhonewsn,newsnowratio=',newsn,rhonewsn,newsnowratio
2210 print *,
'Time-step newsn depth [m], swe [m]',newsn,newsn*rhonewsn
2211 print *,
'Time-step smelt: swe [m]' ,smelt*delt
2212 print *,
'Time-step sublim: swe,[kg m-2]',sublim*delt
2215 snowfallac = snowfallac + newsn * 1.e3_kind_phys
2217 IF (debug_print )
THEN
2219 print *,
'snowfallac,snhei,snwe',snowfallac,snhei,snwe
2230 t3 = stbolt*soilt*soilt*soilt
2232 xinet = emiss*(glw-upflux)
2233 rnet = gswnew + xinet
2234 IF (debug_print )
THEN
2235 print *,
'NO snow on the ground GSWnew -',gswnew,
'RNET=',rnet
2238 if(seaice .LT. 0.5_kind_phys)
then
2240 IF ( add_fire_heat_flux .and. fire_heat_flux>0)
then
2241 IF (debug_print )
THEN
2242 print *,
'RNET no snow, fire_heat_flux, xlat/xlon',rnet, fire_heat_flux,xlat,xlon
2244 rnet = rnet + fire_heat_flux
2247 CALL soil(debug_print,xlat, xlon, testptlat, testptlon,&
2249 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2250 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew,gswin, &
2251 emiss,rnet,qkms,tkms,pc,cst,drip,infwater, &
2252 rho,vegfrac,lai,myj, &
2254 qwrtz,rhocs,dqm,qmin,ref,wilt, &
2255 psis,bclh,ksat,sat,cn, &
2256 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2258 lv,cp,rovcp,g0,cw,stbolt,tabs, &
2261 soilm1d,ts1d,smfrkeep,keepfr, &
2262 dew,soilt,qvg,qsg,qcg,edir1,ec1, &
2263 ett1,eeta,qfx,hfx,s,evapl,prcpl,fltot,runoff1, &
2264 runoff2,mavail,soilice,soiliqw, &
2270 if(alb.ne.albice) gswnew=gsw/(one-alb)*(one-albice)
2272 rnet = gswnew + xinet
2274 CALL sice(debug_print,xlat,xlon, &
2276 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2277 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew, &
2278 emiss,rnet,qkms,tkms,rho,myj, &
2280 tice,rhosice,capice,thdifice, &
2281 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2283 lv,cp,rovcp,cw,stbolt,tabs, &
2285 ts1d,dew,soilt,qvg,qsg,qcg, &
2286 eeta,qfx,hfx,s,evapl,prcpl,fltot &
2288 edir1 = eeta*1.e-3_kind_phys
2308 END SUBROUTINE sfctmp
2316 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN ) :: t
2317 real (kind_phys),
INTENT(IN ) :: tn
2319 real (kind_phys) qsn, r,r1,r2
2322 r=(tn-173.15_kind_dbl_prec)/.05_kind_dbl_prec+one
2327 r = 1._kind_dbl_prec
2328 else if (i .GT. 5000)
then
2330 r = 5001._kind_dbl_prec
2334 qsn=(t(i+1)-r1)*r2 + r1
2342 SUBROUTINE soil (debug_print,xlat,xlon,testptlat,testptlon,&
2343 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,& !--- input variables
2344 PRCPMS,RAINF,PATM,QVATM,QCATM, &
2345 GLW,GSW,GSWin,EMISS,RNET, &
2346 QKMS,TKMS,PC,cst,drip,infwater,rho,vegfrac,lai, &
2348 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & !--- soil fixed fields
2349 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2350 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, & !--- constants
2352 soilmois,tso,smfrkeep,keepfr, & !--- output variables
2353 dew,soilt,qvg,qsg,qcg, &
2354 edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
2355 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
2356 soiliqw,infiltrp,smf)
2420 LOGICAL,
INTENT(IN ) :: debug_print
2421 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
2423 INTEGER,
INTENT(IN ) :: i,j,iland,isoil
2424 real (kind_phys),
INTENT(IN ) :: delt,conflx
2425 real (kind_phys),
INTENT(IN ) :: xlat,xlon,testptlat,testptlon
2426 LOGICAL,
INTENT(IN ) :: myj
2429 INTENT(IN ) :: patm, &
2434 INTENT(IN ) :: glw, &
2448 INTENT(IN ) :: rhocs, &
2458 real (kind_phys),
INTENT(IN ) :: cn, &
2467 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
2471 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
2473 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: tbq
2478 real (kind_phys),
DIMENSION( 1:nzs ) , &
2479 INTENT(INOUT) :: tso, &
2483 real (kind_phys),
DIMENSION( 1:nzs ) , &
2484 INTENT(INOUT) :: keepfr
2488 INTENT(INOUT) :: dew, &
2511 real (kind_phys),
DIMENSION(1:NZS),
INTENT(OUT) :: soilice, &
2516 real (kind_phys) :: infiltrp, transum , &
2518 tabs, t3, upflux, xinet
2519 real (kind_phys) :: cp,rovcp,g0,lv,stbolt,xlmelt,dzstop , &
2520 can,epot,fac,fltot,ft,fq,hft , &
2522 trans,zn,ci,cvw,tln,tavln,pi , &
2523 dd1,cmc2ms,drycan,wetcan , &
2525 real (kind_phys),
DIMENSION(1:NZS) :: transp,cap,diffu,hydro, &
2526 thdif,tranf,tav,soilmoism , &
2527 soilicem,soiliqwm,detal , &
2528 fwsat,lwsat,told,smold
2530 real (kind_phys) :: soiltold,smf
2531 real (kind_phys) :: soilres, alfa, fex, fex_fc, fc, psit
2533 INTEGER :: nzs1,nzs2,k
2573 dzstop=one/(zsmain(2)-zsmain(1))
2574 ras=rho*1.e-3_kind_phys
2575 riw=rhoice*1.e-3_kind_phys
2581 tln=log(tso(k)/tfrz)
2582 if(tln.lt.zero)
then
2583 soiliqw(k)=(dqm+qmin)*(xlmelt* &
2584 (tso(k)-tfrz)/tso(k)/grav/psis) &
2586 soiliqw(k)=max(zero,soiliqw(k))
2587 soiliqw(k)=min(soiliqw(k),soilmois(k))
2588 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2591 if(keepfr(k).eq.one)
then
2592 soilice(k)=min(soilice(k),smfrkeep(k))
2593 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
2598 soiliqw(k)=soilmois(k)
2605 tav(k)=0.5_kind_phys*(tso(k)+tso(k+1))
2606 soilmoism(k)=0.5_kind_phys*(soilmois(k)+soilmois(k+1))
2607 tavln=log(tav(k)/tfrz)
2609 if(tavln.lt.zero)
then
2610 soiliqwm(k)=(dqm+qmin)*(xlmelt* &
2611 (tav(k)-tfrz)/tav(k)/grav/psis) &
2613 fwsat(k)=dqm-soiliqwm(k)
2614 lwsat(k)=soiliqwm(k)+qmin
2615 soiliqwm(k)=max(zero,soiliqwm(k))
2616 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
2617 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
2619 if(keepfr(k).eq.one)
then
2620 soilicem(k)=min(soilicem(k), &
2621 0.5_kind_phys*(smfrkeep(k)+smfrkeep(k+1)))
2622 soiliqwm(k)=max(zero,soilmoism(k)-soilicem(k)*riw)
2623 fwsat(k)=dqm-soiliqwm(k)
2624 lwsat(k)=soiliqwm(k)+qmin
2629 soiliqwm(k)=soilmoism(k)
2637 if(soilice(k).gt.zero)
then
2638 smfrkeep(k)=soilice(k)
2640 smfrkeep(k)=soilmois(k)/riw
2648 CALL soilprop( debug_print, &
2649 xlat, xlon, testptlat, testptlon, &
2651 nzs,fwsat,lwsat,tav,keepfr, &
2652 soilmois,soiliqw,soilice, &
2653 soilmoism,soiliqwm,soilicem, &
2655 qwrtz,rhocs,dqm,qmin,psis,bclh,ksat, &
2657 riw,xlmelt,cp,g0_p,cvw,ci, &
2660 thdif,diffu,hydro,cap)
2667 q1=-qkms*ras*(qvatm - qsg)
2670 IF(qvatm.GE.qsg)
THEN
2678 wetcan=min(0.25_kind_phys,max(zero,(cst/sat))**cn)
2684 CALL transf(debug_print, &
2685 xlat, xlon, testptlat, testptlon, &
2687 nzs,nroot,soiliqw,tabs,lai,gswin, &
2689 dqm,qmin,ref,wilt,zshalf,pc,iland, &
2696 smold(k)=soilmois(k)
2703 fex=min(one,soilmois(1)/dqm)
2704 fex=max(fex,0.01_kind_phys)
2705 psit=psis*fex ** (-bclh)
2706 psit = max(-1.e5_kind_phys, psit)
2707 alfa=min(one,exp(g0_p*psit/r_v/soilt))
2726 if((soilmois(1)+qmin) > fc .or. (qvatm-qvg) > zero)
then
2729 fex_fc=min(one,(soilmois(1)+qmin)/fc)
2730 fex_fc=max(fex_fc,0.01_kind_phys)
2731 soilres=0.25_kind_phys*(one-cos(piconst*fex_fc))**2._kind_phys
2733 IF ( debug_print )
THEN
2734 print *,
'piconst=',piconst
2735 print *,
'fex,psit,psis,bclh,g0_p,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc', &
2736 fex,psit,psis,bclh,g0_p,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc
2743 CALL soiltemp(debug_print,xlat,xlon,testptlat,testptlon,&
2746 delt,ktau,conflx,nzs,nddzs,nroot, &
2748 patm,tabs,qvatm,qcatm,emiss,rnet, &
2749 qkms,tkms,pc,rho,vegfrac, lai, &
2750 thdif,cap,drycan,wetcan, &
2751 transum,dew,mavail,soilres,alfa, &
2753 dqm,qmin,bclh,zsmain,zshalf,dtdzs,tbq, &
2755 xlv,cp,g0_p,cvw,stbolt, &
2757 tso,soilt,qvg,qsg,qcg,x)
2765 IF(qvatm.GE.qsg)
THEN
2766 dew=qkms*(qvatm-qsg)
2774 transp(k)=vegfrac*ras*qkms* &
2776 tranf(k)*drycan/zshalf(nroot+1)
2777 IF(transp(k).GT.zero) transp(k)=zero
2788 tln=log(tso(k)/tfrz)
2789 if(tln.lt.zero)
then
2790 soiliqw(k)=(dqm+qmin)*(xlmelt* &
2791 (tso(k)-tfrz)/tso(k)/grav/psis) &
2793 soiliqw(k)=max(zero,soiliqw(k))
2794 soiliqw(k)=min(soiliqw(k),soilmois(k))
2795 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2797 if(keepfr(k).eq.one)
then
2798 soilice(k)=min(soilice(k),smfrkeep(k))
2799 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
2804 soiliqw(k)=soilmois(k)
2812 CALL soilmoist (debug_print, &
2813 xlat, xlon, testptlat, testptlon, &
2815 delt,nzs,nddzs,dtdzs,dtdzs2,riw, &
2816 zsmain,zshalf,diffu,hydro, &
2817 qsg,qvg,qcg,qcatm,qvatm,-infwater, &
2818 qkms,transp,drip,dew,zero,soilice,vegfrac, &
2821 dqm,qmin,ref,ksat,ras,infmax, &
2823 soilmois,soiliqw,mavail,runoff1, &
2836 if (soilice(k).gt.zero)
then
2837 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k))
then
2847 t3 = stbolt*soiltold*soiltold*soiltold
2848 upflux = t3 * 0.5_kind_phys*(soiltold+soilt)
2849 xinet = emiss*(glw-upflux)
2850 hft=-tkms*cp*rho*(tabs-soilt)
2851 hfx=-tkms*cp*rho*(tabs-soilt) &
2852 *(p1000mb*0.00001_kind_phys/patm)**rovcp
2853 q1=-qkms*ras*(qvatm - qsg)
2856 IF (q1.LE.zero)
THEN
2863 eeta=-qkms*ras*(qvatm/(one+qvatm) - qsg/(one+qsg))*rhowater
2864 cst= cst-eeta*delt*vegfrac
2865 IF (debug_print )
THEN
2867 print *,
'Cond MYJ EETA',eeta,eeta*xlv, i,j
2872 cst=cst+delt*dew*ras * vegfrac
2873 IF (debug_print )
THEN
2875 print *,
'Cond RUC LSM EETA',eeta,eeta*xlv, i,j
2882 edir1 =-soilres*(one-vegfrac)*qkms*ras* &
2885 ec1 = q1 * wetcan * vegfrac
2886 IF (debug_print )
THEN
2887 IF(i.eq.440.and.j.eq.180.or. qfx.gt.1000..or.i.eq.417.and.j.eq.540)
then
2888 print *,
'CST before update=',cst
2889 print *,
'EC1=',ec1,
'CMC2MS=',cmc2ms
2893 cst=max(zero,cst-ec1 * delt)
2897 eeta=-soilres*qkms*ras*(qvatm/(one+qvatm) - qvg/(one+qvg))*rhowater
2899 IF (debug_print )
THEN
2901 print *,
'QKMS,RAS,QVATM/(one+QVATM),QVG/(one+QVG),QSG ', &
2902 qkms,ras,qvatm/(one+qvatm),qvg/(one+qvg),qsg
2903 print *,
'Q1*(1.-vegfrac),EDIR1',q1*(one-vegfrac),edir1
2904 print *,
'CST,WETCAN,DRYCAN',cst,wetcan,drycan
2905 print *,
'EC1=',ec1,
'ETT1=',ett1,
'CMC2MS=',cmc2ms,
'CMC2MS*ras=',cmc2ms*ras
2908 eeta = (edir1 + ec1 + ett1)*rhowater
2909 IF (debug_print )
THEN
2911 print *,
'RUC LSM EETA',eeta,eeta*xlv
2915 eeta = (edir1 + ec1 + ett1)*rhowater
2917 IF (debug_print )
THEN
2918 print *,
'potential temp HFT ',hft
2919 print *,
'abs temp HFX ',hfx
2923 s=thdif(1)*cap(1)*dzstop*(tso(1)-tso(2))
2925 fltot=rnet-hft-xlv*eeta-s-x
2926 IF (debug_print )
THEN
2928 print *,
'SOIL - FLTOT,RNET,HFT,QFX,S,X=',i,j,fltot,rnet,hft,xlv*eeta,s,x
2929 print *,
'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac',&
2930 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac
2932 if(detal(1) .ne. zero)
then
2935 IF (debug_print )
THEN
2936 print *,
'detal(1),xlmelt,soiliqwm(1),delt',detal(1),xlmelt,soiliqwm(1),delt
2937 print *,
'Implicit phase change in the first layer - smf=',smf
2944 1123
FORMAT(i5,8f12.3)
2945 1133
FORMAT(i7,8e12.4)
2946 123
format(i6,f6.2,7f8.1)
2947 122
FORMAT(1x,2i3,6f8.1,f8.3,f8.2)
2957 SUBROUTINE sice ( debug_print,xlat,xlon, &
2958 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & !--- input variables
2959 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSW, &
2960 EMISS,RNET,QKMS,TKMS,rho,myj, &
2961 tice,rhosice,capice,thdifice, & !--- sea ice parameters
2962 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2963 xlv,CP,rovcp,cw,stbolt,tabs, & !--- constants
2964 tso,dew,soilt,qvg,qsg,qcg, & !--- output variables
2965 eeta,qfx,hfx,s,evapl,prcpl,fltot &
2978 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
2980 INTEGER,
INTENT(IN ) :: i,j,iland,isoil
2981 real (kind_phys),
INTENT(IN ) :: delt,conflx,xlat,xlon
2982 LOGICAL,
INTENT(IN ) :: myj, debug_print
2985 INTENT(IN ) :: patm, &
2990 INTENT(IN ) :: glw, &
2997 real (kind_phys),
DIMENSION(1:NZS) , &
3005 real (kind_phys),
INTENT(IN ) :: &
3010 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
3014 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
3016 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: tbq
3021 real (kind_phys),
DIMENSION( 1:nzs ),
INTENT(INOUT) :: tso
3024 INTENT(INOUT) :: dew, &
3038 real (kind_phys) :: x,x1,x2,x4,tn,denom
3039 real (kind_phys) :: rainf, prcpms , &
3040 tabs, t3, upflux, xinet
3042 real (kind_phys) :: cp,rovcp,g0,lv,stbolt,xlmelt,dzstop , &
3043 epot,fltot,ft,fq,hft,ras,cvw
3045 real (kind_phys) :: fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11, &
3046 pi,h,fkq,r210,aa,bb,pp,q1,qs1,ts1,tq2,tx2 , &
3049 real (kind_phys) :: aa1,rhcs, icemelt
3052 real (kind_phys),
DIMENSION(1:NZS) :: cotso,rhtso
3054 INTEGER :: nzs1,nzs2,k,k1,kn,kk
3066 dzstop=1./(zsmain(2)-zsmain(1))
3067 ras=rho*1.e-3_kind_phys
3080 x1=dtdzs(k1)*thdifice(kn-1)
3081 x2=dtdzs(k1+1)*thdifice(kn)
3082 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
3083 -x2*(tso(kn)-tso(kn+1))
3084 denom=1.+x1+x2-x2*cotso(k)
3086 rhtso(k+1)=(ft+x2*rhtso(k))/denom
3097 d9=thdifice(1)*rhcs*dzstop
3099 r211=.5_kind_phys*conflx/delt
3101 r22=.5_kind_phys/(thdifice(1)*delt*dzstop**2)
3102 r6=emiss *stbolt*.5_kind_phys*tn**4
3105 tdenom=d9*(one-d1+r22)+d10+r21+r7 &
3109 aa=xls*(fkq+r210)/tdenom
3110 bb=(d10*tabs+r21*tn+xls*(qvatm*fkq &
3111 +r210*qvg)+d11+d9*(d2+r22*tn) &
3112 +rainf*cvw*prcpms*max(tfrz,tabs) &
3117 IF (debug_print )
THEN
3118 print *,
' VILKA-SEAICE1'
3119 print *,
'D10,TABS,R21,TN,QVATM,FKQ', &
3120 d10,tabs,r21,tn,qvatm,fkq
3121 print *,
'RNET, EMISS, STBOLT, SOILT',rnet, emiss, stbolt, soilt
3122 print *,
'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
3123 r210,qvg,d11,d9,d2,r22,rainf,cvw,prcpms,tdenom
3124 print *,
'tn,aa1,bb,pp,fkq,r210', &
3125 tn,aa1,bb,pp,fkq,r210
3128 CALL vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
3132 tso(1)=min(con_tice,ts1)
3140 tso(k)=min(con_tice,rhtso(kk)+cotso(kk)*tso(k-1))
3146 t3 = stbolt*tn*tn*tn
3147 upflux = t3 *0.5_kind_phys*(tn+soilt)
3148 xinet = emiss*(glw-upflux)
3149 hft=-tkms*cp*rho*(tabs-soilt)
3150 hfx=-tkms*cp*rho*(tabs-soilt) &
3151 *(p1000mb*0.00001_kind_phys/patm)**rovcp
3152 q1=-qkms*ras*(qvatm - qsg)
3153 IF (q1.LE.zero)
THEN
3157 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*rhowater
3158 IF (debug_print )
THEN
3159 print *,
'MYJ EETA',eeta
3163 dew=qkms*(qvatm-qsg)
3165 IF (debug_print )
THEN
3166 print *,
'RUC LSM EETA',eeta
3175 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qvg/(1.+qvg))*rhowater
3176 IF (debug_print )
THEN
3177 print *,
'MYJ EETA',eeta
3183 IF (debug_print )
THEN
3184 print *,
'RUC LSM EETA',eeta
3192 s=thdifice(1)*capice(1)*dzstop*(tso(1)-tso(2))
3196 x= (cp*rho*r211+rhcs*zsmain(2)*0.5_kind_phys/delt)*(soilt-tn) + &
3197 xls*rho*r211*(qsg-qgold)
3200 -rainf*cvw*prcpms*(max(tfrz,tabs)-soilt)
3203 icemelt=rnet-xls*eeta -hft -s -x
3204 IF (debug_print )
THEN
3205 print *,
'icemelt=',icemelt
3208 fltot=rnet-xls*eeta-hft-s-x-icemelt
3209 IF (debug_print )
THEN
3210 print *,
'SICE - FLTOT,RNET,HFT,QFX,S,icemelt,X=', &
3211 fltot,rnet,hft,xls*eeta,s,icemelt,x
3223 SUBROUTINE snowsoil ( debug_print,xlat,xlon, &
3224 testptlat,testptlon, &
3225 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & !--- input variables
3226 isncond_opt,isncovr_opt, &
3227 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
3228 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
3231 GLW,GSW,GSWin,EMISS,RNET,IVGTYP, &
3232 QKMS,TKMS,PC,cst,drip,infwater, &
3233 rho,vegfrac,alb,znt,lai, &
3234 MYJ, & !--- soil fixed fields
3235 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
3236 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
3237 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, & !--- constants
3239 ilnb,snweprint,snheiprint,rsm, & !--- output variables
3240 soilmois,tso,smfrkeep,keepfr, &
3241 dew,soilt,soilt1,tsnav, &
3242 qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, &
3243 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
3244 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
3318 LOGICAL,
INTENT(IN ) :: debug_print
3319 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
3321 INTEGER,
INTENT(IN ) :: i,j,isoil,isncond_opt,isncovr_opt
3323 real (kind_phys),
INTENT(IN ) :: delt,conflx,prcpms, &
3324 rainf,newsnow,rhonewsn, &
3325 testptlat,testptlon, &
3326 snhei_crit,meltfactor,xlat,xlon
3328 LOGICAL,
INTENT(IN ) :: myj
3332 INTENT(IN ) :: patm, &
3336 real (kind_phys) , &
3337 INTENT(IN ) :: glw, &
3348 INTEGER,
INTENT(IN ) :: ivgtyp
3350 real (kind_phys) , &
3351 INTENT(IN ) :: rhocs, &
3362 real (kind_phys),
INTENT(IN ) :: cn, &
3371 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
3375 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
3377 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: tbq
3382 real (kind_phys),
DIMENSION( 1:nzs ) , &
3383 INTENT(INOUT) :: tso, &
3387 real (kind_phys),
DIMENSION( 1:nzs ) , &
3388 INTENT(INOUT) :: keepfr
3391 INTEGER,
INTENT(INOUT) :: iland
3395 real (kind_phys) , &
3396 INTENT(INOUT) :: dew, &
3429 INTEGER,
INTENT(INOUT) :: ilnb
3432 real (kind_phys),
DIMENSION(1:NZS),
INTENT(OUT) :: soilice, &
3435 real (kind_phys),
INTENT(OUT) :: rsm, &
3441 INTEGER :: nzs1,nzs2,k
3443 real (kind_phys) :: infiltrp, transum , &
3445 tabs, t3, upflux, xinet , &
3446 beta, snwepr,epdt,pp
3447 real (kind_phys) :: cp,rovcp,g0,lv,xlvm,stbolt,xlmelt,dzstop, &
3448 can,epot,fac,fltot,ft,fq,hft , &
3450 trans,zn,ci,cvw,tln,tavln,pi , &
3451 dd1,cmc2ms,drycan,wetcan , &
3452 infmax,riw,deltsn,h,umveg
3454 real (kind_phys),
DIMENSION(1:NZS) :: transp,cap,diffu,hydro, &
3455 thdif,tranf,tav,soilmoism , &
3456 soilicem,soiliqwm,detal , &
3457 fwsat,lwsat,told,smold
3458 real (kind_phys) :: soiltold, qgold
3460 real (kind_phys) :: rnet, x
3485 deltsn=0.05_kind_phys*rhowater/rhosn
3486 snth=0.01_kind_phys*rhowater/rhosn
3490 IF(snhei.GE.deltsn+snth)
THEN
3491 if(snhei-deltsn-snth.lt.snth) deltsn=0.5_kind_phys*(snhei-snth)
3492 IF (debug_print )
THEN
3493 print *,
'DELTSN is changed,deltsn,snhei,snth',i,j,deltsn,snhei,snth
3498 ras=rho*1.e-3_kind_dbl_prec
3499 riw=rhoice*1.e-3_kind_dbl_prec
3532 dzstop=one/(zsmain(2)-zsmain(1))
3540 tln=log(tso(k)/tfrz)
3541 if(tln.lt.zero)
then
3542 soiliqw(k)=(dqm+qmin)*(xlmelt* &
3543 (tso(k)-tfrz)/tso(k)/grav/psis) &
3545 soiliqw(k)=max(zero,soiliqw(k))
3546 soiliqw(k)=min(soiliqw(k),soilmois(k))
3547 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3550 if(keepfr(k).eq.1.)
then
3551 soilice(k)=min(soilice(k),smfrkeep(k))
3552 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
3557 soiliqw(k)=soilmois(k)
3564 tav(k)=0.5_kind_phys*(tso(k)+tso(k+1))
3565 soilmoism(k)=0.5_kind_phys*(soilmois(k)+soilmois(k+1))
3566 tavln=log(tav(k)/tfrz)
3568 if(tavln.lt.zero)
then
3569 soiliqwm(k)=(dqm+qmin)*(xlmelt* &
3570 (tav(k)-tfrz)/tav(k)/grav/psis) &
3572 fwsat(k)=dqm-soiliqwm(k)
3573 lwsat(k)=soiliqwm(k)+qmin
3574 soiliqwm(k)=max(zero,soiliqwm(k))
3575 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
3576 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
3578 if(keepfr(k).eq.one)
then
3579 soilicem(k)=min(soilicem(k), &
3580 0.5_kind_phys*(smfrkeep(k)+smfrkeep(k+1)))
3581 soiliqwm(k)=max(zero,soilmoism(k)-soilicem(k)*riw)
3582 fwsat(k)=dqm-soiliqwm(k)
3583 lwsat(k)=soiliqwm(k)+qmin
3588 soiliqwm(k)=soilmoism(k)
3596 if(soilice(k).gt.zero)
then
3597 smfrkeep(k)=soilice(k)
3599 smfrkeep(k)=soilmois(k)/riw
3607 CALL soilprop(debug_print, &
3608 xlat, xlon, testptlat, testptlon, &
3610 nzs,fwsat,lwsat,tav,keepfr, &
3611 soilmois,soiliqw,soilice, &
3612 soilmoism,soiliqwm,soilicem, &
3614 qwrtz,rhocs,dqm,qmin,psis,bclh,ksat, &
3616 riw,xlmelt,cp,g0_p,cvw,ci, &
3619 thdif,diffu,hydro,cap)
3635 epot = -fq*(qvatm-qsg)
3637 IF (debug_print )
THEN
3638 print *,
'SNWE after subtracting intercepted snow - snwe=',snwe,vegfrac,cst
3646 epdt = epot * ras *delt
3647 IF(epdt > zero .and. snwepr.LE.epdt)
THEN
3652 wetcan=min(0.25_kind_phys,max(zero,(cst/sat))**cn)
3658 CALL transf(debug_print, &
3659 xlat, xlon, testptlat, testptlon, &
3661 nzs,nroot,soiliqw,tabs,lai,gswin, &
3663 dqm,qmin,ref,wilt,zshalf,pc,iland, &
3670 smold(k)=soilmois(k)
3677 IF (debug_print )
THEN
3678print *,
'TSO before calling SNOWTEMP: ', tso
3680 CALL snowtemp(debug_print,xlat,xlon,testptlat,testptlon,&
3683 delt,ktau,conflx,nzs,nddzs,nroot, &
3684 isncond_opt,isncovr_opt, &
3685 snwe,snwepr,snhei,newsnow,snowfrac,snhei_crit, &
3686 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, &
3688 patm,tabs,qvatm,qcatm, &
3689 glw,gsw,emiss,rnet, &
3690 qkms,tkms,pc,rho,vegfrac, &
3691 thdif,cap,drycan,wetcan,cst, &
3692 tranf,transum,dew,mavail, &
3694 dqm,qmin,psis,bclh, &
3695 zsmain,zshalf,dtdzs,tbq, &
3697 xlvm,cp,rovcp,g0_p,cvw,stbolt, &
3699 snweprint,snheiprint,rsm, &
3700 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3701 smelt,snoh,snflx,s,ilnb,x)
3708 epot = -fq*(qvatm-qsg)
3709 IF(epot.GT.zero)
THEN
3712 transp(k)=vegfrac*ras*fq*(qvatm-qsg) &
3713 *tranf(k)*drycan/zshalf(nroot+1)
3731 tln=log(tso(k)/tfrz)
3732 if(tln.lt.zero)
then
3733 soiliqw(k)=(dqm+qmin)*(xlmelt* &
3734 (tso(k)-tfrz)/tso(k)/grav/psis) &
3736 soiliqw(k)=max(zero,soiliqw(k))
3737 soiliqw(k)=min(soiliqw(k),soilmois(k))
3738 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3740 if(keepfr(k).eq.one)
then
3741 soilice(k)=min(soilice(k),smfrkeep(k))
3742 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
3747 soiliqw(k)=soilmois(k)
3755 CALL soilmoist (debug_print,xlat,xlon,testptlat,testptlon,&
3757 delt,nzs,nddzs,dtdzs,dtdzs2,riw, &
3758 zsmain,zshalf,diffu,hydro, &
3759 qsg,qvg,qcg,qcatm,qvatm,-infwater, &
3761 zero,smelt,soilice,vegfrac, &
3764 dqm,qmin,ref,ksat,ras,infmax, &
3766 soilmois,soiliqw,mavail,runoff1, &
3772 IF(snhei.EQ.zero)
then
3778 snom=snom+smelt*delt*rhowater
3790 if (soilice(k).gt.zero)
then
3791 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k))
then
3800 t3 = stbolt*soiltold*soiltold*soiltold
3801 upflux = t3 *0.5_kind_phys*(soiltold+soilt)
3802 xinet = emiss*(glw-upflux)
3803 hfx=-tkms*cp*rho*(tabs-soilt) &
3804 *(p1000mb*0.00001_kind_phys/patm)**rovcp
3805 IF (debug_print )
THEN
3806 print *,
'potential temp HFX',hfx
3808 hft=-tkms*cp*rho*(tabs-soilt)
3809 IF (debug_print )
THEN
3810 print *,
'abs temp HFX',hft
3812 q1 = - fq*ras* (qvatm - qsg)
3814 IF (q1.LT.zero)
THEN
3822 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*rhowater
3823 cst= cst-eeta*delt*vegfrac
3824 IF (debug_print )
THEN
3825 print *,
'MYJ EETA cond', eeta
3829 dew=qkms*(qvatm-qsg)
3831 cst=cst+delt*dew*ras * vegfrac
3832 IF (debug_print )
THEN
3833 print *,
'RUC LSM EETA cond',eeta
3840 edir1 = q1*umveg *beta
3842 ec1 = q1 * wetcan * vegfrac
3844 cst=max(zero,cst-ec1 * delt)
3846 IF (debug_print )
THEN
3847 print*,
'Q1,umveg,beta',q1,umveg,beta
3848 print *,
'wetcan,vegfrac',wetcan,vegfrac
3849 print *,
'EC1,CMC2MS',ec1,cmc2ms
3854 eeta=-(qkms*ras*(qvatm/(one+qvatm) - qsg/(one+qsg))*rhowater)*beta
3855 IF (debug_print )
THEN
3856 print *,
'MYJ EETA', eeta*xlvm,eeta
3861 eeta = (edir1 + ec1 + ett1)*rhowater
3862 IF (debug_print )
THEN
3863 print *,
'RUC LSM EETA',eeta*xlvm,eeta
3867 eeta = (edir1 + ec1 + ett1)*rhowater
3872 fltot=rnet-hft-xlvm*eeta-s-snoh-x
3873 IF (debug_print )
THEN
3874 print *,
'SNOWSOIL - FLTOT,RNET,HFT,QFX,S,SNOH,X=',fltot,rnet,hft,xlvm*eeta,s,snoh,x
3875 print *,
'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta',&
3876 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta
3881 1123
FORMAT(i5,8f12.3)
3882 1133
FORMAT(i7,8e12.4)
3883 123
format(i6,f6.2,7f8.1)
3884 122
FORMAT(1x,2i3,6f8.1,f8.3,f8.2)
3887 END SUBROUTINE snowsoil
3895 SUBROUTINE snowseaice( debug_print,xlat,xlon, &
3896 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
3897 isncond_opt,isncovr_opt, &
3898 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
3899 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,snowfrac, &
3900 RHOSN,PATM,QVATM,QCATM, &
3901 GLW,GSW,EMISS,RNET, &
3902 QKMS,TKMS,RHO,myj, &
3903 ALB,ZNT, & !--- sea ice parameters
3904 tice,rhosice,capice,thdifice, &
3905 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
3906 xlv,CP,rovcp,cw,stbolt,tabs, & !--- constants
3907 ilnb,snweprint,snheiprint,rsm,tso, & !--- output variables
3908 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3909 SMELT,SNOH,SNFLX,SNOM,eeta, &
3910 qfx,hfx,s,sublim,prcpl,fltot &
3922 LOGICAL,
INTENT(IN ) :: debug_print
3923 INTEGER,
INTENT(IN ) :: ktau,nzs , &
3925 INTEGER,
INTENT(IN ) :: i,j,isoil,isncond_opt,isncovr_opt
3927 real (kind_phys),
INTENT(IN ) :: delt,conflx,prcpms, &
3928 rainf,newsnow,rhonewsn, &
3929 meltfactor,snhei_crit,xlat,xlon
3932 LOGICAL,
INTENT(IN ) :: myj
3935 INTENT(IN ) :: patm, &
3939 real (kind_phys) , &
3940 INTENT(IN ) :: glw, &
3947 real (kind_phys),
DIMENSION(1:NZS) , &
3954 real (kind_phys),
INTENT(IN ) :: &
3958 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
3962 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
3964 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: tbq
3968 real (kind_phys),
DIMENSION( 1:nzs ) , &
3969 INTENT(INOUT) :: tso
3971 INTEGER,
INTENT(INOUT) :: iland
3975 real (kind_phys) , &
3976 INTENT(INOUT) :: dew, &
4001 INTEGER,
INTENT(INOUT) :: ilnb
4003 real (kind_phys),
INTENT(OUT) :: rsm, &
4009 INTEGER :: nzs1,nzs2,k,k1,kn,kk
4010 real (kind_phys) :: x,x1,x2,dzstop,ft,tn,denom
4012 real (kind_phys) :: snth, newsn , &
4013 tabs, t3, upflux, xinet , &
4014 beta, snwepr,epdt,pp
4015 real (kind_phys) :: cp,rovcp,g0,lv,xlvm,stbolt,xlmelt , &
4016 epot,fltot,fq,hft,q1,ras,ci,cvw , &
4019 real (kind_phys) :: rhocsn,thdifsn, &
4020 xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
4022 real (kind_phys) :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
4023 real (kind_phys) :: fso,fsn, &
4024 fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11, &
4025 fkq,r210,aa,bb,qs1,ts1,tq2,tx2, &
4026 tdenom,aa1,rhcs,h1,tsob, snprim, &
4027 snodif,soh,tnold,qgold,snohgnew
4028 real (kind_phys),
DIMENSION(1:NZS) :: cotso,rhtso
4030 real (kind_phys) :: rnet,rsmfrac,soiltfrac,hsn,icemelt,rr
4033 real (kind_phys) :: keff, fact
4043 keff = 0.265_kind_phys
4054 deltsn=0.05_kind_phys*rhowater/rhosn
4055 snth=0.01_kind_phys*rhowater/rhosn
4059 IF(snhei.GE.deltsn+snth)
THEN
4060 if(snhei-deltsn-snth.lt.snth) deltsn=0.5_kind_phys*(snhei-snth)
4061 IF (debug_print )
THEN
4062 print *,
'DELTSN ICE is changed,deltsn,snhei,snth', &
4063 i,j, deltsn,snhei,snth
4068 ras=rho*1.e-3_kind_dbl_prec
4069 riw=rhoice*1.e-3_kind_dbl_prec
4073 rhocsn=sheatsn * rhosn
4075 rhonewcsn=sheatsn * rhonewsn
4077 if(isncond_opt == 1)
then
4079 thdifsn = 0.265_kind_phys/rhocsn
4084 if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys))
then
4085 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
4087 keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys
4090 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
4096 thdifsn = 4.431718e-7_kind_phys
4098 thdifsn = keff/rhocsn * fact
4102 ras=rho*1.e-3_kind_phys
4122 dzstop=one/(zsmain(2)-zsmain(1))
4137 snhei=snwe*rhowater/rhosn
4142 epot = -fq*(qvatm-qsg)
4143 epdt = epot * ras *delt
4144 IF(epdt.GT.zero .and. snwepr.LE.epdt)
THEN
4145 beta=snwepr/max(1.e-8_kind_phys,epdt)
4158 x1=dtdzs(k1)*thdifice(kn-1)
4159 x2=dtdzs(k1+1)*thdifice(kn)
4160 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
4161 -x2*(tso(kn)-tso(kn+1))
4162 denom=1.+x1+x2-x2*cotso(k)
4164 rhtso(k+1)=(ft+x2*rhtso(k))/denom
4168 IF(snhei.GE.snth)
then
4169 if(snhei.le.deltsn+snth)
then
4172 snprim=max(snth,snhei)
4175 xsn = delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*snprim)
4176 ddzsn = xsn / snprim
4177 x1sn = ddzsn * thdifsn
4178 x2 = dtdzs(1)*thdifice(1)
4179 ft = tso(1)+x1sn*(soilt-tso(1)) &
4181 denom = one + x1sn + x2 -x2*cotso(nzs1)
4182 cotso(nzs)=x1sn/denom
4183 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
4187 tsnav=0.5_kind_phys*(soilt+tso(1)) &
4195 xsn = delt/2._kind_phys/(0.5_kind_phys*snhei)
4196 xsn1= delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*(snhei-deltsn))
4197 ddzsn = xsn / deltsn
4198 ddzsn1 = xsn1 / (snhei-deltsn)
4199 x1sn = ddzsn * thdifsn
4200 x1sn1 = ddzsn1 * thdifsn
4201 x2 = dtdzs(1)*thdifice(1)
4202 ft = tso(1)+x1sn1*(soilt1-tso(1)) &
4204 denom = one + x1sn1 + x2 - x2*cotso(nzs1)
4205 cotso(nzs)=x1sn1/denom
4206 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
4207 ftsnow = soilt1+x1sn*(soilt-soilt1) &
4208 -x1sn1*(soilt1-tso(1))
4209 denomsn = 1. + x1sn + x1sn1 - x1sn1*cotso(nzs)
4211 rhtsn=(ftsnow+x1sn1*rhtso(nzs))/denomsn
4213 tsnav=0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
4214 +(soilt1+tso(1))*(snhei-deltsn)) &
4219 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
4222 snprim=snhei+zsmain(2)
4227 xsn = delt/2._kind_phys/((zshalf(3)-zsmain(2))+0.5_kind_phys*snprim)
4229 x1sn = ddzsn * (fsn*thdifsn+fso*thdifice(1))
4230 x2=dtdzs(2)*thdifice(2)
4231 ft=tso(2)+x1sn*(soilt-tso(2))- &
4233 denom = one + x1sn + x2 - x2*cotso(nzs-2)
4234 cotso(nzs1) = x1sn/denom
4235 rhtso(nzs1)=(ft+x2*rhtso(nzs-2))/denom
4236 tsnav=0.5_kind_phys*(soilt+tso(1)) &
4238 cotso(nzs)=cotso(nzs1)
4239 rhtso(nzs)=rhtso(nzs1)
4250 epot=-qkms*(qvatm-qsg)
4257 d9=thdifice(1)*rhcs*dzstop
4259 r211=.5_kind_phys*conflx/delt
4261 r22=.5_kind_phys/(thdifice(1)*delt*dzstop**2)
4262 r6=emiss *stbolt*.5_kind_phys*tn**4
4266 IF(snhei.GE.snth)
THEN
4267 if(snhei.le.deltsn+snth)
then
4276 d9sn= thdifsn*rhocsn / snprim
4277 r22sn = snprim*snprim*0.5_kind_phys/(thdifsn*delt)
4280 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
4284 d9sn = (fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)/ &
4286 r22sn = snprim*snprim*0.5_kind_phys &
4287 /((fsn*thdifsn+fso*thdifice(1))*delt)
4290 IF(snhei.eq.zero)
then
4300 tdenom = d9sn*(one-d1sn +r22sn)+d10+r21+r7 &
4302 +rhonewcsn*newsnow/delt
4306 aa=xlvm*(beta*fkq+r210)/tdenom
4307 bb=(d10*tabs+r21*tn+xlvm*(qvatm* &
4309 +r210*qvg)+d11+d9sn*(d2sn+r22sn*tn) &
4310 +rainf*cvw*prcpms*max(tfrz,tabs) &
4311 + rhonewcsn*newsnow/delt*min(tfrz,tabs) &
4314 pp=patm*1.e3_kind_phys
4319 IF (debug_print )
THEN
4320 print *,
'VILKA-SNOW on SEAICE'
4321 print *,
'tn,aa1,bb,pp,fkq,r210', &
4322 tn,aa1,bb,pp,fkq,r210
4323 print *,
'TABS,QVATM,TN,QVG=',tabs,qvatm,tn,qvg
4326 CALL vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
4334 if(nmelt==1 .and. snowfrac==one)
then
4335 soilt = min(tfrz,soilt)
4338 IF (debug_print )
THEN
4339 print *,
' AFTER VILKA-SNOW on SEAICE'
4340 print *,
' TS1,QS1: ', ts1,qs1
4343 IF(snhei.GE.snth)
THEN
4344 if(snhei.gt.deltsn+snth)
then
4346 soilt1=min(tfrz,rhtsn+cotsn*soilt)
4347 tso(1)=min(con_tice,(rhtso(nzs)+cotso(nzs)*soilt1))
4351 tso(1)=min(con_tice,(rhtso(nzs)+cotso(nzs)*soilt))
4355 ELSEIF (snhei > zero .and. snhei < snth)
THEN
4357 tso(2)=min(con_tice,(rhtso(nzs1)+cotso(nzs1)*soilt))
4358 tso(1)=min(con_tice,(tso(2)+(soilt-tso(2))*fso))
4363 tso(1)=min(con_tice,soilt)
4364 soilt1=min(con_tice,soilt)
4368 IF (snhei > zero .and. snhei < snth)
THEN
4372 tso(k)=min(con_tice,rhtso(kk)+cotso(kk)*tso(k-1))
4377 tso(k)=min(con_tice,rhtso(kk)+cotso(kk)*tso(k-1))
4389 if(nmelt.eq.1)
go to 220
4393 IF(soilt>tfrz .AND. beta==one .AND. snhei>zero)
THEN
4396 soiltfrac=snowfrac*tfrz+(1.-snowfrac)*min(con_tice,soilt)
4398 qsg= qsn(soiltfrac,tbq)/pp
4399 t3 = stbolt*tnold*tnold*tnold
4400 upflux = t3 * 0.5_kind_phys*(tnold+soiltfrac)
4401 xinet = emiss*(glw-upflux)
4402 epot = -qkms*(qvatm-qsg)
4405 IF (q1.LE.zero)
THEN
4413 eeta = q1 * beta * rhowater
4418 hfx=d10*(tabs-soiltfrac)
4420 IF(snhei.GE.snth)
then
4421 soh=thdifsn*rhocsn*(soiltfrac-tsob)/snprim
4424 soh=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4425 (soiltfrac-tsob)/snprim
4428 x= (r21+d9sn*r22sn)*(soiltfrac-tnold) + &
4429 xlvm*r210*(qsg-qgold)
4431 snoh=rnet+qfx +hfx &
4432 +rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soiltfrac) &
4433 -soh-x+rainf*cvw*prcpms* &
4434 (max(tfrz,tabs)-soiltfrac)
4436 IF (debug_print )
THEN
4437 print *,
'SNOWSEAICE melt I,J,SNOH,RNET,QFX,HFX,SOH,X',i,j,snoh,rnet,qfx,hfx,soh,x
4438 print *,
'RHOnewCSN*NEWSNOW/DELT*(min(tfrz,TABS)-soiltfrac)', &
4439 rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soiltfrac)
4440 print *,
'RAINF*CVW*PRCPMS*(max(tfrz,TABS)-soiltfrac)', &
4441 rainf*cvw*prcpms*(max(tfrz,tabs)-soiltfrac)
4443 snoh=amax1(zero,snoh)
4445 smelt= snoh /xlmelt*1.e-3_kind_phys
4446 smelt=amin1(smelt,snwepr/delt-beta*epot*ras)
4447 smelt=amax1(zero,smelt)
4449 IF (debug_print )
THEN
4450 print *,
'1-SMELT i,j',smelt,i,j
4453 smelt= amin1(smelt,delt/60._kind_phys* 5.6e-8_kind_phys*meltfactor*max(one,(soilt-tfrz)))
4454 IF (debug_print )
THEN
4455 print *,
'2-SMELT i,j',smelt,i,j
4459 rr=snwepr/delt-beta*epot*ras
4461 IF (debug_print )
THEN
4462 print *,
'3- SMELT i,j,smelt,rr',i,j,smelt,rr
4464 snohgnew=smelt*xlmelt*1.e3
4465 snodif=amax1(zero,(snoh-snohgnew))
4469 IF (debug_print )
THEN
4470 print*,
'soiltfrac,soilt,SNOHGNEW,SNODIF=', &
4471 i,j,soiltfrac,soilt,snohgnew,snodif
4472 print *,
'SNOH,SNODIF',snoh,snodif
4476 rsmfrac=min(0.18_kind_phys,(max(0.08_kind_phys,snwepr/0.10_kind_phys*0.13_kind_phys)))
4477 if(snhei > 0.01_kind_phys)
then
4478 rsm=rsmfrac*smelt*delt
4484 smelt=amax1(zero,smelt-rsm/delt)
4485 IF (debug_print )
THEN
4486 print *,
'4-SMELT i,j,smelt,rsm,snwepr,rsmfrac', &
4487 i,j,smelt,rsm,snwepr,rsmfrac
4492 snwe = amax1(zero,(snwepr- &
4493 (smelt+beta*epot*ras)*delt &
4499 if(snhei > zero.and. beta == one)
then
4500 epot=-qkms*(qvatm-qsg)
4501 snwe = amax1(zero,(snwepr- &
4502 beta*epot*ras*delt))
4514 if(smelt > zero .and. rsm > zero)
then
4515 if(snwe.le.rsm)
then
4516 IF (debug_print )
THEN
4517 print *,
'SEAICE SNWE<RSM snwe,rsm,smelt*delt,epot*ras*delt,beta', &
4518 snwe,rsm,smelt*delt,epot*ras*delt,beta
4526 xsn=(rhosn*(snwe-rsm)+rhowater*rsm)/ &
4528 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
4530 rhocsn=sheatsn* rhosn
4531 if(isncond_opt == 1)
then
4533 thdifsn = 0.265_kind_phys/rhocsn
4538 if(rhosn < 156._kind_phys .or. (newsn > zero .and. rhonewsn < 156._kind_phys))
then
4539 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
4541 keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys
4544 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
4550 thdifsn = 4.431718e-7_kind_phys
4552 thdifsn = keff/rhocsn * fact
4563 snheiprint=snweprint*rhowater / rhosn
4565 IF (debug_print )
THEN
4566print *,
'snweprint : ',snweprint
4567print *,
'D9SN,SOILT,TSOB : ', d9sn,soilt,tsob
4569 IF(snhei.GT.zero)
THEN
4571 tsnav=0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
4572 +(soilt1+tso(1))*(snhei-deltsn)) &
4575 tsnav=0.5_kind_phys*(soilt+tso(1)) - tfrz
4580 pp=patm*1.e3_kind_phys
4581 qsg= qsn(soilt,tbq)/pp
4582 epot = -fq*(qvatm-qsg)
4583 IF(epot.LT.zero)
THEN
4588 snom=snom+smelt*delt*rhowater
4592 t3 = stbolt*tnold*tnold*tnold
4593 upflux = t3 *0.5_kind_phys*(soilt+tnold)
4594 xinet = emiss*(glw-upflux)
4595 hft=-tkms*cp*rho*(tabs-soilt)
4596 hfx=-tkms*cp*rho*(tabs-soilt) &
4597 *(p1000mb*0.00001_kind_phys/patm)**rovcp
4598 q1 = - fq*ras* (qvatm - qsg)
4599 IF (q1.LT.zero)
THEN
4603 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*rhowater
4606 dew=qkms*(qvatm-qsg)
4616 eeta=-qkms*ras*beta*(qvatm/(1.+qvatm) - qvg/(1.+qvg))*rhowater
4620 eeta = q1*beta*rhowater
4623 eeta = q1*beta*rhowater
4628 IF(snhei.GE.snth)
then
4629 s=thdifsn*rhocsn*(soilt-tsob)/snprim
4631 ELSEIF(snhei.lt.snth.and.snhei.GT.zero)
then
4632 s=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4635 IF (debug_print )
THEN
4636 print *,
'SNOW is thin, snflx',i,j,snflx
4639 snflx=d9sn*(soilt-tsob)
4640 IF (debug_print )
THEN
4641 print *,
'SNOW is GONE, snflx',i,j,snflx
4645 snhei=snwe *rhowater / rhosn
4647 IF (debug_print )
THEN
4648 print *,
'SNHEI,SNOH',i,j,snhei,snoh
4651 x= (r21+d9sn*r22sn)*(soilt-tnold) + &
4652 xlvm*r210*(qsg-qgold)
4653 IF (debug_print )
THEN
4654 print *,
'SNOWSEAICE storage ',i,j,x
4655 print *,
'R21,D9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim', &
4656 r21,d9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim
4659 -rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soilt) &
4660 -rainf*cvw*prcpms*(max(tfrz,tabs)-soilt)
4663 icemelt = rnet-hft-xlvm*eeta-s-snoh-x
4664 IF (debug_print )
THEN
4665 print *,
'SNOWSEAICE icemelt=',icemelt
4668 fltot=rnet-hft-xlvm*eeta-s-snoh-x-icemelt
4669 IF (debug_print )
THEN
4670 print *,
'i,j,snhei,qsg,soilt,soilt1,tso,TABS,QVATM', &
4671 i,j,snhei,qsg,soilt,soilt1,tso,tabs,qvatm
4672 print *,
'SNOWSEAICE - FLTOT,RNET,HFT,QFX,S,SNOH,icemelt,snodif,X,SOILT=' &
4673 ,fltot,rnet,hft,xlvm*eeta,s,snoh,icemelt,snodif,x,soilt
4676 IF(snhei.EQ.zero)
then
4678 emiss=0.98_kind_phys
4685 END SUBROUTINE snowseaice
4691 SUBROUTINE soiltemp( debug_print,xlat,xlon,testptlat,testptlon,&
4692 i,j,iland,isoil, & !--- input variables
4693 delt,ktau,conflx,nzs,nddzs,nroot, &
4694 PRCPMS,RAINF,PATM,TABS,QVATM,QCATM, &
4696 QKMS,TKMS,PC,RHO,VEGFRAC,lai, &
4697 THDIF,CAP,DRYCAN,WETCAN, &
4698 TRANSUM,DEW,MAVAIL,soilres,alfa, &
4699 DQM,QMIN,BCLH, & !---soil fixed fields
4700 ZSMAIN,ZSHALF,DTDZS,TBQ, &
4701 XLV,CP,G0_P,CVW,STBOLT, & !--- constants
4702 TSO,SOILT,QVG,QSG,QCG,X)
4757 LOGICAL,
INTENT(IN ) :: debug_print
4758 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
4760 INTEGER,
INTENT(IN ) :: i,j,iland,isoil
4761 real (kind_phys),
INTENT(IN ) :: delt,conflx,prcpms, rainf
4762 real (kind_phys),
INTENT(IN ) :: xlat, xlon, testptlat, testptlon
4763 real (kind_phys),
INTENT(INOUT) :: drycan,wetcan,transum
4766 INTENT(IN ) :: patm, &
4770 real (kind_phys) , &
4783 real (kind_phys) , &
4788 real (kind_phys) , &
4793 real (kind_phys),
INTENT(IN ) :: cp, &
4801 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
4806 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
4808 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: tbq
4813 real (kind_phys),
DIMENSION( 1:nzs ) , &
4814 INTENT(INOUT) :: tso
4817 real (kind_phys) , &
4828 real (kind_phys) :: x,x1,x2,x4,dzstop,can,ft,sph , &
4829 tn,trans,umveg,denom,fex
4831 real (kind_phys) :: fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11, &
4832 pi,h,fkq,r210,aa,bb,pp,q1,qs1,ts1,tq2,tx2 , &
4835 real (kind_phys) :: c,cc,aa1,rhcs,h1, qgold
4837 real (kind_phys),
DIMENSION(1:NZS) :: cotso,rhtso
4839 INTEGER :: nzs1,nzs2,k,k1,kn,kk, iter
4848 dzstop=1./(zsmain(2)-zsmain(1))
4864 x1=dtdzs(k1)*thdif(kn-1)
4865 x2=dtdzs(k1+1)*thdif(kn)
4866 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
4867 -x2*(tso(kn)-tso(kn+1))
4868 denom=1.+x1+x2-x2*cotso(k)
4870 rhtso(k+1)=(ft+x2*rhtso(k))/denom
4880 trans=transum*drycan/zshalf(nroot+1)
4882 umveg=(1.-vegfrac) * soilres
4888 d9=thdif(1)*rhcs*dzstop
4890 r211=.5_kind_phys*conflx/delt
4892 r22=.5_kind_phys/(thdif(1)*delt*dzstop**2)
4893 r6=emiss *stbolt*.5_kind_phys*tn**4
4896 tdenom=d9*(one-d1+r22)+d10+r21+r7 &
4902 aa=xlv*(fkq*umveg+r210)/tdenom
4903 bb=(d10*tabs+r21*tn+xlv*(qvatm* &
4905 +r210*qvg)+d11+d9*(d2+r22*tn) &
4906 +rainf*cvw*prcpms*max(tfrz,tabs) &
4909 pp=patm*1.e3_kind_phys
4911 CALL vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
4915 IF (debug_print )
THEN
4916 print *,
'VILKA1 - TS1,QS1,TQ2,H,TX2,Q1',ts1,qs1,tq2,h,tx2,q1
4918 IF(q1.LT.qs1)
GOTO 100
4924 qcg=max(zero,q1-qs1)
4925 IF (debug_print )
THEN
4926 print *,
'90 QVG,QSG,QCG,TSO(1)',qvg,qsg,qcg,tso(1)
4933 CALL vilka(tn,aa,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
4935 IF (debug_print )
THEN
4937 print *,
'VILKA2 - TS1,QS1,TQ2,H,TX2,Q1',ts1,qs1,tq2,h,tx2,q1
4939 IF(q1.GE.qs1)
GOTO 90
4950 IF (debug_print )
THEN
4951 print *,
'200 QVG,QSG,QCG,TSO(1)',qvg,qsg,qcg,tso(1)
4953 IF (debug_print )
THEN
4955 print *,
'QVATM,QVG,QSG,QCG,TS1',qvatm,qvg,qsg,qcg,ts1
4965 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
4968 x= (cp*rho*r211+rhcs*zsmain(2)*0.5_kind_phys/delt)*(soilt-tn) + &
4969 xlv*rho*r211*(qvg-qgold)
4971 IF (debug_print )
THEN
4972 print*,
'SOILTEMP storage, i,j,x,soilt,tn,qvg,qvgold', &
4973 i,j,x,soilt,tn,qvg,qgold
4974 print *,
'TEMP term (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN)',&
4975 (cp*rho*r211+rhcs*zsmain(2)*0.5_kind_phys/delt)*(soilt-tn)
4976 print *,
'QV term XLV*rho*r211*(QVG-QGOLD)',xlv*rho*r211*(qvg-qgold)
4980 -rainf*cvw*prcpms*(max(tfrz,tabs)-soilt)
4982 IF (debug_print )
THEN
4987 END SUBROUTINE soiltemp
4993 SUBROUTINE snowtemp( debug_print,xlat,xlon, &
4994 testptlat,testptlon,i,j,iland,isoil, & !--- input variables
4995 delt,ktau,conflx,nzs,nddzs,nroot, &
4996 isncond_opt,isncovr_opt, &
4997 snwe,snwepr,snhei,newsnow,snowfrac,snhei_crit, &
4998 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
5000 PATM,TABS,QVATM,QCATM, &
5001 GLW,GSW,EMISS,RNET, &
5002 QKMS,TKMS,PC,RHO,VEGFRAC, &
5003 THDIF,CAP,DRYCAN,WETCAN,CST, &
5004 TRANF,TRANSUM,DEW,MAVAIL, &
5005 DQM,QMIN,PSIS,BCLH, & !--- soil fixed fields
5006 ZSMAIN,ZSHALF,DTDZS,TBQ, &
5007 XLVM,CP,rovcp,G0_P,CVW,STBOLT, & !--- constants
5008 SNWEPRINT,SNHEIPRINT,RSM, & !--- output variables
5009 TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, &
5010 SMELT,SNOH,SNFLX,S,ILNB,X)
5064 LOGICAL,
INTENT(IN ) :: debug_print
5065 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
5068 INTEGER,
INTENT(IN ) :: i,j,iland,isoil,isncond_opt,isncovr_opt
5069 real (kind_phys),
INTENT(IN ) :: delt,conflx,prcpms , &
5070 rainf,newsnow,deltsn,snth , &
5071 tabs,transum,snwepr , &
5072 testptlat,testptlon , &
5073 rhonewsn,meltfactor,xlat,xlon,snhei_crit
5078 INTENT(IN ) :: patm, &
5082 real (kind_phys) , &
5083 INTENT(IN ) :: glw, &
5092 real (kind_phys) , &
5099 real (kind_phys),
INTENT(IN ) :: cp, &
5107 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
5113 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
5115 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: tbq
5120 real (kind_phys),
DIMENSION( 1:nzs ) , &
5121 INTENT(INOUT) :: tso
5125 real (kind_phys) , &
5126 INTENT(INOUT) :: dew, &
5145 real (kind_phys),
INTENT(INOUT) :: drycan, wetcan
5147 real (kind_phys),
INTENT(OUT) :: rsm, &
5150 INTEGER,
INTENT(OUT) :: ilnb
5154 INTEGER :: nzs1,nzs2,k,k1,kn,kk
5156 real (kind_phys) :: x,x1,x2,x4,dzstop,can,ft,sph, &
5157 tn,trans,umveg,denom
5159 real (kind_phys) :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
5161 real (kind_phys) :: t3,upflux,xinet,ras, &
5162 xlmelt,rhocsn,thdifsn, &
5163 beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
5165 real (kind_phys) :: fso,fsn, &
5166 fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11, &
5167 pi,h,fkq,r210,aa,bb,pp,q1,qs1,ts1,tq2,tx2, &
5168 tdenom,c,cc,aa1,rhcs,h1, &
5169 tsob, snprim, sh1, sh2, &
5170 smeltg,snohg,snodif,soh, &
5171 cmc2ms,tnold,qgold,snohgnew
5173 real (kind_phys),
DIMENSION(1:NZS) :: transp,cotso,rhtso
5174 real (kind_phys) :: edir1, &
5181 real (kind_phys) :: rnet,rsmfrac,soiltfrac,hsn,rr,keff,fact
5182 integer :: nmelt, iter
5191 keff = 0.265_kind_phys
5199 IF (debug_print )
THEN
5200print *,
'SNOWTEMP: SNHEI,SNTH,SOILT1: ',snhei,snth,soilt1,soilt
5203 rhocsn=sheatsn* rhosn
5204 rhonewcsn=sheatsn* rhonewsn
5205 if(isncond_opt == 1)
then
5207 thdifsn = 0.265_kind_phys/rhocsn
5212 if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys))
then
5213 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
5215 keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys
5216 if(debug_print)
then
5217 print *,
'SnowTemp xlat,xlon,rhosn,keff', xlat,xlon,rhosn,keff,keff/rhocsn*fact
5218 print *,
'SNOWTEMP - 0.265/rhocsn',0.265_kind_phys/rhocsn
5221 if ( debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5222 print *,
'SNOWTEMP - xlat,xlon,newsnow,rhonewsn,rhosn,fact,keff',xlat,xlon,newsnow, rhonewsn,rhosn,fact,keff
5225 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
5231 thdifsn = 4.431718e-7_kind_phys
5233 thdifsn = keff/rhocsn * fact
5235 if (debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5236 print *,
'SNOWTEMP - thdifsn',xlat,xlon,thdifsn
5237 print *,
'SNOWTEMP - 0.265/rhocsn',0.265_kind_phys/rhocsn
5242 ras=rho*1.e-3_kind_phys
5260 dzstop=one/(zsmain(2)-zsmain(1))
5270 x1=dtdzs(k1)*thdif(kn-1)
5271 x2=dtdzs(k1+1)*thdif(kn)
5272 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
5273 -x2*(tso(kn)-tso(kn+1))
5274 denom=1.+x1+x2-x2*cotso(k)
5276 rhtso(k+1)=(ft+x2*rhtso(k))/denom
5280 IF(snhei.GE.snth)
then
5281 if(snhei.le.deltsn+snth)
then
5283 IF (debug_print )
THEN
5284 print *,
'1-layer - snth,snhei,deltsn',snth,snhei,deltsn
5287 snprim=max(snth,snhei)
5290 xsn = delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*snprim)
5291 ddzsn = xsn / snprim
5292 x1sn = ddzsn * thdifsn
5293 x2 = dtdzs(1)*thdif(1)
5294 ft = tso(1)+x1sn*(soilt-tso(1)) &
5296 denom = one + x1sn + x2 -x2*cotso(nzs1)
5297 cotso(nzs)=x1sn/denom
5298 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
5302 tsnav=min(zero,0.5_kind_phys*(soilt+tso(1))-tfrz)
5306 IF (debug_print )
THEN
5307 print *,
'2-layer - snth,snhei,deltsn',snth,snhei,deltsn
5312 xsn = delt/2._kind_phys/(0.5_kind_phys*deltsn)
5313 xsn1= delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*(snhei-deltsn))
5314 ddzsn = xsn / deltsn
5315 ddzsn1 = xsn1 / (snhei-deltsn)
5316 x1sn = ddzsn * thdifsn
5317 x1sn1 = ddzsn1 * thdifsn
5318 x2 = dtdzs(1)*thdif(1)
5319 ft = tso(1)+x1sn1*(soilt1-tso(1)) &
5321 denom = 1. + x1sn1 + x2 - x2*cotso(nzs1)
5322 cotso(nzs)=x1sn1/denom
5323 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
5324 ftsnow = soilt1+x1sn*(soilt-soilt1) &
5325 -x1sn1*(soilt1-tso(1))
5326 denomsn = one + x1sn + x1sn1 - x1sn1*cotso(nzs)
5328 rhtsn=(ftsnow+x1sn1*rhtso(nzs))/denomsn
5330 tsnav=min(zero,0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
5331 +(soilt1+tso(1))*(snhei-deltsn)) &
5335 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
5338 snprim=snhei+zsmain(2)
5343 xsn = delt/2._kind_phys/((zshalf(3)-zsmain(2))+0.5_kind_phys*snprim)
5345 x1sn = ddzsn * (fsn*thdifsn+fso*thdif(1))
5346 x2=dtdzs(2)*thdif(2)
5347 ft=tso(2)+x1sn*(soilt-tso(2))- &
5349 denom = one + x1sn + x2 - x2*cotso(nzs-2)
5350 cotso(nzs1) = x1sn/denom
5351 rhtso(nzs1)=(ft+x2*rhtso(nzs-2))/denom
5352 tsnav=min(zero,0.5_kind_phys*(soilt+tso(1)) &
5354 cotso(nzs)=cotso(nzs1)
5355 rhtso(nzs)=rhtso(nzs1)
5368 epot=-qkms*(qvatm-qgold)
5371 trans=transum*drycan/zshalf(nroot+1)
5378 d9=thdif(1)*rhcs*dzstop
5380 r211=.5_kind_phys*conflx/delt
5382 r22=.5_kind_phys/(thdif(1)*delt*dzstop**2)
5383 r6=emiss *stbolt*.5_kind_phys*tn**4
5387 IF(snhei.GE.snth)
THEN
5388 if(snhei.le.deltsn+snth)
then
5392 IF (debug_print )
THEN
5393 print *,
'1 layer d1sn,d2sn',i,j,d1sn,d2sn
5399 IF (debug_print )
THEN
5400 print *,
'2 layers d1sn,d2sn',i,j,d1sn,d2sn
5403 d9sn= thdifsn*rhocsn / snprim
5404 r22sn = snprim*snprim*0.5_kind_phys/(thdifsn*delt)
5405 IF (debug_print )
THEN
5406 print *,
'1 or 2 layers D9sn,R22sn',d9sn,r22sn
5410 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
5414 d9sn = (fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)/ &
5416 r22sn = snprim*snprim*0.5_kind_phys &
5417 /((fsn*thdifsn+fso*thdif(1))*delt)
5418 IF (debug_print )
THEN
5419 print *,
' Combined D9SN,R22SN,D1SN,D2SN: ',d9sn,r22sn,d1sn,d2sn
5422 IF(snhei.eq.zero)
then
5428 IF (debug_print )
THEN
5429 print *,
' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',d9sn,r22sn,d1sn,d2sn
5439 tdenom = d9sn*(one-d1sn +r22sn)+d10+r21+r7 &
5441 +rhonewcsn*newsnow/delt
5447 aa=xlvm*(beta*fkq*umveg+r210)/tdenom
5448 bb=(d10*tabs+r21*tn+xlvm*(qvatm* &
5449 (beta*fkq*umveg+c) &
5450 +r210*qgold)+d11+d9sn*(d2sn+r22sn*tn) &
5451 +rainf*cvw*prcpms*max(tfrz,tabs) &
5452 + rhonewcsn*newsnow/delt*min(tfrz,tabs) &
5455 pp=patm*1.e3_kind_phys
5459 IF (debug_print )
THEN
5460 if (abs(xlat-33.35).lt.0.2 .and. abs(xlon-272.55).lt.0.2)
then
5461 print *,
'1-', i,rnet,tabs,tn,aa1,bb,pp,ktau,newsnow,snwepr,snwe,snhei,snowfrac,soilt,soilt1,tso,rhosn
5462 print *,
'2-', i,tdenom,fkq,vegfrac,can,r210,d10,r21,d9sn,d1sn,r22sn,r7,prcpms
5465 CALL vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
5469 IF (debug_print )
THEN
5471 print *,
'VILKA1 - TS1,QS1,TQ2,H,TX2,Q1',ts1,qs1,tq2,h,tx2,q1,xlat,xlon
5473 IF(q1.LT.qs1)
GOTO 100
5478 qcg=max(zero,q1-qs1)
5479 IF (debug_print )
THEN
5480 print *,
'90 QVG,QSG,QCG,TSO(1)',qvg,qsg,qcg,tso(1)
5485 CALL vilka(tn,aa,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
5487 IF (debug_print )
THEN
5489 print *,
'VILKA2 - TS1,QS1,H,TX2,Q1',ts1,qs1,tq2,h,tx2,q1
5491 IF(q1.GT.qs1)
GOTO 90
5495 IF (debug_print )
THEN
5496 print *,
'No Saturation QVG,QSG,QCG,TSO(1)',qvg,qsg,qcg,tso(1)
5502 if(nmelt==1 .and. snowfrac==one .and. snwe > zero .and. soilt > tfrz)
then
5507 if (debug_print )
then
5509 print *,
'soilt is too high =',soilt,xlat,xlon
5510 soilt = min(tfrz,soilt)
5514 IF (debug_print )
THEN
5516 print *,
'snwe,snwepr,snhei,snowfr,soilt,soilt1,tso',i,j,snwe,snwepr,snhei,snowfrac,soilt,soilt1,tso
5519 IF(snhei.GE.snth)
THEN
5520 if(snhei.gt.deltsn+snth)
then
5522 soilt1=rhtsn+cotsn*soilt
5523 tso(1)=rhtso(nzs)+cotso(nzs)*soilt1
5527 tso(1)=rhtso(nzs)+cotso(nzs)*soilt
5531 ELSEIF (snhei > zero .and. snhei < snth)
THEN
5533 tso(2)=rhtso(nzs1)+cotso(nzs1)*soilt
5534 tso(1)=(tso(2)+(soilt-tso(2))*fso)
5544 if(nmelt==1.and.snowfrac==one)
then
5546 soilt1= min(tfrz,soilt1)
5547 tso(1)= min(tfrz,tso(1))
5548 tsob = min(tfrz,tsob)
5552 IF (snhei > zero .and. snhei < snth)
THEN
5556 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
5562 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
5575 IF (debug_print )
THEN
5577 print *,
'Final SOILT,SOILT1,tso,TSOB,QSG',xlat,xlon,soilt,soilt1,tso,tsob,qsg,
'nmelt=',nmelt
5578 print *,
'SNWEPR-BETA*EPOT*RAS*DELT',snwepr-beta*epot*ras*delt,beta,snwepr,epot
5581 if(nmelt.eq.1)
go to 220
5585 IF(soilt > tfrz.AND.beta==one.AND.snhei>zero)
THEN
5588 soiltfrac=snowfrac*tfrz+(one-snowfrac)*soilt
5589 qsg=min(qsg, qsn(soiltfrac,tbq)/pp)
5590 qvg=snowfrac*qsg+(one-snowfrac)*qvg
5591 t3 = stbolt*tn*tn*tn
5592 upflux = t3 * 0.5_kind_phys*(tn + soiltfrac)
5593 xinet = emiss*(glw-upflux)
5594 epot = -qkms*(qvatm-qsg)
5598 IF (q1.LE.0..or.iter==1)
THEN
5610 transp(k)=-vegfrac*q1 &
5611 *tranf(k)*drycan/zshalf(nroot+1)
5618 edir1 = q1*umveg * beta
5619 ec1 = q1 * wetcan * vegfrac
5621 eeta = (edir1 + ec1 + ett1)*rhowater
5626 hfx=-d10*(tabs-soiltfrac)
5628 IF(snhei.GE.snth)
then
5629 soh=thdifsn*rhocsn*(soiltfrac-tsob)/snprim
5632 soh=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5633 (soiltfrac-tsob)/snprim
5638 x= (r21+d9sn*r22sn)*(soiltfrac-tn) + &
5639 xlvm*r210*(qvg-qgold)
5640 IF (debug_print )
THEN
5641 print *,
'SNOWTEMP storage ',i,j,x
5642 print *,
'R21,D9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim', &
5643 r21,d9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim
5647 snoh=rnet-qfx -hfx - soh - x &
5648 +rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soiltfrac) &
5649 +rainf*cvw*prcpms*(max(tfrz,tabs)-soiltfrac)
5652 smelt= snoh /xlmelt*1.e-3_kind_phys
5653 IF (debug_print )
THEN
5655 print *,
'1- SMELT',smelt,snoh,xlat,xlon
5658 IF(epot.gt.zero .and. snwepr.LE.epot*ras*delt)
THEN
5660 beta=snwepr/(epot*ras*delt)
5661 smelt=amax1(zero,amin1(smelt,snwepr/delt-beta*epot*ras))
5663 IF (debug_print )
THEN
5665 print *,
'2- SMELT',xlat,xlon,snwe,smelt,rhonewsn,xlat,xlon
5673 if( (rhosn < 350._kind_phys .or. (newsnow > zero .and. rhonewsn < 450._kind_phys)) .and. soilt < 283._kind_phys )
then
5674 smelt= amin1(smelt, delt/60._kind_phys*5.6e-8_kind_phys*meltfactor*max(one,(soilt-tfrz)))
5675 IF (debug_print )
THEN
5677 print *,
'3- SMELT',xlat,xlon,smelt,rhosn,rhonewsn,xlat,xlon
5682 rr=max(zero,snwepr/delt-beta*epot*ras)
5684 smelt = min(smelt,rr)
5686 IF (debug_print )
THEN
5688 print *,
'4- SMELT i,j,smelt,rr',xlat,xlon,smelt,rr
5692 snohgnew=smelt*xlmelt*rhowater
5693 snodif=amax1(zero,(snoh-snohgnew))
5696 IF (debug_print )
THEN
5698 print *,
'SNOH,SNODIF',snoh,snodif
5699 print *,
' xlat, xlon', xlat, xlon
5702 IF( smelt > zero)
then
5704 rsmfrac=min(0.18_kind_phys,(max(0.08_kind_phys,snwepr/0.10_kind_phys*0.13_kind_phys)))
5705 if(snhei > 0.01_kind_phys .and. rhosn < 350._kind_phys)
then
5706 rsm=min(snwe,rsmfrac*smelt*delt)
5714 smelt=max(zero,smelt-rsm/delt)
5715 IF (debug_print )
THEN
5717 print *,
'5- SMELT i,j,smelt,rsm,snwepr,rsmfrac', &
5718 i,j,smelt,rsm,snwepr,rsmfrac
5719 print *,
' xlat, xlon', xlat, xlon
5727 if(snwe > zero)
then
5728 snwe = amax1(zero,(snwepr- &
5729 (smelt+beta*epot*ras)*delt &
5731 IF (debug_print )
THEN
5733 print *,
' Snow is melting and sublimating, snwe', xlat, xlon, snwe
5737 IF (debug_print )
THEN
5739 print *,
' all snwe is sublimated or melted', xlat, xlon, snwe
5746 if(snhei.ne.zero .and. beta == one)
then
5747 epot=-qkms*(qvatm-qsg)
5748 snwe = amax1(zero,(snwepr- &
5749 beta*epot*ras*delt))
5759 if(nmelt.eq.1)
goto 212
5762 if(smelt > zero .and. rsm > zero)
then
5763 if(snwe.le.rsm)
then
5764 IF ( debug_print )
THEN
5765 print *,
'SNWE<RSM snwe,rsm,smelt*delt,epot*ras*delt,beta', &
5766 snwe,rsm,smelt*delt,epot*ras*delt,beta
5773 xsn=(rhosn*(snwe-rsm)+rhowater*rsm)/ &
5775 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
5777 rhocsn=sheatsn* rhosn
5778 if(isncond_opt == 1)
then
5780 thdifsn = 0.265_kind_phys/rhocsn
5785 if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys))
then
5786 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
5788 keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys
5789 if(debug_print)
then
5790 print *,
'End SNOWTEMP - xlat,xlon,rhosn,keff',xlat,xlon,rhosn,keff
5791 print *,
'End SNOWTEMP - 0.265/rhocsn',0.265/rhocsn
5794 if (debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5795 print *,
'END SNOWTEMP - newsnow, rhonewsn,rhosn,fact,keff', &
5796 xlat,xlon,newsnow, rhonewsn,rhosn,fact,keff,keff/rhocsn*fact
5799 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
5805 thdifsn = 4.431718e-7_kind_phys
5807 thdifsn = keff/rhocsn * fact
5811 if (debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5812 print *,
'END SNOWTEMP - thdifsn',xlat,xlon,thdifsn
5813 print *,
'END SNOWTEMP - 0.265/rhocsn',0.265/rhocsn
5819 IF(snhei.GE.snth)
then
5820 s=thdifsn*rhocsn*(soilt-tsob)/snprim
5822 s=d9*(tso(1)-tso(2))
5823 ELSEIF(snhei.lt.snth.and.snhei.GT.zero)
then
5824 s=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5827 s=d9*(tso(1)-tso(2))
5831 s=d9*(tso(1)-tso(2))
5835 snhei=snwe * rhowater / rhosn
5841 IF (debug_print )
THEN
5843 print *,
'snhei,snwe,rhosn,snowfr',snhei,snwe,rhosn,snowfrac,xlat,xlon
5846 IF(tso(1).GT.tfrz .and. snhei > zero)
THEN
5848 if (snhei.GT.deltsn+snth)
then
5849 hsn = snhei - deltsn
5850 IF (debug_print )
THEN
5851 print*,
'2 layer snow - snhei,hsn',snhei,hsn
5854 IF (debug_print )
THEN
5855 print*,
'1 layer snow or blended - snhei',snhei
5860 soiltfrac=snowfrac*tfrz+(one-snowfrac)*tso(1)
5862 snohg=(tso(1)-soiltfrac)*(cap(1)*zshalf(2)+ &
5863 rhocsn*0.5_kind_phys*hsn) / delt
5864 snohg=amax1(zero,snohg)
5866 smeltg=snohg/xlmelt*1.e-3_kind_phys
5867 IF (debug_print )
THEN
5869 print *,
' SMELTG =',smeltg,xlat,xlon
5873 if( (rhosn < 350._kind_phys .or. (newsnow > zero .and. rhonewsn < 450._kind_phys)) .and. soilt < 283._kind_phys )
then
5874 smelt=amin1(smeltg, 5.8e-9_kind_phys)
5879 smeltg=amin1(smeltg, rr)
5881 snohgnew=smeltg*xlmelt*rhowater
5882 snodif=amax1(zero,(snohg-snohgnew))
5883 IF (debug_print )
THEN
5885 print *,
'TSO(1),soiltfrac,snowfrac,smeltg,SNODIF',tso(1),soiltfrac,snowfrac,smeltg,snodif
5886 print *,
' xlat, xlon', xlat, xlon
5889 snwe=max(zero,snwe-smeltg*delt)
5890 snhei=snwe * rhowater / rhosn
5892 smelt = smelt + smeltg
5894 if(snhei > zero) tso(1) = soiltfrac
5896 IF (debug_print )
THEN
5898 print *,
'Melt from the bottom snwe,snhei',snwe,snhei
5899 print *,
' xlat, xlon', xlat, xlon
5900 print *,
'TSO(1),soiltfrac,snowfrac,smeltg,SNODIF',tso(1),soiltfrac,snowfrac,smeltg,snodif
5901 print *,
'Melt from the bottom snwe,snhei,snoh',snwe,snhei,snoh
5902 print *,
' Final TSO ',tso
5904 print *,
'Snow is all melted on the warm ground'
5910 snheiprint=snweprint*rhowater / rhosn
5912 x= (r21+d9sn*r22sn)*(soilt-tn) + &
5913 xlvm*r210*(qsg-qgold)
5914 IF (debug_print )
THEN
5916 print *,
'end SNOWTEMP storage ',xlat,xlon,x
5917 print *,
'R21,D9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim', &
5918 r21,d9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim
5919 print *,
'snwe, snhei ',snwe,snhei
5924 -rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soilt) &
5925 -rainf*cvw*prcpms*(max(tfrz,tabs)-soilt)
5926 IF (debug_print )
THEN
5928 print *,
'SNHEI=',snhei
5929 print *,
'SNFLX=',snflx
5932 IF(snhei.GT.zero)
THEN
5934 tsnav=min(zero,0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
5935 +(soilt1+tso(1))*(snhei-deltsn)) &
5938 tsnav=min(zero,0.5_kind_phys*(soilt+tso(1)) - tfrz)
5941 tsnav= min(zero,soilt - tfrz)
5945 END SUBROUTINE snowtemp
5951 SUBROUTINE soilmoist ( debug_print, &
5952 xlat, xlon, testptlat, testptlon, &
5953 DELT,NZS,NDDZS,DTDZS,DTDZS2,RIW, & !--- input parameters
5954 ZSMAIN,ZSHALF,DIFFU,HYDRO, &
5955 QSG,QVG,QCG,QCATM,QVATM,PRCP, &
5957 DEW,SMELT,SOILICE,VEGFRAC,SNOWFRAC,soilres, &
5958 DQM,QMIN,REF,KSAT,RAS,INFMAX, & !--- soil properties
5959 SOILMOIS,SOILIQW,MAVAIL,RUNOFF,RUNOFF2,INFILTRP)
6004 LOGICAL,
INTENT(IN ) :: debug_print
6005 real (kind_phys),
INTENT(IN ) :: delt
6006 real (kind_phys),
INTENT(IN ) :: xlat, xlon, testptlat, testptlon
6007 INTEGER,
INTENT(IN ) :: nzs,nddzs
6011 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: zsmain, &
6019 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: dtdzs
6021 real (kind_phys),
INTENT(IN ) :: qsg,qvg,qcg,qcatm,qvatm, &
6022 qkms,vegfrac,drip,prcp , &
6023 dew,smelt,snowfrac , &
6024 dqm,qmin,ref,ksat,ras,riw,soilres
6028 real (kind_phys),
DIMENSION( 1:nzs ) , &
6029 INTENT(INOUT) :: soilmois,soiliqw
6031 real (kind_phys),
INTENT(INOUT) :: mavail,runoff,runoff2,infiltrp, &
6036 real (kind_phys),
DIMENSION( 1:nzs ) :: cosmc,rhsmc
6038 real (kind_phys) :: dzs,r1,r2,r3,r4,r5,r6,r7,r8,r9,r10
6039 real (kind_phys) :: refkdt,refdk,delt1,f1max,f2max
6040 real (kind_phys) :: f1,f2,fd,kdt,val,ddt,px,fk,fkmax
6041 real (kind_phys) :: qq,umveg,infmax1,trans
6042 real (kind_phys) :: totliq,flx,flxsat,qtot
6043 real (kind_phys) :: did,x1,x2,x4,denom,q2,q4
6044 real (kind_phys) :: dice,fcr,acrt,frzx,sum,cvfrz
6046 INTEGER :: nzs1,nzs2,k,kk,k1,kn,ialp1,jj,jk
6054 118
format(6(10pf23.19))
6061 did=(zsmain(nzs)-zshalf(nzs))
6062 x1=zsmain(nzs)-zsmain(nzs1)
6064 denom=(one+diffu(nzs1)/x1/did*delt+hydro(nzs)/(2._kind_phys*did)*delt)
6065 cosmc(1)=delt*(diffu(nzs1)/did/x1 &
6066 +hydro(nzs1)/2._kind_phys/did)/denom
6067 rhsmc(1)=(soilmois(nzs)+transp(nzs)*delt/ &
6074 denom=1.+diffu(nzs1)/x1/did*delt
6075 cosmc(1)=delt*(diffu(nzs1)/did/x1 &
6076 +hydro(nzs1)/did)/denom
6078 rhsmc(1)=(soilmois(nzs)-hydro(nzs)*delt/did*soilmois(nzs) &
6079 +transp(nzs)*delt/did)/denom
6082 rhsmc(1)=soilmois(nzs)
6087 x4=2.*dtdzs(k1)*diffu(kn-1)
6088 x2=2.*dtdzs(k1+1)*diffu(kn)
6089 q4=x4+hydro(kn-1)*dtdzs2(kn-1)
6090 q2=x2-hydro(kn+1)*dtdzs2(kn-1)
6091 denom=one+x2+x4-q2*cosmc(k)
6093 IF (debug_print )
THEN
6094 if (abs(xlat-testptlat).lt.0.05 .and. &
6095 abs(xlon-testptlon).lt.0.05)
then
6096 print *,
'xlat,xlon=',xlat,xlon
6097 print *,
'q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k' &
6098 ,q2,soilmois(kn),diffu(kn),x2,hydro(kn+1),dtdzs2(kn-1),kn,k
6101 rhsmc(k+1)=(soilmois(kn)+q2*rhsmc(k) &
6103 /(zshalf(kn+1)-zshalf(kn)) &
6110 umveg=(one-vegfrac)*soilres
6118 r4=r3+hydro(1)*.5_kind_phys
6119 r5=r3-hydro(2)*.5_kind_phys
6128 totliq=prcp-drip/delt-(one-vegfrac)*dew*ras-smelt
6129 IF (debug_print )
THEN
6130 if (abs(xlat-testptlat).lt.0.05 .and. &
6131 abs(xlon-testptlon).lt.0.05)
then
6132 print *,
'xlat,xlon=',xlat,xlon
6133 print *,
'UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT', &
6134 umveg*prcp,drip/delt,umveg*dew*ras,smelt
6150 cvfrz = 3._kind_phys
6154 refdk=3.4341e-6_kind_phys
6155 delt1=delt/86400._kind_phys
6157 f2max=dqm*(zshalf(3)-zshalf(2))
6158 f1=f1max*(one-soilmois(1)/dqm)
6159 dice=soilice(1)*zshalf(2)
6162 dice=dice+(zshalf(k+1)-zshalf(k))*soilice(k)
6163 fkmax=dqm*(zshalf(k+1)-zshalf(k))
6164 fk=fkmax*(one-soilmois(k)/dqm)
6167 kdt=refkdt*ksat/refdk
6168 val=(1.-exp(-kdt*delt1))
6171 IF(px < zero) px = zero
6173 infmax1 = (px*(ddt/(px+ddt)))/delt
6177 IF (debug_print )
THEN
6178 print *,
'INFMAX1 before frozen part',infmax1
6186 frzx= 0.15_kind_phys*((dqm+qmin)/ref) * (0.412_kind_phys / 0.468_kind_phys)
6189 IF ( dice .GT. 1.e-2_kind_phys)
THEN
6190 acrt = cvfrz * frzx / dice
6198 sum = sum + (acrt ** ( cvfrz-jk)) / float(k)
6200 fcr = one - exp(-acrt) * sum
6202 IF (debug_print )
THEN
6203 print *,
'FCR--------',fcr
6204 print *,
'DICE=',dice
6206 infmax1 = infmax1* fcr
6209 infmax = max(infmax1,hydro(1)*soilmois(1))
6210 infmax = min(infmax, -totliq)
6211 IF (debug_print )
THEN
6212 print *,
'INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ', &
6213 infmax,infmax1,hydro(1)*soiliqw(1),-totliq
6216 IF (-totliq.GT.infmax)
THEN
6217 runoff=-totliq-infmax
6219 IF (debug_print )
THEN
6220 print *,
'FLX,RUNOFF1=',flx,runoff
6226 r7=.5_kind_phys*dzs/delt
6228 flx=flx-soilmois(1)*r7
6231 r8=umveg*r6*(one-snowfrac)
6237 IF(r10.LE.zero)
THEN
6238 qq=(r5*r2-flx+r9)/(r4-r5*r1-r10*r8/(ref-qmin))
6239 flxsat=-dqm*(r4-r5*r1-r10*r8/(ref-qmin)) &
6243 qq=(r2*r5-flx+r8*(qtot-qcg-qvg)+r9)/(r4-r1*r5)
6244 flxsat=-dqm*(r4-r1*r5)+r2*r5+r8*(qtot-qvg-qcg)+r9
6249 soilmois(1)=1.e-8_kind_phys
6251 ELSE IF(qq.GT.dqm)
THEN
6254 IF (debug_print )
THEN
6255 print *,
'FLXSAT,FLX,DELT',flxsat,flx,delt,runoff2
6257 runoff=runoff+(flxsat-flx)
6259 soilmois(1)=min(dqm,max(1.e-8_kind_phys,qq))
6262 IF (debug_print )
THEN
6263 if (abs(xlat-testptlat).lt.0.05 .and. &
6264 abs(xlon-testptlon).lt.0.05)
then
6265 print *,
'xlat,xlon=',xlat,xlon
6266 print *,
'SOILMOIS,SOILIQW, soilice',soilmois,soiliqw,soilice*riw
6267 print *,
'COSMC,RHSMC',cosmc,rhsmc
6274 qq=cosmc(kk)*soilmois(k-1)+rhsmc(kk)
6276 IF (qq.LT.zero)
THEN
6278 ELSE IF(qq.GT.dqm)
THEN
6282 IF (debug_print )
THEN
6283 print *,
'hydro(k),QQ,DQM,k',hydro(k),qq,dqm,k
6285 runoff2=runoff2+((qq-dqm)*(zsmain(k)-zshalf(k)))/delt
6287 runoff2=runoff2+((qq-dqm)*(zshalf(k+1)-zshalf(k)))/delt
6290 soilmois(k)=min(dqm,max(1.e-8_kind_phys,qq))
6293 IF (debug_print )
THEN
6294 if (abs(xlat-testptlat).lt.0.05 .and. &
6295 abs(xlon-testptlon).lt.0.05)
then
6296 print *,
'xlat,xlon=',xlat,xlon
6297 print *,
'END soilmois,soiliqw,soilice',soilmois,soiliqw,soilice*riw
6301 mavail=max(.00001_kind_phys,min(one,(soilmois(1)/(ref-qmin)*(one-snowfrac)+one*snowfrac)))
6303 END SUBROUTINE soilmoist
6309 SUBROUTINE soilprop( debug_print, &
6310 xlat, xlon, testptlat, testptlon, &
6311 nzs,fwsat,lwsat,tav,keepfr, & !--- input variables
6312 soilmois,soiliqw,soilice, &
6313 soilmoism,soiliqwm,soilicem, &
6314 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, & !--- soil fixed fields
6315 riw,xlmelt,CP,G0_P,cvw,ci, & !--- constants
6317 thdif,diffu,hydro,cap)
6342 LOGICAL,
INTENT(IN ) :: debug_print
6343 INTEGER,
INTENT(IN ) :: nzs
6344 real (kind_phys),
INTENT(IN ) :: xlat, xlon, testptlat, testptlon
6346 real (kind_phys) , &
6347 INTENT(IN ) :: rhocs, &
6355 real (kind_phys),
DIMENSION( 1:nzs ) , &
6356 INTENT(IN ) :: soilmois, &
6360 real (kind_phys),
INTENT(IN ) :: cp, &
6372 real (kind_phys),
DIMENSION(1:NZS) , &
6373 INTENT(INOUT) :: cap,diffu,hydro , &
6377 soilicem,soiliqwm , &
6381 real (kind_phys),
DIMENSION(1:NZS) :: hk,detal,kasat,kjpl
6383 real (kind_phys) :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci
6384 real (kind_phys) :: tln,tavln,tn,pf,a,am,ame,h
6388 real (kind_phys) :: kzero,gamd,kdry,kas,x5,sr,ke
6405 x1=xlmelt/(g0_p*psis)
6409 gamd=(one-ws)*2700._kind_phys
6410 kdry=(0.135_kind_phys*gamd+64.7_kind_phys)/(2700._kind_phys-0.947_kind_phys*gamd)
6412 if(qwrtz > 0.2_kind_phys)
then
6413 kas=kqwrtz**qwrtz*kzero**(1.-qwrtz)
6415 kas=kqwrtz**qwrtz*3._kind_phys**(one-qwrtz)
6420 wd=ws - riw*soilicem(k)
6421 psif=psis*100._kind_phys*(wd/(soiliqwm(k)+qmin))**bclh &
6422 * (ws/wd)**3._kind_phys
6425 fact=one+riw*soilicem(k)
6427 IF(pf.LE.5.2_kind_phys)
THEN
6428 hk(k)=420._kind_phys*exp(-(pf+2.7_kind_phys))*fact
6430 hk(k)=.1744_kind_phys*fact
6433 IF(soilicem(k).NE.zero.AND.tn.LT.zero)
then
6437 detal(k)=tfrz*x2/(tav(k)*tav(k))* &
6438 (tav(k)/(x1*tn))**x4
6440 if(keepfr(k).eq.one)
then
6447 kasat(k)=kas**(one-ws)*kice**fwsat(k) &
6450 x5=(soilmoism(k)+qmin)/ws
6451 if(soilicem(k).eq.zero)
then
6452 sr=max(0.101_kind_phys,x5)
6458 kjpl(k)=ke*(kasat(k)-kdry)+kdry
6461 cap(k)=(one-ws)*rhocs &
6462 + (soiliqwm(k)+qmin)*cvw &
6464 + (dqm-soilmoism(k))*cp*1.2_kind_phys &
6465 - detal(k)*rhowater*xlmelt
6469 if((ws-a).lt.0.12_kind_phys)
then
6472 h=max(zero,(soilmoism(k)+qmin-a)/(max(1.e-8_kind_phys,(ws-a))))
6474 if(a.ne.zero)facd=one-a/max(1.e-8_kind_phys,soilmoism(k))
6475 ame=max(1.e-8_kind_phys,ws-riw*soilicem(k))
6477 diffu(k)=-bclh*ksat*psis/ame* &
6478 (ws/ame)**3._kind_phys &
6479 *h**(bclh+2._kind_phys)*facd
6485 thdif(k)=kjpl(k)/cap(k)
6489 IF (debug_print )
THEN
6490 print *,
'soilice*riw,soiliqw,soilmois,ws',soilice*riw,soiliqw,soilmois,ws
6494 if((ws-riw*soilice(k)).lt.0.12_kind_phys)
then
6498 if(soilice(k).ne.zero) &
6499 fach=one-riw*soilice(k)/max(1.e-8_kind_phys,soilmois(k))
6500 am=max(1.e-8_kind_phys,ws-riw*soilice(k))
6502 hydro(k)=min(ksat,ksat/am* &
6504 **(2._kind_phys*bclh+2._kind_phys) &
6506 if(hydro(k)<1.e-10_kind_phys)hydro(k)=zero
6510 IF (debug_print )
THEN
6511 print *,
'hydro=',hydro
6515 END SUBROUTINE soilprop
6521 SUBROUTINE transf( debug_print, &
6522 xlat,xlon,testptlat,testptlon, &
6523 nzs,nroot,soiliqw,tabs,lai,gswin, & !--- input variables
6524 dqm,qmin,ref,wilt,zshalf,pc,iland, & !--- soil fixed fields
6541 LOGICAL,
INTENT(IN ) :: debug_print
6542 INTEGER,
INTENT(IN ) :: nroot,nzs,iland
6543 real (kind_phys),
INTENT(IN ) :: xlat,xlon,testptlat,testptlon
6545 real (kind_phys) , &
6546 INTENT(IN ) :: gswin, tabs, lai
6548 real (kind_phys) , &
6549 INTENT(IN ) :: dqm, &
6555 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: soiliqw, &
6559 real (kind_phys),
DIMENSION(1:NZS),
INTENT(OUT) :: tranf
6560 real (kind_phys),
INTENT(OUT) :: transum
6563 real (kind_phys) :: totliq, did
6567 real (kind_phys) :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4
6568 real (kind_phys) :: ftem, pctot, fsol, f1, cmin, cmax, totcnd
6569 real (kind_phys),
DIMENSION(1:NZS) :: part
6578 totliq=soiliqw(1)+qmin
6584 ap1=-8.152_kind_phys
6585 ap2=61.653_kind_phys
6586 ap3=-115.876_kind_phys
6587 ap4=59.656_kind_phys
6588 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
6589 if(totliq.ge.ref) gx=one
6590 if(totliq.le.wilt) gx=zero
6591 if(gx.gt.one) gx=one
6592 if(gx.lt.zero) gx=zero
6595 IF(totliq.GT.ref)
THEN
6597 ELSE IF(totliq.LE.wilt)
THEN
6600 tranf(1)=(totliq-wilt)/(ref-wilt)*did
6606 totliq=soiliqw(k)+qmin
6611 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
6612 if(totliq.ge.ref) gx=one
6613 if(totliq.le.wilt) gx=zero
6614 if(gx.gt.one) gx=one
6615 if(gx.lt.zero) gx=zero
6616 did=zshalf(k+1)-zshalf(k)
6618 IF(totliq.GE.ref)
THEN
6620 ELSE IF(totliq.LE.wilt)
THEN
6623 tranf(k)=(totliq-wilt) &
6629 IF (debug_print )
THEN
6630 if (abs(xlat-testptlat).lt.0.05 .and. &
6631 abs(xlon-testptlon).lt.0.05)
then
6632 print *,
'xlat,xlon=',xlat,xlon
6633 print *,
'soiliqw =',soiliqw,
'wilt=',wilt,
'qmin= ',qmin
6634 print *,
'tranf = ',tranf
6639 if(lai > 4._kind_phys)
then
6647 IF ( debug_print )
THEN
6648 if (abs(xlat-testptlat).lt.0.05 .and. &
6649 abs(xlon-testptlon).lt.0.05)
then
6650 print *,
'xlat,xlon=',xlat,xlon
6651 print *,
'pctot,lai,pc',pctot,lai,pc
6657 IF (tabs .LE. 302.15_kind_phys)
THEN
6658 ftem = one / (one + exp(-0.41_kind_phys * (tabs - 282.05_kind_phys)))
6660 ftem = one / (one + exp(0.5_kind_phys * (tabs - 314.0_kind_phys)))
6663 cmin = one/rsmax_data
6664 cmax = one/rstbl(iland)
6666 cmax = lai/rstbl(iland)
6679 if (gswin < rgltbl(iland))
then
6680 fsol = one / (one + exp(-0.034_kind_phys * (gswin - 3.5_kind_phys)))
6685 totcnd =(cmin + (cmax - cmin)*pctot*ftem*fsol)/cmax
6687 IF ( debug_print )
THEN
6688 if (abs(xlat-testptlat).lt.0.05 .and. &
6689 abs(xlon-testptlon).lt.0.05)
then
6690 print *,
'xlat,xlon=',xlat,xlon
6691 print *,
'GSWin,Tabs,lai,f1,cmax,cmin,pc,pctot,ftem,fsol',gswin,tabs,lai,f1,cmax,cmin,pc,pctot,ftem,fsol
6692 print *,
'iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd' &
6693 ,iland,rgltbl(iland),rstbl(iland),rsmax_data,totcnd
6701 tranf(k)=max(zero,tranf(k)*totcnd)
6702 transum=transum+tranf(k)
6704 IF ( debug_print )
THEN
6705 if (abs(xlat-testptlat).lt.0.05 .and. &
6706 abs(xlon-testptlon).lt.0.05)
then
6707 print *,
'xlat,xlon=',xlat,xlon
6708 print *,
'transum,TRANF',transum,tranf
6713 END SUBROUTINE transf
6720 SUBROUTINE vilka(TN,D1,D2,PP,QS,TS,TT,NSTEP,ii,j,iland,isoil,xlat,xlon)
6725 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN ) :: tt
6726 real (kind_phys),
INTENT(IN ) :: tn,d1,d2,pp,xlat,xlon
6727 INTEGER,
INTENT(IN ) :: nstep,ii,j,iland,isoil
6729 real (kind_phys),
INTENT(OUT ) :: qs, ts
6731 real (kind_phys) :: f1,t1,t2,rn
6734 i=(tn-1.7315e2_kind_dbl_prec)/.05_kind_dbl_prec+1
6735 t1=173.1_kind_dbl_prec+float(i)*.05_kind_dbl_prec
6737 i1=i-f1/(.05_kind_dbl_prec+d1*(tt(i+1)-tt(i)))
6739 IF(i.GT.5000.OR.i.LT.1)
GOTO 1
6741 t1=173.1_kind_dbl_prec+float(i)*.05_kind_dbl_prec
6743 rn=f1/(.05_kind_dbl_prec+d1*(tt(i+1)-tt(i)))
6745 IF(i.GT.5000.OR.i.LT.1)
GOTO 1
6747 ts=t1-.05_kind_dbl_prec*rn
6748 qs=(tt(i)+(tt(i)-tt(i+1))*rn)/pp
6750 1 print *,
' AVOST IN VILKA Table index= ',i
6751 print *,
'I,J=',ii,j,
'LU_index = ',iland,
'Psfc[hPa] = ',pp,
'Tsfc = ',tn
6752 print *,
'AVOST point at xlat/xlon=',xlat,xlon
6755 END SUBROUTINE vilka
6762 SUBROUTINE soilvegin ( debug_print,mosaic_lu,mosaic_soil, &
6763 soilfrac,nscat,shdmin, shdmax, &
6764 NLCAT,IVGTYP,ISLTYP,iswater,MYJ, &
6765 IFOREST,lufrac,vegfrac,EMISS,PC, &
6766 MSNF,FACSNF,ZNT,LAI,RDLAI2D, &
6767 QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,I,J, &
6792 integer,
parameter :: nsoilclas=19
6793 integer,
parameter :: nvegclas=24+3
6794 integer,
parameter :: ilsnow=99
6796 LOGICAL,
INTENT(IN ) :: debug_print
6797 INTEGER,
INTENT(IN ) :: mosaic_lu, mosaic_soil
6798 INTEGER,
INTENT(IN ) :: nlcat, nscat, iswater, i, j
6824 real (kind_phys) lqma(nsoilclas),lrhc(nsoilclas), &
6825 lpsi(nsoilclas),lqmi(nsoilclas), &
6826 lbcl(nsoilclas),lkas(nsoilclas), &
6827 lwil(nsoilclas),lref(nsoilclas), &
6837 DATA lqma /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
6838 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
6839 0.20, 0.435, 0.468, 0.200, 0.339/
6846 DATA lref /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
6847 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
6848 0.1, 0.249, 0.454, 0.17, 0.236/
6855 DATA lwil/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
6856 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
6857 0.006, 0.114, 0.030, 0.006, 0.01/
6863 DATA lqmi/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
6864 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
6865 0.004, 0.065, 0.020, 0.004, 0.008/
6873 DATA lpsi/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
6874 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
6875 0.121, 0.218, 0.468, 0.069, 0.069/
6883 DATA lkas/1.76e-4, 1.56e-4, 3.47e-5, 7.20e-6, 7.20e-6, &
6884 6.95e-6, 6.30e-6, 1.70e-6, 2.45e-6, 2.17e-6, &
6885 1.03e-6, 1.28e-6, 6.95e-6, 0.0, 1.41e-4, &
6886 3.47e-5, 1.28e-6, 1.41e-4, 1.76e-4/
6893 DATA lbcl/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
6894 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
6895 4.05, 4.90, 11.55, 2.79, 2.79/
6897 DATA lrhc /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
6898 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
6900 DATA datqtz/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
6901 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
6961 real (kind_phys) lalb(nvegclas),lmoi(nvegclas),lemi(nvegclas), &
6962 lrou(nvegclas),lthi(nvegclas),lsig(nvegclas), &
6970 DATA lalb/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
6971 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
6973 DATA lemi/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
6974 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
6977 DATA lrou/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
6978 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
6981 DATA lmoi/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
6982 .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95,.40,.50,.40/
6986 DATA lpc /0.4,0.3,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,5*0.55,0.,0.55,0.55, &
6987 0.3,0.3,0.4,0.4,0.3,0.,.3,0.,0./
6995 LOGICAL,
INTENT(IN ) :: myj
6996 real (kind_phys),
INTENT(IN ) :: shdmax
6997 real (kind_phys),
INTENT(IN ) :: shdmin
6998 real (kind_phys),
INTENT(IN ) :: vegfrac
6999 real (kind_phys),
DIMENSION( 1:NLCAT ),
INTENT(IN):: lufrac
7000 real (kind_phys),
DIMENSION( 1:NSCAT ),
INTENT(IN):: soilfrac
7002 real (kind_phys) , &
7003 INTENT ( OUT) :: pc, &
7007 real (kind_phys) , &
7008 INTENT (INOUT ) :: emiss, &
7011 LOGICAL,
intent(in) :: rdlai2d
7013 real (kind_phys) , &
7014 INTENT( OUT) :: rhocs, &
7023 INTEGER,
INTENT ( OUT) :: iforest
7024 character(len=*),
intent(out) :: errmsg
7025 integer,
intent(out) :: errflg
7026 INTEGER :: kstart, kfin, lstart, lfin
7028 real (kind_phys) :: area, factor, znt1, lb
7029 real (kind_phys),
DIMENSION( 1:NLCAT ) :: znttoday, laitoday, deltalai
7040 iforest = ifortbl(ivgtyp)
7042 IF (debug_print )
THEN
7043 print *,
'ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)', &
7044 ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)
7053 if((shdmax - shdmin) .lt. one)
then
7056 factor = one - max(zero,min(one,(vegfrac - shdmin)/max(one,(shdmax-shdmin))))
7061 if(ifortbl(k) == 1) deltalai(k)=min(0.2_kind_phys,0.8_kind_phys*laitbl(k))
7062 if(ifortbl(k) == 2 .or. ifortbl(k) == 7) deltalai(k)=min(0.5_kind_phys,0.8_kind_phys*laitbl(k))
7063 if(ifortbl(k) == 3) deltalai(k)=min(0.45_kind_phys,0.8_kind_phys*laitbl(k))
7064 if(ifortbl(k) == 4) deltalai(k)=min(0.75_kind_phys,0.8_kind_phys*laitbl(k))
7065 if(ifortbl(k) == 5) deltalai(k)=min(0.86_kind_phys,0.8_kind_phys*laitbl(k))
7067 if(k.ne.iswater)
then
7069 laitoday(k) = laitbl(k) - deltalai(k) * factor
7071 if(ifortbl(k) == 7)
then
7073 znttoday(k) = z0tbl(k) - 0.125_kind_phys * factor
7075 znttoday(k) = z0tbl(k)
7078 laitoday(k) = laitbl(k)
7083 IF (debug_print )
THEN
7084 print *,
'ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai,laitoday(ivgtyp),znttoday(ivgtyp)', &
7085 i,j,ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai(ivgtyp),laitoday(ivgtyp),znttoday(ivgtyp)
7094 if(.not.rdlai2d) lai = zero
7101 if(mosaic_lu == 1)
then
7103 area = area + lufrac(k)
7104 emiss = emiss+ lemitbl(k)*lufrac(k)
7105 znt = znt + lufrac(k)/alog(lb/znttoday(k))**2._kind_phys
7107 znt1 = znt1 + lufrac(k)*znttoday(k)
7108 if(.not.rdlai2d) lai = lai + laitoday(k)*lufrac(k)
7109 pc = pc + pctbl(k)*lufrac(k)
7110 msnf = msnf + mfsno(k)*lufrac(k)
7111 facsnf= facsnf + sncovfac(k)*lufrac(k)
7114 if (area.gt.one) area=one
7115 if (area <= zero)
then
7116 print *,
'Bad area of grid box', area
7118 errmsg =
'ERROR(SOILVEGIN): Bad area of grid box'
7122 IF (debug_print )
THEN
7123 print *,
'area=',area,i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),emiss,znt,znt1,lai,pc
7128 znt = lb/exp(sqrt(one/znt))
7129 if(.not.rdlai2d) lai = lai/area
7132 facsnf= facsnf /area
7134 IF (debug_print )
THEN
7135 print *,
'mosaic=',j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),emiss,znt,znt1,lai,pc
7140 emiss = lemitbl(ivgtyp)
7141 znt = znttoday(ivgtyp)
7143 msnf = mfsno(ivgtyp)
7144 facsnf= sncovfac(ivgtyp)
7145 if(.not.rdlai2d) lai = laitoday(ivgtyp)
7160 if(mosaic_soil == 1 )
then
7164 area = area + soilfrac(k)
7165 rhocs = rhocs + hc(k)*1.e6_kind_phys*soilfrac(k)
7166 bclh = bclh + bb(k)*soilfrac(k)
7167 dqm = dqm + (maxsmc(k)- &
7168 drysmc(k))*soilfrac(k)
7169 ksat = ksat + satdk(k)*soilfrac(k)
7170 psis = psis - satpsi(k)*soilfrac(k)
7171 qmin = qmin + drysmc(k)*soilfrac(k)
7172 ref = ref + refsmc(k)*soilfrac(k)
7173 wilt = wilt + wltsmc(k)*soilfrac(k)
7174 qwrtz = qwrtz + qtz(k)*soilfrac(k)
7177 if (area.gt.one) area=one
7178 if (area <= zero)
then
7181 rhocs = hc(isltyp)*1.e6_kind_phys
7183 dqm = maxsmc(isltyp)- &
7185 ksat = satdk(isltyp)
7186 psis = - satpsi(isltyp)
7187 qmin = drysmc(isltyp)
7188 ref = refsmc(isltyp)
7189 wilt = wltsmc(isltyp)
7205 if(isltyp.ne.14)
then
7206 rhocs = hc(isltyp)*1.e6_kind_phys
7208 dqm = maxsmc(isltyp)- &
7210 ksat = satdk(isltyp)
7211 psis = - satpsi(isltyp)
7212 qmin = drysmc(isltyp)
7213 ref = refsmc(isltyp)
7214 wilt = wltsmc(isltyp)
7220 END SUBROUTINE soilvegin
7227 SUBROUTINE ruclsminit( debug_print, landfrac, fice, min_seaice, &
7228 nzs, isltyp, ivgtyp, mavail, &
7229 sh2o, smfr3d, tslb, smois, &
7230 ims,ime, jms,jme, kms,kme, &
7231 its,ite, jts,jte, kts,kte )
7235#if ( WRF_CHEM == 1 )
7236 USE module_data_gocart_dust
7239 LOGICAL,
INTENT(IN ) :: debug_print
7240 real (kind_phys),
DIMENSION( ims:ime),
INTENT(IN ) :: landfrac, fice
7241 real (kind_phys),
INTENT(IN ) :: min_seaice
7243 INTEGER,
INTENT(IN ) :: &
7244 ims,ime, jms,jme, kms,kme, &
7245 its,ite, jts,jte, kts,kte, &
7248 real (kind_phys),
DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
7249 INTENT(IN) :: tslb, &
7252 INTEGER,
DIMENSION( ims:ime, jms:jme ) , &
7253 INTENT(INOUT) :: isltyp,ivgtyp
7255 real (kind_phys),
DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
7256 INTENT(OUT) :: smfr3d, &
7259 real (kind_phys),
DIMENSION( ims:ime, jms:jme ) , &
7260 INTENT(OUT) :: mavail
7263 real (kind_phys),
DIMENSION ( 1:nzs ) :: soiliqw
7265 INTEGER :: i,j,l,itf,jtf
7266 real (kind_phys) :: riw,xlmelt,tln,dqm,ref,psis,qmin,bclh
7270 riw=rhoice*1.e-3_kind_phys
7281 IF ( isltyp( i,j ) .LT. 0 )
THEN
7284 "module_sf_ruclsm.F: lsminit: out of range ISLTYP ",i,j,isltyp( i,j )
7288 IF ( errflag .EQ. 1 )
THEN
7290 "module_sf_ruclsm.F: lsminit: out of range value "// &
7291 "of ISLTYP. Is this field in the input?"
7299 if (isltyp(i,j) == 0) isltyp(i,j)=14
7301 if(landfrac(i) > zero )
then
7305 dqm = maxsmc(isltyp(i,j)) - &
7307 ref = refsmc(isltyp(i,j))
7308 psis = - satpsi(isltyp(i,j))
7309 qmin = drysmc(isltyp(i,j))
7310 bclh = bb(isltyp(i,j))
7312 mavail(i,j) = max(0.00001_kind_phys,min(one,(smois(i,1,j)-qmin)/(ref-qmin)))
7316 tln=log(tslb(i,l,j)/tfrz)
7318 if(tln.lt.zero)
then
7319 soiliqw(l)=(dqm+qmin)*(xlmelt* &
7320 (tslb(i,l,j)-tfrz)/tslb(i,l,j)/grav/psis) &
7322 soiliqw(l)=max(zero,soiliqw(l))
7323 soiliqw(l)=min(soiliqw(l),smois(i,l,j))
7324 sh2o(i,l,j)=soiliqw(l)
7325 smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/riw
7329 sh2o(i,l,j)=smois(i,l,j)
7333 elseif( fice(i) > min_seaice)
then
7360 SUBROUTINE ruclsm_soilvegparm( debug_print,MMINLURUC, MMINSL)
7364 LOGICAL,
INTENT(IN ) :: debug_print
7366 integer :: lumatch, iindex, lc, num_slope
7368 INTEGER ,
PARAMETER :: open_ok = 0
7370 character*8 :: mminluruc, mminsl
7371 character*128 :: vege_parm_string
7403 OPEN(19, file=
'VEGPARM.TBL',form=
'FORMATTED',status=
'OLD',iostat=ierr)
7404 IF(ierr .NE. open_ok )
THEN
7406 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
7410 'INPUT VEGPARM FOR ',mminluruc
7416 READ (19,
'(A)') vege_parm_string
7420 READ (19,2000,
END=2002)lutype
7421 READ (19,*)lucats,iindex
7425 'VEGPARM FOR ',lutype,
' FOUND', lucats,
' CATEGORIES'
7427 IF(lutype.NE.mminluruc)
THEN
7429 'Skipping ', lutype,
' table'
7437 READ (19,
'(A)',
END=2002) vege_parm_string
7439 IF (trim(vege_parm_string) .EQ.
"Vegetation Parameters")
THEN
7446 'Found ', lutype,
' table'
7452 IF (lumatch == 1)
then
7454 'Reading ',lutype,
' table'
7457 READ (19,*)iindex,albtbl(lc),z0tbl(lc),lemitbl(lc),pctbl(lc), &
7458 shdtbl(lc),ifortbl(lc),rstbl(lc),rgltbl(lc), &
7459 hstbl(lc),snuptbl(lc),laitbl(lc),maxalb(lc)
7465 READ (19,*)topt_data
7467 READ (19,*)cmcmax_data
7469 READ (19,*)cfactr_data
7471 READ (19,*)rsmax_data
7481 READ (19,*,iostat=ierr)urban
7483 if ( ierr /= 0 ) print *,
"-------- VEGPARM.TBL READ ERROR --------"
7484 if ( ierr /= 0 ) print *,
"Problem read URBAN from VEGPARM.TBL"
7485 if ( ierr /= 0 ) print *,
" -- Use updated version of VEGPARM.TBL "
7486 if ( ierr /= 0 ) print *,
"Problem read URBAN from VEGPARM.TBL"
7493 IF (debug_print )
THEN
7494 print *,
' LEMITBL, PCTBL, Z0TBL, LAITBL --->', lemitbl, pctbl, z0tbl, laitbl
7498 IF (lumatch == 0)
then
7505 OPEN(19, file=
'SOILPARM.TBL',form=
'FORMATTED',status=
'OLD',iostat=ierr)
7506 IF(ierr .NE. open_ok )
THEN
7508 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
7511 print *,
'INPUT SOIL TEXTURE CLASSIFICATION = ',mminsl
7516 READ (19,
'(A)') vege_parm_string
7520 READ (19,2000,
END=2003)sltype
7521 READ (19,*)slcats,iindex
7525 'SOILPARM FOR ',sltype,
' FOUND', slcats,
' CATEGORIES'
7527 IF(sltype.NE.mminsl)
THEN
7529 'Skipping ', sltype,
' table'
7537 READ (19,
'(A)',
END=2002) vege_parm_string
7539 IF (trim(vege_parm_string) .EQ.
"Soil Parameters")
THEN
7546 'Found ', sltype,
' table'
7552 IF (lumatch == 1)
then
7553 print *,
'SLCATS=',slcats
7556 READ (19,*) iindex,bb(lc),drysmc(lc),hc(lc),maxsmc(lc),&
7557 refsmc(lc),satpsi(lc),satdk(lc), satdw(lc), &
7567 IF(lumatch.EQ.0)
THEN
7568 print *,
'SOIl TEXTURE IN INPUT FILE DOES NOT '
7569 print *,
'MATCH SOILPARM TABLE'
7570 print *,
'INCONSISTENT OR MISSING SOILPARM FILE'
7576 OPEN(19, file=
'GENPARM.TBL',form=
'FORMATTED',status=
'OLD',iostat=ierr)
7577 IF(ierr .NE. open_ok )
THEN
7579 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
7585 READ (19,*) num_slope
7592 READ (19,*)slope_data(lc)
7598 READ (19,*)sbeta_data
7600 READ (19,*)fxexp_data
7602 READ (19,*)csoil_data
7604 READ (19,*)salp_data
7606 READ (19,*)refdk_data
7608 READ (19,*)refkdt_data
7610 READ (19,*)frzk_data
7612 READ (19,*)zbot_data
7614 READ (19,*)czil_data
7616 READ (19,*)smlow_data
7618 READ (19,*)smhigh_data
7623 END SUBROUTINE ruclsm_soilvegparm
7629 SUBROUTINE soilin (ISLTYP, DQM, REF, PSIS, QMIN, BCLH )
7656 integer,
parameter :: nsoilclas=19
7658 integer,
intent ( in) :: isltyp
7659 real,
intent ( out) :: dqm,ref,qmin,psis,bclh
7661 real (kind_phys) lqma(nsoilclas),lref(nsoilclas),lbcl(nsoilclas), &
7662 lpsi(nsoilclas),lqmi(nsoilclas)
7671 DATA lqma /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
7672 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
7673 0.20, 0.435, 0.468, 0.200, 0.339/
7676 DATA lref /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
7677 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
7678 0.1, 0.249, 0.454, 0.17, 0.236/
7681 DATA lqmi/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
7682 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
7683 0.004, 0.065, 0.020, 0.004, 0.008/
7686 DATA lpsi/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
7687 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
7688 0.121, 0.218, 0.468, 0.069, 0.069/
7691 DATA lbcl/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
7692 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
7693 4.05, 4.90, 11.55, 2.79, 2.79/
7696 dqm = lqma(isltyp)- &
7699 psis = - lpsi(isltyp)
7703 END SUBROUTINE soilin
7709 real (kind_phys) function
rslf(p,t)
7712 real (kind_phys),
INTENT(IN):: p, t
7713 real (kind_phys):: esl,x
7714 real (kind_phys),
PARAMETER:: c0= .611583699e03
7715 real (kind_phys),
PARAMETER:: c1= .444606896e02
7716 real (kind_phys),
PARAMETER:: c2= .143177157e01
7717 real (kind_phys),
PARAMETER:: c3= .264224321e-1
7718 real (kind_phys),
PARAMETER:: c4= .299291081e-3
7719 real (kind_phys),
PARAMETER:: c5= .203154182e-5
7720 real (kind_phys),
PARAMETER:: c6= .702620698e-8
7721 real (kind_phys),
PARAMETER:: c7= .379534310e-11
7722 real (kind_phys),
PARAMETER:: c8=-.321582393e-13
7724 x=max(-80._kind_dbl_prec,t-273.16_kind_dbl_prec)
7726 esl=c0+x*(c1+x*(c2+x*(c3+x*(c4+x*(c5+x*(c6+x*(c7+x*c8)))))))
7727 esl=min(esl, p*0.15_kind_dbl_prec)
7728 rslf=.622_kind_dbl_prec*esl/max(1.e-4_kind_dbl_prec,(p-esl))
subroutine, public ruclsminit(debug_print, landfrac, fice, min_seaice, nzs, isltyp, ivgtyp, mavail, sh2o, smfr3d, tslb, smois, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
This subroutine computes liquid and forezen soil moisture from the total soil moisture,...
real(kind_phys) function, public rslf(p, t)
This function calculates the liquid saturation vapor mixing ratio as a function of temperature and pr...
subroutine, public lsmruc(xlat, xlon, dt, init, lsm_cold_start, ktau, iter, nsl, graupelncv, snowncv, rainncv, raincv, zs, rainbl, snow, snowh, snowc, frzfrac, frpcpn, rhosnf, precipfr, exticeden, hgt, stdev, z3d, p8w, t3d, qv3d, qc3d, rho3d, emisbck, glw, gswdn, gsw, emiss, chklowq, chs, flqc, flhc, rhonewsn_ex, mosaic_lu, mosaic_soil, isncond_opt, isncovr_opt, mavail, canwat, vegfra, alb, znt, z0, snoalb, albbck, lai, landusef, nlcat, soilctop, nscat, smcwlt, smcref, qsfc, qsg, qvg, qcg, dew, soilt1, tsnav, tbot, ivgtyp, isltyp, xland, iswater, isice, xice, xice_threshold, cp, rv, rd, g0, pi, lv, stbolt, soilmois, sh2o, smavail, smmax, tso, soilt, edir, ec, ett, sublim, snoh, hfx, qfx, lh, infiltr, runoff1, runoff2, acrunoff, sfcexc, sfcevp, grdflx, snowfallac, acsnow, snom, smfr3d, keepfr3dflag, add_fire_heat_flux, fire_heat_flux, myj, shdmin, shdmax, rdlai2d, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte, errmsg, errflg)
The RUN LSM model is described in Smirnova et al.(1997) and Smirnova et al.(2000) .
This module contains the entity of the RUC LSM model, which is a soil/veg/snowpack and ice/snowpack...
This module contains the namelist options of soil/vegetation in RUC.
This module contains the entity of GFS Noah LSM Model(Version 2.7).