SUBROUTINE GBPHYS(IM,IX,levs,lsoil,lsm,ntrac,ncld,
& ntoz,ntcw,nmtvr,lonf,latg,jcap,ras,nlons,xkt2,nrcm,pre_rad,
& UGRS,VGRS,PGR,TGRS,QGRS,vvel,
& GT0,GQ0,GU0,GV0,sinlat,coslat,rcs2,
& prsi,prsl,prslk,prsik,phii,phil,dpshc,fhour,lssav,solhr,
! & prsi,prsl,prslk,prsik,phii,phil,prsshc,fhour,lssav,solhr,
& lsfwd,clstp,dtp,dtf,poz,prdout,ko3,pl_coeff,
& HICE,FICE,TISFC,SFCDSW, ! FOR SEA-ICE - XW Nov04
Clu [+2L]: add (tprcp,srflag),(slc,snwdph,slope,shdmin,shdmax,snoalb),sfalb
+ TPRCP, SRFLAG,
+ SLC ,SNWDPH,SLOPE ,SHDMIN,SHDMAX,SNOALB,SFALB ,
Cwei added 10/24/2006
+ CHH,CMM,EPI,DLWSFCI,ULWSFCI,USWSFCI,DSWSFCI,DTSFCI,
+ DQSFCI,GFLUXI,SRUNOFF,T1,Q1,U1,V1,ZLVL,EVBSA,EVCWA,
+ TRANSA,SBSNOA,SNOWCA,SOILM,SNOHFA,
& TSEA ,SHELEG,SNCOVR, TG3 ,ZORL ,CV ,CVB ,CVT ,
& SLMSK ,VFRAC ,CANOPY,F10M ,VTYPE ,STYPE ,UUSTAR,FFMM ,FFHH ,
& TMPMIN,TMPMAX,
!jwang add spfhmax/spfhmin
& SPFHMIN,SPFHMAX,
& GESHEM,DUSFC ,DVSFC ,DTSFC ,DQSFC ,DLWSFC,ULWSFC,
& GFLUX ,RUNOFF,EP ,CLDWRK,DUGWD ,DVGWD ,PSMEAN,BENGSH,XLON ,
& COSZEN,SFCNSW,XLAT ,
& SFCDLW,TSFLW ,PSURF ,U10M ,V10M ,T2M ,Q2M ,
! & COSZEN,SFCNSW,SFCDLW,TSFLW ,PSURF ,U10M ,V10M ,T2M ,Q2M ,
& HPBL ,PWAT ,SWH,HLW,SMC,STC,HPRIME,slag,sdec,cdec,
& acv,acvb,acvt,
& phy_f3d, phy_f2d, num_p3d, num_p2d, flgmin,
& DT3DT, DQ3DT, DU3DT, DV3DT, upd_mf, dwn_mf, det_mf, LDIAG3D,
! & DT3DT, DQ3DT, DU3DT, DV3DT, LDIAG3D,
!
! Coupling deletion->
! & flipv, me,kdt,lat,oro, crtrh, ncw, old_monin)
!<-Coupling deletion
! Coupling insertion->
& flipv, me,kdt,lat,oro, crtrh, ncw, old_monin,cnvgwd,ccwf,ctei_rm,
& sashal,newsas,mom4ice,mstrat,trans_trac,
> lssav_cc_dummy,DLWSFC_cc_dummy,ULWSFC_cc_dummy,SWSFC_cc_dummy,
> XMU_cc_dummy,
> DLW_cc_dummy,DSW_cc_dummy,SNW_cc_dummy,LPREC_cc_dummy, !cpl insertion
> DUSFC_cc_dummy,DVSFC_cc_dummy,DTSFC_cc_dummy,DQSFC_cc_dummy,
> PRECR_cc_dummy)
USE SURFACE_cc, ONLY: JCAL_cc
!<-Coupling insertion
!
USE MACHINE , ONLY : kind_phys
USE PHYSCONS, ROCP => con_rocp, CP => con_cp, FV => con_fvirt
&, grav => con_g, RD => con_RD
&, RV => con_RV, HVAP => con_HVAP
&, HFUS => con_HFUS
&, rerth => con_rerth, pi => con_pi
implicit none
!
integer levs,lsoil,lsm,ix,im,ntrac,ncld,ntoz,ntcw,nmtvr,lonf,
& latg,jcap,nlons(im),num_p3d,num_p2d,nrcm,lat
&, pl_coeff
!
real, parameter :: hocp=hvap/cp
LOGICAL lssav,lsfwd, old_monin, cnvgwd
integer levshc(im), levshcm ! Needed for pry version
integer ncw(2)
! real(kind=kind_phys) dtp,dtf,FHOUR,solhr, prsshc
real(kind=kind_phys) dtp,dtf,FHOUR,solhr, dpshc(im), crtrh(3)
real(kind=kind_phys) flgmin(2), ccwf, ctei_rm
real, parameter :: fhourpr=0.0
! Coupling insertion->
logical lssav_cc_dummy
real(kind=kind_phys),dimension(IM)::
> DLWSFC_cc_dummy,ULWSFC_cc_dummy,SWSFC_cc_dummy,XMU_cc_dummy,
> DLW_cc_dummy,DSW_cc_dummy,SNW_cc_dummy,LPREC_cc_dummy,
> DUSFC_cc_dummy,DVSFC_cc_dummy,DTSFC_cc_dummy,DQSFC_cc_dummy,
> PRECR_cc_dummy
real(kind=kind_phys),parameter:: CONVRAD_cc=JCAL_cc*1.E4/60.
!<-Coupling insertion
real(kind=kind_phys) UGRS(IX,LEVS), VGRS(IX,LEVS),
& TGRS(IX,LEVS), qgrs(IX,levs,ntrac),
& VVEL(IX,LEVS),
!
& GT0(IX,LEVS), GU0(IX,LEVS),
& GV0(IX,LEVS), gq0(IX,levs,ntrac),
!
& DEL(IX,LEVS), PRSI(IX,LEVS+1),
& PRSL(IX,LEVS), PRSLK(IX,LEVS),
& PRSIK(IX,LEVS+1), PHII(IX,LEVS+1),
& PHIL(IX,LEVS),
& PGR(IM), XKT2(IX,nrcm)
&, ccwfac(im)
real(kind=kind_phys) RCS2(IM), SINLAT(IM), COSLAT(IM),clstp
real(kind=kind_phys) SMC(IX,LSOIL), STC(IX,LSOIL), SWH(IX,LEVS)
&, HICE(IM), FICE(IM), SFCDSW(IM) ! SEA-ICE
&, TISFC(IM)
&, HLW(IX,LEVS), HPRIME(IX,NMTVR), TSEA(IM)
&, SHELEG(IM), TG3(IM), ZORL(IM)
&, SNCOVR(IM)
&, CV(IM), CVB(IM), CVT(IM)
&, COSZEN(IM), PWAT(IM), SLMSK(IM)
&, VFRAC(IM), CANOPY(IM), F10M(IM)
&, VTYPE(IM), STYPE(IM), UUSTAR(IM)
&, FFMM(IM), FFHH(IM), TMPMIN(IM)
&, TMPMAX(IM), GESHEM(IM), DUSFC(IM)
&, DVSFC(IM), DTSFC(IM), DQSFC(IM)
&, DLWSFC(IM), ULWSFC(IM), GFLUX(IM)
&, RUNOFF(IM), EP(IM), CLDWRK(IM)
&, DUGWD(IM), DVGWD(IM), PSMEAN(IM)
&, BENGSH(IM), XLON(IM), SFCNSW(IM)
&, SFCDLW(IM), TSFLW(IM), PSURF(IM)
&, U10M(IM), V10M(IM), T2M(IM)
&, Q2M(IM), HPBL(IM), xlat(IM)
!jwang add spfhmax/spfhmin
&, SPFHMIN(IM), SPFHMAX(IM)
Clu [+5L]: add (tprcp,srflag),(slc,snwdph,shdmin,shdmax,snoalb,slope),sfalb,(fm10,fh2)
&, TPRCP(IM), SRFLAG(IM)
&, SLC(IX,LSOIL)
&, SNWDPH(IM), SHDMIN(IM), SHDMAX(IM)
&, SNOALB(IM), SLOPE(IM), SFALB(IM)
Cwei added 10/24/2006
&, CHH(IM),CMM(IM),EPI(IM),DLWSFCI(IM),ULWSFCI(IM),USWSFCI(IM)
&, DSWSFCI(IM),DTSFCI(IM),DQSFCI(IM),GFLUXI(IM),SRUNOFF(IM),T1(IM)
&, Q1(IM),U1(IM),V1(IM),ZLVL(IM)
&, EVBSA(IM),EVCWA(IM),TRANSA(IM),SBSNOA(IM),SNOWCA(IM),SOILM(IM)
&, SNOHFA(IM)
!
&, phy_f3d(IX,LEVS,num_p3d), phy_f2d(IX,num_p2d)
&, acv(IM), acvb(IM), acvt(IM)
&, oro(im)
real(kind=kind_phys) slag,sdec,cdec
!
real(kind=kind_phys) dt3dt(IX,levs,6), dq3dt(IX,levs,5+pl_coeff),
& du3dt(IX,levs,4), dv3dt(IX,levs,4)
! &, cumflx(ix,levs)
&, upd_mf(ix,levs), dwn_mf(ix,levs)
&, det_mf(ix,levs)
!
integer me, kdt
logical RAS,LDIAG3D,pre_rad,sashal,newsas,mom4ice,mstrat
&, trans_trac
!
! In CLW, the first two varaibles are cloud water and ice.
! From third to ntrac are convective transportable tracers,
! third being the ozone, when ntrac=3 (valid only with RAS)
!
real(kind=kind_phys), allocatable :: CLW(:,:,:)
real(kind=kind_phys) garea(im), dlength(im)
! real(kind=kind_phys) CLW(IX,LEVS,ntrac)
! &, garea(im), dlength(im)
!
integer KO3
real(kind=kind_phys) poz(KO3), prdout(IX,ko3,pl_coeff)
!
real(kind=kind_phys) RHC(IM,LEVS), SR(IM,levs)
&, xncw(IM), rhbbot, rhbtop, rhpbl
&, dxmax, dxmin, dxinv
&, ud_mf(ix,levs), dd_mf(ix,levs)
&, dt_mf(ix,levs)
logical flipv
!
real(kind=kind_phys), parameter :: rhc_max=0.9999 ! 20060512
! real(kind=kind_phys), parameter :: rhc_max=0.999 ! for pry
!
real (kind=kind_phys), parameter :: qmin=1.0e-10, p850=85.0
!
!
! parameter (dxmax=log(1.0/7200.0), dxmin=log(1.0/192.0)
! parameter (dxmax=-8.8818363, dxmin=-5.2574954
! &, dxinv=1.0/(dxmax-dxmin))
! parameter (dxmax=ln(1.0/14000.0), dxmin=ln(1.0/192.0)
! parameter (dxmax=-9.5468126, dxmin=-5.2574954
! &, dxinv=1.0/(dxmax-dxmin))
!
! parameter (dxmax=ln(1.0/(14000.0*7000.0)), dxmin=ln(1.0/(192.0*94.0))
!
! parameter (dxmax=-18.40047804, dxmin=-9.800790154
! &, dxinv=1.0/(dxmax-dxmin))
!
!
! parameter (dxmax=ln(1.0/(4000.0*2000.0)), dxmin=ln(1.0/(192.0*94.0))
!
! parameter (dxmax=-15.8949521, dxmin=-9.800790154
! &, dxinv=1.0/(dxmax-dxmin))
!
! parameter (dxmax=ln(1.0/(3000.0*1500.0)), dxmin=ln(1.0/(192.0*94.0))
!
! parameter (dxmax=-15.31958795, dxmin=-9.800790154
! &, dxinv=1.0/(dxmax-dxmin))
!
! parameter (dxmax=ln(1.0/(2500.0*1250.0)), dxmin=ln(1.0/(192.0*94.0))
!
! parameter (dxmax=-14.95494484, dxmin=-9.800790154
! &, dxinv=1.0/(dxmax-dxmin))
!
! parameter (dxmax=ln(1.0/(2000.0*1000.0)), dxmin=ln(1.0/(192.0*94.0))
! parameter (dxmax=-14.50865774, dxmin=-9.800790154
! &, dxinv=1.0/(dxmax-dxmin))
!
! parameter (dxmax=ln(1.0/(5000.0*2500.0)), dxmin=ln(1.0/(192.0*94.0))
parameter (dxmax=-16.118095651, dxmin=-9.800790154
&, dxinv=1.0/(dxmax-dxmin))
!
real (kind=kind_phys), parameter :: cb2mb=10.0
c
c Local variables
c ---------------
real(kind=kind_phys) DTDT(IM,LEVS), DQDT(IM,LEVS,ntrac),
& DUDT(IM,LEVS), DVDT(IM,LEVS),
& GWDCU(IM,LEVS), GWDCV(IM,LEVS),
& DIAGN1(IM,LEVS), DIAGN2(IM,LEVS)
&, cuhr(im,levs)
real(kind=kind_phys) cumchr(IM,levs),cumabs(IM)
real(kind=kind_phys) qmax(IM)
!
Clu [-3L/+1L]: add slsoil; comment out ai, bi, cci, rhsmc, zsoil
Cwei [+1L]: uncomment and add slsoil because of adding OSU LSM
real(kind=kind_phys) SMSOIL(IM,LSOIL), STSOIL(IM,LSOIL),
& AI(IM,LSOIL), BI(IM,LSOIL),
& CCI(IM,LSOIL), RHSMC(IM,LSOIL),
& ZSOIL(IM,LSOIL),
+ SLSOIL(IM,LSOIL)
!c-- XW: FOR SEA-ICE Nov04
real(kind=kind_phys) CICE(IM), DSWSFC(IM), ZICE(IM)
&, TICE(IM)
!c-- XW: END SEA-ICE
!
real(kind=kind_phys) GFLX(IM), RAIN(IM), RAINC(IM),
& RAINL(IM), RAIN1(IM), EVAPC(IM),
& SNOWMT(IM), CD(IM), CDQ(IM),
& QSS(IM), radsl(IM), DUSFCG(IM),
& DVSFCG(IM), DUSFC1(IM), DVSFC1(IM),
& DTSFC1(IM), DQSFC1(IM), DLWSF1(IM),
& ULWSF1(IM), RB(IM), RHSCNPY(IM),
& DRAIN(IM), CLD1D(IM), RAINCS(IM),
& EVAP(IM), HFLX(IM), STRESS(IM),
! & EVAP(IM), HFLX(IM), RNET(IM),
& T850(IM),
& EP1D(IM), GAMT(IM), GAMQ(IM),
& sigmaf(IM),
& rcl(IM), rcs(IM),
& oc(IM), oa4(IM,4), clx(IM,4),
& theta(IM),gamma(IM),sigma(IM),elvmax(IM),
& wind(IM), work1(IM), work2(IM),
& runof(IM), xmu(IM)
&, qr_col(im,levs), fc_ice(im,levs)
&, oro_land(im)
Cwei added 10/24/2006
&, EVBS(im),EVCW(im),TRANS(im),SBSNO(im),SNOWC(im),SNOHF(im)
!
Clu [+1L]: add (fm10,fh2)
&, FM10(IM), FH2(IM)
Clu_q2m_iter [+1L]: add tsurf
&, tsurf(im)
Clu_q2m_iter [+1L]: add flag_iter and flag_guess
logical flag_iter(im), flag_guess(im)
Clu_q2m [+1L]: add wrk1, wrk2
real(kind=kind_phys) wrk1, wrk2
!
!
REAL(kind=kind_phys), PARAMETER :: EPSQ=1.E-20, HSUB=HVAP+HFUS
integer KBOT(IM), KTOP(IM), kcnv(im),
& soiltyp(IM), vegtype(IM), kpbl(IM)
Clu [+1L]: add slopetyp
&, SLOPETYP(IM)
integer i, nvdiff, kk, ic, k, n, kinver(im), ipr, lv, k1
&, lmh(im), tracers, trc_shft, tottracer
Clu_q2m_iter [+1L]: add iter
+, iter
!
logical lprnt, invrsn(im)
real(kind=kind_phys) frain, tem, tem1, qi, qw, qr, wc
&, f_rain, f_ice, tem0d, tx1(im)
&, tem2, ctei_r(im), flgmin_l(im)
! &, sumq, sumr, tem2
! &, pwato(im), raino(im), evapo(im), sume
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! lprnt = .true.
lprnt = .false.
! ipr = 1
! lprnt = kdt .gt. 0 .and. ilon .eq. 1
! do i=1,im
! work1(1) = xlon(i) * 57.29578
! if (work1(1) .ge. 180.0) work1(1) = work1(1) - 360.0
! work2(1) = xlat(i) * 57.29578
! print *,' me=',me,' work1=',work1(1),' work2=',work2(1),' i=',i
! lprnt = kdt .gt. 4320
! lprnt = kdt .gt. 0
! & .and. abs(work1(1)-103.0) .lt. 1.0
! & .and. abs(work2(1)-20.0) .lt. 0.5
! lprnt = kdt .ge. 14 .and. lat .eq. 43
! lprnt = kdt .ge. 1
! & .and. abs(xlon(i)*57.29578-114.325) .lt. 0.101
! & .and. abs(xlat(i)*57.29578+6.66) .lt. 0.101
! print *,' i=',i,' xlon=',xlon(i)*57.29578
! &, ' xlat=',xlat(i)*57.29578
! &,' i=',i,' me=',me
! if (lprnt) then
! ipr = i
! exit
! endif
! enddo
! lprnt = .false.
! if(lprnt) then
! print *,' IM=',IM,' IX=',IX,' levs=',levs,' lsoil=',lsoil
! *,' ntrac=',ntrac,' ntoz=',ntoz,' ntcw=',ntcw,' me=',me
! *,' xlat=',xlat(ipr),' kdt=',kdt,' slmsk=',slmsk(ipr)
! *,' nrcm=',nrcm
! &,' xlon=',xlon(ipr),' sfalb=',sfalb(ipr),' kdt=',kdt
! print *,' pgr=',pgr(ipr),' kdt=',kdt,' ipr=',ipr
! print *,' ipr=',ipr,' phy_f2d=',phy_f2d(ipr,1:num_p2d)
! print *,' ugrs=',ugrs(ipr,:)
! print *,' vgrs=',vgrs(ipr,:)
! print *,' tgrs=',tgrs(ipr,:),' kdt=',kdt,' ipr=',ipr
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! print *,' qgrs=',qgrs(ipr,:,1)
! print *,' ozg=',qgrs(ipr,:,2)
! print *,' clw=',qgrs(ipr,:,3)
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! endif
!
!
lmh(:) = levs
NVDIFF = ntrac ! vertical diffusion of all tracers!
!
! Figure out how many extra tracers are there
!
if (trans_trac) then
if (ntcw > 0) then
if (ntoz < ntcw) then
trc_shft = ntcw + ncld
else
trc_shft = ntoz
endif
elseif (ntoz > 0) then
trc_shft = ntoz
else
trc_shft = 1
endif
tracers = ntrac - trc_shft
if (ntoz > 0) tottracer = tracers + 1 ! Ozone is added separately
else
tottracer = 0 ! No convective transport of tracers
endif
allocate (clw(ix,levs,tottracer+2))
! if (lprnt) print *,' trans_trac=',trans_trac,' tottracer=',
! & tottracer,' trc_shft=',trc_shft,' kdt=',kdt
!
if (ras) then
if (ccwf >= 0.0) then
ccwfac = ccwf
else
ccwfac = -999.0
endif
endif
!
!
CALL GET_PRS(im,ix,levs,ntrac,TGRS,QGRS
&, prsi,prsik,prsl,prslk,phii,phil,del)
!..........................................................................
!
rhbbot = crtrh(1)
rhpbl = crtrh(2)
rhbtop = crtrh(3)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! do i=1,im
! PWATo(i) = 0.
! evapo(i) = dqsfc(i)
! raino(i) = geshem(i)
! enddo
! DO K=1,LEVS
! do i=1,im
! work1(i) = 0.0
! enddo
! if (ncld .gt. 0) then
! do ic=ntcw, ntcw+ncld-1
! do i=1,im
! work1(i) = work1(i) + qgrs(i,k,ic)
! enddo
! enddo
! endif
! do i=1,im
! PWATo(i) = PWATo(i) + DEL(i,K)*(qgrs(i,K,1)+work1(i))
! enddo
! ENDDO
! do i=1,im
! PWATo(i) = PWATo(i)*(1.E3/grav)
! enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
! print *,' prsi=',prsi(ipr,:)
! print *,' prsl=',prsl(ipr,:)
! print *,' del=',del(ipr,:)
! print *,' phii=',phii(ipr,:)
! print *,' phil=',phil(ipr,:)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! For pry version
do i=1,im
levshc(i) = 0
enddo
do k=2,levs
do i=1,im
if (prsi(i,1)-prsi(i,k) .le. dpshc(i)) levshc(i) = k
enddo
enddo
levshcm = 1
do i=1,im
levshcm = max(levshcm, levshc(i))
enddo
if (mstrat) then
levshcm = max(levshcm, levs/2)
else
levshcm = min(levshcm, levs/2)
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! FRAIN=FACTOR FOR CENTERED DIFFERENCE SCHEME CORRECTION OF RAIN AMOUNT.
!
FRAIN = DTF/DTP
do i=1,im
SOILTYP(i) = INT(STYPE(i)+.5)
SIGMAF(i) = MAX(VFRAC(i),0.01)
! SIGMAF(i) = MAX(VFRAC(i),.3)
if (lsm == 0) SIGMAF(i) = 0.5 + VFRAC(i) * 0.5
VEGTYPE(i) = INT(VTYPE(i)+.5)
SLOPETYP(i) = INT(SLOPE(i)+.5) !! Clu [+1L]: slope -> slopetyp
! fm1(i) = ffmm(i)
! fh1(i) = ffhh(i)
! ustar1(i) = uustar(i)
IF(SLMSK(i).EQ.2.) THEN
SOILTYP(i) = 9
VEGTYPE(i) = 13
SLOPETYP(i) = 9 !! Clu [+1L]: QA(slopetyp)
ENDIF
enddo
!sela if(vegtype.eq.0 )print*,' vegtyp ilon ilat =',vegtype, ilon,ilat
!
! TRANSFER SOIL MOISTURE AND TEMPERATURE FROM GLOBAL TO LOCAL VARIABLES
!
DO K = 1, LSOIL
do i=1,im
SMSOIL(i,K) = SMC(i,k)
STSOIL(i,K) = STC(i,k)
SLSOIL(i,K) = SLC(i,k) !! Clu [+1L]: slc -> slsoil
enddo
ENDDO
!
!c-- XW: FOR SEA-ICE Nov04
!
! TRANSFER ICE THICKNESS & CONCENTRATION FROM GLOBAL TO LOCAL VARIABLES
!
do i=1,im
ZICE(i) = HICE(i)
CICE(i) = FICE(i)
TICE(i) = TISFC(i)
enddo
!c-- XW: END SEA-ICE
!
do k=1,levs
do i=1,im
dudt(i,k) = 0.
dvdt(i,k) = 0.
dtdt(i,k) = 0.
enddo
enddo
do n=1,ntrac
do k=1,levs
do i=1,im
dqdt(i,k,n) = 0. ! dqdt may be dimensioned (levs,ntrac)
enddo
enddo
enddo
do i=1,im
RCL(i) = RCS2(i)
RCS(i) = SQRT(RCL(i))
psurf(i) = pgr(i)
work1(i) = prsik(i,1) / prslk(i,1)
! garea(i) = rerth * rerth * (pi+pi)*pi*coslat(i)/(nlons(i)*latg)
tem1 = rerth * (pi+pi)*coslat(i)/nlons(i)
tem2 = rerth * pi/latg
garea(i) = tem1 * tem2
dlength(i) = sqrt(tem1*tem1+tem2*tem2)
enddo
if(lssav) then
do i=1,im
PSMEAN(i) = PSMEAN(i) + PGR(i)*DTF
enddo
endif
!
! INITIALIZE DTDT WITH HEATING RATE FROM DCYC2
!
! if(lprnt) then
! do ipr=1,im
! print *,' before DCYC2: IM=',IM,' LSOIL=',LSOIL,' levs=',levs
! &,' sde=',sdec,' cdec=',cdec,' tsea=',tsea(ipr),' ipr=',ipr
! &,' lat=',lat,' me=',me,' kdt=',kdt
! &,' sfcdlw=',sfcdlw(ipr),' sfcnsw=',sfcnsw(ipr)
! print *,' hlw=',hlw(ipr,:),' me=',me,' lat=',lat,xlon(ipr)
! print *,' swh=',swh(ipr,:),' me=',me,' lat=',lat,xlon(ipr)
! enddo
! endif
!
!
if(pre_rad)then
CALL DCYC2t3_pre_rad(IX,IM,LEVS,SOLHR,SLAG,
& SINLAT,COSLAT,SDEC,CDEC,
& XLON,COSZEN,SFCDLW,SFCNSW,TGRS(1,1),
& SFCDSW,DSWSFC, ! FOR SEA-ICE - XW Nov04
& TSEA,TGRS(1,1),SWH,HLW,
& DLWSF1,ULWSF1,radsl,DTDT,xmu)
else
CALL DCYC2t3(IX,IM,LEVS,SOLHR,SLAG,
& SINLAT,COSLAT,SDEC,CDEC,
& XLON,COSZEN,SFCDLW,SFCNSW,TGRS(1,1),
& SFCDSW,DSWSFC, ! FOR SEA-ICE - XW Nov04
& TSEA,TSFLW,SWH,HLW,
& DLWSF1,ULWSF1,radsl,DTDT,xmu)
endif
!-> Coupling insertion
if(lssav_cc_dummy) then
XMU_cc_dummy(1:IM)=MAX(0.,MIN(1.0,XMU(1:IM)*COSZEN(1:IM)))
DSW_cc_dummy(1:IM)=DSWSFC(1:IM)
DLW_cc_dummy(1:IM)=DLWSF1(1:IM)
DLWSFC_cc_dummy(1:IM)=DLWSFC_cc_dummy(1:IM)+DLWSF1(1:IM) !*DTF (?)
ULWSFC_cc_dummy(1:IM)=ULWSFC_cc_dummy(1:IM)+ULWSF1(1:IM) !*DTF (?)
SWSFC_cc_dummy(1:IM)= !*DTF (?)
! > SWSFC_cc_dummy(1:IM)+radsl(1:IM)*CONVRAD_cc+DLWSF1(1:IM)
! <- see progtmr.f, subr. progtm (called below), and
! SUBROUTINE DCYC2T3 just called
! ** check signs here, esp. for SWR.
! - there are many pecularities here, the whole rad.
! stuff must be kept an eye on
! Signs for SWSFC_cc: sign for SWSFC_cc is intended positive
! upward, as for other fluxes. Downward LWR needs to be subtracted
! from radsl*const, where the downward LWR is considered positive
! upward. Since DLWSF1 is downward LWR positive downward,
! DLWSF1 must be added rather than subtracted. Similarly (see
! progtmr.f), radsl must be taken with +.
> SWSFC_cc_dummy(1:IM)+radsl(1:IM)+DLWSF1(1:IM)
end if
!<- Coupling insertion
if(lssav)then
do i=1,im
DLWSFC(i) = DLWSFC(i) + DLWSF1(i)*DTF
ULWSFC(i) = ULWSFC(i) + ULWSF1(i)*DTF
enddo
IF (LDIAG3D) THEN
DO K=1,LEVS
do i=1,im
DT3DT(i,k,1) = DT3DT(i,k,1) + HLW(i,K)*DTF
DT3DT(i,k,2) = DT3DT(i,k,2) + SWH(i,K)*DTF*xmu(i)
enddo
ENDDO
ENDIF
endif
!
!
do i=1,im
kcnv(i) = 0
kinver(i) = levs
invrsn(i) = .false.
tx1(i) = 0.0
ctei_r(i) = 10.0
enddo
! ctei_rm = 0.20
! ctei_rm = 0.25
! ctei_rm = 0.50
! ctei_rm = 0.70
do k=1,levs/2
do i=1,im
if (prsi(i,1)-prsi(i,k+1) .lt. 0.35*prsi(i,1)
& .and. (.not. invrsn(i))) then
tem = (TGRS(i,K+1) - TGRS(i,K)) / (prsl(i,k) - prsl(i,k+1))
! if (tem .gt. 0.02 .and. tx1(i) .lt. 0.0) then
! if (tem .gt. 0.10 .and. tx1(i) .lt. 0.0) then
! if (tem .gt. 0.05 .and. tx1(i) .lt. 0.0) then
if (tem .gt. 0.025 .and. tx1(i) .lt. 0.0) then
invrsn(i) = .true.
!
if (qgrs(i,k,1) .ne. qgrs(i,k+1,1) .and.
& .not. sashal) then
tem1 = (1.0 + hocp*max(qgrs(i,k+1,1),qmin)/tgrs(i,k+1))
tem2 = (1.0 + hocp*max(qgrs(i,k,1),qmin)/tgrs(i,k))
tem1 = tem1 * tgrs(i,k+1) / prslk(i,k+1)
& - tem2 * tgrs(i,k) / prslk(i,k)
! (Cp/L)(delthetae)/(deltotwater) < 0.7
ctei_r(i) = (1.0/hocp)*tem1/(qgrs(i,k+1,1)-qgrs(i,k,1)
& + qgrs(i,k+1,ntcw)-qgrs(i,k,ntcw))
else
ctei_r(i) = 10
endif
kinver(i) = k
endif
tx1(i) = tem
endif
enddo
enddo
! ipr = 1
! if(lprnt) then
! print *,' before PROGTM: IM=',IM,' LSOIL=',LSOIL
! &,' nvdiff=',nvdiff,' radsl=',radsl(ipr),' dlwsf1=',dlwsf1(ipr)
! &,' dlwsf1=',dlwsf1(ipr),' tsea2=',tsea(ipr)
! &,' ipr=',ipr,' me=',me,' lat=',lat,' xlon=',xlon(ipr)
! &,' kdt=',kdt
! print *,' dtdth=',dtdt(ipr,:),' kdt=',kdt
! endif
! phy_f2d(:,num_p2d) = 0.0 ! disable downdraft effect on evap
!
Clu [-11L/+34L]: Revision starts here ...................................
Clu CALL PROGTM(IM,LSOIL,PGR,UGRS,VGRS,TGRS,QGRS,
Clu & SHELEG,TSEA,QSS,
Clu & SMSOIL,STSOIL,EVAPC,SOILTYP,SIGMAF,
Clu & vegtype,CANOPY,
Clu & dlwsf1,
Clu & radsl,SNOWMT,dtf ,ZORL,TG3,
Clu & GFLX,F10M,U10M,V10M,T2M,Q2M,ZSOIL,
Clu & CD,CDQ,RB,RHSCNPY,RHSMC,AI,BI,CCI,
Clu & RCL,PRSL(1,1),work1,SLMSK,
Clu & DRAIN,EVAP,HFLX,STRESS,EP1D,ffmm,ffhh,
Clu & uustar,WIND,phy_f2d(1,num_p2d))
Clu_q2m [+22L]: print selected variables
! do i=1,im
! wrk1 = xlon(i) * 57.29578
! if (wrk1 .ge. 180.0) wrk1 = wrk1 - 360.0
! wrk2 = xlat(i) * 57.29578
! lprnt = abs(wrk1+96.5) .lt. 0.5
! & .and. abs(wrk2-39.1) .lt. 0.5
! & .and. me .eq. 12
! if(lprnt) then
! write(221) LSOIL,PGR(i),UGRS(i,1),VGRS(i,1),
! & TGRS(i,1),QGRS(i,1,1),
! & SHELEG(i),TSEA(i),QSS(i),
! & (SMSOIL(i,k),k=1,LSOIL),(STSOIL(i,k),k=1,LSOIL),
! & (SLSOIL(i,k),k=1,LSOIL),EVAPC(i),SOILTYP(i),SIGMAF(I),
! & vegtype(i),CANOPY(i),dlwsf1(i),
! & radsl(i),dtf ,ZORL(i),
! & TG3(i),GFLX(i),U10M(i),
! & V10M(i),T2M(i),Q2M(i),CD(i),CDQ(i),
! & RCL(i),PRSL(1,1),work1(i),SLMSK(i),
! + DRAIN(i),EVAP(i),HFLX(i),EP1D(i),ffmm(i),ffhh(i),
! & uustar(i),WIND(i)
! endif
! enddo
Clu_q2m_iter [+5L]: initialize flag_guess, flag_iter, tsurf
do i=1, im
tsurf(i) = tsea(i)
flag_guess(i) = .False.
flag_iter(i) = .True.
drain(i) = 0.0
ep1d(i) = 0.0
runof(i) = 0.0
hflx(i) = 0.0
evap(i) = 0.0
!
evbs(i) = 0.0
evcw(i) = 0.0
trans(i) = 0.0
sbsno(i) = 0.0
snowc(i) = 0.0
snohf(i) = 0.0
enddo
Clu_q2m_iter [+2L]: add iter-loop over (sfc_diff,sfc_drv,sfc_ocean,sfc_sice)
do iter = 1, 2 !!!!! <---- Clu_q2m_iter
!
!** surface exchange coefficients
!
CALL SFC_DIFF(IM,PGR,UGRS,VGRS,TGRS,QGRS,
& TSEA,ZORL,CD,CDQ,RB,
& RCL,PRSL(1,1),work1,SLMSK,
& STRESS,ffmm,ffhh,
Clu_q2m_iter [-1L/+2L]: add tsurf, flag_iter
!* & uustar,WIND,phy_f2d(1,num_p2d),fm10,fh2)
+ uustar,WIND,phy_f2d(1,num_p2d),fm10,fh2,
+ tsurf, flag_iter)
! if (lprnt) print *,' cdq=',cdq(ipr),' iter=',iter
! &,' wind=',wind(ipr),'phy_f2d=',phy_f2d(ipr,num_p2d),' ugrs='
! &,ugrs(ipr,1),' vgrs=',vgrs(ipr,1)
Clu_q2m_iter [+4L]: update flag_guess
do i=1, im
if((iter.eq.1) .and. (wind(i).lt.2.))
+ flag_guess(i) = .True.
enddo
!
!** surface energy balance over ocean
!
CALL SFC_OCEAN(IM,LSOIL,PGR,UGRS,VGRS,TGRS,QGRS,
& TSEA,QSS,EVAPC,GFLX,CD,CDQ,
& RCL,PRSL(1,1),work1,SLMSK,
Cwei added 10/24/2006
& CMM,CHH,
Clu_q2m_iter [-1L/+1L]: add flag_iter
!* & EVAP,HFLX,EP1D,phy_f2d(1,num_p2d))
+ EVAP,HFLX,EP1D,phy_f2d(1,num_p2d),flag_iter)
! if (lprnt) print *,' sfalb=',sfalb(ipr),' ipr=',ipr
! &,' sheleg=',sheleg(ipr),' snwdph=',snwdph(ipr)
! &,' tprcp=',tprcp(ipr),' kdt=',kdt,' iter=',iter
!
!** surface energy balance over land
!
if (lsm == 1) then ! NOAH LSM CALL
CALL SFC_DRV(IM,LSOIL,PGR,UGRS,VGRS,TGRS,QGRS,
& SHELEG,SNCOVR,SNWDPH,TSEA,QSS,TPRCP,SRFLAG,
& SMSOIL,STSOIL,SLSOIL,EVAPC,SOILTYP,SIGMAF,
& vegtype,CANOPY,dlwsf1,DSWSFC,
& radsl,dtf,TG3,GFLX,CD,CDQ,
& RCL,PRSL(1,1),work1,SLMSK,
& DRAIN,EVAP,HFLX,EP1D,phy_f2d(1,num_p2d),
Clu_q2m_iter [-1L/+2L]: add tsurf, flag_iter, flag_guess
!* + RUNOF,SLOPETYP,SHDMIN,SHDMAX,SNOALB,SFALB)
+ RUNOF,SLOPETYP,SHDMIN,SHDMAX,SNOALB,SFALB,
Cwei added 10/24/2006
+ CMM,CHH,ZLVL,EVBS,EVCW,TRANS,SBSNO,
+ SNOWC,SOILM,SNOHF,
+ tsurf, flag_iter, flag_guess)
else ! OSU LSM CALL
CALL SFC_LAND(IM,LSOIL,PGR,UGRS,VGRS,TGRS,QGRS,
& SHELEG,TSEA,QSS,TPRCP,SRFLAG,
& SMSOIL,STSOIL,EVAPC,SOILTYP,SIGMAF,
& vegtype,CANOPY,dlwsf1,
& radsl,SNOWMT,dtf,ZORL,TG3,
& GFLX,ZSOIL,
& CD,CDQ,RHSCNPY,RHSMC,AI,BI,CCI,
& RCL,PRSL(1,1),work1,SLMSK,
& DRAIN,EVAP,HFLX,EP1D,phy_f2d(1,num_p2d),
Cwei added 10/24/2006
+ CMM,CHH,ZLVL,EVBS,EVCW,TRANS,SBSNO,
+ SNOWC,SOILM,SNOHF,
& tsurf,flag_iter, flag_guess)
endif
!
! if (lprnt) print *,' tseabeficemodel =',tsea(ipr),' me=',me
! &,' kdt=',kdt
!
!** surface energy balance over seaice
!
!<-- cpl insertion
CALL SFC_SICE(IM,LSOIL,PGR,UGRS,VGRS,TGRS,QGRS,
& ZICE,CICE,TICE,DSWSFC, ! FOR SEA-ICE - XW Nov04
& SHELEG,SNWDPH,TSEA,QSS,TPRCP,SRFLAG,STSOIL,EVAPC,
& dlwsf1,radsl,SNOWMT,dtf,GFLX,CD,CDQ,
& RCL,PRSL(1,1),work1,SLMSK,
Cwei added 10/24/2006
+ CMM,CHH,ZLVL,
Clu_q2m_iter [-1L/+1L]: add flag_iter
!* & EVAP,HFLX,EP1D,phy_f2d(1,num_p2d))
+ EVAP,HFLX,EP1D,phy_f2d(1,num_p2d),flag_iter,
& mom4ice,lsm)
Clu_q2m_iter [+7L]: update flag_iter and flag_guess
do i=1, im
flag_iter(i) = .False.
flag_guess(i) = .False.
if((slmsk(i) .eq. 1.) .and. (iter .eq. 1)) then
if(wind(i).lt.2.) flag_iter(i) = .True.
endif
enddo
enddo !!!!! <---- Clu_q2m_iter
Cwei added 10/24/2006
do i=1,im
epi(i) = ep1d(i)
dlwsfci(i) = dlwsf1(i)
ulwsfci(i) = ulwsf1(i)
uswsfci(i) = sfcnsw(i)*xmu(i) + dswsfc(i)
dswsfci(i) = dswsfc(i)
gfluxi(i) = gflx(i)
t1(i) = tgrs(i,1)
q1(i) = qgrs(i,1,1)
u1(i) = ugrs(i,1)
v1(i) = vgrs(i,1)
enddo
if (lsm == 0) then ! OSU LSM CALL
do i=1,im
SNCOVR(i) = 0.0
if (SHELEG(i) > 0.0) SNCOVR(i) = 1.0
enddo
endif
!
!** update near surface fields
!
CALL SFC_DIAG(IM,LSOIL,PGR,UGRS,VGRS,TGRS,QGRS,
& TSEA,QSS,F10M,U10M,V10M,T2M,Q2M,RCL,work1,SLMSK,
& EVAP,ffmm,ffhh,fm10,fh2)
!
Clu_q2m [+22L]: print selected variables
! do i=1,im
! wrk1 = xlon(i) * 57.29578
! if (wrk1 .ge. 180.0) wrk1 = wrk1 - 360.0
! wrk2 = xlat(i) * 57.29578
! lprnt = abs(wrk1+96.5) .lt. 0.5
! & .and. abs(wrk2-39.1) .lt. 0.5
! & .and. me .eq. 12
! if(lprnt) then
! write(222) LSOIL,PGR(i),UGRS(i,1),VGRS(i,1),
! & TGRS(i,1),QGRS(i,1,1),
! & SHELEG(i),TSEA(i),QSS(i),
! & (SMSOIL(i,k),k=1,LSOIL),(STSOIL(i,k),k=1,LSOIL),
! & (SLSOIL(i,k),k=1,LSOIL),EVAPC(i),SOILTYP(i),SIGMAF(I),
! & vegtype(i),CANOPY(i),dlwsf1(i),
! & radsl(i),dtf ,ZORL(i),
! & TG3(i),GFLX(i),U10M(i),
! & V10M(i),T2M(i),Q2M(i),CD(i),CDQ(i),
! & RCL(i),PRSL(1,1),work1(i),SLMSK(i),
! + DRAIN(i),EVAP(i),HFLX(i),EP1D(i),ffmm(i),ffhh(i),
! & uustar(i),WIND(i)
! endif
! enddo
Clu [-11L/+34L]: Revision ends here .....................................
!
! if(lprnt) then
! print *,' hflx=',hflx(ipr)
! print *,' evap=',evap(ipr)
! print *,' stress=',stress(ipr)
! endif
!
do i=1,im
phy_f2d(i,num_p2d) = 0.0
enddo
! if (lprnt) print *,' tseaim=',tsea(ipr),' me=',me,' kdt=',kdt
!
if(lssav)then
do i=1,im
GFLUX(i) = GFLUX(i) + GFLX(i)*DTF
Cwei added 10/24/2006
EVBSA(I) = EVBSA(I) + EVBS(i)*DTF
EVCWA(I) = EVCWA(I) + EVCW(i)*DTF
TRANSA(I) = TRANSA(I) + TRANS(i)*DTF
SBSNOA(I) = SBSNOA(I) + SBSNO(i)*DTF
SNOWCA(I) = SNOWCA(I) + SNOWC(i)*DTF
SNOHFA(I) = SNOHFA(I) + SNOHF(i)*DTF
TMPMAX(i) = MAX(TMPMAX(i),T2M(i))
TMPMIN(i) = MIN(TMPMIN(i),T2M(i))
!jwang [+2L] add SPFHMIN & SPFHMAX
SPFHMAX(i) = MAX(SPFHMAX(i),Q2M(i))
SPFHMIN(i) = MIN(SPFHMIN(i),Q2M(i))
!
EP(i) = EP(i) + EP1D(i) * DTF
CWei move to the place after calling progt2
Clu [+3L]: compute total runoff
! TOTAL RUNOFF IS COMPOSED OF DRAINAGE INTO WATER TABLE AND
! RUNOFF AT THE SURFACE AND IS ACCUMULATED IN UNIT OF METERS
! RUNOFF(i) = RUNOFF(i) + (DRAIN(i)+RUNOF(i)) * DTF / 1000.0
! uustar(i) = ustar1(i) ! should be under lsfwd
! ffmm(i) = fm1(i) ! should be under lsfwd
! ffhh(i) = fh1(i) ! should be under lsfwd
enddo
endif
csela the following is the correct code, the above duplicates oper. code.
csela if(lsfwd)then
csela uustar(ilon,ilat)=ustar1
csela ffmm(ilon,ilat)=fm1
csela ffhh(ilon,ilat)=fh1
csela endif
!
Clu [+11L]: update smc, stc, slc
!
! RETURN UPDATED SMSOIL AND STSOIL TO GLOBAL ARRAYS
!
!Wei move to the place after calling progt2
! DO K = 1, LSOIL
! do i=1,im
! SMC(i,k) = SMSOIL(i,K)
! STC(i,k) = STSOIL(i,K)
! SLC(i,k) = SLSOIL(i,K)
! enddo
! ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Commented by Moorthi on 20050505
!
do i=1,im
!
! COMPUTE COEFFICIENT OF EVAPORATION IN EVAPC
!
IF (EVAPC(i) .GT. 1.0E0) EVAPC(i) = 1.0E0
!
! OVER SNOW COVER OR ICE OR SEA, COEF OF EVAP =1.0E0
!
IF((SHELEG(i).GT.0.) .OR. (SLMSK(i).NE.1.0)) EVAPC(i) = 1.0E0
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! VERTICAL DIFFUSION
!
! if (lprnt) print *,' tsea3=',tsea(ipr),' slmsk=',slmsk(ipr)
! &,' kdt=',kdt,' evap=',evap(ipr)
! if (lprnt) print *,' dtdtb=',dtdt(ipr,:)
!
do i=1,im
if (slmsk(i) .eq. 0) then
oro_land(i) = 0.0
else
oro_land(i) = oro(i)
endif
enddo
if (OLD_MONIN) then
if (mstrat) then
CALL MONINP1(IX,IM,LEVS,nvdiff,DVDT,DUDT,DTDT,DQDT,
& UGRS,VGRS,TGRS,QGRS,
& prsik(1,1),RB,ffmm,ffhh,TSEA,QSS,HFLX,EVAP,STRESS,WIND,KPBL,
& PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,dtp,
& DUSFC1,DVSFC1,DTSFC1,DQSFC1,HPBL,GAMT,GAMQ,
& kinver, oro_land)
else
CALL MONINP(IX,IM,LEVS,nvdiff,DVDT,DUDT,DTDT,DQDT,
& UGRS,VGRS,TGRS,QGRS,
& prsik(1,1),RB,ffmm,ffhh,TSEA,QSS,HFLX,EVAP,STRESS,WIND,KPBL,
& PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,dtp,
& DUSFC1,DVSFC1,DTSFC1,DQSFC1,HPBL,GAMT,GAMQ)
endif
ELSE
if (mstrat) then
CALL MONINQ1(IX,IM,LEVS,nvdiff,DVDT,DUDT,DTDT,DQDT,
& UGRS,VGRS,TGRS,QGRS,swh,hlw,xmu,
& prsik(1,1),RB,ffmm,ffhh,TSEA,QSS,HFLX,EVAP,STRESS,WIND,KPBL,
& PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,dtp,
& DUSFC1,DVSFC1,DTSFC1,DQSFC1,HPBL,GAMT,GAMQ,
& kinver)
else
CALL MONINQ(IX,IM,LEVS,nvdiff,DVDT,DUDT,DTDT,DQDT,
& UGRS,VGRS,TGRS,QGRS,swh,hlw,xmu,
& prsik(1,1),RB,ffmm,ffhh,TSEA,QSS,HFLX,EVAP,STRESS,WIND,KPBL,
& PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,dtp,
& DUSFC1,DVSFC1,DTSFC1,DQSFC1,HPBL,GAMT,GAMQ)
endif
ENDIF
!
! if (ntoz .gt. 0) dqdt(:,:,ntoz) = 0.0
!
! if (lprnt) then
! print *,' dusfc1=',dusfc1(ipr)
! print *,' dtsfc1=',dtsfc1(ipr)
! print *,' dqsfc1=',dqsfc1(ipr)
! print *,' dtdt=',dtdt(ipr,:)
! print *,' dudtm=',dudt(ipr,:)
! endif
!-> Coupling insertion
if(lssav_cc_dummy) then
DUSFC_cc_dummy(1:IM)=DUSFC_cc_dummy(1:IM)+DUSFC1 !*DTF <-na h.(?)
DVSFC_cc_dummy(1:IM)=DVSFC_cc_dummy(1:IM)+DVSFC1 !*DTF <-na h.(?)
DTSFC_cc_dummy(1:IM)=DTSFC_cc_dummy(1:IM)+DTSFC1 !*DTF <-na h.(?)
DQSFC_cc_dummy(1:IM)=DQSFC_cc_dummy(1:IM)+DQSFC1 !*DTF <-na h.(?)
end if
!<- Coupling insertion
if(lssav) then
do i=1,im
DUSFC(i) = DUSFC(i) + DUSFC1(i)*DTF
DVSFC(i) = DVSFC(i) + DVSFC1(i)*DTF
DTSFC(i) = DTSFC(i) + DTSFC1(i)*DTF
DQSFC(i) = DQSFC(i) + DQSFC1(i)*DTF
dtsfci(i) = dtsfc1(i)
dqsfci(i) = dqsfc1(i)
enddo
IF (LDIAG3D) THEN
DO K=1,LEVS
do i=1,im
tem1 = rcs(i) * dtf
tem = dtdt(i,k) - (hlw(i,k)+swh(i,k)*xmu(i))
DT3DT(i,k,3) = DT3DT(i,k,3) + tem*DTF
DQ3DT(i,k,1) = DQ3DT(i,k,1) + DQDT(i,K,1)*DTF
DU3DT(i,k,1) = DU3DT(i,k,1) + DUDT(i,K)*tem1
DU3DT(i,k,2) = DU3DT(i,k,2) - DUDT(i,K)*tem1
DV3DT(i,k,1) = DV3DT(i,k,1) + DVDT(i,K)*tem1
DV3DT(i,k,2) = DV3DT(i,k,2) - DVDT(i,K)*tem1
enddo
ENDDO
if (ntoz .gt. 0) then
DO K=1,LEVS
do i=1,im
DQ3DT(i,k,5) = DQ3DT(i,k,5) + DQDT(i,K,ntoz)*DTF
enddo
enddo
endif
ENDIF
endif
!
if (NMTVR .eq. 6) then
do i=1,im
oc(i) = hprime(i,2)
enddo
do k = 1, 4
do i=1,im
oa4(i,k) = hprime(i,k+2)
clx(i,k) = 0.0
enddo
enddo
elseif(NMTVR .eq. 10) then
do i=1,im
oc(i) = hprime(i,2)
enddo
do k = 1, 4
do i=1,im
oa4(i,k) = hprime(i,k+2)
clx(i,k) = hprime(i,k+6)
enddo
enddo
! --- lm mb (*j*)
elseif(NMTVR .eq. 14) then
! --- get the kim fields (until this is changed)
do i=1,im
oc(i) = hprime(i,2)
enddo
do k = 1, 4
do i=1,im
oa4(i,k) = hprime(i,k+2)
clx(i,k) = hprime(i,k+6)
enddo
enddo
do i=1,im
theta(i) = hprime(i,11)
gamma(i) = hprime(i,12)
sigma(i) = hprime(i,13)
elvmax(i) = hprime(i,14)
enddo
else
oc = 0
oa4 = 0
clx = 0
theta = 0
gamma = 0
sigma = 0
elvmax = 0
endif
!sela also replace lonf2 with cleff in this routine and compute
!sela cleff for all columns once before the loops on lat and lon
!sela do this in gloopb once and for all, better yet, in step1
!
!
! Call (old) operational gravity-wave drag
!
! CALL GWDPS(lonf2,LEVS,DVDT,DUDT,UGRS,VGRS,TGRS,QGRS,
! & PGR,SI,DEL,CL,SL,slk,RCL,dtp,LAT,HPRIME(ilon,ilat,1),
! & oc,oa4,DUSFCG,DVSFCG)
!
! if(lssav) then
! DUGWD(ilon,ilat)=DUGWD(ilon,ilat)+DUSFCG*DTF
! DVGWD(ilon,ilat)=DVGWD(ilon,ilat)+DVSFCG*DTF
! IF (LDIAG3D) THEN
! work1 = rcs * dtf
! DO K=1,LEVS
! DU3DT(k,2,ilon,ilat) = DU3DT(k,2,ilon,ilat) + DUDT(K)*work1
! DV3DT(k,2,ilon,ilat) = DV3DT(k,2,ilon,ilat) + DVDT(K)*work1
! ENDDO
! ENDIF
! endif
!#ifdef NEWGWD
!
! Call New operational gravity-wave drag
!
! if(lprnt) then
! print *,' kdt=',kdt,' rcl=',rcl(ipr),' dtp=',dtp,
! &' pgr=',pgr(ipr),' kpbl=',kpbl,grav, CP, RD, RV, LONF
! print *,' ugrs=',ugrs
! print *,' vgrs=',vgrs
! print *,' tgrs=',tgrs(ipr,:)
! print *,' qgrs=',qgrs
! print *,' hprime=',hprime(ipr,:)
! print *,' oa4=',oa4,' oc=',oc,' clx=',clx
! print *,' si=',si
! print *,' prsl=',prsl(1,:),' me=',me
! print *,' del=',del(1,:),' me=',me
! print *,' prslk=',prslk(1,:),' me=',me
! print *,' prsik=',prsik(1,:),' me=',me
! endif
!
CALL GWDPS(IM, IX, IM, LEVS, DVDT, DUDT,
& UGRS, VGRS, TGRS, QGRS,
& KPBL, PRSI, DEL, PRSL, PRSLK,
& PHII, PHIL, RCL, DTP,
!! & PGR, KPBL, PRSI, DEL, PRSL, PRSLK, RCL, DTP,
& KDT, hprime(1,1), oc, oa4, clx,
& theta,sigma,gamma,elvmax,DUSFCG, DVSFCG,
& grav, CP, RD, RV, LONF, nmtvr, me, lprnt,ipr)
!
if (lprnt) print *,' dudtg=',dudt(ipr,:)
if(lssav) then
do i=1,im
DUGWD(i) = DUGWD(i) + DUSFCG(i)*DTF
DVGWD(i) = DVGWD(i) + DVSFCG(i)*DTF
enddo
! if (lprnt) print *,' dugwd=',dugwd(ipr),' dusfcg=',dusfcg(ipr)
! if (lprnt) print *,' dvgwd=',dvgwd(ipr),' dvsfcg=',dvsfcg(ipr)
IF (LDIAG3D) THEN
DO K=1,LEVS
do i=1,im
tem = rcs(i) * dtf
DU3DT(i,k,2) = DU3DT(i,k,2) + DUDT(i,K)*tem
DV3DT(i,k,2) = DV3DT(i,k,2) + DVDT(i,K)*tem
enddo
ENDDO
! if(lprnt) then
! print *,' gwdu=',du3dt(:,:,2)
! endif
ENDIF
endif
!#endif
!
DO K=1,LEVS
do i=1,im
GT0(i,K) = TGRS(i,K) + DTDT(i,K) * DTp
GU0(i,K) = UGRS(i,K) + DUDT(i,K) * DTp
GV0(i,K) = VGRS(i,K) + DVDT(i,K) * DTp
enddo
enddo
!
DO N=1,ntrac
DO K=1,LEVS
do i=1,im
GQ0(i,K,N) = QGRS(i,K,N) + dqdt(i,k,n) * dtp
enddo
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! if (me .eq. 0) then
! sumq = 0.0
! DO K=1,LEVS
! do i=1,im
! sumq = sumq + (dqdt(i,k,1)+dqdt(i,k,ntcw)) * del(i,k)
! enddo
! enddo
! sume = 0.0
! do i=1,im
! sume = sume + dqsfc1(i)
! enddo
! sumq = sumq * 1000.0 / grav
! sume = sume / hvap
! print *,' after MON: sumq=',sumq,' sume=',sume, ' kdt=',kdt
! endif
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
! OZONE PHYSICS
!
if(ntoz .gt. 0 .and. ntrac .ge. ntoz)then
CALL OZPHYS(IX, IM, LEVS, KO3, DTP, GQ0(1,1,ntoz), GQ0(1,1,ntoz)
&, GT0, poz, prsl, prdout, pl_coeff, DEL, LDIAG3D
&, dq3dt(1,1,6), me)
endif
!
! to side-step the ozone physics
! if(NTRAC .GE. 2) then
! do k = 1, LEVS
! GQ0(k,ntoz) = qgrs(k,ntoz)
! enddo
! endif
!
! if (lprnt) then
! print *,' levs=',levs,' jcap=',jcap,' dtp',dtp
! *,' slmsk=',slmsk(ilon,ilat),' rcs=',rcs,' kdt=',kdt
! print *,' xkt2=',xkt2,' ncld=',ncld,' iq=',iq,' lat=',lat
! print *,' pgr=',pgr
! print *,' del=',del(1,:)
! print *,' prsl=',prsl(1,:)
! print *,' prslk=',prslk(1,:)
! print *,' xkt2=',xkt2(ipr,1)
! print *,' GT0=',GT0(ipr,:)
! &,' kdt=',kdt,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! print *,' dtdt=',dtdt(ipr,:)
! print *,' gu0=',gu0(ipr,:)
! print *,' gv0=',gv0(ipr,:)
! print *,' gq0=',(gq0(ipr,k,1),k=1,levs)
! print *,' gq1=',(gq0(ipr,k,ntcw),k=1,levs)
! print *,' vvel=',vvel
! endif
!
IF (LDIAG3D) THEN
do k=1,levs
do i=1,im
dtdt(i,k) = GT0(i,k)
dqdt(i,k,1) = gq0(i,k,1)
dudt(i,k) = GU0(i,k)
dvdt(i,k) = GV0(i,k)
enddo
enddo
elseif (cnvgwd) then
do k=1,levs
do i=1,im
dtdt(i,k) = GT0(i,k)
enddo
enddo
ENDIF
!
call GET_PHI(im,ix,levs,ntrac,gt0,gq0,
& prsi,prsik,prsl,prslk,phii,phil)
!
!
! if (lprnt) then
! print *,' phii2=',phii(ipr,:)
! print *,' phil2=',phil(ipr,:)
! endif
!
do k=1,levs
do i=1,im
CLW(i,k,1) = 0.0
CLW(i,k,2) = -999.9
enddo
enddo
! For convective tracer transport (while using RAS)
if (ras) then
if (tottracer > 0) then
if (ntoz > 0) then
CLW(:,:,3) = GQ0(:,:,ntoz)
if (tracers > 0) then
do n=1,tracers
CLW(:,:,3+n) = GQ0(:,:,n+trc_shft)
enddo
endif
else
do n=1,tracers
CLW(:,:,2+n) = GQ0(:,:,n+trc_shft)
enddo
endif
endif
endif
!
do i=1,im
ktop(i) = 1
kbot(i) = levs
enddo
!
! Calling precipitation processes
!
do i=1,im
!pry work1(i) = (log(1.0 / (rcs(i)*nlons(i))) - dxmin) * dxinv
! work1(i) = (log(1.0 / (rcs(i)*nlons(i)*latg)) - dxmin) * dxinv
work1(i) = (log(coslat(i) / (nlons(i)*latg)) - dxmin) * dxinv
work1(i) = max(0.0, min(1.0,work1(i)))
work2(i) = 1.0 - work1(i)
enddo
!
! Calling convective parameterization
!
if (ntcw .gt. 0) then
!
do k=1,levs
do i=1,im
rhc(i,k) = rhbbot - (rhbbot-rhbtop) * (1.0-prslk(i,k))
rhc(i,k) = rhc_max * work1(i) + rhc(i,k) * work2(i)
rhc(i,k) = max(0.0, min(1.0,rhc(i,k)))
enddo
enddo
if (num_p3d .eq. 3) then ! Call Brad Ferrier's Microphysics
do i=1,im
flgmin_l(i) = flgmin(1)*work1(i) + flgmin(2)*work2(i)
enddo
!
!*** ALGORITHM TO SEPARATE DIFFERENT HYDROMETEOR SPECIES
!
DO K=1,LEVS
do i=1,im
WC = GQ0(i,K,ntcw)
QI = 0.
QR = 0.
QW = 0.
F_ice = max(0.0, min(1.0, phy_f3d(I,K,1)))
F_rain = max(0.0, min(1.0, phy_f3d(I,K,2)))
!
QI = F_ice*WC
QW = WC-QI
IF (QW .GT. 0.0) THEN
QR = F_rain*QW
QW = QW-QR
ENDIF
!
! IF (F_ice .GE. 1.) THEN
! QI = WC
! ELSE IF (F_ice .LE. 0.) THEN
! QW = WC
! ELSE
! QI = F_ice*WC
! QW = WC-QI
! ENDIF
!
! IF (QW.GT.0. .AND. F_rain.GT.0.) THEN
! IF (F_rain .GE. 1.) THEN
! QR = QW
! QW = 0.
! ELSE
! QR = F_rain*QW
! QW = QW-QR
! ENDIF
! ENDIF
!
QR_col(I,K) = QR
! CLW(I,K) = QI + QW
CLW(I,K,1) = QI
CLW(I,K,2) = QW
!
!
!*** ARRAY TO TRACK FRACTION OF "CLOUD" IN THE FORM OF ICE
!
! IF (QI+QW .GT. EPSQ) THEN
! FC_ice(I,K) = QI / (QI+QW)
! ELSE
! FC_ice(I,K) = 0.
! ENDIF
enddo
ENDDO
else
DO K=1,LEVS
do i=1,im
clw(i,K,1) = GQ0(i,K,ntcw)
enddo
ENDDO
endif
else
rhc(:,:) = 1.0
endif
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (me .eq. 0) then
! sumq = 0.0
! DO K=1,LEVS
! do i=1,im
! sumq = sumq - (gq0(i,k,1)+clw(i,k,1)+clw(i,k,2)) * del(i,k)
! enddo
! enddo
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IF (.NOT. RAS) THEN
if (.not. newsas) then
CALL SASCNV(IM,IX,LEVS,JCAP,dtp,DEL,PRSL,PGR,PHIL,
& CLW,GQ0,GT0,gu0,gv0,rcs,CLD1D,
& RAIN1,KBOT,KTOP,kcnv,SLMSK,
& VVEL,xkt2,ncld)
else
CALL SASCNVN(IM,IX,LEVS,JCAP,dtp,DEL,PRSL,PGR,PHIL,
& CLW,GQ0,GT0,gu0,gv0,rcs,CLD1D,
& RAIN1,KBOT,KTOP,kcnv,SLMSK,
& VVEL,xkt2,ncld)
endif
! if(lprnt) print *,' rain1=',rain1(ipr),' xkt2=',xkt2(ipr,1)
ELSE
! if(lprnt) print *,' calling RAS for kdt=',kdt,' me=',me
! &,' lprnt=',lprnt
CALL RASCNV(IM, IX, LEVS, DTP, DTF, xkt2
! &, GT0, GQ0, GU0, GV0, clw
&, GT0, GQ0, GU0, GV0, clw, tottracer
&, prsi, prsl, prsik, prslk, phil, phii
&, KPBL, CD, RAIN1, KBOT, KTOP, kcnv
&, phy_f2d(1,num_p2d), flipv, cb2mb
&, me, garea, lmh, ccwfac, nrcm, rhc
&, ud_mf, dd_mf, dt_mf, lprnt, ipr)
! &, me, 1, 1, garea, lprnt, ipr)
!
! do i=1,im
! if (tsea(i) .gt. 380.0 .or. tsea(i) .lt. 10) then
! print *,' tsea=', tsea(i),' i=',i,' lat=',lat,
! &' kdt=',kdt,' xlon=',xlon(i),' xlat=',xlat(i),' slmsk=',
! &slmsk(i),' me=',me
! stop
! endif
! enddo
! do k=1,levs
! do i=1,im
! if (gt0(i,k) .gt. 330.0 .or. gt0(i,k) .lt. 80) then
! print *,' gt0=', gt0(i,k),' i=',i,' k=',k,' lat=',lat,
! &' kdt=',kdt,' xlon=',xlon(i),' xlat=',xlat(i)
! stop
! endif
! if (gq0(i,k,1) .gt. 1.0 ) then
! print *,' gq0=', gq0(i,k,1),' i=',i,' k=',k,' lat=',lat,
! &' kdt=',kdt
! stop
! endif
! enddo
! enddo
! if(lprnt) print *,' returning from RAS for kdt=', kdt,' me=',me
! &,' lat=',lat
!
CLD1D = 0
!
! Update the tracers due to convective transport
!
if (tottracer > 0) then
if (ntoz > 0) then ! For ozone
GQ0(:,:,ntoz) = clw(:,:,3)
if (tracers > 0) then ! For other tracers
do n=1,tracers
GQ0(:,:,n+trc_shft) = CLW(:,:,3+n)
enddo
endif
else
do n=1,tracers
GQ0(:,:,n+trc_shft) = CLW(:,:,2+n)
enddo
endif
endif
ENDIF
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (me .eq. 0) then
! DO K=1,LEVS
! do i=1,im
! sumq = sumq + (gq0(i,k,1)+clw(i,k,1)+clw(i,k,2)) * del(i,k)
! enddo
! enddo
! sumr = 0.0
! do i=1,im
! sumr = sumr + rain1(i)
! enddo
! sumq = sumq * 1000.0 / grav
! sumr = sumr *1000
! print *,' after RAS: sumq=',sumq,' sumr=',sumr, ' kdt=',kdt
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
if(lssav) then
do i=1,im
CLDWRK(i) = CLDWRK(i) + CLD1D(i) * DTF
enddo
IF (LDIAG3D) THEN
DO K=1,LEVS
do i=1,im
tem = rcs(i) * FRAIN
DT3DT(i,k,4) = DT3DT(i,k,4) + (GT0(i,k)-dtdt(i,k)) * FRAIN
DQ3DT(i,k,2) = DQ3DT(i,k,2) + (gq0(i,k,1)-dqdt(i,k,1))
& * FRAIN
DU3DT(i,k,3) = DU3DT(i,k,3) + (GU0(i,k)-dudt(i,k)) * tem
DV3DT(i,k,3) = DV3DT(i,k,3) + (GV0(i,k)-dvdt(i,k)) * tem
! cumflx(i,k) = cumflx(i,k) + (ud_mf(i,k)-dd_mf(i,k))
! & * (grav*FRAIN)
upd_mf(i,k) = upd_mf(i,k) + ud_mf(i,k) * (grav*FRAIN)
dwn_mf(i,k) = dwn_mf(i,k) + dd_mf(i,k) * (grav*FRAIN)
det_mf(i,k) = det_mf(i,k) + dt_mf(i,k) * (grav*FRAIN)
enddo
ENDDO
ENDIF
endif
!
do i=1,im
RAINC(i) = FRAIN * RAIN1(i)
enddo
if(lssav) then
do i=1,im
BENGSH(i) = BENGSH(i) + RAINC(i)
enddo
endif
!
!************************************************************************
if (cnvgwd) then ! call CONVECTIVE GRAVITY WAVE DRAG
!
!-----------------------------------------------------------------------
! Calculate Maximum Convective Heating Rate qmax [K/s]
! cuhr = Temperature change due to deep convection
!-----------------------------------------------------------------------
do i=1,IM
qmax(i) = 0.
cumabs(i) = 0.
enddo
do k=1,LEVS
do i=1,IM
! cuhr(i,k) = (GT0(i,k)-DTDT(i,k)) / DTF
cuhr(i,k) = (GT0(i,k)-DTDT(i,k)) / DTP ! Moorthi
cumchr(i,k) = 0.
GWDCU(i,k) = 0.
GWDCV(i,k) = 0.
DIAGN1(i,k) = 0.
DIAGN2(i,k) = 0.
!
if (k >= kbot(i) .and. k <= ktop(i)) then
qmax(i) = max(qmax(i),cuhr(i,k))
cumabs(i) = cuhr(i,k) + cumabs(i)
endif
enddo
enddo
do i=1,IM
do k=KBOT(i),KTOP(i)
do k1=KBOT(i),k
cumchr(i,k) = cuhr(i,k1) + cumchr(i,k)
enddo
cumchr(i,k) = cumchr(i,k) / cumabs(i)
enddo
enddo
!
if (lprnt) then
if (KBOT(ipr).le.KTOP(ipr)) then
write(*,*) 'KBOT <= KTOP for (lat,lon) = ',
& xlon(ipr)*57.29578,xlat(ipr)*57.29578
write(*,*) 'kcnv KBOT KTOP QMAX DLENGTH ',
+ kcnv(ipr),kbot(ipr),ktop(ipr),(86400.*qmax(ipr)),dlength(ipr)
write(*,9000) kdt
do k=KTOP(ipr),KBOT(ipr),-1
write(*,9010) k,(86400.*cuhr(ipr,k)),(100.*cumchr(ipr,k))
enddo
endif
9000 format(/,3x,'K',5x,'CUHR(K)',4x,'CUMCHR(K)',5x,'at KDT = ',i4,/)
9010 format(2x,i2,2x,f8.2,5x,f6.0)
!
print *,' Before GWDC in GBPHYS fhour ',fhour
if (fhour >= fhourpr) then
print *,' Before GWDC in GBPHYS start print'
write(*,*) 'FHOUR IX IM LEVS = ',fhour,ix,im,levs
print *,'dtp dtf RCS = ',dtp,dtf,RCS(ipr)
write(*,9100)
k=LEVS+1
write(*,9110) k,(10.*prsi(ipr,k))
do k=LEVS,1,-1
write(*,9120) k,(10.*prsl(ipr,k)),(10.*del(ipr,k))
write(*,9110) k,(10.*prsi(ipr,k))
enddo
9100 format(//,14x,'PRESSURE LEVELS',//,
+' ILEV',7x,'PRSI',8x,'PRSL',8x,'DELP',/)
9110 format(i4,2x,f10.3)
9120 format(i4,12x,2(2x,f10.3))
write(*,9130)
do k=LEVS,1,-1
write(*,9140) k,UGRS(ipr,k),GU0(ipr,k),
+ VGRS(ipr,k),GV0(ipr,k),
+ TGRS(ipr,k),GT0(ipr,k),DTDT(ipr,k),
+ dudt(ipr,k),dvdt(ipr,k)
enddo
9130 format(//,10x,'Before GWDC in GBPHYS',//,' ILEV',6x,
+'UGRS',9x,'GU0',8x,'VGRS',9x,'GV0',8x,
+'TGRS',9x,'GT0',8x,'GT0B',8x,'DUDT',8x,'DVDT',/)
9140 format(i4,9(2x,f10.3))
print *,' Before GWDC in GBPHYS end print'
endif
endif
CALL GWDC(IM, IX, IM, LEVS, LAT, UGRS, VGRS, TGRS, QGRS,
& RCS, PRSL, PRSI, DEL, QMAX, CUMCHR, KTOP, KBOT, kcnv,
& GWDCU, GWDCV, grav, CP, RD, fv, dlength,
& lprnt, ipr, fhour,
& DUSFCG,DVSFCG,DIAGN1,DIAGN2)
!
if (lprnt) then
if (fhour.ge.fhourpr) then
print *,' After GWDC in GBPHYS start print'
write(*,9131)
do k=LEVS,1,-1
write(*,9141) k,UGRS(ipr,k),GU0(ipr,k),
+ VGRS(ipr,k),GV0(ipr,k),
+ TGRS(ipr,k),GT0(ipr,k),DTDT(ipr,k),
+ GWDCU(ipr,k),GWDCV(ipr,k)
enddo
9131 format(//,10x,'After GWDC in GBPHYS',//,' ILEV',6x,
+'UGRS',9x,'GU0',8x,'VGRS',9x,'GV0',8x,
+'TGRS',9x,'GT0',8x,'GT0B',7x,'GWDCU',7x,'GWDCV',/)
9141 format(i4,9(2x,f10.3))
print *,' After GWDC in GBPHYS end print'
endif
endif
!
!-----------------------------------------------------------------------
! Write out cloud top stress and wind tendencies
!-----------------------------------------------------------------------
if(lssav) then
do i=1,im
DUGWD(i) = DUGWD(i) + DUSFCG(i)*DTF
DVGWD(i) = DVGWD(i) + DVSFCG(i)*DTF
enddo
IF (LDIAG3D) THEN
DO K=1,LEVS
do i=1,im
tem = rcs(i) * dtf
DU3DT(i,k,4) = DU3DT(i,k,4) + GWDCU(i,K)*tem
DV3DT(i,k,4) = DV3DT(i,k,4) + GWDCV(i,K)*tem
! DU3DT(i,k,2) = DU3DT(i,k,2) + DIAGN1(i,K)*tem
! DV3DT(i,k,2) = DV3DT(i,k,2) + DIAGN2(i,K)*tem
enddo
ENDDO
ENDIF
endif
!-----------------------------------------------------------------------
! Update the wind components with GWDC tendencies
!-----------------------------------------------------------------------
DO K=1,LEVS
do i=1,im
GU0(i,K) = GU0(i,K) + GWDCU(i,K) * DTP
GV0(i,K) = GV0(i,K) + GWDCV(i,K) * DTP
enddo
enddo
!
if (lprnt) then
if (fhour.ge.fhourpr) then
print *,' After Tendency GWDC in GBPHYS start print'
write(*,9132)
do k=LEVS,1,-1
write(*,9142) k,UGRS(ipr,k),GU0(ipr,k),VGRS(ipr,k),
& GV0(ipr,k),TGRS(ipr,k),GT0(ipr,k),DTDT(ipr,k),
+ GWDCU(ipr,k),GWDCV(ipr,k)
enddo
9132 format(//,10x,'After Tendency GWDC in GBPHYS',//,' ILEV',6x,
+'UGRS',9x,'GU0',8x,'VGRS',9x,'GV0',8x,
+'TGRS',9x,'GT0',8x,'GT0B',7x,'GWDCU',7x,'GWDCV',/)
9142 format(i4,9(2x,f10.3))
print *,' After Tendency GWDC in GBPHYS end print'
endif
endif
!
endif ! END CONVECTIVE GRAVITY WAVE DRAG
!************************************************************************
!
!
IF (LDIAG3D) THEN
do k=1,levs
do i=1,im
dtdt(i,k) = GT0(i,k)
dqdt(i,k,1) = gq0(i,k,1)
enddo
enddo
ENDIF
!
if (.not. sashal) then
! if (lprnt) print *,' levshcm=',levshcm,' gt0sh=',gt0(ipr,:)
if (.not. mstrat) then
CALL SHALCVt3(IM,IX,LEVSHCM,dtp,DEL,PRSI,PRSL,PRSLK,
& kcnv,GQ0,GT0) ! for pry
else
CALL SHALCV(IM,IX,LEVSHCM,dtp,DEL,PRSI,PRSL,PRSLK,kcnv,GQ0,GT0
&, levshc,phil,kinver,ctei_r,ctei_rm)
! &, DPSHC,phil,kinver)
! &, DPSHC,kinver)
! if (lprnt) print *,' levshcm=',levshcm,' gt0sha=',gt0(ipr,:)
endif
!
else
CALL SHALCNV(IM,IX,LEVS,JCAP,dtp,DEL,PRSL,PGR,PHIL,
& CLW,GQ0,GT0,gu0,gv0,rcs,CLD1D,
& RAIN1,KBOT,KTOP,kcnv,SLMSK,
& VVEL,ncld,hpbl)
!
do i=1,im
RAINCS(i) = FRAIN * RAIN1(i)
enddo
if(lssav) then
do i=1,im
BENGSH(i) = BENGSH(i) + RAINCS(i)
enddo
endif
do i=1,im
RAINC(i) = RAINC(i) + RAINCS(i)
enddo
if(lssav) then
do i=1,im
CLDWRK(i) = CLDWRK(i) + CLD1D(i) * DTF
enddo
endif
endif
!
!
if(lssav) then
IF (LDIAG3D) THEN
DO K=1,LEVS
do i=1,im
DT3DT(i,k,5) = DT3DT(i,k,5) + (GT0(i,k)-dtdt(i,k)) * FRAIN
DQ3DT(i,k,3) = DQ3DT(i,k,3) + (gq0(i,k,1)-dqdt(i,k,1))
& * FRAIN
enddo
ENDDO
do k=1,levs
do i=1,im
dtdt(i,k) = GT0(i,k)
dqdt(i,k,1) = gq0(i,k,1)
enddo
enddo
ENDIF
endif
!
!
DO K=1,LEVS
do i=1,im
if (CLW(I,K,2) .le. -999.0) CLW(I,K,2) = 0.0
enddo
ENDDO
if (ntcw .gt. 0) then
if (num_p3d .eq. 3) then ! Call Brad Ferrier's Microphysics
!
!*** EXTRACT CLOUD WATER & ICE FROM FC_ice
!
DO K=1,LEVS
do i=1,im
! QI = CLW(I,K)*FC_ice(I,K)
! QW = CLW(I,K) - QI
!
QI = CLW(I,K,1)
QW = CLW(I,K,2)
!
!*** ALGORITHM TO COMBINE DIFFERENT HYDROMETEOR SPECIES
!
! GQ0(i,K,ntcw) = MAX(EPSQ, QI+QW+QR_col(i,K))
GQ0(i,K,ntcw) = QI + QW + QR_col(i,K)
IF (QI .LE. EPSQ) THEN
phy_f3d(I,K,1) = 0.
ELSE
phy_f3d(I,K,1) = QI/GQ0(i,K,ntcw)
ENDIF
IF (QR_col(I,K) .LE. EPSQ) THEN
phy_f3d(I,K,2) = 0.
ELSE
phy_f3d(I,K,2) = QR_col(I,K) / (QW+QR_col(I,K))
ENDIF
enddo
ENDDO
else
DO K=1,LEVS
do i=1,im
! GQ0(i,K,ntcw) = CLW(i,K) + GQ0(i,K,ntcw)
GQ0(i,K,ntcw) = CLW(i,K,1) + clw(i,K,2)
! GQ0(i,K,ntcw) = CLW(i,K,1) ! for pry
enddo
ENDDO
endif
else
DO K=1,LEVS
do i=1,im
clw(i,K,1) = CLW(i,K,1) + clw(i,K,2)
enddo
ENDDO
endif
!
CALL CNVC90(CLSTP, IM, IX, RAINC, KBOT, KTOP, LEVS, PRSI,
& aCV, aCVB, aCVT, CV, CVB, CVT)
!
!
IF (NCLD .EQ. 0) THEN
CALL LRGSCL(IX,IM,LEVS,dtp,GT0,GQ0,PRSL,DEL,PRSLK,RAIN1,CLW)
ELSEIF (NCLD .EQ. 1) THEN
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (me .eq. 0) then
! sumq = 0.0
! DO K=1,LEVS
! do i=1,im
! sumq = sumq - (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k)
! enddo
! enddo
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! To call moist convective adjustment
!
! CALL MSTCNV(IM,IX,LEVS,DTP,GT0,GQ0,PRSL,DEL,PRSLK,RAIN1
! &, lprnt,ipr)
! do i=1,im
! RAINC(i) = RAINC(i) + FRAIN * RAIN1(i)
! enddo
! if(lssav) then
! do i=1,im
! BENGSH(i) = BENGSH(i) + RAIN1(i) * FRAIN
! enddo
! IF (LDIAG3D) THEN
! DO K=1,LEVS
! do i=1,im
! DT3DT(i,k,4) = DT3DT(i,k,4) + (GT0(i,k)-dtdt(i,k))
! & * FRAIN
! DQ3DT(i,k,2) = DQ3DT(i,k,2) + (gq0(i,k,1)-dqdt(i,k,1))
! & * FRAIN
! enddo
! ENDDO
! ENDIF
! endif
!
! IF (LDIAG3D) THEN
! do k=1,levs
! do i=1,im
! dtdt(i,k) = GT0(i,k)
! dqdt(i,k,1) = gq0(i,k,1)
! enddo
! enddo
! ENDIF
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (me .eq. 0) then
! DO K=1,LEVS
! do i=1,im
! sumq = sumq + (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k)
! enddo
! enddo
! sumr = 0.0
! do i=1,im
! sumr = sumr + rain1(i)
! enddo
! sumq = sumq * 1000.0 / grav
! sumr = sumr *1000
! print *,' after MCN: sumq=',sumq,' sumr=',sumr, ' kdt=',kdt
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! moist convective adjustment over
!
!
if (num_p3d .eq. 3) then ! Call Brad Ferrier's Microphysics
!
do i=1,im
xncw(i) = ncw(2) * work1(i) + ncw(1) * work2(i)
enddo
if (kdt .eq. 1 .and. abs(xlon(1)) .lt. 0.0001)
& print *,' xncw=',xncw(1),' rhc=',rhc(1,1)
&, ' work1=',work1(1),' work2=',work2(1),' flgmin=',flgmin_l(1)
&, ' lon=',xlon(1) * 57.29578,' lat=',lat,' me=',me
! &, ' lon=',xlon(1) * 57.29578,' lat=',xlat(1) * 57.29578
! &,' kinver=',kinver(1)
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (me .eq. 0) then
! sumq = 0.0
! DO K=1,LEVS
! do i=1,im
! sumq = sumq - (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k)
! enddo
! enddo
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! if (lprnt) print *,' ipr=',ipr,' gt0_gsmb=',gt0(ipr,:)
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
call GSMDRIVE(IM, IX, LEVS, DTP, PRSL, DEL
&, GT0, GQ0(1,1,1), GQ0(1,1,ntcw), slmsk
&, phy_f3d(1,1,1), phy_f3d(1,1,2)
&, phy_f3d(1,1,3), RAIN1, SR, grav
&, hvap, hsub, cp, rhc, xncw, flgmin_l
&, me, lprnt, ipr)
! &, hvap, hsub, cp, rhc, xncw, me, lprnt, ipr)
!
! if (lprnt) print *,' ipr=',ipr,' gt0_gsma=',gt0(ipr,:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (me .eq. 0) then
! DO K=1,LEVS
! do i=1,im
! sumq = sumq + (gq0(i,k,1)+gq0(i,k,ntcw)) * del(i,k)
! enddo
! enddo
! sumr = 0.0
! do i=1,im
! sumr = sumr + rain1(i)
! enddo
! sumq = sumq * 1000.0 / grav
! sumr = sumr *1000
! print *,' after GSM: sumq=',sumq,' sumr=',sumr, ' kdt=',kdt
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
elseif (num_p3d .eq. 4) then ! Call Zhao/Carr/Sundqvist Microphysics
!
! The following commented by Moorthi on 11/01/2002
! do k=1,levs/2
! do i=1,im
!! if (k .le. kinver(i) .and. slmsk(i) .eq. 0.0) then
! if (k .le. kinver(i)) then
! rhc(i,k) = 0.85
! endif
! enddo
! enddo
! if (me .eq. 0 .and. kdt .eq. 1 .and. ilon .eq. 1)
! & print *,' rhc=',rhc
! if (me .eq. 0 .and. kdt .eq. 1)
! & print *,' rhc=',rhc(1,1)
! &, ' lon=',xlon(1) * 57.29578,' lat=',xlat(1) * 57.29578
! &,' kinver=',kinver(1)
! if (kdt .eq. 1 .and. abs(xlon(1)) .lt. 0.0001)
! & print *,' rhc=',rhc(1,1)
! &, ' work1=',work1(1),' work2=',work2(1)
! &, ' lon=',xlon(1) * 57.29578,' lat=',lat,' me=',me
!
! if (lprnt) print *,' ipr=',ipr,' gt0_gsc=',gt0(ipr,:)
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! if (lprnt) print *,' ipr=',ipr,' gq0_gsc=',gq0(ipr,:,1)
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! if (lprnt) print *,' ipr=',ipr,' gc0_gsc=',gq0(ipr,:,ntcw)
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! if (lprnt) print *,' ipr=',ipr,' phy_f3d1=',phy_f3d(ipr,:,1)
! if (lprnt) print *,' ipr=',ipr,' phy_f3d2=',phy_f3d(ipr,:,2)
! if (lprnt) print *,' ipr=',ipr,' phy_f3d3=',phy_f3d(ipr,:,3)
! if (lprnt) print *,' ipr=',ipr,' phy_f2d=',
! &phy_f2d(ipr,1:num_p2d),' num_p2d=',num_p2d
! &, phy_f2d(ipr,2)
CALL GSCOND(IM, IX, LEVS, DTP, PRSL, PGR,
& GQ0(1,1,1), GQ0(1,1,ntcw), GT0,
& phy_f3d(1,1,1), phy_f3d(1,1,2), phy_f2d(1,1),
& phy_f3d(1,1,3), phy_f3d(1,1,4), phy_f2d(1,2),
! & rhc,lprnt)
& rhc,lprnt, ipr)
CALL PRECPD(IM, IX, LEVS, DTP, DEL, PRSL, PGR,
& GQ0(1,1,1), GQ0(1,1,ntcw), GT0, RAIN1,
& rhc, lprnt, ipr)
! & GQ0(1,1,1), GQ0(1,1,ntcw), GT0, RAIN1, rhc, lprnt)
! if (lprnt) print *,' rain1=',rain1(ipr)
endif
!
! ELSEIF (NCLD .EQ. 2) THEN
! CALL CLOUD3(1, 1, LEVS,DTP,PSEXP,
! 1 GT0(1,1),GQ0(1,1,1),GQ0(1,1,ntcw),NCLD,SL,DEL,SLK,
! 2 RAIN1,LAT,VVEL,KDT,FHOUR)
! ELSEIF (NCLD .EQ. 4) THEN
! CALL CLOUD5(1, 1, LEVS,DTP,PSEXP,
! 1 GT0,GQ0(1,1),GQ0(1,ntcw),NCLD,SL,DEL,SLK,
! 2 RAIN1,LAT,VVEL,KDT,FHOUR)
! ELSEIF (NCLD .EQ. 5) THEN
! CALL CLOUD6(1, 1, LEVS,DTP,PSEXP,
! 1 GT0,GQ0(1,1),GQ0(1,ntcw),NCLD,SL,DEL,SLK,
! 2 RAIN1,LAT,VVEL,KDT,FHOUR)
ENDIF
!
! if (lprnt) print *,' rain1=',rain1(ipr),' rainc=',rainc(ipr)
do i=1,im
RAINL(i) = FRAIN * RAIN1(i)
RAIN(i) = RAINC(i) + RAINL(i)
enddo
!!-> Coupling insertion
if(lssav_cc_dummy) then
PRECR_cc_dummy(1:IM)=PRECR_cc_dummy(1:IM)+RAIN(1:IM)
end if
!!<- Coupling insertion
if(lssav) then
do i=1,im
GESHEM(i) = GESHEM(i) + RAIN(i)
enddo
IF (LDIAG3D) THEN
DO K=1,LEVS
do i=1,im
DT3DT(i,k,6) = DT3DT(i,k,6) + (GT0(i,k)-dtdt(i,k)) * FRAIN
DQ3DT(i,k,4) = DQ3DT(i,k,4) + (gq0(i,k,1)-dqdt(i,k,1))
& * FRAIN
enddo
ENDDO
ENDIF
endif
!
! ESTIMATE T850 FOR RAIN-SNOW DECISION
!
do i=1,im
T850(i) = GT0(I,1)
enddo
DO K = 1, LEVS - 1
do i=1,im
IF(PRSL(I,K) .GT. P850 .AND. PRSL(I,K+1) .LE. P850) THEN
T850(i) = GT0(i,K) - (PRSL(i,k)-P850)
& / (PRSL(I,K)-PRSL(I,K+1)) * (GT0(I,K)-GT0(I,K+1))
ENDIF
enddo
ENDDO
!lu [-4L/+3L]: snow-rain detection is performed in land/sice module
!
! FACTOR=WEIGHTED MEAN TEP.
do i=1,im
TPRCP(i) = RAIN(i) !! Clu [+1L]: rain -> tprcp
SRFLAG(i) = 0. !! Clu [+1L]: default srflag as 'rain'
IF(T850(i).LE.273.16) THEN
SRFLAG(i) = 1. !! Clu [+1L]: set srflag to 'snow'
CWei when call OSU LSM
if(lsm.eq.0 .and. SLMSK(i).NE.0.) THEN
SHELEG(i) = SHELEG(i) + 1.E3*RAIN(i) !! neutral impact, for code
consistency
TPRCP(i) = 0.
endif
ENDIF
enddo
! if (lprnt) print *,' TPRCP=',tprcp(ipr),' rain=',rain(ipr)
!-->cpl insertion
! SNW_cc_dummy=0.
do i=1,im
! if(T850(i).LE.273.16.AND.SLMSK(i).NE.0.) THEN
if(T850(i).LE.273.16) THEN
LPREC_cc_dummy(i)=0.0
SNW_cc_dummy(i)=RAIN(i)
else
LPREC_cc_dummy(i)=RAIN(i)
SNW_cc_dummy(i)=0.0
endif
enddo
!<-- cpl insertion
CWei [+2] when call OSU LSM
!
! UPDATE SOIL MOISTURE AND CANOPY WATER AFTER PRECIPITATION computaion
!
if (lsm == 0) then
CALL PROGT2(IM,LSOIL,RHSCNPY,RHSMC,AI,BI,CCI,SMSOIL,
& SLMSK,CANOPY,TPRCP,RUNOF,SNOWMT,
& ZSOIL,SOILTYP,SIGMAF,dtf,me)
CWei [+5]: let soil liquid water equal to soil total water
DO K = 1, LSOIL
do i=1,im
IF(SLMSK(i).EQ.1.) THEN
SLSOIL(i,K) = SMSOIL(i,k)
ENDIF
enddo
ENDDO
endif
!
! TOTAL RUNOFF IS COMPOSED OF DRAINAGE INTO WATER TABLE AND
! RUNOFF AT THE SURFACE AND IS ACCUMULATED IN UNIT OF METERS
!
if(lssav) then
do i=1,im
RUNOFF(i) = RUNOFF(i) + (DRAIN(i)+RUNOF(i)) * DTF / 1000.0
Cwei added 10/24/2006
SRUNOFF(i) = SRUNOFF(i) + RUNOF(i) * DTF / 1000.0
enddo
endif
!c-- XW: FOR SEA-ICE Nov04
!
! RETURN UPDATED ICE THICKNESS & CONCENTRATION TO GLOBAL ARRAYS
!
do i=1,im
IF(SLMSK(i).EQ.2.) THEN
HICE(i) = ZICE(i)
FICE(i) = CICE(i)
TISFC(i) = TICE(i)
else
HICE(i) = 0.0
FICE(i) = 0.0
TISFC(i) = TSEA(i)
endif
enddo
!c-- XW: END SEA-ICE
Clu [-10L]: comment out smc/stc update
!
! RETURN UPDATED SMSOIL AND STSOIL TO GLOBAL ARRAYS
!
DO K = 1, LSOIL
do i=1,im
SMC(i,k) = SMSOIL(i,K)
STC(i,k) = STSOIL(i,K)
SLC(i,k) = SLSOIL(i,K)
enddo
ENDDO
!
! calc. integral of moistue in pwat
!
do i=1,im
PWAT(i) = 0.
enddo
DO K=1,LEVS
do i=1,im
work1(i) = 0.0
enddo
if (ncld .gt. 0) then
do ic=ntcw, ntcw+ncld-1
do i=1,im
! work1(i) = work1(i) + max(gq0(i,k,ic), qmin)
work1(i) = work1(i) + gq0(i,k,ic)
enddo
enddo
endif
do i=1,im
! inside this routine, let t as dry temperature only ! hmhj
! work2(i) = 1.0 + fv * max(gq0(i,k,1),qmin) ! hmhj
! GT0(i,K) = GT0(i,K) * work2(i) ! hmhj
PWAT(i) = PWAT(i) + DEL(i,K)*(GQ0(i,K,1)+work1(i))
enddo
ENDDO
490 CONTINUE
do i=1,im
PWAT(i) = PWAT(i)*(1.E3/grav)
enddo
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (me .eq. 0) then
! do i=1,im
! tem1 = dqsfc(i) - evapo(i)
! tem2 = geshem(i) - raino(i)
! print *,' pwatdif=',pwat(i)-pwato(i),' Edif=',tem1
! &,' Pdif=',tem2,' E-P=',(tem1/hvap-tem2*1000)/frain
! &,' pwato=',pwato(i),' pwat=',pwat(i)
! enddo
! endif
! if (lprnt) then
! do i=1,im
! print *,' i=',i,' gt0=',gt0(i,:),' kdt=',kdt
! &,' xlon=',xlon(i)*57.296,' xlat(i)=',xlat(i)*57.296
! print *,' ipr=',ipr,' gt0=',gt0(ipr,:),' kdt=',kdt,' ipr='
! print *,' gt0E=',gt0(ipr,:)
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! print *,' gu0=',gu0(ipr,:)
! print *,' gv0=',gv0(ipr,:)
! print *,' gq0=',gq0(ipr,:,3)
! print *,' gq0=',gq0(ipr,14,3)
! &,' xlon=',xlon(ipr),' xlat=',xlat(ipr)
! &,ipr
! enddo
! endif
! do i=1,im
! print *,' i=',i,' me=',me,' lat=',lat,' gt0=',gt0(i,:)
! enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
! CALL dscal(LEVS*IX,rcl,GU0,1) ! Commented by Moorthi -- Moved to gbphys_call.f
! CALL dscal(LEVS*IX,rcl,GV0,1)
!
csela write(66,103)gq0(1,2),gq0(levs,2)
103 format(1h ,' at end ozphys gq0(1,2) gq0(levs,2)',2(2x,e12.3))
!
deallocate (clw)
! if(lprnt) call mpi_quit(7)
! if (kdt .gt. 1) call mpi_quit(7)
!
RETURN
END