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
820 IF (debug_print )
THEN
821 if (abs(xlat-testptlat).lt.0.2 .and. &
822 abs(xlon-testptlon).lt.0.2)
then
823 print*,
' lat,lon=',xlat,xlon
824 print *,
' ZNT, LAI, VEGFRA, SAT, EMIS, PC --->', &
825 znt(i,j),lai(i,j),vegfra(i,j),sat,emissl(i,j),pc(i,j)
826 print *,
' ZS, ZSMAIN, ZSHALF, CONFLX, CN, SAT, --->', zs,zsmain,zshalf,conflx,cn,sat
827 print *,
'NROOT, meltfactor, iforest, ivgtyp, i,j ', nroot,meltfactor,iforest,ivgtyp(i,j),i,j
831 IF((xland(i,j)-1.5).GE.0._kind_phys)
THEN
840 acsnow(i,j)=acsnow(i,j)+precipfr(i,j)
845 patmb=p8w(i,1,j)*1.e-2_kind_phys
846 qvg(i,j) =
qsn(soilt(i,j),tbq)/patmb
847 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
849 q2sat=
qsn(tabs,tbq)/patmb
854 tso(i,k,j)= soilt(i,j)
857 IF (debug_print )
THEN
858 if (abs(xlat-testptlat).lt.0.2 .and. &
859 abs(xlon-testptlon).lt.0.2)
then
860 print*,
' water point'
861 print*,
' lat,lon=',xlat,xlon,
'SOILT=', soilt(i,j)
868 if(xice(i,j).ge.xice_threshold)
then
874 IF(seaice(i,j).GT.0.5_kind_phys)
THEN
876 IF (debug_print )
THEN
877 if (abs(xlat-testptlat).lt.0.2 .and. &
878 abs(xlon-testptlon).lt.0.2)
then
879 print*,
' sea-ice at water point'
880 print*,
' lat,lon=',xlat,xlon
889 znt(i,j) = 0.011_kind_phys
891 emissl(i,j) = emisbck(i,j)
897 patmb=p8w(i,1,j)*1.e-2_kind_phys
898 qvg(i,j) =
qsn(soilt(i,j),tbq)/patmb
900 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
903 soilmois(i,k,j) = one
906 keepfr3dflag(i,k,j) = zero
907 tso(i,k,j) = min(con_tice,tso(i,k,j))
916 soilm1d(k) = min(max(zero,soilmois(i,k,j)-qmin),dqm)
917 tso1d(k) = tso(i,k,j)
918 soiliqw(k) = min(max(zero,sh2o(i,k,j)-qmin),soilm1d(k))
919 soilice(k) =(soilm1d(k) - soiliqw(k))/0.9_kind_phys
923 smfrkeep(k) = smfr3d(i,k,j)
924 keepfr(k) = keepfr3dflag(i,k,j)
927 lmavail(i,j)=max(0.00001_kind_phys,min(one,soilm1d(1)/(ref-qmin)))
929 IF (debug_print )
THEN
930 if (abs(xlat-testptlat).lt.0.2 .and. &
931 abs(xlon-testptlon).lt.0.2)
then
932 print*,
' lat,lon=',xlat,xlon
933 print *,
'LAND, i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO', &
934 i,j,tso1d,soilm1d,patm,tabs,qvatm,qcatm,rho
935 print *,
'CONFLX =',conflx
936 print *,
'SMFRKEEP,KEEPFR ',smfrkeep,keepfr
943 smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(k))* &
944 (zshalf(k+1)-zshalf(k))
947 if (debug_print .and. abs(xlat-testptlat).lt.0.2 &
948 .and. abs(xlon-testptlon).lt.0.2)
then
949 print *,
'Old soilm1d ',i,soilm1d
952 canwatold(i,j) = canwatr
954 CALL sfctmp (debug_print, dt,ktau,conflx,i,j, &
955 xlat, xlon, testptlat, testptlon, &
957 nzs,nddzs,nroot,meltfactor, &
958 isncond_opt,isncovr_opt, &
959 iland,isoil,ivgtyp(i,j),isltyp(i,j), &
960 prcpms, newsnms,snwe,snhei,snowfrac, &
961 exticeden,rhosn,rhonewsn_ex(i),rhonewsn, &
962 rhosnfall,snowrat,grauprat,icerat,curat, &
963 patm,tabs,qvatm,qcatm,rho, &
964 glw(i,j),gswdn(i,j),gsw(i,j), &
965 emissl(i,j),emisbck(i,j), &
966 msnf(i,j), facsnf(i,j), &
967 qkms,tkms,pc(i,j),lmavail(i,j), &
968 canwatr,vegfra(i,j),alb(i,j),znt(i,j), &
969 snoalb(i,j),albbck(i,j),lai(i,j), &
970 hgt(i,j),stdev(i,j), &
971 myj,seaice(i,j),isice, &
972 add_fire_heat_flux,fire_heat_flux(i,j), &
975 rhocs,dqm,qmin,ref, &
976 wilt,psis,bclh,ksat, &
977 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
979 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
982 snweprint,snheiprint,rsm, &
983 soilm1d,tso1d,smfrkeep,keepfr, &
984 soilt(i,j),soilt1(i,j),tsnav(i,j),dew(i,j), &
985 qvg(i,j),qsg(i,j),qcg(i,j),smelt(i,j), &
986 snoh(i,j),snflx(i,j),snom(i,j),snowfallac(i,j), &
987 acsnow(i,j),edir(i,j),ec(i,j),ett(i,j),qfx(i,j), &
988 lh(i,j),hfx(i,j),
sflx(i,j),sublim(i,j), &
989 evapl(i,j),prcpl(i,j),budget(i,j),runoff1(i,j), &
990 runoff2(i,j),soilice,soiliqw,infiltrp,smf(i,j))
1000 if(mosaic_lu == 1)
then
1002 factor = max(zero,min(one,(vegfra(i,j)-shdmin(i,j))/max(one,(shdmax(i,j)-shdmin(i,j)))))
1003 if (debug_print )
then
1004 if (abs(xlat-testptlat).lt.0.1 .and. &
1005 abs(xlon-testptlon).lt.0.1)
then
1006 print *,
' lat,lon=',xlat,xlon,
' factor=',factor
1010 if((ivgtyp(i,j) == natural .or. ivgtyp(i,j) == crop) .and. factor > 0.75)
then
1015 cropsm=1.05_kind_phys*wilt - qmin
1016 cropfr = min(one,lufrac(crop) + 0.4*lufrac(natural))
1017 newsm = cropsm*cropfr + (1.-cropfr)*soilm1d(k)
1018 if(soilm1d(k) < newsm)
then
1019 IF (debug_print )
THEN
1020 if (abs(xlat-testptlat).lt.0.1 .and. &
1021 abs(xlon-testptlon).lt.0.1)
then
1022 print * ,
'Soil moisture is below wilting in cropland areas at time step',ktau
1023 print * ,
' lat,lon=',xlat,xlon
1024 print * ,
' lufrac=',lufrac,
'factor=',factor &
1025 ,
'lai,ivgtyp,lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm,', &
1026 lai(i,j),ivgtyp(i,j),lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm
1030 IF (debug_print )
THEN
1031 if (abs(xlat-testptlat).lt.0.1 .and. &
1032 abs(xlon-testptlon).lt.0.1)
then
1033 print*,
' lat,lon=',xlat,xlon
1034 print * ,
'Added soil water to cropland areas, k,soilm1d(k)',k,soilm1d(k)
1052 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* &
1053 (zshalf(k+1)-zshalf(k))
1054 smmax(i,j) =smmax(i,j)+(qmin+dqm)* &
1055 (zshalf(k+1)-zshalf(k))
1058 if (debug_print)
then
1059 if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
1060 print 100,
'(RUC runoff) i=',i,
' lat,lon=',xlat,xlon, &
1061 'RUNOFF1', runoff1(i,j),
'RUNOFF2 ',runoff2(i,j), &
1062 'edir ',edir(i,j),
'ec ',ec(i,j),
'ett ',ett(i,j)
1068 acrunoff(i,j) = (runoff1(i,j)+runoff2(i,j))*dt*rhowater
1069 smavail(i,j) = smavail(i,j) * rhowater
1070 smmax(i,j) = smmax(i,j) * rhowater
1071 smtotold(i,j) = smtotold(i,j) * rhowater
1075 soilmois(i,k,j) = soilm1d(k) + qmin
1076 sh2o(i,k,j) = min(soiliqw(k) + qmin,soilmois(i,k,j))
1077 tso(i,k,j) = tso1d(k)
1080 tso(i,nzs,j) = tbot(i,j)
1083 smfr3d(i,k,j) = smfrkeep(k)
1084 keepfr3dflag(i,k,j) = keepfr(k)
1092 patmb=p8w(i,1,j)*1.e-2_kind_phys
1093 q2sat=
qsn(tabs,tbq)/patmb
1094 qsfc(i,j) = qvg(i,j)/(one+qvg(i,j))
1098 IF((qvatm.GE.q2sat*0.95_kind_phys).AND.qvatm.LT.qvg(i,j))
THEN
1104 if(snow(i,j)==zero) emissl(i,j) = emisbck(i,j)
1105 emiss(i,j) = emissl(i,j)
1107 snow(i,j) = snwe*1000._kind_phys
1109 canwat(i,j) = canwatr*1000._kind_phys
1111 if (debug_print)
then
1112 if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
1113 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)
1116 infiltr(i,j) = infiltrp
1118 mavail(i,j) = lmavail(i,j)
1119 IF (debug_print )
THEN
1120 if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
1121 print *,
' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
1124 sfcevp(i,j) = sfcevp(i,j) + qfx(i,j) * dt
1125 grdflx(i,j) = -one *
sflx(i,j)
1135 rhosnf(i,j)=rhosnfall
1138 sfcevp(i,j) = sfcevp(i,j) + qfx(i,j) * dt
1146 if (debug_print )
then
1147 if (abs(xlat-testptlat).lt.0.2 .and. &
1148 abs(xlon-testptlon).lt.0.2)
then
1154 ac=canwat(i,j)-canwatold(i,j)*rhowater
1155 as=snwe-snowold(i,j)
1156 wb = smavail(i,j)-smtotold(i,j)
1157 waterbudget(i,j)=rainbl(i,j)+smelt(i,j)*dt*rhowater &
1159 -runoff1(i,j)*dt*rhowater-runoff2(i,j)*dt*rhowater &
1162 print *,
'soilm1d ',i,soilm1d
1163 print 100,
'(RUC budgets) i=',i,
' lat,lon=',xlat,xlon, &
1164 'budget ',budget(i,j),
'waterbudget',waterbudget(i,j), &
1165 'rainbl ',rainbl(i,j),
'runoff1 ',runoff1(i,j), &
1166 'smelt ',smelt(i,j)*dt*1.e3_kind_phys,
'smc change ',wb, &
1167 'snwe change ',as,
'canw change ',ac,
'runoff2 ',runoff2(i,j), &
1168 'qfx*dt ',qfx(i,j)*dt,
'smavail ',smavail(i,j),
'smcold',smtotold(i,j)
1170 print *,
'Smf=',smf(i,j),i,j
1171 print *,
'SNOW,SNOWold',i,j,snwe,snowold(i,j)
1172 print *,
'SNOW-SNOWold',i,j,max(zero,snwe-snowold(i,j))
1173 print *,
'CANWATold, canwat ',i,j,canwatold(i,j),canwat(i,j)
1174 print *,
'canwat(i,j)-canwatold(i,j)',max(zero,canwat(i,j)-canwatold(i,j))
1178 100
format (
";;; ",a,i4,a,2f14.7/(4(a10,
'='es14.7)))
1181 IF (debug_print )
THEN
1182 if (abs(xlat-testptlat).lt.0.2 .and. &
1183 abs(xlon-testptlon).lt.0.2)
then
1184 print *,
'LAND, i,tso1d,soilm1d,soilt - end of time step', &
1185 i,tso1d,soilm1d,soilt(i,j)
1186 print *,
'LAND, QFX, HFX after SFCTMP', i,lh(i,j),hfx(i,j)
1211 SUBROUTINE sfctmp (debug_print, delt,ktau,conflx,i,j, & !--- input variables
1212 xlat,xlon,testptlat,testptlon, &
1213 nzs,nddzs,nroot,meltfactor, &
1214 isncond_opt,isncovr_opt, &
1215 ILAND,ISOIL,IVGTYP,ISLTYP,PRCPMS, &
1216 NEWSNMS,SNWE,SNHEI,SNOWFRAC, &
1217 exticeden,RHOSN,RHONEWSN_ex,RHONEWSN,RHOSNFALL, &
1218 snowrat,grauprat,icerat,curat, &
1219 PATM,TABS,QVATM,QCATM,rho, &
1220 GLW,GSWdn,GSW,EMISS,EMISBCK,msnf,facsnf, &
1221 QKMS,TKMS,PC,MAVAIL,CST,VEGFRA,ALB,ZNT, &
1222 ALB_SNOW,ALB_SNOW_FREE,lai,hgt,stdev, &
1224 add_fire_heat_flux,fire_heat_flux, &
1225 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & !--- soil fixed fields
1226 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1227 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, & !--- constants
1229 snweprint,snheiprint,rsm, & !---output variables
1230 soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
1231 tsnav,dew,qvg,qsg,qcg, &
1232 SMELT,SNOH,SNFLX,SNOM,SNOWFALLAC,ACSNOW, &
1233 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
1234 evapl,prcpl,fltot,runoff1,runoff2,soilice, &
1235 soiliqw,infiltr,smf)
1242 INTEGER,
INTENT(IN ) :: isice,i,j,nroot,ktau,nzs , &
1244 integer,
intent(in ) :: isncond_opt,isncovr_opt
1246 real (kind_phys),
INTENT(IN ) :: DELT,CONFLX,meltfactor,xlat,xlon
1247 real (kind_phys),
INTENT(IN ) :: testptlat,testptlon
1248 real (kind_phys),
INTENT(IN ) :: C1SN,C2SN,RHONEWSN_ex
1249 LOGICAL,
INTENT(IN ) :: myj, debug_print, exticeden
1251 real (kind_phys) , &
1252 INTENT(IN ) :: PATM, &
1256 real (kind_phys) , &
1257 INTENT(IN ) :: GLW, &
1271 LOGICAL,
INTENT(IN ) :: add_fire_heat_flux
1273 INTEGER,
INTENT(IN ) :: IVGTYP, ISLTYP
1275 real (kind_phys) , &
1276 INTENT(INOUT) :: EMISS, &
1285 real (kind_phys) :: &
1297 real (kind_phys),
INTENT(IN ) :: CN, &
1308 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: ZSMAIN, &
1313 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: DTDZS
1315 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: TBQ
1320 real (kind_phys),
DIMENSION( 1:nzs ) , &
1321 INTENT(INOUT) :: TS1D, &
1324 real (kind_phys),
DIMENSION( 1:nzs ) , &
1325 INTENT(INOUT) :: KEEPFR
1327 real (kind_phys),
DIMENSION(1:NZS),
INTENT(INOUT) :: SOILICE, &
1331 INTEGER,
INTENT(INOUT) :: ILAND,ISOIL
1335 real (kind_phys) , &
1336 INTENT(INOUT) :: DEW, &
1375 real (kind_phys),
DIMENSION(1:NZS) :: &
1387 real (kind_phys) :: &
1412 real (kind_phys),
INTENT(INOUT) :: RSM, &
1419 real (kind_phys) :: BSN, XSN , &
1420 RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , &
1421 T3, UPFLUX, XINET, snowfrac2, m
1422 real (kind_phys) :: snhei_crit, snhei_crit_newsn, keep_snow_albedo, SNOWFRACnewsn
1423 real (kind_phys) :: newsnowratio, dd1
1425 real (kind_phys) :: rhonewgr,rhonewice
1427 real (kind_phys) :: RNET,GSWNEW,GSWIN,EMISSN,ZNTSN,EMISS_snowfree
1428 real (kind_phys) :: VEGFRAC, snow_mosaic, snfr, vgfr
1429 real :: cice, albice, albsn, drip, dripsn, dripliq
1430 real :: interw, intersn, infwater, intwratio
1433 integer,
parameter :: ilsnow=99
1435 IF (debug_print )
THEN
1436 print *,
' in SFCTMP',i,j,nzs,nddzs,nroot, &
1437 snwe,rhosn,snom,smelt,ts1d
1448 snhei_crit=0.01601_kind_phys*rhowater/rhosn
1449 snhei_crit_newsn=0.0005_kind_phys*rhowater/rhosn
1451 zntsn = z0tbl(isice)
1458 rhonewsn = 100._kind_phys
1459 if(snhei == zero) snowfrac=zero
1472 vegfrac=0.01_kind_phys*vegfra
1492 albice=alb_snow_free
1494 emissn = 0.99_kind_phys
1495 emiss_snowfree = emisbck
1501 if(seaice.ge.0.5_kind_phys)
then
1503 tice(k) = ts1d(k) - tfrz
1504 rhosice(k) = 917.6_kind_phys/(one-0.000165_kind_phys*tice(k))
1505 cice = 2115.85_kind_phys +7.7948_kind_phys*tice(k)
1506 capice(k) = cice*rhosice(k)
1507 thdifice(k) = 2.260872_kind_phys/capice(k)
1513 albice = min(alb_snow_free,max(alb_snow_free - 0.05_kind_phys, &
1514 alb_snow_free - 0.1_kind_phys*(tice(1)+10._kind_phys)/10._kind_phys ))
1517 IF (debug_print )
THEN
1518 print *,
'alb_snow_free',alb_snow_free
1519 print *,
'GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE',&
1520 gsw,gswnew,glw,soilt,emiss,alb,albice,snwe
1523 if(snhei.gt.0.0081_kind_phys*rhowater/rhosn)
then
1525 bsn=delt/3600._kind_phys*c1sn*exp(0.08_kind_phys*min(zero,tsnav)-c2sn*rhosn*1.e-3_kind_phys)
1526 if(bsn*snwe*100._kind_phys.lt.1.e-4_kind_phys)
goto 777
1527 xsn=rhosn*(exp(bsn*snwe*100._kind_phys)-one)/(bsn*snwe*100._kind_phys)
1528 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
1533 if(snowfrac < 0.75_kind_phys) snow_mosaic = one
1537 acsnow=acsnow+newsn*rhowater
1539 IF(newsn.GT.zero)
THEN
1541 IF (debug_print )
THEN
1542 print *,
'THERE IS NEW SNOW, newsn', newsn
1545 newsnowratio = min(one,newsn/(snwe+newsn))
1558 rhonewsn = rhonewsn_ex
1560 rhonewsn=min(125._kind_phys,rhowater/max(8._kind_phys,(17._kind_phys*tanh((276.65_kind_phys-tabs)*0.15_kind_phys))))
1561 rhonewgr=min(500._kind_phys,rhowater/max(2._kind_phys,(3.5_kind_phys*tanh((274.15_kind_phys-tabs)*0.3333_kind_phys))))
1567 rhosnfall = min(500._kind_phys,max(58.8_kind_phys,(rhonewsn*snowrat + &
1568 rhonewgr*grauprat + rhonewice*icerat + rhonewgr*curat)))
1570 if (debug_print)
then
1572 print *,
' xlat, xlon', xlat, xlon
1573 print *,
'snow_mosaic = ',snow_mosaic
1574 print *,
'new snow,newsnowratio,rhosnfall =',newsn,newsnowratio,rhosnfall
1575 print *,
'snowrat,grauprat,icerat,curat,rhonewgr,rhonewsn,rhonewice',snowrat,grauprat,icerat,curat,rhonewgr,rhonewsn,rhonewice
1584 xsn=(rhosn*snwe+rhonewsn*newsn)/ &
1586 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
1589 IF(prcpms > zero)
THEN
1602 if(vegfrac > 0.01_kind_phys)
then
1605 interw=0.25_kind_phys*delt*prcpms*(one-exp(-0.5_kind_phys*lai))*vegfrac
1606 intersn=0.25_kind_phys*newsn*(one-exp(-0.5_kind_phys*lai))*vegfrac
1607 infwater=prcpms - interw/delt
1608 if((interw+intersn) > zero)
then
1609 intwratio=interw/(interw+intersn)
1613 dd1=cst + interw + intersn
1628 IF(newsn.GT.zero)
THEN
1630 snwe=max(zero,snwe+newsn-intersn)
1632 if(drip > zero)
then
1633 if (snow_mosaic==one)
then
1634 dripliq=drip*intwratio
1635 dripsn = drip - dripliq
1637 infwater=infwater+dripliq
1644 snhei=snwe*rhowater/rhosn
1645 newsn=newsn*rhowater/rhonewsn
1648 IF(snhei.GT.zero)
THEN
1665 if(isncovr_opt == 1)
then
1666 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
1667 elseif(isncovr_opt == 2)
then
1668 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
1669 if(ivgtyp == glacier .or. ivgtyp == bare)
then
1671 snowfrac2 = tanh( snhei/(2.5_kind_phys * 0.2_kind_phys *(rhosn/rhonewsn)**m))
1676 snowfrac2 = tanh( snhei/(2.5_kind_phys *min(0.2_kind_phys,znt) *(rhosn/rhonewsn)**m))
1679 snowfrac = 0.5_kind_phys*(snowfrac+snowfrac2)
1686 snowfrac = tanh( snhei/(10._kind_phys * facsnf *(rhosn/rhonewsn)**m))
1689 if(newsn > zero )
then
1690 snowfracnewsn=min(one,snowfallac*1.e-3_kind_phys/snhei_crit_newsn)
1695 if(hgt > 2500._kind_phys .and. ivgtyp == glacier) snowfrac=min(0.85_kind_phys,snowfrac)
1698 if(ivgtyp == urban) snowfrac=min(0.75_kind_phys,snowfrac)
1700 if(snowfrac < 0.75_kind_phys) snow_mosaic = one
1702 keep_snow_albedo = zero
1703 IF (snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys)
THEN
1705 keep_snow_albedo = one
1710 IF (debug_print )
THEN
1711 print *,
'SNHEI_CRIT,SNOWFRAC,SNHEI_CRIT_newsn,SNOWFRACnewsn', &
1712 snhei_crit,snowfrac,snhei_crit_newsn,snowfracnewsn
1717 IF(newsn.eq.zero .and. znt.le.0.2_kind_phys .and. ivgtyp.ne.isice)
then
1718 if( snhei .le. 2._kind_phys*znt)
then
1720 znt=0.55_kind_phys*znt+0.45_kind_phys*z0tbl(iland)
1721 elseif( snhei .gt. 2._kind_phys*znt .and. snhei .le. 4._kind_phys*znt)
then
1722 znt=0.2_kind_phys*znt+0.8_kind_phys*z0tbl(iland)
1723 elseif(snhei > 4._kind_phys*znt)
then
1730 IF(seaice .LT. 0.5_kind_phys)
THEN
1736 if( snow_mosaic == one)
then
1738 if(keep_snow_albedo > 0.9_kind_phys .and. albsn < 0.4_kind_phys)
then
1748 albsn = max(keep_snow_albedo*alb_snow, &
1749 min((alb_snow_free + &
1750 (alb_snow - alb_snow_free) * snowfrac), alb_snow))
1751 if(newsn > zero .and. keep_snow_albedo > 0.9_kind_phys .and. albsn < 0.4_kind_phys)
then
1760 emiss = max(keep_snow_albedo*emissn, &
1761 min((emiss_snowfree + &
1762 (emissn - emiss_snowfree) * snowfrac), emissn))
1765 IF (debug_print )
THEN
1767 print *,
'Snow on soil ALBsn,emiss,snow_mosaic',i,j,albsn,emiss,snow_mosaic
1778 if(albsn.lt.0.4_kind_phys .or. keep_snow_albedo==1)
then
1782 alb = min(albsn,max(albsn - 0.1_kind_phys*(soilt - 263.15_kind_phys)/ &
1783 (tfrz-263.15_kind_phys)*albsn, albsn - 0.05_kind_phys))
1787 if( snow_mosaic == one)
then
1791 albsn = max(keep_snow_albedo*alb_snow, &
1792 min((albice + (alb_snow - albice) * snowfrac), alb_snow))
1793 emiss = max(keep_snow_albedo*emissn, &
1795 min((emiss_snowfree + &
1796 (emissn - emiss_snowfree) * snowfrac), emissn))
1799 IF (debug_print )
THEN
1800 print *,
'Snow on ice snow_mosaic,ALBsn,emiss',i,j,albsn,emiss,snow_mosaic
1805 if(albsn.lt.alb_snow .or. keep_snow_albedo .eq.one)
then
1809 alb = min(albsn,max(albsn - 0.15_kind_phys*albsn*(soilt - 263.15_kind_phys)/ &
1810 (tfrz-263.15_kind_phys), albsn - 0.1_kind_phys))
1815 if (snow_mosaic==one)
then
1818 if(seaice .LT. 0.5_kind_phys)
then
1823 gswnew=gswin*(one-alb_snow_free)
1825 t3 = stbolt*soilt*soilt*soilt
1827 xinet = emiss_snowfree*(glw-upflux)
1828 rnet = gswnew + xinet
1829 IF ( add_fire_heat_flux .and. fire_heat_flux >0 )
then
1830 IF (debug_print )
THEN
1831 print *,
'RNET snow-free, fire_heat_flux, xlat/xlon',rnet, fire_heat_flux,xlat,xlon
1833 rnet = rnet + fire_heat_flux
1836 IF (debug_print )
THEN
1837 print *,
'Fractional snow - snowfrac=',snowfrac
1838 print *,
'Snowfrac<1 GSWin,GSWnew -',gswin,gswnew,
'SOILT, RNET',soilt,rnet
1841 soilm1ds(k) = soilm1d(k)
1843 smfrkeeps(k) = smfrkeep(k)
1844 keepfrs(k) = keepfr(k)
1845 soilices(k) = soilice(k)
1846 soiliqws(k) = soiliqw(k)
1860 CALL soil(debug_print,xlat, xlon, testptlat, testptlon,&
1862 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1863 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew,gswin, &
1864 emiss_snowfree,rnet,qkms,tkms,pc,csts,dripliq, &
1865 infwater,rho,vegfrac,lai,myj, &
1867 qwrtz,rhocs,dqm,qmin,ref,wilt, &
1868 psis,bclh,ksat,sat,cn, &
1869 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1871 lv,cp,rovcp,g0,cw,stbolt,tabs, &
1874 soilm1ds,ts1ds,smfrkeeps,keepfrs, &
1875 dews,soilts,qvgs,qsgs,qcgs,edir1s,ec1s, &
1876 ett1s,eetas,qfxs,hfxs,ss,evapls,prcpls,fltots,runoff1s, &
1877 runoff2s,mavails,soilices,soiliqws, &
1884 gswnew=gswin*(one-albice)
1886 t3 = stbolt*soilt*soilt*soilt
1888 xinet = emiss_snowfree*(glw-upflux)
1889 rnet = gswnew + xinet
1890 IF (debug_print )
THEN
1891 print *,
'Fractional snow - snowfrac=',snowfrac
1892 print *,
'Snowfrac<1 GSWin,GSWnew -',gswin,gswnew,
'SOILT, RNET',soilt,rnet
1905 CALL sice(debug_print,xlat,xlon, &
1907 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1908 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew, &
1909 0.98_kind_phys,rnet,qkms,tkms,rho,myj, &
1911 tice,rhosice,capice,thdifice, &
1912 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1914 lv,cp,rovcp,cw,stbolt,tabs, &
1916 ts1ds,dews,soilts,qvgs,qsgs,qcgs, &
1917 eetas,qfxs,hfxs,ss,evapls,prcpls,fltots &
1919 edir1 = eeta*1.e-3_kind_phys
1940 gswnew=gswin*(one-alb)
1942 t3 = stbolt*soilt*soilt*soilt
1944 xinet = emiss*(glw-upflux)
1945 rnet = gswnew + xinet
1946 IF (debug_print )
THEN
1948 print *,
'RNET=',rnet
1949 print *,
'SNOW - I,J,newsn,snwe,snhei,GSW,GSWnew,GLW,UPFLUX,ALB',&
1950 i,j,newsn,snwe,snhei,gsw,gswnew,glw,upflux,alb
1951 print *,
'GSWnew',gswnew,
'alb=',alb
1954 if (seaice .LT. 0.5_kind_phys)
then
1956 IF ( add_fire_heat_flux .and. fire_heat_flux>0 )
then
1957 IF (debug_print )
THEN
1958 print *,
'RNET snow, fire_heat_flux, xlat/xlon',rnet, fire_heat_flux,xlat,xlon
1960 rnet = rnet + fire_heat_flux
1962 if(snow_mosaic==one)
then
1967 CALL snowsoil (debug_print,xlat,xlon,testptlat,testptlon, &
1968 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1969 isncond_opt,isncovr_opt, &
1970 meltfactor,rhonewsn,snhei_crit, &
1971 iland,prcpms,rainf,newsn,snhei,snwe,snfr, &
1972 rhosn,patm,qvatm,qcatm, &
1973 glw,gswnew,gswin,emiss,rnet,ivgtyp, &
1974 qkms,tkms,pc,cst,dripsn,infwater, &
1975 rho,vegfrac,alb,znt,lai, &
1978 qwrtz,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1979 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1981 lv,cp,rovcp,g0,cw,stbolt,tabs, &
1984 ilnb,snweprint,snheiprint,rsm, &
1985 soilm1d,ts1d,smfrkeep,keepfr, &
1986 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1987 smelt,snoh,snflx,snom,edir1,ec1,ett1,eeta, &
1988 qfx,hfx,s,sublim,prcpl,fltot,runoff1,runoff2, &
1989 mavail,soilice,soiliqw,infiltr )
1992 if(snow_mosaic==one)
then
1999 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
2000 isncond_opt,isncovr_opt, &
2001 meltfactor,rhonewsn,snhei_crit, &
2002 iland,prcpms,rainf,newsn,snhei,snwe,snfr, &
2003 rhosn,patm,qvatm,qcatm, &
2004 glw,gswnew,emiss,rnet, &
2005 qkms,tkms,rho,myj, &
2008 tice,rhosice,capice,thdifice, &
2009 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2011 lv,cp,rovcp,cw,stbolt,tabs, &
2013 ilnb,snweprint,snheiprint,rsm,ts1d, &
2014 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
2015 smelt,snoh,snflx,snom,eeta, &
2016 qfx,hfx,s,sublim,prcpl,fltot &
2018 edir1 = eeta*1.e-3_kind_phys
2036 if (snow_mosaic==one)
then
2039 if(seaice .LT. 0.5_kind_phys)
then
2041 IF (debug_print )
THEN
2043 print *,
' xlat, xlon', xlat, xlon
2044 print *,
' snowfrac = ',snowfrac
2045 print *,
' SOILT snow on land', ktau, i,j,soilt
2046 print *,
' SOILT on snow-free land', i,j,soilts
2047 print *,
' ts1d,ts1ds',i,j,ts1d,ts1ds
2048 print *,
' SNOW flux',i,j, snflx
2049 print *,
' Ground flux on snow-covered land',i,j, s
2050 print *,
' Ground flux on snow-free land', i,j,ss
2051 print *,
' CSTS, CST', i,j,csts,cst
2055 soilm1d(k) = soilm1ds(k)*(1.-snowfrac) + soilm1d(k)*snowfrac
2056 ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac
2057 smfrkeep(k) = smfrkeeps(k)*(1.-snowfrac) + smfrkeep(k)*snowfrac
2058 if(snowfrac > 0.5_kind_phys)
then
2059 keepfr(k) = keepfr(k)
2061 keepfr(k) = keepfrs(k)
2063 soilice(k) = soilices(k)*(1.-snowfrac) + soilice(k)*snowfrac
2064 soiliqw(k) = soiliqws(k)*(1.-snowfrac) + soiliqw(k)*snowfrac
2066 dew = dews*(one-snowfrac) + dew*snowfrac
2067 soilt = soilts*(one-snowfrac) + soilt*snowfrac
2068 qvg = qvgs*(one-snowfrac) + qvg*snowfrac
2069 qsg = qsgs*(one-snowfrac) + qsg*snowfrac
2070 qcg = qcgs*(one-snowfrac) + qcg*snowfrac
2071 edir1 = edir1s*(one-snowfrac) + edir1*snowfrac
2072 ec1 = ec1s*(one-snowfrac) + ec1*snowfrac
2073 cst = csts*(one-snowfrac) + cst*snowfrac
2074 ett1 = ett1s*(one-snowfrac) + ett1*snowfrac
2075 eeta = eetas*(one-snowfrac) + eeta*snowfrac
2076 qfx = qfxs*(one-snowfrac) + qfx*snowfrac
2077 hfx = hfxs*(one-snowfrac) + hfx*snowfrac
2078 s = ss*(one-snowfrac) + s*snowfrac
2079 evapl = evapls*(one-snowfrac)
2080 prcpl = prcpls*(one-snowfrac) + prcpl*snowfrac
2081 fltot = fltots*(one-snowfrac) + fltot*snowfrac
2082 alb = max(keep_snow_albedo*alb, &
2083 min((alb_snow_free + (alb - alb_snow_free) * snowfrac), alb))
2085 emiss = max(keep_snow_albedo*emissn, &
2086 min((emiss_snowfree + &
2087 (emissn - emiss_snowfree) * snowfrac), emissn))
2089 runoff1 = runoff1s*(one-snowfrac) + runoff1*snowfrac
2090 runoff2 = runoff2s*(one-snowfrac) + runoff2*snowfrac
2091 mavail = mavails*(one-snowfrac) + one*snowfrac
2092 infiltr = infiltrs*(one-snowfrac) + infiltr*snowfrac
2094 IF (debug_print )
THEN
2096 print *,
' Ground flux combined', xlat,xlon, s
2097 print *,
' SOILT combined on land', soilt
2098 print *,
' TS combined on land', ts1d
2103 IF (debug_print )
THEN
2104 print *,
'SOILT snow on ice', soilt
2107 ts1d(k) = ts1ds(k)*(one-snowfrac) + ts1d(k)*snowfrac
2109 dew = dews*(one-snowfrac) + dew*snowfrac
2110 soilt = soilts*(one-snowfrac) + soilt*snowfrac
2111 qvg = qvgs*(one-snowfrac) + qvg*snowfrac
2112 qsg = qsgs*(one-snowfrac) + qsg*snowfrac
2113 qcg = qcgs*(one-snowfrac) + qcg*snowfrac
2115 eeta = eetas*(one-snowfrac) + eeta*snowfrac
2116 qfx = qfxs*(one-snowfrac) + qfx*snowfrac
2117 hfx = hfxs*(one-snowfrac) + hfx*snowfrac
2118 s = ss*(one-snowfrac) + s*snowfrac
2119 prcpl = prcpls*(one-snowfrac) + prcpl*snowfrac
2120 fltot = fltots*(one-snowfrac) + fltot*snowfrac
2121 alb = max(keep_snow_albedo*alb, &
2122 min((albice + (alb - alb_snow_free) * snowfrac), alb))
2123 emiss = max(keep_snow_albedo*emissn, &
2124 min((emiss_snowfree + &
2125 (emissn - emiss_snowfree) * snowfrac), emissn))
2126 runoff1 = runoff1s*(one-snowfrac) + runoff1*snowfrac
2127 runoff2 = runoff2s*(one-snowfrac) + runoff2*snowfrac
2128 IF (debug_print )
THEN
2129 print *,
'SOILT combined on ice', soilt
2165 IF(snhei == zero)
then
2170 emiss = emiss_snowfree
2174 if(isncovr_opt == 1)
then
2175 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
2176 elseif(isncovr_opt == 2)
then
2178 snowfrac=min(one,snhei/(2._kind_phys*snhei_crit))
2179 if(ivgtyp == glacier .or. ivgtyp == bare)
then
2181 snowfrac2 = tanh( snhei/(2.5_kind_phys * 0.2_kind_phys *(rhosn/rhonewsn)**m))
2186 snowfrac2 = tanh( snhei/(2.5 *min(0.2,znt) *(rhosn/rhonewsn)**m))
2189 snowfrac = 0.5_kind_phys*(snowfrac+snowfrac2)
2196 snowfrac = tanh( snhei/(10._kind_phys * facsnf *(rhosn/rhonewsn)**m))
2201 if(hgt > 2500._kind_phys .and. ivgtyp == glacier) snowfrac=min(0.85_kind_phys,snowfrac)
2203 if(ivgtyp == urban) snowfrac=min(0.75_kind_phys,snowfrac)
2207 IF (debug_print )
then
2209 print *,
'Snowfallac xlat, xlon',xlat,xlon
2210 print *,
'newsn [m],rhonewsn,newsnowratio=',newsn,rhonewsn,newsnowratio
2211 print *,
'Time-step newsn depth [m], swe [m]',newsn,newsn*rhonewsn
2212 print *,
'Time-step smelt: swe [m]' ,smelt*delt
2213 print *,
'Time-step sublim: swe,[kg m-2]',sublim*delt
2216 snowfallac = snowfallac + newsn * 1.e3_kind_phys
2218 IF (debug_print )
THEN
2220 print *,
'snowfallac,snhei,snwe',snowfallac,snhei,snwe
2231 t3 = stbolt*soilt*soilt*soilt
2233 xinet = emiss*(glw-upflux)
2234 rnet = gswnew + xinet
2235 IF (debug_print )
THEN
2236 print *,
'NO snow on the ground GSWnew -',gswnew,
'RNET=',rnet
2239 if(seaice .LT. 0.5_kind_phys)
then
2241 IF ( add_fire_heat_flux .and. fire_heat_flux>0)
then
2242 IF (debug_print )
THEN
2243 print *,
'RNET no snow, fire_heat_flux, xlat/xlon',rnet, fire_heat_flux,xlat,xlon
2245 rnet = rnet + fire_heat_flux
2248 CALL soil(debug_print,xlat, xlon, testptlat, testptlon,&
2250 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2251 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew,gswin, &
2252 emiss,rnet,qkms,tkms,pc,cst,drip,infwater, &
2253 rho,vegfrac,lai,myj, &
2255 qwrtz,rhocs,dqm,qmin,ref,wilt, &
2256 psis,bclh,ksat,sat,cn, &
2257 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2259 lv,cp,rovcp,g0,cw,stbolt,tabs, &
2262 soilm1d,ts1d,smfrkeep,keepfr, &
2263 dew,soilt,qvg,qsg,qcg,edir1,ec1, &
2264 ett1,eeta,qfx,hfx,s,evapl,prcpl,fltot,runoff1, &
2265 runoff2,mavail,soilice,soiliqw, &
2271 if(alb.ne.albice) gswnew=gsw/(one-alb)*(one-albice)
2273 rnet = gswnew + xinet
2275 CALL sice(debug_print,xlat,xlon, &
2277 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2278 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew, &
2279 emiss,rnet,qkms,tkms,rho,myj, &
2281 tice,rhosice,capice,thdifice, &
2282 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2284 lv,cp,rovcp,cw,stbolt,tabs, &
2286 ts1d,dew,soilt,qvg,qsg,qcg, &
2287 eeta,qfx,hfx,s,evapl,prcpl,fltot &
2289 edir1 = eeta*1.e-3_kind_phys
2341 SUBROUTINE soil (debug_print,xlat,xlon,testptlat,testptlon,&
2342 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,& !--- input variables
2343 PRCPMS,RAINF,PATM,QVATM,QCATM, &
2344 GLW,GSW,GSWin,EMISS,RNET, &
2345 QKMS,TKMS,PC,cst,drip,infwater,rho,vegfrac,lai, &
2347 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & !--- soil fixed fields
2348 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2349 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, & !--- constants
2351 soilmois,tso,smfrkeep,keepfr, & !--- output variables
2352 dew,soilt,qvg,qsg,qcg, &
2353 edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
2354 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
2355 soiliqw,infiltrp,smf)
2419 LOGICAL,
INTENT(IN ) :: debug_print
2420 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
2422 INTEGER,
INTENT(IN ) :: i,j,iland,isoil
2423 real (kind_phys),
INTENT(IN ) :: DELT,CONFLX
2424 real (kind_phys),
INTENT(IN ) :: xlat,xlon,testptlat,testptlon
2425 LOGICAL,
INTENT(IN ) :: myj
2428 INTENT(IN ) :: PATM, &
2433 INTENT(IN ) :: GLW, &
2447 INTENT(IN ) :: RHOCS, &
2457 real (kind_phys),
INTENT(IN ) :: CN, &
2466 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: ZSMAIN, &
2470 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: DTDZS
2472 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: TBQ
2477 real (kind_phys),
DIMENSION( 1:nzs ) , &
2478 INTENT(INOUT) :: TSO, &
2482 real (kind_phys),
DIMENSION( 1:nzs ) , &
2483 INTENT(INOUT) :: KEEPFR
2487 INTENT(INOUT) :: DEW, &
2510 real (kind_phys),
DIMENSION(1:NZS),
INTENT(OUT) :: SOILICE, &
2515 real (kind_phys) :: INFILTRP, transum , &
2517 TABS, T3, UPFLUX, XINET
2518 real (kind_phys) :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , &
2519 can,epot,fac,fltot,ft,fq,hft , &
2521 trans,zn,ci,cvw,tln,tavln,pi , &
2522 DD1,CMC2MS,DRYCAN,WETCAN , &
2524 real (kind_phys),
DIMENSION(1:NZS) :: transp,cap,diffu,hydro, &
2525 thdif,tranf,tav,soilmoism , &
2526 soilicem,soiliqwm,detal , &
2527 fwsat,lwsat,told,smold
2529 real (kind_phys) :: soiltold,smf
2530 real (kind_phys) :: soilres, alfa, fex, fex_fc, fc, psit
2532 INTEGER :: nzs1,nzs2,k
2572 dzstop=one/(zsmain(2)-zsmain(1))
2573 ras=rho*1.e-3_kind_phys
2574 riw=rhoice*1.e-3_kind_phys
2580 tln=log(tso(k)/tfrz)
2581 if(tln.lt.zero)
then
2582 soiliqw(k)=(dqm+qmin)*(xlmelt* &
2583 (tso(k)-tfrz)/tso(k)/grav/psis) &
2585 soiliqw(k)=max(zero,soiliqw(k))
2586 soiliqw(k)=min(soiliqw(k),soilmois(k))
2587 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2590 if(keepfr(k).eq.one)
then
2591 soilice(k)=min(soilice(k),smfrkeep(k))
2592 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
2597 soiliqw(k)=soilmois(k)
2604 tav(k)=0.5_kind_phys*(tso(k)+tso(k+1))
2605 soilmoism(k)=0.5_kind_phys*(soilmois(k)+soilmois(k+1))
2606 tavln=log(tav(k)/tfrz)
2608 if(tavln.lt.zero)
then
2609 soiliqwm(k)=(dqm+qmin)*(xlmelt* &
2610 (tav(k)-tfrz)/tav(k)/grav/psis) &
2612 fwsat(k)=dqm-soiliqwm(k)
2613 lwsat(k)=soiliqwm(k)+qmin
2614 soiliqwm(k)=max(zero,soiliqwm(k))
2615 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
2616 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
2618 if(keepfr(k).eq.one)
then
2619 soilicem(k)=min(soilicem(k), &
2620 0.5_kind_phys*(smfrkeep(k)+smfrkeep(k+1)))
2621 soiliqwm(k)=max(zero,soilmoism(k)-soilicem(k)*riw)
2622 fwsat(k)=dqm-soiliqwm(k)
2623 lwsat(k)=soiliqwm(k)+qmin
2628 soiliqwm(k)=soilmoism(k)
2636 if(soilice(k).gt.zero)
then
2637 smfrkeep(k)=soilice(k)
2639 smfrkeep(k)=soilmois(k)/riw
2648 xlat, xlon, testptlat, testptlon, &
2650 nzs,fwsat,lwsat,tav,keepfr, &
2651 soilmois,soiliqw,soilice, &
2652 soilmoism,soiliqwm,soilicem, &
2654 qwrtz,rhocs,dqm,qmin,psis,bclh,ksat, &
2656 riw,xlmelt,cp,g0_p,cvw,ci, &
2659 thdif,diffu,hydro,cap)
2666 q1=-qkms*ras*(qvatm - qsg)
2669 IF(qvatm.GE.qsg)
THEN
2677 wetcan=min(0.25_kind_phys,max(zero,(cst/sat))**cn)
2683 CALL transf(debug_print, &
2684 xlat, xlon, testptlat, testptlon, &
2686 nzs,nroot,soiliqw,tabs,lai,gswin, &
2688 dqm,qmin,ref,wilt,zshalf,pc,iland, &
2695 smold(k)=soilmois(k)
2702 fex=min(one,soilmois(1)/dqm)
2703 fex=max(fex,0.01_kind_phys)
2704 psit=psis*fex ** (-bclh)
2705 psit = max(-1.e5_kind_phys, psit)
2706 alfa=min(one,exp(g0_p*psit/r_v/soilt))
2725 if((soilmois(1)+qmin) > fc .or. (qvatm-qvg) > zero)
then
2728 fex_fc=min(one,(soilmois(1)+qmin)/fc)
2729 fex_fc=max(fex_fc,0.01_kind_phys)
2730 soilres=0.25_kind_phys*(one-cos(piconst*fex_fc))**2._kind_phys
2732 IF ( debug_print )
THEN
2733 print *,
'piconst=',piconst
2734 print *,
'fex,psit,psis,bclh,g0_p,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc', &
2735 fex,psit,psis,bclh,g0_p,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc
2742 CALL soiltemp(debug_print,xlat,xlon,testptlat,testptlon,&
2745 delt,ktau,conflx,nzs,nddzs,nroot, &
2747 patm,tabs,qvatm,qcatm,emiss,rnet, &
2748 qkms,tkms,pc,rho,vegfrac, lai, &
2749 thdif,cap,drycan,wetcan, &
2750 transum,dew,mavail,soilres,alfa, &
2752 dqm,qmin,bclh,zsmain,zshalf,dtdzs,tbq, &
2754 xlv,cp,g0_p,cvw,stbolt, &
2756 tso,soilt,qvg,qsg,qcg,x)
2764 IF(qvatm.GE.qsg)
THEN
2765 dew=qkms*(qvatm-qsg)
2773 transp(k)=vegfrac*ras*qkms* &
2775 tranf(k)*drycan/zshalf(nroot+1)
2776 IF(transp(k).GT.zero) transp(k)=zero
2787 tln=log(tso(k)/tfrz)
2788 if(tln.lt.zero)
then
2789 soiliqw(k)=(dqm+qmin)*(xlmelt* &
2790 (tso(k)-tfrz)/tso(k)/grav/psis) &
2792 soiliqw(k)=max(zero,soiliqw(k))
2793 soiliqw(k)=min(soiliqw(k),soilmois(k))
2794 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2796 if(keepfr(k).eq.one)
then
2797 soilice(k)=min(soilice(k),smfrkeep(k))
2798 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
2803 soiliqw(k)=soilmois(k)
2812 xlat, xlon, testptlat, testptlon, &
2814 delt,nzs,nddzs,dtdzs,dtdzs2,riw, &
2815 zsmain,zshalf,diffu,hydro, &
2816 qsg,qvg,qcg,qcatm,qvatm,-infwater, &
2817 qkms,transp,drip,dew,zero,soilice,vegfrac, &
2820 dqm,qmin,ref,ksat,ras,infmax, &
2822 soilmois,soiliqw,mavail,runoff1, &
2835 if (soilice(k).gt.zero)
then
2836 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k))
then
2846 t3 = stbolt*soiltold*soiltold*soiltold
2847 upflux = t3 * 0.5_kind_phys*(soiltold+soilt)
2848 xinet = emiss*(glw-upflux)
2849 hft=-tkms*cp*rho*(tabs-soilt)
2850 hfx=-tkms*cp*rho*(tabs-soilt) &
2851 *(p1000mb*0.00001_kind_phys/patm)**rovcp
2852 q1=-qkms*ras*(qvatm - qsg)
2855 IF (q1.LE.zero)
THEN
2862 eeta=-qkms*ras*(qvatm/(one+qvatm) - qsg/(one+qsg))*rhowater
2863 cst= cst-eeta*delt*vegfrac
2864 IF (debug_print )
THEN
2866 print *,
'Cond MYJ EETA',eeta,eeta*xlv, i,j
2871 cst=cst+delt*dew*ras * vegfrac
2872 IF (debug_print )
THEN
2874 print *,
'Cond RUC LSM EETA',eeta,eeta*xlv, i,j
2881 edir1 =-soilres*(one-vegfrac)*qkms*ras* &
2884 ec1 = q1 * wetcan * vegfrac
2885 IF (debug_print )
THEN
2886 IF(i.eq.440.and.j.eq.180.or. qfx.gt.1000..or.i.eq.417.and.j.eq.540)
then
2887 print *,
'CST before update=',cst
2888 print *,
'EC1=',ec1,
'CMC2MS=',cmc2ms
2892 cst=max(zero,cst-ec1 * delt)
2896 eeta=-soilres*qkms*ras*(qvatm/(one+qvatm) - qvg/(one+qvg))*rhowater
2898 IF (debug_print )
THEN
2900 print *,
'QKMS,RAS,QVATM/(one+QVATM),QVG/(one+QVG),QSG ', &
2901 qkms,ras,qvatm/(one+qvatm),qvg/(one+qvg),qsg
2902 print *,
'Q1*(1.-vegfrac),EDIR1',q1*(one-vegfrac),edir1
2903 print *,
'CST,WETCAN,DRYCAN',cst,wetcan,drycan
2904 print *,
'EC1=',ec1,
'ETT1=',ett1,
'CMC2MS=',cmc2ms,
'CMC2MS*ras=',cmc2ms*ras
2907 eeta = (edir1 + ec1 + ett1)*rhowater
2908 IF (debug_print )
THEN
2910 print *,
'RUC LSM EETA',eeta,eeta*xlv
2914 eeta = (edir1 + ec1 + ett1)*rhowater
2916 IF (debug_print )
THEN
2917 print *,
'potential temp HFT ',hft
2918 print *,
'abs temp HFX ',hfx
2922 s=thdif(1)*cap(1)*dzstop*(tso(1)-tso(2))
2924 fltot=rnet-hft-xlv*eeta-s-x
2925 IF (debug_print )
THEN
2927 print *,
'SOIL - FLTOT,RNET,HFT,QFX,S,X=',i,j,fltot,rnet,hft,xlv*eeta,s,x
2928 print *,
'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac',&
2929 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac
2931 if(detal(1) .ne. zero)
then
2934 IF (debug_print )
THEN
2935 print *,
'detal(1),xlmelt,soiliqwm(1),delt',detal(1),xlmelt,soiliqwm(1),delt
2936 print *,
'Implicit phase change in the first layer - smf=',smf
2943 1123
FORMAT(i5,8f12.3)
2944 1133
FORMAT(i7,8e12.4)
2945 123
format(i6,f6.2,7f8.1)
2946 122
FORMAT(1x,2i3,6f8.1,f8.3,f8.2)
3223 testptlat,testptlon, &
3224 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & !--- input variables
3225 isncond_opt,isncovr_opt, &
3226 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
3227 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
3230 GLW,GSW,GSWin,EMISS,RNET,IVGTYP, &
3231 QKMS,TKMS,PC,cst,drip,infwater, &
3232 rho,vegfrac,alb,znt,lai, &
3233 MYJ, & !--- soil fixed fields
3234 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
3235 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
3236 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, & !--- constants
3238 ilnb,snweprint,snheiprint,rsm, & !--- output variables
3239 soilmois,tso,smfrkeep,keepfr, &
3240 dew,soilt,soilt1,tsnav, &
3241 qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, &
3242 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
3243 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
3317 LOGICAL,
INTENT(IN ) :: debug_print
3318 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
3320 INTEGER,
INTENT(IN ) :: i,j,isoil,isncond_opt,isncovr_opt
3322 real (kind_phys),
INTENT(IN ) :: DELT,CONFLX,PRCPMS, &
3323 RAINF,NEWSNOW,RHONEWSN, &
3324 testptlat,testptlon, &
3325 SNHEI_CRIT,meltfactor,xlat,xlon
3327 LOGICAL,
INTENT(IN ) :: myj
3331 INTENT(IN ) :: PATM, &
3335 real (kind_phys) , &
3336 INTENT(IN ) :: GLW, &
3347 INTEGER,
INTENT(IN ) :: IVGTYP
3349 real (kind_phys) , &
3350 INTENT(IN ) :: RHOCS, &
3361 real (kind_phys),
INTENT(IN ) :: CN, &
3370 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: ZSMAIN, &
3374 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: DTDZS
3376 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: TBQ
3381 real (kind_phys),
DIMENSION( 1:nzs ) , &
3382 INTENT(INOUT) :: TSO, &
3386 real (kind_phys),
DIMENSION( 1:nzs ) , &
3387 INTENT(INOUT) :: KEEPFR
3390 INTEGER,
INTENT(INOUT) :: ILAND
3394 real (kind_phys) , &
3395 INTENT(INOUT) :: DEW, &
3428 INTEGER,
INTENT(INOUT) :: ILNB
3431 real (kind_phys),
DIMENSION(1:NZS),
INTENT(OUT) :: SOILICE, &
3434 real (kind_phys),
INTENT(OUT) :: RSM, &
3440 INTEGER :: nzs1,nzs2,k
3442 real (kind_phys) :: INFILTRP, TRANSUM , &
3444 TABS, T3, UPFLUX, XINET , &
3445 BETA, SNWEPR,EPDT,PP
3446 real (kind_phys) :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt,dzstop, &
3447 can,epot,fac,fltot,ft,fq,hft , &
3449 trans,zn,ci,cvw,tln,tavln,pi , &
3450 DD1,CMC2MS,DRYCAN,WETCAN , &
3451 INFMAX,RIW,DELTSN,H,UMVEG
3453 real (kind_phys),
DIMENSION(1:NZS) :: transp,cap,diffu,hydro, &
3454 thdif,tranf,tav,soilmoism , &
3455 soilicem,soiliqwm,detal , &
3456 fwsat,lwsat,told,smold
3457 real (kind_phys) :: soiltold, qgold
3459 real (kind_phys) :: RNET, X
3484 deltsn=0.05_kind_phys*rhowater/rhosn
3485 snth=0.01_kind_phys*rhowater/rhosn
3489 IF(snhei.GE.deltsn+snth)
THEN
3490 if(snhei-deltsn-snth.lt.snth) deltsn=0.5_kind_phys*(snhei-snth)
3491 IF (debug_print )
THEN
3492 print *,
'DELTSN is changed,deltsn,snhei,snth',i,j,deltsn,snhei,snth
3497 ras=rho*1.e-3_kind_dbl_prec
3498 riw=rhoice*1.e-3_kind_dbl_prec
3531 dzstop=one/(zsmain(2)-zsmain(1))
3539 tln=log(tso(k)/tfrz)
3540 if(tln.lt.zero)
then
3541 soiliqw(k)=(dqm+qmin)*(xlmelt* &
3542 (tso(k)-tfrz)/tso(k)/grav/psis) &
3544 soiliqw(k)=max(zero,soiliqw(k))
3545 soiliqw(k)=min(soiliqw(k),soilmois(k))
3546 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3549 if(keepfr(k).eq.1.)
then
3550 soilice(k)=min(soilice(k),smfrkeep(k))
3551 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
3556 soiliqw(k)=soilmois(k)
3563 tav(k)=0.5_kind_phys*(tso(k)+tso(k+1))
3564 soilmoism(k)=0.5_kind_phys*(soilmois(k)+soilmois(k+1))
3565 tavln=log(tav(k)/tfrz)
3567 if(tavln.lt.zero)
then
3568 soiliqwm(k)=(dqm+qmin)*(xlmelt* &
3569 (tav(k)-tfrz)/tav(k)/grav/psis) &
3571 fwsat(k)=dqm-soiliqwm(k)
3572 lwsat(k)=soiliqwm(k)+qmin
3573 soiliqwm(k)=max(zero,soiliqwm(k))
3574 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
3575 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
3577 if(keepfr(k).eq.one)
then
3578 soilicem(k)=min(soilicem(k), &
3579 0.5_kind_phys*(smfrkeep(k)+smfrkeep(k+1)))
3580 soiliqwm(k)=max(zero,soilmoism(k)-soilicem(k)*riw)
3581 fwsat(k)=dqm-soiliqwm(k)
3582 lwsat(k)=soiliqwm(k)+qmin
3587 soiliqwm(k)=soilmoism(k)
3595 if(soilice(k).gt.zero)
then
3596 smfrkeep(k)=soilice(k)
3598 smfrkeep(k)=soilmois(k)/riw
3607 xlat, xlon, testptlat, testptlon, &
3609 nzs,fwsat,lwsat,tav,keepfr, &
3610 soilmois,soiliqw,soilice, &
3611 soilmoism,soiliqwm,soilicem, &
3613 qwrtz,rhocs,dqm,qmin,psis,bclh,ksat, &
3615 riw,xlmelt,cp,g0_p,cvw,ci, &
3618 thdif,diffu,hydro,cap)
3634 epot = -fq*(qvatm-qsg)
3636 IF (debug_print )
THEN
3637 print *,
'SNWE after subtracting intercepted snow - snwe=',snwe,vegfrac,cst
3645 epdt = epot * ras *delt
3646 IF(epdt > zero .and. snwepr.LE.epdt)
THEN
3651 wetcan=min(0.25_kind_phys,max(zero,(cst/sat))**cn)
3657 CALL transf(debug_print, &
3658 xlat, xlon, testptlat, testptlon, &
3660 nzs,nroot,soiliqw,tabs,lai,gswin, &
3662 dqm,qmin,ref,wilt,zshalf,pc,iland, &
3669 smold(k)=soilmois(k)
3676 IF (debug_print )
THEN
3677print *,
'TSO before calling SNOWTEMP: ', tso
3679 CALL snowtemp(debug_print,xlat,xlon,testptlat,testptlon,&
3682 delt,ktau,conflx,nzs,nddzs,nroot, &
3683 isncond_opt,isncovr_opt, &
3684 snwe,snwepr,snhei,newsnow,snowfrac,snhei_crit, &
3685 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, &
3687 patm,tabs,qvatm,qcatm, &
3688 glw,gsw,emiss,rnet, &
3689 qkms,tkms,pc,rho,vegfrac, &
3690 thdif,cap,drycan,wetcan,cst, &
3691 tranf,transum,dew,mavail, &
3693 dqm,qmin,psis,bclh, &
3694 zsmain,zshalf,dtdzs,tbq, &
3696 xlvm,cp,rovcp,g0_p,cvw,stbolt, &
3698 snweprint,snheiprint,rsm, &
3699 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3700 smelt,snoh,snflx,s,ilnb,x)
3707 epot = -fq*(qvatm-qsg)
3708 IF(epot.GT.zero)
THEN
3711 transp(k)=vegfrac*ras*fq*(qvatm-qsg) &
3712 *tranf(k)*drycan/zshalf(nroot+1)
3730 tln=log(tso(k)/tfrz)
3731 if(tln.lt.zero)
then
3732 soiliqw(k)=(dqm+qmin)*(xlmelt* &
3733 (tso(k)-tfrz)/tso(k)/grav/psis) &
3735 soiliqw(k)=max(zero,soiliqw(k))
3736 soiliqw(k)=min(soiliqw(k),soilmois(k))
3737 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3739 if(keepfr(k).eq.one)
then
3740 soilice(k)=min(soilice(k),smfrkeep(k))
3741 soiliqw(k)=max(zero,soilmois(k)-soilice(k)*riw)
3746 soiliqw(k)=soilmois(k)
3754 CALL soilmoist (debug_print,xlat,xlon,testptlat,testptlon,&
3756 delt,nzs,nddzs,dtdzs,dtdzs2,riw, &
3757 zsmain,zshalf,diffu,hydro, &
3758 qsg,qvg,qcg,qcatm,qvatm,-infwater, &
3760 zero,smelt,soilice,vegfrac, &
3763 dqm,qmin,ref,ksat,ras,infmax, &
3765 soilmois,soiliqw,mavail,runoff1, &
3771 IF(snhei.EQ.zero)
then
3777 snom=snom+smelt*delt*rhowater
3789 if (soilice(k).gt.zero)
then
3790 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k))
then
3799 t3 = stbolt*soiltold*soiltold*soiltold
3800 upflux = t3 *0.5_kind_phys*(soiltold+soilt)
3801 xinet = emiss*(glw-upflux)
3802 hfx=-tkms*cp*rho*(tabs-soilt) &
3803 *(p1000mb*0.00001_kind_phys/patm)**rovcp
3804 IF (debug_print )
THEN
3805 print *,
'potential temp HFX',hfx
3807 hft=-tkms*cp*rho*(tabs-soilt)
3808 IF (debug_print )
THEN
3809 print *,
'abs temp HFX',hft
3811 q1 = - fq*ras* (qvatm - qsg)
3813 IF (q1.LT.zero)
THEN
3821 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*rhowater
3822 cst= cst-eeta*delt*vegfrac
3823 IF (debug_print )
THEN
3824 print *,
'MYJ EETA cond', eeta
3828 dew=qkms*(qvatm-qsg)
3830 cst=cst+delt*dew*ras * vegfrac
3831 IF (debug_print )
THEN
3832 print *,
'RUC LSM EETA cond',eeta
3839 edir1 = q1*umveg *beta
3841 ec1 = q1 * wetcan * vegfrac
3843 cst=max(zero,cst-ec1 * delt)
3845 IF (debug_print )
THEN
3846 print*,
'Q1,umveg,beta',q1,umveg,beta
3847 print *,
'wetcan,vegfrac',wetcan,vegfrac
3848 print *,
'EC1,CMC2MS',ec1,cmc2ms
3853 eeta=-(qkms*ras*(qvatm/(one+qvatm) - qsg/(one+qsg))*rhowater)*beta
3854 IF (debug_print )
THEN
3855 print *,
'MYJ EETA', eeta*xlvm,eeta
3860 eeta = (edir1 + ec1 + ett1)*rhowater
3861 IF (debug_print )
THEN
3862 print *,
'RUC LSM EETA',eeta*xlvm,eeta
3866 eeta = (edir1 + ec1 + ett1)*rhowater
3871 fltot=rnet-hft-xlvm*eeta-s-snoh-x
3872 IF (debug_print )
THEN
3873 print *,
'SNOWSOIL - FLTOT,RNET,HFT,QFX,S,SNOH,X=',fltot,rnet,hft,xlvm*eeta,s,snoh,x
3874 print *,
'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta',&
3875 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta
3880 1123
FORMAT(i5,8f12.3)
3881 1133
FORMAT(i7,8e12.4)
3882 123
format(i6,f6.2,7f8.1)
3883 122
FORMAT(1x,2i3,6f8.1,f8.3,f8.2)
3895 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
3896 isncond_opt,isncovr_opt, &
3897 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
3898 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,snowfrac, &
3899 RHOSN,PATM,QVATM,QCATM, &
3900 GLW,GSW,EMISS,RNET, &
3901 QKMS,TKMS,RHO,myj, &
3902 ALB,ZNT, & !--- sea ice parameters
3903 tice,rhosice,capice,thdifice, &
3904 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
3905 xlv,CP,rovcp,cw,stbolt,tabs, & !--- constants
3906 ilnb,snweprint,snheiprint,rsm,tso, & !--- output variables
3907 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3908 SMELT,SNOH,SNFLX,SNOM,eeta, &
3909 qfx,hfx,s,sublim,prcpl,fltot &
3921 LOGICAL,
INTENT(IN ) :: debug_print
3922 INTEGER,
INTENT(IN ) :: ktau,nzs , &
3924 INTEGER,
INTENT(IN ) :: i,j,isoil,isncond_opt,isncovr_opt
3926 real (kind_phys),
INTENT(IN ) :: DELT,CONFLX,PRCPMS, &
3927 RAINF,NEWSNOW,RHONEWSN, &
3928 meltfactor,snhei_crit,xlat,xlon
3931 LOGICAL,
INTENT(IN ) :: myj
3934 INTENT(IN ) :: PATM, &
3938 real (kind_phys) , &
3939 INTENT(IN ) :: GLW, &
3946 real (kind_phys),
DIMENSION(1:NZS) , &
3953 real (kind_phys),
INTENT(IN ) :: &
3957 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: ZSMAIN, &
3961 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: DTDZS
3963 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: TBQ
3967 real (kind_phys),
DIMENSION( 1:nzs ) , &
3968 INTENT(INOUT) :: TSO
3970 INTEGER,
INTENT(INOUT) :: ILAND
3974 real (kind_phys) , &
3975 INTENT(INOUT) :: DEW, &
4000 INTEGER,
INTENT(INOUT) :: ILNB
4002 real (kind_phys),
INTENT(OUT) :: RSM, &
4008 INTEGER :: nzs1,nzs2,k,k1,kn,kk
4009 real (kind_phys) :: x,x1,x2,dzstop,ft,tn,denom
4011 real (kind_phys) :: SNTH, NEWSN , &
4012 TABS, T3, UPFLUX, XINET , &
4013 BETA, SNWEPR,EPDT,PP
4014 real (kind_phys) :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt , &
4015 epot,fltot,fq,hft,q1,ras,ci,cvw , &
4018 real (kind_phys) :: rhocsn,thdifsn, &
4019 xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
4021 real (kind_phys) :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
4022 real (kind_phys) :: fso,fsn, &
4023 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
4024 FKQ,R210,AA,BB,QS1,TS1,TQ2,TX2, &
4025 TDENOM,AA1,RHCS,H1,TSOB, SNPRIM, &
4026 SNODIF,SOH,TNOLD,QGOLD,SNOHGNEW
4027 real (kind_phys),
DIMENSION(1:NZS) :: cotso,rhtso
4029 real (kind_phys) :: RNET,rsmfrac,soiltfrac,hsn,icemelt,rr
4032 real (kind_phys) :: keff, fact
4042 keff = 0.265_kind_phys
4053 deltsn=0.05_kind_phys*rhowater/rhosn
4054 snth=0.01_kind_phys*rhowater/rhosn
4058 IF(snhei.GE.deltsn+snth)
THEN
4059 if(snhei-deltsn-snth.lt.snth) deltsn=0.5_kind_phys*(snhei-snth)
4060 IF (debug_print )
THEN
4061 print *,
'DELTSN ICE is changed,deltsn,snhei,snth', &
4062 i,j, deltsn,snhei,snth
4067 ras=rho*1.e-3_kind_dbl_prec
4068 riw=rhoice*1.e-3_kind_dbl_prec
4072 rhocsn=sheatsn * rhosn
4074 rhonewcsn=sheatsn * rhonewsn
4076 if(isncond_opt == 1)
then
4078 thdifsn = 0.265_kind_phys/rhocsn
4083 if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys))
then
4084 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
4086 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
4089 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
4095 thdifsn = 4.431718e-7_kind_phys
4097 thdifsn = keff/rhocsn * fact
4101 ras=rho*1.e-3_kind_phys
4121 dzstop=one/(zsmain(2)-zsmain(1))
4136 snhei=snwe*rhowater/rhosn
4141 epot = -fq*(qvatm-qsg)
4142 epdt = epot * ras *delt
4143 IF(epdt.GT.zero .and. snwepr.LE.epdt)
THEN
4144 beta=snwepr/max(1.e-8_kind_phys,epdt)
4157 x1=dtdzs(k1)*thdifice(kn-1)
4158 x2=dtdzs(k1+1)*thdifice(kn)
4159 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
4160 -x2*(tso(kn)-tso(kn+1))
4161 denom=1.+x1+x2-x2*cotso(k)
4163 rhtso(k+1)=(ft+x2*rhtso(k))/denom
4167 IF(snhei.GE.snth)
then
4168 if(snhei.le.deltsn+snth)
then
4171 snprim=max(snth,snhei)
4174 xsn = delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*snprim)
4175 ddzsn = xsn / snprim
4176 x1sn = ddzsn * thdifsn
4177 x2 = dtdzs(1)*thdifice(1)
4178 ft = tso(1)+x1sn*(soilt-tso(1)) &
4180 denom = one + x1sn + x2 -x2*cotso(nzs1)
4181 cotso(nzs)=x1sn/denom
4182 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
4186 tsnav=0.5_kind_phys*(soilt+tso(1)) &
4194 xsn = delt/2._kind_phys/(0.5_kind_phys*snhei)
4195 xsn1= delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*(snhei-deltsn))
4196 ddzsn = xsn / deltsn
4197 ddzsn1 = xsn1 / (snhei-deltsn)
4198 x1sn = ddzsn * thdifsn
4199 x1sn1 = ddzsn1 * thdifsn
4200 x2 = dtdzs(1)*thdifice(1)
4201 ft = tso(1)+x1sn1*(soilt1-tso(1)) &
4203 denom = one + x1sn1 + x2 - x2*cotso(nzs1)
4204 cotso(nzs)=x1sn1/denom
4205 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
4206 ftsnow = soilt1+x1sn*(soilt-soilt1) &
4207 -x1sn1*(soilt1-tso(1))
4208 denomsn = 1. + x1sn + x1sn1 - x1sn1*cotso(nzs)
4210 rhtsn=(ftsnow+x1sn1*rhtso(nzs))/denomsn
4212 tsnav=0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
4213 +(soilt1+tso(1))*(snhei-deltsn)) &
4218 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
4221 snprim=snhei+zsmain(2)
4226 xsn = delt/2._kind_phys/((zshalf(3)-zsmain(2))+0.5_kind_phys*snprim)
4228 x1sn = ddzsn * (fsn*thdifsn+fso*thdifice(1))
4229 x2=dtdzs(2)*thdifice(2)
4230 ft=tso(2)+x1sn*(soilt-tso(2))- &
4232 denom = one + x1sn + x2 - x2*cotso(nzs-2)
4233 cotso(nzs1) = x1sn/denom
4234 rhtso(nzs1)=(ft+x2*rhtso(nzs-2))/denom
4235 tsnav=0.5_kind_phys*(soilt+tso(1)) &
4237 cotso(nzs)=cotso(nzs1)
4238 rhtso(nzs)=rhtso(nzs1)
4249 epot=-qkms*(qvatm-qsg)
4256 d9=thdifice(1)*rhcs*dzstop
4258 r211=.5_kind_phys*conflx/delt
4260 r22=.5_kind_phys/(thdifice(1)*delt*dzstop**2)
4261 r6=emiss *stbolt*.5_kind_phys*tn**4
4265 IF(snhei.GE.snth)
THEN
4266 if(snhei.le.deltsn+snth)
then
4275 d9sn= thdifsn*rhocsn / snprim
4276 r22sn = snprim*snprim*0.5_kind_phys/(thdifsn*delt)
4279 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
4283 d9sn = (fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)/ &
4285 r22sn = snprim*snprim*0.5_kind_phys &
4286 /((fsn*thdifsn+fso*thdifice(1))*delt)
4289 IF(snhei.eq.zero)
then
4299 tdenom = d9sn*(one-d1sn +r22sn)+d10+r21+r7 &
4301 +rhonewcsn*newsnow/delt
4305 aa=xlvm*(beta*fkq+r210)/tdenom
4306 bb=(d10*tabs+r21*tn+xlvm*(qvatm* &
4308 +r210*qvg)+d11+d9sn*(d2sn+r22sn*tn) &
4309 +rainf*cvw*prcpms*max(tfrz,tabs) &
4310 + rhonewcsn*newsnow/delt*min(tfrz,tabs) &
4313 pp=patm*1.e3_kind_phys
4318 IF (debug_print )
THEN
4319 print *,
'VILKA-SNOW on SEAICE'
4320 print *,
'tn,aa1,bb,pp,fkq,r210', &
4321 tn,aa1,bb,pp,fkq,r210
4322 print *,
'TABS,QVATM,TN,QVG=',tabs,qvatm,tn,qvg
4325 CALL vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
4333 if(nmelt==1 .and. snowfrac==one)
then
4334 soilt = min(tfrz,soilt)
4337 IF (debug_print )
THEN
4338 print *,
' AFTER VILKA-SNOW on SEAICE'
4339 print *,
' TS1,QS1: ', ts1,qs1
4342 IF(snhei.GE.snth)
THEN
4343 if(snhei.gt.deltsn+snth)
then
4345 soilt1=min(tfrz,rhtsn+cotsn*soilt)
4346 tso(1)=min(con_tice,(rhtso(nzs)+cotso(nzs)*soilt1))
4350 tso(1)=min(con_tice,(rhtso(nzs)+cotso(nzs)*soilt))
4354 ELSEIF (snhei > zero .and. snhei < snth)
THEN
4356 tso(2)=min(con_tice,(rhtso(nzs1)+cotso(nzs1)*soilt))
4357 tso(1)=min(con_tice,(tso(2)+(soilt-tso(2))*fso))
4362 tso(1)=min(con_tice,soilt)
4363 soilt1=min(con_tice,soilt)
4367 IF (snhei > zero .and. snhei < snth)
THEN
4371 tso(k)=min(con_tice,rhtso(kk)+cotso(kk)*tso(k-1))
4376 tso(k)=min(con_tice,rhtso(kk)+cotso(kk)*tso(k-1))
4388 if(nmelt.eq.1)
go to 220
4392 IF(soilt>tfrz .AND. beta==one .AND. snhei>zero)
THEN
4395 soiltfrac=snowfrac*tfrz+(1.-snowfrac)*min(con_tice,soilt)
4397 qsg=
qsn(soiltfrac,tbq)/pp
4398 t3 = stbolt*tnold*tnold*tnold
4399 upflux = t3 * 0.5_kind_phys*(tnold+soiltfrac)
4400 xinet = emiss*(glw-upflux)
4401 epot = -qkms*(qvatm-qsg)
4404 IF (q1.LE.zero)
THEN
4412 eeta = q1 * beta * rhowater
4417 hfx=d10*(tabs-soiltfrac)
4419 IF(snhei.GE.snth)
then
4420 soh=thdifsn*rhocsn*(soiltfrac-tsob)/snprim
4423 soh=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4424 (soiltfrac-tsob)/snprim
4427 x= (r21+d9sn*r22sn)*(soiltfrac-tnold) + &
4428 xlvm*r210*(qsg-qgold)
4430 snoh=rnet+qfx +hfx &
4431 +rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soiltfrac) &
4432 -soh-x+rainf*cvw*prcpms* &
4433 (max(tfrz,tabs)-soiltfrac)
4435 IF (debug_print )
THEN
4436 print *,
'SNOWSEAICE melt I,J,SNOH,RNET,QFX,HFX,SOH,X',i,j,snoh,rnet,qfx,hfx,soh,x
4437 print *,
'RHOnewCSN*NEWSNOW/DELT*(min(tfrz,TABS)-soiltfrac)', &
4438 rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soiltfrac)
4439 print *,
'RAINF*CVW*PRCPMS*(max(tfrz,TABS)-soiltfrac)', &
4440 rainf*cvw*prcpms*(max(tfrz,tabs)-soiltfrac)
4442 snoh=amax1(zero,snoh)
4444 smelt= snoh /xlmelt*1.e-3_kind_phys
4445 smelt=amin1(smelt,snwepr/delt-beta*epot*ras)
4446 smelt=amax1(zero,smelt)
4448 IF (debug_print )
THEN
4449 print *,
'1-SMELT i,j',smelt,i,j
4452 smelt= amin1(smelt,delt/60._kind_phys* 5.6e-8_kind_phys*meltfactor*max(one,(soilt-tfrz)))
4453 IF (debug_print )
THEN
4454 print *,
'2-SMELT i,j',smelt,i,j
4458 rr=snwepr/delt-beta*epot*ras
4460 IF (debug_print )
THEN
4461 print *,
'3- SMELT i,j,smelt,rr',i,j,smelt,rr
4463 snohgnew=smelt*xlmelt*1.e3
4464 snodif=amax1(zero,(snoh-snohgnew))
4468 IF (debug_print )
THEN
4469 print*,
'soiltfrac,soilt,SNOHGNEW,SNODIF=', &
4470 i,j,soiltfrac,soilt,snohgnew,snodif
4471 print *,
'SNOH,SNODIF',snoh,snodif
4475 rsmfrac=min(0.18_kind_phys,(max(0.08_kind_phys,snwepr/0.10_kind_phys*0.13_kind_phys)))
4476 if(snhei > 0.01_kind_phys)
then
4477 rsm=rsmfrac*smelt*delt
4483 smelt=amax1(zero,smelt-rsm/delt)
4484 IF (debug_print )
THEN
4485 print *,
'4-SMELT i,j,smelt,rsm,snwepr,rsmfrac', &
4486 i,j,smelt,rsm,snwepr,rsmfrac
4491 snwe = amax1(zero,(snwepr- &
4492 (smelt+beta*epot*ras)*delt &
4498 if(snhei > zero.and. beta == one)
then
4499 epot=-qkms*(qvatm-qsg)
4500 snwe = amax1(zero,(snwepr- &
4501 beta*epot*ras*delt))
4513 if(smelt > zero .and. rsm > zero)
then
4514 if(snwe.le.rsm)
then
4515 IF (debug_print )
THEN
4516 print *,
'SEAICE SNWE<RSM snwe,rsm,smelt*delt,epot*ras*delt,beta', &
4517 snwe,rsm,smelt*delt,epot*ras*delt,beta
4525 xsn=(rhosn*(snwe-rsm)+rhowater*rsm)/ &
4527 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
4529 rhocsn=sheatsn* rhosn
4530 if(isncond_opt == 1)
then
4532 thdifsn = 0.265_kind_phys/rhocsn
4537 if(rhosn < 156._kind_phys .or. (newsn > zero .and. rhonewsn < 156._kind_phys))
then
4538 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
4540 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
4543 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
4549 thdifsn = 4.431718e-7_kind_phys
4551 thdifsn = keff/rhocsn * fact
4562 snheiprint=snweprint*rhowater / rhosn
4564 IF (debug_print )
THEN
4565print *,
'snweprint : ',snweprint
4566print *,
'D9SN,SOILT,TSOB : ', d9sn,soilt,tsob
4568 IF(snhei.GT.zero)
THEN
4570 tsnav=0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
4571 +(soilt1+tso(1))*(snhei-deltsn)) &
4574 tsnav=0.5_kind_phys*(soilt+tso(1)) - tfrz
4579 pp=patm*1.e3_kind_phys
4580 qsg=
qsn(soilt,tbq)/pp
4581 epot = -fq*(qvatm-qsg)
4582 IF(epot.LT.zero)
THEN
4587 snom=snom+smelt*delt*rhowater
4591 t3 = stbolt*tnold*tnold*tnold
4592 upflux = t3 *0.5_kind_phys*(soilt+tnold)
4593 xinet = emiss*(glw-upflux)
4594 hft=-tkms*cp*rho*(tabs-soilt)
4595 hfx=-tkms*cp*rho*(tabs-soilt) &
4596 *(p1000mb*0.00001_kind_phys/patm)**rovcp
4597 q1 = - fq*ras* (qvatm - qsg)
4598 IF (q1.LT.zero)
THEN
4602 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*rhowater
4605 dew=qkms*(qvatm-qsg)
4615 eeta=-qkms*ras*beta*(qvatm/(1.+qvatm) - qvg/(1.+qvg))*rhowater
4619 eeta = q1*beta*rhowater
4622 eeta = q1*beta*rhowater
4627 IF(snhei.GE.snth)
then
4628 s=thdifsn*rhocsn*(soilt-tsob)/snprim
4630 ELSEIF(snhei.lt.snth.and.snhei.GT.zero)
then
4631 s=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4634 IF (debug_print )
THEN
4635 print *,
'SNOW is thin, snflx',i,j,snflx
4638 snflx=d9sn*(soilt-tsob)
4639 IF (debug_print )
THEN
4640 print *,
'SNOW is GONE, snflx',i,j,snflx
4644 snhei=snwe *rhowater / rhosn
4646 IF (debug_print )
THEN
4647 print *,
'SNHEI,SNOH',i,j,snhei,snoh
4650 x= (r21+d9sn*r22sn)*(soilt-tnold) + &
4651 xlvm*r210*(qsg-qgold)
4652 IF (debug_print )
THEN
4653 print *,
'SNOWSEAICE storage ',i,j,x
4654 print *,
'R21,D9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim', &
4655 r21,d9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim
4658 -rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soilt) &
4659 -rainf*cvw*prcpms*(max(tfrz,tabs)-soilt)
4662 icemelt = rnet-hft-xlvm*eeta-s-snoh-x
4663 IF (debug_print )
THEN
4664 print *,
'SNOWSEAICE icemelt=',icemelt
4667 fltot=rnet-hft-xlvm*eeta-s-snoh-x-icemelt
4668 IF (debug_print )
THEN
4669 print *,
'i,j,snhei,qsg,soilt,soilt1,tso,TABS,QVATM', &
4670 i,j,snhei,qsg,soilt,soilt1,tso,tabs,qvatm
4671 print *,
'SNOWSEAICE - FLTOT,RNET,HFT,QFX,S,SNOH,icemelt,snodif,X,SOILT=' &
4672 ,fltot,rnet,hft,xlvm*eeta,s,snoh,icemelt,snodif,x,soilt
4675 IF(snhei.EQ.zero)
then
4677 emiss=0.98_kind_phys
4993 testptlat,testptlon,i,j,iland,isoil, & !--- input variables
4994 delt,ktau,conflx,nzs,nddzs,nroot, &
4995 isncond_opt,isncovr_opt, &
4996 snwe,snwepr,snhei,newsnow,snowfrac,snhei_crit, &
4997 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
4999 PATM,TABS,QVATM,QCATM, &
5000 GLW,GSW,EMISS,RNET, &
5001 QKMS,TKMS,PC,RHO,VEGFRAC, &
5002 THDIF,CAP,DRYCAN,WETCAN,CST, &
5003 TRANF,TRANSUM,DEW,MAVAIL, &
5004 DQM,QMIN,PSIS,BCLH, & !--- soil fixed fields
5005 ZSMAIN,ZSHALF,DTDZS,TBQ, &
5006 XLVM,CP,rovcp,G0_P,CVW,STBOLT, & !--- constants
5007 SNWEPRINT,SNHEIPRINT,RSM, & !--- output variables
5008 TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, &
5009 SMELT,SNOH,SNFLX,S,ILNB,X)
5063 LOGICAL,
INTENT(IN ) :: debug_print
5064 INTEGER,
INTENT(IN ) :: nroot,ktau,nzs , &
5067 INTEGER,
INTENT(IN ) :: i,j,iland,isoil,isncond_opt,isncovr_opt
5068 real (kind_phys),
INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
5069 RAINF,NEWSNOW,DELTSN,SNTH , &
5070 TABS,TRANSUM,SNWEPR , &
5071 testptlat,testptlon , &
5072 rhonewsn,meltfactor,xlat,xlon,snhei_crit
5077 INTENT(IN ) :: PATM, &
5081 real (kind_phys) , &
5082 INTENT(IN ) :: GLW, &
5091 real (kind_phys) , &
5098 real (kind_phys),
INTENT(IN ) :: CP, &
5106 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: ZSMAIN, &
5112 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: DTDZS
5114 real (kind_phys),
DIMENSION(1:5001),
INTENT(IN) :: TBQ
5119 real (kind_phys),
DIMENSION( 1:nzs ) , &
5120 INTENT(INOUT) :: TSO
5124 real (kind_phys) , &
5125 INTENT(INOUT) :: DEW, &
5144 real (kind_phys),
INTENT(INOUT) :: DRYCAN, WETCAN
5146 real (kind_phys),
INTENT(OUT) :: RSM, &
5149 INTEGER,
INTENT(OUT) :: ilnb
5153 INTEGER :: nzs1,nzs2,k,k1,kn,kk
5155 real (kind_phys) :: x,x1,x2,x4,dzstop,can,ft,sph, &
5156 tn,trans,umveg,denom
5158 real (kind_phys) :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
5160 real (kind_phys) :: t3,upflux,xinet,ras, &
5161 xlmelt,rhocsn,thdifsn, &
5162 beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
5164 real (kind_phys) :: fso,fsn, &
5165 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
5166 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2, &
5167 TDENOM,C,CC,AA1,RHCS,H1, &
5168 tsob, snprim, sh1, sh2, &
5169 smeltg,snohg,snodif,soh, &
5170 CMC2MS,TNOLD,QGOLD,SNOHGNEW
5172 real (kind_phys),
DIMENSION(1:NZS) :: transp,cotso,rhtso
5173 real (kind_phys) :: edir1, &
5180 real (kind_phys) :: RNET,rsmfrac,soiltfrac,hsn,rr,keff,fact
5181 integer :: nmelt, iter
5190 keff = 0.265_kind_phys
5198 IF (debug_print )
THEN
5199print *,
'SNOWTEMP: SNHEI,SNTH,SOILT1: ',snhei,snth,soilt1,soilt
5202 rhocsn=sheatsn* rhosn
5203 rhonewcsn=sheatsn* rhonewsn
5204 if(isncond_opt == 1)
then
5206 thdifsn = 0.265_kind_phys/rhocsn
5211 if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys))
then
5212 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
5214 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
5215 if(debug_print)
then
5216 print *,
'SnowTemp xlat,xlon,rhosn,keff', xlat,xlon,rhosn,keff,keff/rhocsn*fact
5217 print *,
'SNOWTEMP - 0.265/rhocsn',0.265_kind_phys/rhocsn
5220 if ( debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5221 print *,
'SNOWTEMP - xlat,xlon,newsnow,rhonewsn,rhosn,fact,keff',xlat,xlon,newsnow, rhonewsn,rhosn,fact,keff
5224 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
5230 thdifsn = 4.431718e-7_kind_phys
5232 thdifsn = keff/rhocsn * fact
5234 if (debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5235 print *,
'SNOWTEMP - thdifsn',xlat,xlon,thdifsn
5236 print *,
'SNOWTEMP - 0.265/rhocsn',0.265_kind_phys/rhocsn
5241 ras=rho*1.e-3_kind_phys
5259 dzstop=one/(zsmain(2)-zsmain(1))
5269 x1=dtdzs(k1)*thdif(kn-1)
5270 x2=dtdzs(k1+1)*thdif(kn)
5271 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
5272 -x2*(tso(kn)-tso(kn+1))
5273 denom=1.+x1+x2-x2*cotso(k)
5275 rhtso(k+1)=(ft+x2*rhtso(k))/denom
5279 IF(snhei.GE.snth)
then
5280 if(snhei.le.deltsn+snth)
then
5282 IF (debug_print )
THEN
5283 print *,
'1-layer - snth,snhei,deltsn',snth,snhei,deltsn
5286 snprim=max(snth,snhei)
5289 xsn = delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*snprim)
5290 ddzsn = xsn / snprim
5291 x1sn = ddzsn * thdifsn
5292 x2 = dtdzs(1)*thdif(1)
5293 ft = tso(1)+x1sn*(soilt-tso(1)) &
5295 denom = one + x1sn + x2 -x2*cotso(nzs1)
5296 cotso(nzs)=x1sn/denom
5297 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
5301 tsnav=min(zero,0.5_kind_phys*(soilt+tso(1))-tfrz)
5305 IF (debug_print )
THEN
5306 print *,
'2-layer - snth,snhei,deltsn',snth,snhei,deltsn
5311 xsn = delt/2._kind_phys/(0.5_kind_phys*deltsn)
5312 xsn1= delt/2._kind_phys/(zshalf(2)+0.5_kind_phys*(snhei-deltsn))
5313 ddzsn = xsn / deltsn
5314 ddzsn1 = xsn1 / (snhei-deltsn)
5315 x1sn = ddzsn * thdifsn
5316 x1sn1 = ddzsn1 * thdifsn
5317 x2 = dtdzs(1)*thdif(1)
5318 ft = tso(1)+x1sn1*(soilt1-tso(1)) &
5320 denom = 1. + x1sn1 + x2 - x2*cotso(nzs1)
5321 cotso(nzs)=x1sn1/denom
5322 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
5323 ftsnow = soilt1+x1sn*(soilt-soilt1) &
5324 -x1sn1*(soilt1-tso(1))
5325 denomsn = one + x1sn + x1sn1 - x1sn1*cotso(nzs)
5327 rhtsn=(ftsnow+x1sn1*rhtso(nzs))/denomsn
5329 tsnav=min(zero,0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
5330 +(soilt1+tso(1))*(snhei-deltsn)) &
5334 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
5337 snprim=snhei+zsmain(2)
5342 xsn = delt/2._kind_phys/((zshalf(3)-zsmain(2))+0.5_kind_phys*snprim)
5344 x1sn = ddzsn * (fsn*thdifsn+fso*thdif(1))
5345 x2=dtdzs(2)*thdif(2)
5346 ft=tso(2)+x1sn*(soilt-tso(2))- &
5348 denom = one + x1sn + x2 - x2*cotso(nzs-2)
5349 cotso(nzs1) = x1sn/denom
5350 rhtso(nzs1)=(ft+x2*rhtso(nzs-2))/denom
5351 tsnav=min(zero,0.5_kind_phys*(soilt+tso(1)) &
5353 cotso(nzs)=cotso(nzs1)
5354 rhtso(nzs)=rhtso(nzs1)
5367 epot=-qkms*(qvatm-qgold)
5370 trans=transum*drycan/zshalf(nroot+1)
5377 d9=thdif(1)*rhcs*dzstop
5379 r211=.5_kind_phys*conflx/delt
5381 r22=.5_kind_phys/(thdif(1)*delt*dzstop**2)
5382 r6=emiss *stbolt*.5_kind_phys*tn**4
5386 IF(snhei.GE.snth)
THEN
5387 if(snhei.le.deltsn+snth)
then
5391 IF (debug_print )
THEN
5392 print *,
'1 layer d1sn,d2sn',i,j,d1sn,d2sn
5398 IF (debug_print )
THEN
5399 print *,
'2 layers d1sn,d2sn',i,j,d1sn,d2sn
5402 d9sn= thdifsn*rhocsn / snprim
5403 r22sn = snprim*snprim*0.5_kind_phys/(thdifsn*delt)
5404 IF (debug_print )
THEN
5405 print *,
'1 or 2 layers D9sn,R22sn',d9sn,r22sn
5409 IF(snhei.LT.snth.AND.snhei.GT.zero)
then
5413 d9sn = (fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)/ &
5415 r22sn = snprim*snprim*0.5_kind_phys &
5416 /((fsn*thdifsn+fso*thdif(1))*delt)
5417 IF (debug_print )
THEN
5418 print *,
' Combined D9SN,R22SN,D1SN,D2SN: ',d9sn,r22sn,d1sn,d2sn
5421 IF(snhei.eq.zero)
then
5427 IF (debug_print )
THEN
5428 print *,
' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',d9sn,r22sn,d1sn,d2sn
5438 tdenom = d9sn*(one-d1sn +r22sn)+d10+r21+r7 &
5440 +rhonewcsn*newsnow/delt
5446 aa=xlvm*(beta*fkq*umveg+r210)/tdenom
5447 bb=(d10*tabs+r21*tn+xlvm*(qvatm* &
5448 (beta*fkq*umveg+c) &
5449 +r210*qgold)+d11+d9sn*(d2sn+r22sn*tn) &
5450 +rainf*cvw*prcpms*max(tfrz,tabs) &
5451 + rhonewcsn*newsnow/delt*min(tfrz,tabs) &
5454 pp=patm*1.e3_kind_phys
5458 IF (debug_print )
THEN
5459 if (abs(xlat-33.35).lt.0.2 .and. abs(xlon-272.55).lt.0.2)
then
5460 print *,
'1-', i,rnet,tabs,tn,aa1,bb,pp,ktau,newsnow,snwepr,snwe,snhei,snowfrac,soilt,soilt1,tso,rhosn
5461 print *,
'2-', i,tdenom,fkq,vegfrac,can,r210,d10,r21,d9sn,d1sn,r22sn,r7,prcpms
5464 CALL vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
5468 IF (debug_print )
THEN
5470 print *,
'VILKA1 - TS1,QS1,TQ2,H,TX2,Q1',ts1,qs1,tq2,h,tx2,q1,xlat,xlon
5472 IF(q1.LT.qs1)
GOTO 100
5477 qcg=max(zero,q1-qs1)
5478 IF (debug_print )
THEN
5479 print *,
'90 QVG,QSG,QCG,TSO(1)',qvg,qsg,qcg,tso(1)
5484 CALL vilka(tn,aa,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil,xlat,xlon)
5486 IF (debug_print )
THEN
5488 print *,
'VILKA2 - TS1,QS1,H,TX2,Q1',ts1,qs1,tq2,h,tx2,q1
5490 IF(q1.GT.qs1)
GOTO 90
5494 IF (debug_print )
THEN
5495 print *,
'No Saturation QVG,QSG,QCG,TSO(1)',qvg,qsg,qcg,tso(1)
5501 if(nmelt==1 .and. snowfrac==one .and. snwe > zero .and. soilt > tfrz)
then
5506 if (debug_print )
then
5508 print *,
'soilt is too high =',soilt,xlat,xlon
5509 soilt = min(tfrz,soilt)
5513 IF (debug_print )
THEN
5515 print *,
'snwe,snwepr,snhei,snowfr,soilt,soilt1,tso',i,j,snwe,snwepr,snhei,snowfrac,soilt,soilt1,tso
5518 IF(snhei.GE.snth)
THEN
5519 if(snhei.gt.deltsn+snth)
then
5521 soilt1=rhtsn+cotsn*soilt
5522 tso(1)=rhtso(nzs)+cotso(nzs)*soilt1
5526 tso(1)=rhtso(nzs)+cotso(nzs)*soilt
5530 ELSEIF (snhei > zero .and. snhei < snth)
THEN
5532 tso(2)=rhtso(nzs1)+cotso(nzs1)*soilt
5533 tso(1)=(tso(2)+(soilt-tso(2))*fso)
5543 if(nmelt==1.and.snowfrac==one)
then
5545 soilt1= min(tfrz,soilt1)
5546 tso(1)= min(tfrz,tso(1))
5547 tsob = min(tfrz,tsob)
5551 IF (snhei > zero .and. snhei < snth)
THEN
5555 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
5561 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
5574 IF (debug_print )
THEN
5576 print *,
'Final SOILT,SOILT1,tso,TSOB,QSG',xlat,xlon,soilt,soilt1,tso,tsob,qsg,
'nmelt=',nmelt
5577 print *,
'SNWEPR-BETA*EPOT*RAS*DELT',snwepr-beta*epot*ras*delt,beta,snwepr,epot
5580 if(nmelt.eq.1)
go to 220
5584 IF(soilt > tfrz.AND.beta==one.AND.snhei>zero)
THEN
5587 soiltfrac=snowfrac*tfrz+(one-snowfrac)*soilt
5588 qsg=min(qsg,
qsn(soiltfrac,tbq)/pp)
5589 qvg=snowfrac*qsg+(one-snowfrac)*qvg
5590 t3 = stbolt*tn*tn*tn
5591 upflux = t3 * 0.5_kind_phys*(tn + soiltfrac)
5592 xinet = emiss*(glw-upflux)
5593 epot = -qkms*(qvatm-qsg)
5597 IF (q1.LE.0..or.iter==1)
THEN
5609 transp(k)=-vegfrac*q1 &
5610 *tranf(k)*drycan/zshalf(nroot+1)
5617 edir1 = q1*umveg * beta
5618 ec1 = q1 * wetcan * vegfrac
5620 eeta = (edir1 + ec1 + ett1)*rhowater
5625 hfx=-d10*(tabs-soiltfrac)
5627 IF(snhei.GE.snth)
then
5628 soh=thdifsn*rhocsn*(soiltfrac-tsob)/snprim
5631 soh=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5632 (soiltfrac-tsob)/snprim
5637 x= (r21+d9sn*r22sn)*(soiltfrac-tn) + &
5638 xlvm*r210*(qvg-qgold)
5639 IF (debug_print )
THEN
5640 print *,
'SNOWTEMP storage ',i,j,x
5641 print *,
'R21,D9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim', &
5642 r21,d9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim
5646 snoh=rnet-qfx -hfx - soh - x &
5647 +rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soiltfrac) &
5648 +rainf*cvw*prcpms*(max(tfrz,tabs)-soiltfrac)
5651 smelt= snoh /xlmelt*1.e-3_kind_phys
5652 IF (debug_print )
THEN
5654 print *,
'1- SMELT',smelt,snoh,xlat,xlon
5657 IF(epot.gt.zero .and. snwepr.LE.epot*ras*delt)
THEN
5659 beta=snwepr/(epot*ras*delt)
5660 smelt=amax1(zero,amin1(smelt,snwepr/delt-beta*epot*ras))
5662 IF (debug_print )
THEN
5664 print *,
'2- SMELT',xlat,xlon,snwe,smelt,rhonewsn,xlat,xlon
5672 if( (rhosn < 350._kind_phys .or. (newsnow > zero .and. rhonewsn < 450._kind_phys)) .and. soilt < 283._kind_phys )
then
5673 smelt= amin1(smelt, delt/60._kind_phys*5.6e-8_kind_phys*meltfactor*max(one,(soilt-tfrz)))
5674 IF (debug_print )
THEN
5676 print *,
'3- SMELT',xlat,xlon,smelt,rhosn,rhonewsn,xlat,xlon
5681 rr=max(zero,snwepr/delt-beta*epot*ras)
5683 smelt = min(smelt,rr)
5685 IF (debug_print )
THEN
5687 print *,
'4- SMELT i,j,smelt,rr',xlat,xlon,smelt,rr
5691 snohgnew=smelt*xlmelt*rhowater
5692 snodif=amax1(zero,(snoh-snohgnew))
5695 IF (debug_print )
THEN
5697 print *,
'SNOH,SNODIF',snoh,snodif
5698 print *,
' xlat, xlon', xlat, xlon
5701 IF( smelt > zero)
then
5703 rsmfrac=min(0.18_kind_phys,(max(0.08_kind_phys,snwepr/0.10_kind_phys*0.13_kind_phys)))
5704 if(snhei > 0.01_kind_phys .and. rhosn < 350._kind_phys)
then
5705 rsm=min(snwe,rsmfrac*smelt*delt)
5713 smelt=max(zero,smelt-rsm/delt)
5714 IF (debug_print )
THEN
5716 print *,
'5- SMELT i,j,smelt,rsm,snwepr,rsmfrac', &
5717 i,j,smelt,rsm,snwepr,rsmfrac
5718 print *,
' xlat, xlon', xlat, xlon
5726 if(snwe > zero)
then
5727 snwe = amax1(zero,(snwepr- &
5728 (smelt+beta*epot*ras)*delt &
5730 IF (debug_print )
THEN
5732 print *,
' Snow is melting and sublimating, snwe', xlat, xlon, snwe
5736 IF (debug_print )
THEN
5738 print *,
' all snwe is sublimated or melted', xlat, xlon, snwe
5745 if(snhei.ne.zero .and. beta == one)
then
5746 epot=-qkms*(qvatm-qsg)
5747 snwe = amax1(zero,(snwepr- &
5748 beta*epot*ras*delt))
5758 if(nmelt.eq.1)
goto 212
5761 if(smelt > zero .and. rsm > zero)
then
5762 if(snwe.le.rsm)
then
5763 IF ( debug_print )
THEN
5764 print *,
'SNWE<RSM snwe,rsm,smelt*delt,epot*ras*delt,beta', &
5765 snwe,rsm,smelt*delt,epot*ras*delt,beta
5772 xsn=(rhosn*(snwe-rsm)+rhowater*rsm)/ &
5774 rhosn=min(max(58.8_kind_phys,xsn),500._kind_phys)
5776 rhocsn=sheatsn* rhosn
5777 if(isncond_opt == 1)
then
5779 thdifsn = 0.265_kind_phys/rhocsn
5784 if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys))
then
5785 keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys
5787 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
5788 if(debug_print)
then
5789 print *,
'End SNOWTEMP - xlat,xlon,rhosn,keff',xlat,xlon,rhosn,keff
5790 print *,
'End SNOWTEMP - 0.265/rhocsn',0.265/rhocsn
5793 if (debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5794 print *,
'END SNOWTEMP - newsnow, rhonewsn,rhosn,fact,keff', &
5795 xlat,xlon,newsnow, rhonewsn,rhosn,fact,keff,keff/rhocsn*fact
5798 if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys)
then
5804 thdifsn = 4.431718e-7_kind_phys
5806 thdifsn = keff/rhocsn * fact
5810 if (debug_print .and. abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)
then
5811 print *,
'END SNOWTEMP - thdifsn',xlat,xlon,thdifsn
5812 print *,
'END SNOWTEMP - 0.265/rhocsn',0.265/rhocsn
5818 IF(snhei.GE.snth)
then
5819 s=thdifsn*rhocsn*(soilt-tsob)/snprim
5821 s=d9*(tso(1)-tso(2))
5822 ELSEIF(snhei.lt.snth.and.snhei.GT.zero)
then
5823 s=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5826 s=d9*(tso(1)-tso(2))
5830 s=d9*(tso(1)-tso(2))
5834 snhei=snwe * rhowater / rhosn
5840 IF (debug_print )
THEN
5842 print *,
'snhei,snwe,rhosn,snowfr',snhei,snwe,rhosn,snowfrac,xlat,xlon
5845 IF(tso(1).GT.tfrz .and. snhei > zero)
THEN
5847 if (snhei.GT.deltsn+snth)
then
5848 hsn = snhei - deltsn
5849 IF (debug_print )
THEN
5850 print*,
'2 layer snow - snhei,hsn',snhei,hsn
5853 IF (debug_print )
THEN
5854 print*,
'1 layer snow or blended - snhei',snhei
5859 soiltfrac=snowfrac*tfrz+(one-snowfrac)*tso(1)
5861 snohg=(tso(1)-soiltfrac)*(cap(1)*zshalf(2)+ &
5862 rhocsn*0.5_kind_phys*hsn) / delt
5863 snohg=amax1(zero,snohg)
5865 smeltg=snohg/xlmelt*1.e-3_kind_phys
5866 IF (debug_print )
THEN
5868 print *,
' SMELTG =',smeltg,xlat,xlon
5872 if( (rhosn < 350._kind_phys .or. (newsnow > zero .and. rhonewsn < 450._kind_phys)) .and. soilt < 283._kind_phys )
then
5873 smelt=amin1(smeltg, 5.8e-9_kind_phys)
5878 smeltg=amin1(smeltg, rr)
5880 snohgnew=smeltg*xlmelt*rhowater
5881 snodif=amax1(zero,(snohg-snohgnew))
5882 IF (debug_print )
THEN
5884 print *,
'TSO(1),soiltfrac,snowfrac,smeltg,SNODIF',tso(1),soiltfrac,snowfrac,smeltg,snodif
5885 print *,
' xlat, xlon', xlat, xlon
5888 snwe=max(zero,snwe-smeltg*delt)
5889 snhei=snwe * rhowater / rhosn
5891 smelt = smelt + smeltg
5893 if(snhei > zero) tso(1) = soiltfrac
5895 IF (debug_print )
THEN
5897 print *,
'Melt from the bottom snwe,snhei',snwe,snhei
5898 print *,
' xlat, xlon', xlat, xlon
5899 print *,
'TSO(1),soiltfrac,snowfrac,smeltg,SNODIF',tso(1),soiltfrac,snowfrac,smeltg,snodif
5900 print *,
'Melt from the bottom snwe,snhei,snoh',snwe,snhei,snoh
5901 print *,
' Final TSO ',tso
5903 print *,
'Snow is all melted on the warm ground'
5909 snheiprint=snweprint*rhowater / rhosn
5911 x= (r21+d9sn*r22sn)*(soilt-tn) + &
5912 xlvm*r210*(qsg-qgold)
5913 IF (debug_print )
THEN
5915 print *,
'end SNOWTEMP storage ',xlat,xlon,x
5916 print *,
'R21,D9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim', &
5917 r21,d9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim
5918 print *,
'snwe, snhei ',snwe,snhei
5923 -rhonewcsn*newsnow/delt*(min(tfrz,tabs)-soilt) &
5924 -rainf*cvw*prcpms*(max(tfrz,tabs)-soilt)
5925 IF (debug_print )
THEN
5927 print *,
'SNHEI=',snhei
5928 print *,
'SNFLX=',snflx
5931 IF(snhei.GT.zero)
THEN
5933 tsnav=min(zero,0.5_kind_phys/snhei*((soilt+soilt1)*deltsn &
5934 +(soilt1+tso(1))*(snhei-deltsn)) &
5937 tsnav=min(zero,0.5_kind_phys*(soilt+tso(1)) - tfrz)
5940 tsnav= min(zero,soilt - tfrz)
5951 xlat, xlon, testptlat, testptlon, &
5952 DELT,NZS,NDDZS,DTDZS,DTDZS2,RIW, & !--- input parameters
5953 ZSMAIN,ZSHALF,DIFFU,HYDRO, &
5954 QSG,QVG,QCG,QCATM,QVATM,PRCP, &
5956 DEW,SMELT,SOILICE,VEGFRAC,SNOWFRAC,soilres, &
5957 DQM,QMIN,REF,KSAT,RAS,INFMAX, & !--- soil properties
5958 SOILMOIS,SOILIQW,MAVAIL,RUNOFF,RUNOFF2,INFILTRP)
6003 LOGICAL,
INTENT(IN ) :: debug_print
6004 real (kind_phys),
INTENT(IN ) :: DELT
6005 real (kind_phys),
INTENT(IN ) :: xlat, xlon, testptlat, testptlon
6006 INTEGER,
INTENT(IN ) :: NZS,NDDZS
6010 real (kind_phys),
DIMENSION(1:NZS),
INTENT(IN) :: ZSMAIN, &
6018 real (kind_phys),
DIMENSION(1:NDDZS),
INTENT(IN) :: DTDZS
6020 real (kind_phys),
INTENT(IN ) :: QSG,QVG,QCG,QCATM,QVATM, &
6021 QKMS,VEGFRAC,DRIP,PRCP , &
6022 DEW,SMELT,SNOWFRAC , &
6023 DQM,QMIN,REF,KSAT,RAS,RIW,SOILRES
6027 real (kind_phys),
DIMENSION( 1:nzs ) , &
6028 INTENT(INOUT) :: SOILMOIS,SOILIQW
6030 real (kind_phys),
INTENT(INOUT) :: MAVAIL,RUNOFF,RUNOFF2,INFILTRP, &
6035 real (kind_phys),
DIMENSION( 1:nzs ) :: COSMC,RHSMC
6037 real (kind_phys) :: DZS,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10
6038 real (kind_phys) :: REFKDT,REFDK,DELT1,F1MAX,F2MAX
6039 real (kind_phys) :: F1,F2,FD,KDT,VAL,DDT,PX,FK,FKMAX
6040 real (kind_phys) :: QQ,UMVEG,INFMAX1,TRANS
6041 real (kind_phys) :: TOTLIQ,FLX,FLXSAT,QTOT
6042 real (kind_phys) :: DID,X1,X2,X4,DENOM,Q2,Q4
6043 real (kind_phys) :: dice,fcr,acrt,frzx,sum,cvfrz
6045 INTEGER :: NZS1,NZS2,K,KK,K1,KN,ialp1,jj,jk
6053 118
format(6(10pf23.19))
6060 did=(zsmain(nzs)-zshalf(nzs))
6061 x1=zsmain(nzs)-zsmain(nzs1)
6063 denom=(one+diffu(nzs1)/x1/did*delt+hydro(nzs)/(2._kind_phys*did)*delt)
6064 cosmc(1)=delt*(diffu(nzs1)/did/x1 &
6065 +hydro(nzs1)/2._kind_phys/did)/denom
6066 rhsmc(1)=(soilmois(nzs)+transp(nzs)*delt/ &
6073 denom=1.+diffu(nzs1)/x1/did*delt
6074 cosmc(1)=delt*(diffu(nzs1)/did/x1 &
6075 +hydro(nzs1)/did)/denom
6077 rhsmc(1)=(soilmois(nzs)-hydro(nzs)*delt/did*soilmois(nzs) &
6078 +transp(nzs)*delt/did)/denom
6081 rhsmc(1)=soilmois(nzs)
6086 x4=2.*dtdzs(k1)*diffu(kn-1)
6087 x2=2.*dtdzs(k1+1)*diffu(kn)
6088 q4=x4+hydro(kn-1)*dtdzs2(kn-1)
6089 q2=x2-hydro(kn+1)*dtdzs2(kn-1)
6090 denom=one+x2+x4-q2*cosmc(k)
6092 IF (debug_print )
THEN
6093 if (abs(xlat-testptlat).lt.0.05 .and. &
6094 abs(xlon-testptlon).lt.0.05)
then
6095 print *,
'xlat,xlon=',xlat,xlon
6096 print *,
'q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k' &
6097 ,q2,soilmois(kn),diffu(kn),x2,hydro(kn+1),dtdzs2(kn-1),kn,k
6100 rhsmc(k+1)=(soilmois(kn)+q2*rhsmc(k) &
6102 /(zshalf(kn+1)-zshalf(kn)) &
6109 umveg=(one-vegfrac)*soilres
6117 r4=r3+hydro(1)*.5_kind_phys
6118 r5=r3-hydro(2)*.5_kind_phys
6127 totliq=prcp-drip/delt-(one-vegfrac)*dew*ras-smelt
6128 IF (debug_print )
THEN
6129 if (abs(xlat-testptlat).lt.0.05 .and. &
6130 abs(xlon-testptlon).lt.0.05)
then
6131 print *,
'xlat,xlon=',xlat,xlon
6132 print *,
'UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT', &
6133 umveg*prcp,drip/delt,umveg*dew*ras,smelt
6149 cvfrz = 3._kind_phys
6153 refdk=3.4341e-6_kind_phys
6154 delt1=delt/86400._kind_phys
6156 f2max=dqm*(zshalf(3)-zshalf(2))
6157 f1=f1max*(one-soilmois(1)/dqm)
6158 dice=soilice(1)*zshalf(2)
6161 dice=dice+(zshalf(k+1)-zshalf(k))*soilice(k)
6162 fkmax=dqm*(zshalf(k+1)-zshalf(k))
6163 fk=fkmax*(one-soilmois(k)/dqm)
6166 kdt=refkdt*ksat/refdk
6167 val=(1.-exp(-kdt*delt1))
6170 IF(px < zero) px = zero
6172 infmax1 = (px*(ddt/(px+ddt)))/delt
6176 IF (debug_print )
THEN
6177 print *,
'INFMAX1 before frozen part',infmax1
6185 frzx= 0.15_kind_phys*((dqm+qmin)/ref) * (0.412_kind_phys / 0.468_kind_phys)
6188 IF ( dice .GT. 1.e-2_kind_phys)
THEN
6189 acrt = cvfrz * frzx / dice
6197 sum = sum + (acrt ** ( cvfrz-jk)) / float(k)
6199 fcr = one - exp(-acrt) * sum
6201 IF (debug_print )
THEN
6202 print *,
'FCR--------',fcr
6203 print *,
'DICE=',dice
6205 infmax1 = infmax1* fcr
6208 infmax = max(infmax1,hydro(1)*soilmois(1))
6209 infmax = min(infmax, -totliq)
6210 IF (debug_print )
THEN
6211 print *,
'INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ', &
6212 infmax,infmax1,hydro(1)*soiliqw(1),-totliq
6215 IF (-totliq.GT.infmax)
THEN
6216 runoff=-totliq-infmax
6218 IF (debug_print )
THEN
6219 print *,
'FLX,RUNOFF1=',flx,runoff
6225 r7=.5_kind_phys*dzs/delt
6227 flx=flx-soilmois(1)*r7
6230 r8=umveg*r6*(one-snowfrac)
6236 IF(r10.LE.zero)
THEN
6237 qq=(r5*r2-flx+r9)/(r4-r5*r1-r10*r8/(ref-qmin))
6238 flxsat=-dqm*(r4-r5*r1-r10*r8/(ref-qmin)) &
6242 qq=(r2*r5-flx+r8*(qtot-qcg-qvg)+r9)/(r4-r1*r5)
6243 flxsat=-dqm*(r4-r1*r5)+r2*r5+r8*(qtot-qvg-qcg)+r9
6248 soilmois(1)=1.e-8_kind_phys
6250 ELSE IF(qq.GT.dqm)
THEN
6253 IF (debug_print )
THEN
6254 print *,
'FLXSAT,FLX,DELT',flxsat,flx,delt,runoff2
6256 runoff=runoff+(flxsat-flx)
6258 soilmois(1)=min(dqm,max(1.e-8_kind_phys,qq))
6261 IF (debug_print )
THEN
6262 if (abs(xlat-testptlat).lt.0.05 .and. &
6263 abs(xlon-testptlon).lt.0.05)
then
6264 print *,
'xlat,xlon=',xlat,xlon
6265 print *,
'SOILMOIS,SOILIQW, soilice',soilmois,soiliqw,soilice*riw
6266 print *,
'COSMC,RHSMC',cosmc,rhsmc
6273 qq=cosmc(kk)*soilmois(k-1)+rhsmc(kk)
6275 IF (qq.LT.zero)
THEN
6277 ELSE IF(qq.GT.dqm)
THEN
6281 IF (debug_print )
THEN
6282 print *,
'hydro(k),QQ,DQM,k',hydro(k),qq,dqm,k
6284 runoff2=runoff2+((qq-dqm)*(zsmain(k)-zshalf(k)))/delt
6286 runoff2=runoff2+((qq-dqm)*(zshalf(k+1)-zshalf(k)))/delt
6289 soilmois(k)=min(dqm,max(1.e-8_kind_phys,qq))
6292 IF (debug_print )
THEN
6293 if (abs(xlat-testptlat).lt.0.05 .and. &
6294 abs(xlon-testptlon).lt.0.05)
then
6295 print *,
'xlat,xlon=',xlat,xlon
6296 print *,
'END soilmois,soiliqw,soilice',soilmois,soiliqw,soilice*riw
6300 mavail=max(.00001_kind_phys,min(one,(soilmois(1)/(ref-qmin)*(one-snowfrac)+one*snowfrac)))
6762 soilfrac,nscat,shdmin, shdmax, &
6763 NLCAT,IVGTYP,ISLTYP,iswater,MYJ, &
6764 IFOREST,lufrac,vegfrac,EMISS,PC, &
6765 MSNF,FACSNF,ZNT,LAI,RDLAI2D, &
6766 QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,I,J, &
6791 integer,
parameter :: nsoilclas=19
6792 integer,
parameter :: nvegclas=24+3
6793 integer,
parameter :: ilsnow=99
6795 LOGICAL,
INTENT(IN ) :: debug_print
6796 INTEGER,
INTENT(IN ) :: mosaic_lu, mosaic_soil
6797 INTEGER,
INTENT(IN ) :: nlcat, nscat, iswater, i, j
6823 real (kind_phys) LQMA(nsoilclas),LRHC(nsoilclas), &
6824 LPSI(nsoilclas),LQMI(nsoilclas), &
6825 LBCL(nsoilclas),LKAS(nsoilclas), &
6826 LWIL(nsoilclas),LREF(nsoilclas), &
6836 DATA lqma /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
6837 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
6838 0.20, 0.435, 0.468, 0.200, 0.339/
6845 DATA lref /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
6846 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
6847 0.1, 0.249, 0.454, 0.17, 0.236/
6854 DATA lwil/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
6855 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
6856 0.006, 0.114, 0.030, 0.006, 0.01/
6862 DATA lqmi/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
6863 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
6864 0.004, 0.065, 0.020, 0.004, 0.008/
6872 DATA lpsi/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
6873 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
6874 0.121, 0.218, 0.468, 0.069, 0.069/
6882 DATA lkas/1.76e-4, 1.56e-4, 3.47e-5, 7.20e-6, 7.20e-6, &
6883 6.95e-6, 6.30e-6, 1.70e-6, 2.45e-6, 2.17e-6, &
6884 1.03e-6, 1.28e-6, 6.95e-6, 0.0, 1.41e-4, &
6885 3.47e-5, 1.28e-6, 1.41e-4, 1.76e-4/
6892 DATA lbcl/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
6893 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
6894 4.05, 4.90, 11.55, 2.79, 2.79/
6896 DATA lrhc /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
6897 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
6899 DATA datqtz/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
6900 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
6960 real (kind_phys) LALB(nvegclas),LMOI(nvegclas),LEMI(nvegclas), &
6961 LROU(nvegclas),LTHI(nvegclas),LSIG(nvegclas), &
6969 DATA lalb/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
6970 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
6972 DATA lemi/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
6973 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
6976 DATA lrou/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
6977 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
6980 DATA lmoi/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
6981 .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95,.40,.50,.40/
6985 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, &
6986 0.3,0.3,0.4,0.4,0.3,0.,.3,0.,0./
6994 LOGICAL,
INTENT(IN ) :: myj
6995 real (kind_phys),
INTENT(IN ) :: SHDMAX
6996 real (kind_phys),
INTENT(IN ) :: SHDMIN
6997 real (kind_phys),
INTENT(IN ) :: VEGFRAC
6998 real (kind_phys),
DIMENSION( 1:NLCAT ),
INTENT(IN):: LUFRAC
6999 real (kind_phys),
DIMENSION( 1:NSCAT ),
INTENT(IN):: SOILFRAC
7001 real (kind_phys) , &
7002 INTENT ( OUT) :: pc, &
7006 real (kind_phys) , &
7007 INTENT (INOUT ) :: emiss, &
7010 LOGICAL,
intent(in) :: rdlai2d
7012 real (kind_phys) , &
7013 INTENT( OUT) :: RHOCS, &
7022 INTEGER,
INTENT ( OUT) :: iforest
7023 character(len=*),
intent(out) :: errmsg
7024 integer,
intent(out) :: errflg
7025 INTEGER :: kstart, kfin, lstart, lfin
7027 real (kind_phys) :: area, factor, znt1, lb
7028 real (kind_phys),
DIMENSION( 1:NLCAT ) :: ZNTtoday, LAItoday, deltalai
7039 iforest = ifortbl(ivgtyp)
7041 IF (debug_print )
THEN
7042 print *,
'ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)', &
7043 ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)
7052 if((shdmax - shdmin) .lt. one)
then
7055 factor = one - max(zero,min(one,(vegfrac - shdmin)/max(one,(shdmax-shdmin))))
7060 if(ifortbl(k) == 1) deltalai(k)=min(0.2_kind_phys,0.8_kind_phys*laitbl(k))
7061 if(ifortbl(k) == 2 .or. ifortbl(k) == 7) deltalai(k)=min(0.5_kind_phys,0.8_kind_phys*laitbl(k))
7062 if(ifortbl(k) == 3) deltalai(k)=min(0.45_kind_phys,0.8_kind_phys*laitbl(k))
7063 if(ifortbl(k) == 4) deltalai(k)=min(0.75_kind_phys,0.8_kind_phys*laitbl(k))
7064 if(ifortbl(k) == 5) deltalai(k)=min(0.86_kind_phys,0.8_kind_phys*laitbl(k))
7066 if(k.ne.iswater)
then
7068 laitoday(k) = laitbl(k) - deltalai(k) * factor
7070 if(ifortbl(k) == 7)
then
7072 znttoday(k) = z0tbl(k) - 0.125_kind_phys * factor
7074 znttoday(k) = z0tbl(k)
7077 laitoday(k) = laitbl(k)
7082 IF (debug_print )
THEN
7083 print *,
'ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai,laitoday(ivgtyp),znttoday(ivgtyp)', &
7084 i,j,ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai(ivgtyp),laitoday(ivgtyp),znttoday(ivgtyp)
7093 if(.not.rdlai2d) lai = zero
7100 if(mosaic_lu == 1)
then
7102 area = area + lufrac(k)
7103 emiss = emiss+ lemitbl(k)*lufrac(k)
7104 znt = znt + lufrac(k)/alog(lb/znttoday(k))**2._kind_phys
7106 znt1 = znt1 + lufrac(k)*znttoday(k)
7107 if(.not.rdlai2d) lai = lai + laitoday(k)*lufrac(k)
7108 pc = pc + pctbl(k)*lufrac(k)
7109 msnf = msnf + mfsno(k)*lufrac(k)
7110 facsnf= facsnf + sncovfac(k)*lufrac(k)
7113 if (area.gt.one) area=one
7114 if (area <= zero)
then
7115 print *,
'Bad area of grid box', area
7117 errmsg =
'ERROR(SOILVEGIN): Bad area of grid box'
7121 IF (debug_print )
THEN
7122 print *,
'area=',area,i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),emiss,znt,znt1,lai,pc
7127 znt = lb/exp(sqrt(one/znt))
7128 if(.not.rdlai2d) lai = lai/area
7131 facsnf= facsnf /area
7133 IF (debug_print )
THEN
7134 print *,
'mosaic=',j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),emiss,znt,znt1,lai,pc
7139 emiss = lemitbl(ivgtyp)
7140 znt = znttoday(ivgtyp)
7142 msnf = mfsno(ivgtyp)
7143 facsnf= sncovfac(ivgtyp)
7144 if(.not.rdlai2d) lai = laitoday(ivgtyp)
7159 if(mosaic_soil == 1 )
then
7163 area = area + soilfrac(k)
7164 rhocs = rhocs + hc(k)*1.e6_kind_phys*soilfrac(k)
7165 bclh = bclh + bb(k)*soilfrac(k)
7166 dqm = dqm + (maxsmc(k)- &
7167 drysmc(k))*soilfrac(k)
7168 ksat = ksat + satdk(k)*soilfrac(k)
7169 psis = psis - satpsi(k)*soilfrac(k)
7170 qmin = qmin + drysmc(k)*soilfrac(k)
7171 ref = ref + refsmc(k)*soilfrac(k)
7172 wilt = wilt + wltsmc(k)*soilfrac(k)
7173 qwrtz = qwrtz + qtz(k)*soilfrac(k)
7176 if (area.gt.one) area=one
7177 if (area <= zero)
then
7180 rhocs = hc(isltyp)*1.e6_kind_phys
7182 dqm = maxsmc(isltyp)- &
7184 ksat = satdk(isltyp)
7185 psis = - satpsi(isltyp)
7186 qmin = drysmc(isltyp)
7187 ref = refsmc(isltyp)
7188 wilt = wltsmc(isltyp)
7204 if(isltyp.ne.14)
then
7205 rhocs = hc(isltyp)*1.e6_kind_phys
7207 dqm = maxsmc(isltyp)- &
7209 ksat = satdk(isltyp)
7210 psis = - satpsi(isltyp)
7211 qmin = drysmc(isltyp)
7212 ref = refsmc(isltyp)
7213 wilt = wltsmc(isltyp)