Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sm jun142020 #466

Closed
wants to merge 36 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
36 commits
Select commit Hold shift + click to select a range
262d1d9
Merge remote-tracking branch 'ncarphy/master'
SMoorthi-emc Jan 10, 2020
99ac1a0
updating cycle to distingush lakes from ocean
SMoorthi-emc Mar 4, 2020
f8eb82c
minot changes in sfcsub
SMoorthi-emc Mar 5, 2020
1e0a583
Merge remote-tracking branch 'upstream/master'
SMoorthi-emc Mar 9, 2020
fd0bc5a
Merge branch 'master' of https://github.com/SMoorthi-EMC/ccpp-physics…
SMoorthi-emc Mar 9, 2020
ab6cb5a
correcting errors in ccpp when levr < levs; this fix is essential for…
SMoorthi-emc Mar 10, 2020
3ba810b
a bug fix
SMoorthi-emc Mar 30, 2020
465f019
after merging with ncar master and some fix to RAS
SMoorthi-emc Apr 10, 2020
74de553
Merge remote-tracking branch 'upstream/master'
SMoorthi-emc Apr 19, 2020
861e4eb
Merge branch 'master' into SM_Apr182020
SMoorthi-emc Apr 19, 2020
f1c24fb
after merging with NCAR/ccpp-physics/master and some bug fixes in som…
SMoorthi-emc Apr 20, 2020
c316f79
removing some unneeded do loops in sfc_drv.f - results reproduce
SMoorthi-emc Apr 21, 2020
5953c52
updating precision of constants in several physics routines
SMoorthi-emc Apr 22, 2020
e19953d
adding _kind_phys to constants in some physics routines
SMoorthi-emc Apr 22, 2020
f85730d
updating to fix a potential snow bug
SMoorthi-emc Apr 27, 2020
0241449
fixing some comment lines
SMoorthi-emc Apr 27, 2020
4694c00
some minor change with same result
SMoorthi-emc Apr 27, 2020
e216116
update consistent with ipd
SMoorthi-emc Apr 28, 2020
13c3ee3
Merge remote-tracking branch 'upstream/master' into SM_Apr282020
SMoorthi-emc Apr 28, 2020
cdac822
testing an alternate option in ccpp
SMoorthi-emc Apr 29, 2020
d61ecbe
some additional updates to the code
SMoorthi-emc Apr 30, 2020
f245b7a
some updates for several routines
SMoorthi-emc May 2, 2020
8ff5aa5
Merge remote-tracking branch 'origin/SM_Apr182020' into SM_Apr282020
SMoorthi-emc May 2, 2020
9d0808d
Merge remote-tracking branch 'upstream/master'
SMoorthi-emc May 3, 2020
2b34488
after merging with ccpp master
SMoorthi-emc May 4, 2020
e146375
removed tisfcin_cpl and tseain_cpl as they are not needed
SMoorthi-emc May 4, 2020
513cb29
minor update to surface_composites
SMoorthi-emc May 5, 2020
bbf56bb
merged with ccpp-physics, updated some code, tested coupled model wit…
SMoorthi-emc May 10, 2020
3cdcdaa
change 633.0 to 622.0_r8
SMoorthi-emc May 11, 2020
4c08f73
updating nst model coupled with ocean model
SMoorthi-emc May 19, 2020
760e9ea
after merging with ccpp-physics-master
SMoorthi-emc May 24, 2020
d810799
some fix related to ice in surface cycling
SMoorthi-emc May 28, 2020
8b77f36
updating sfc_diff.f to compute z0 overocean when ww3 sends z0 values …
SMoorthi-emc Jun 14, 2020
e450d81
after merging with ccpp_physics/master on Jun 14, 2020
SMoorthi-emc Jun 15, 2020
37444dc
updating sfc_diff.f to recompute z0 over ocean when coupled to ww3 an…
SMoorthi-emc Jun 29, 2020
3af3d7f
fixing errors/logic with fractional grid option to reproduce a contin…
SMoorthi-emc Jul 8, 2020
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
134 changes: 70 additions & 64 deletions physics/GFS_MP_generic.F90
Original file line number Diff line number Diff line change
Expand Up @@ -97,8 +97,8 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
logical, intent(in) :: cal_pre, lssav, ldiag3d, cplflx, cplchm

real(kind=kind_phys), intent(in) :: dtf, frain, con_g
real(kind=kind_phys), dimension(im), intent(in) :: rainc, rain1, xlat, xlon, tsfc
real(kind=kind_phys), dimension(im), intent(inout) :: ice, snow, graupel
real(kind=kind_phys), dimension(im), intent(in) :: rain1, xlat, xlon, tsfc
real(kind=kind_phys), dimension(im), intent(inout) :: ice, snow, graupel, rainc
real(kind=kind_phys), dimension(im), intent(in) :: rain0, ice0, snow0, graupel0
real(kind=kind_phys), dimension(ix,nrcm), intent(in) :: rann
real(kind=kind_phys), dimension(im,levs), intent(in) :: gt0, prsl, save_t, save_qv, del
Expand Down Expand Up @@ -140,24 +140,24 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
integer, intent(out) :: errflg

! DH* TODO: CLEANUP, all of these should be coming in through the argument list
real(kind=kind_phys), parameter :: con_p001= 0.001d0
real(kind=kind_phys), parameter :: con_day = 86400.0d0
real(kind=kind_phys), parameter :: rainmin = 1.0d-13
real(kind=kind_phys), parameter :: p850 = 85000.0d0
real(kind=kind_phys), parameter :: con_p001= 0.001_kind_phys
real(kind=kind_phys), parameter :: con_day = 86400.0_kind_phys
real(kind=kind_phys), parameter :: rainmin = 1.0e-13_kind_phys
real(kind=kind_phys), parameter :: p850 = 85000.0_kind_phys
! *DH

integer :: i, k, ic

real(kind=kind_phys), parameter :: zero = 0.0d0, one = 1.0d0
real(kind=kind_phys) :: crain, csnow, onebg, tem, total_precip
real(kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys
real(kind=kind_phys) :: crain, csnow, onebg, tem, total_precip, tem1, tem2
real(kind=kind_phys), dimension(im) :: domr, domzr, domip, doms, t850, work1

! Initialize CCPP error handling variables
errmsg = ''
errflg = 0

onebg = one/con_g

do i = 1, im
rain(i) = rainc(i) + frain * rain1(i) ! time-step convective plus explicit
enddo
Expand All @@ -171,7 +171,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
! physics timestep, while Diag%{rain,rainc} and all totprecip etc
! are on the dynamics timestep. Confusing, but works if frain=1. *DH
if (imp_physics == imp_physics_gfdl) then
tprcp = max(0., rain) ! clu: rain -> tprcp
tprcp = max(zero, rain) ! clu: rain -> tprcp
!graupel = frain*graupel0
!ice = frain*ice0
!snow = frain*snow0
Expand All @@ -180,13 +180,13 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
snow = snow0
! Do it right from the beginning for Thompson
else if (imp_physics == imp_physics_thompson) then
tprcp = max (0.,rainc + frain * rain1) ! time-step convective and explicit precip
tprcp = max (zero, rainc + frain * rain1) ! time-step convective and explicit precip
graupel = frain*graupel0 ! time-step graupel
ice = frain*ice0 ! time-step ice
snow = frain*snow0 ! time-step snow

else if (imp_physics == imp_physics_fer_hires) then
tprcp = max (0.,rain) ! time-step convective and explicit precip
tprcp = max (zero, rain) ! time-step convective and explicit precip
ice = frain*rain1*sr ! time-step ice
end if

Expand All @@ -200,7 +200,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
!Note (GJF): Precipitation LWE thicknesses are multiplied by the frain factor, and are thus on the dynamics time step, but the conversion as written
! (with dtp in the denominator) assumes the rate is calculated on the physics time step. This only works as expected when dtf=dtp (i.e. when frain=1).
if (lsm == lsm_noahmp) then
tem = 1.0 / (dtp*con_p001) !GJF: This conversion was taken from GFS_physics_driver.F90, but should denominator also have the frain factor?
tem = one / (dtp*con_p001) !GJF: This conversion was taken from GFS_physics_driver.F90, but should denominator also have the frain factor?
draincprv(:) = tem * raincprv(:)
drainncprv(:) = tem * rainncprv(:)
dsnowprv(:) = tem * snowprv(:)
Expand All @@ -221,11 +221,11 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt

if (imp_physics /= imp_physics_gfdl .and. imp_physics /= imp_physics_thompson) then
do i=1,im
tprcp(i) = max(0.0, rain(i) )
if(doms(i) > 0.0 .or. domip(i) > 0.0) then
srflag(i) = 1.
tprcp(i) = max(zero, rain(i) )
if(doms(i) > zero .or. domip(i) > zero) then
srflag(i) = one
else
srflag(i) = 0.
srflag(i) = zero
end if
enddo
endif
Expand All @@ -240,34 +240,6 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt

endif

if (lssav) then
! if (Model%me == 0) print *,'in phys drive, kdt=',Model%kdt, &
! 'totprcpb=', Diag%totprcpb(1),'totprcp=',Diag%totprcp(1), &
! 'rain=',Diag%rain(1)
do i=1,im
cnvprcp (i) = cnvprcp (i) + rainc(i)
totprcp (i) = totprcp (i) + rain(i)
totice (i) = totice (i) + ice(i)
totsnw (i) = totsnw (i) + snow(i)
totgrp (i) = totgrp (i) + graupel(i)

cnvprcpb(i) = cnvprcpb(i) + rainc(i)
totprcpb(i) = totprcpb(i) + rain(i)
toticeb (i) = toticeb (i) + ice(i)
totsnwb (i) = totsnwb (i) + snow(i)
totgrpb (i) = totgrpb (i) + graupel(i)
enddo

if (ldiag3d) then
do k=1,levs
do i=1,im
dt3dt(i,k) = dt3dt(i,k) + (gt0(i,k)-save_t(i,k)) * frain
! dq3dt(i,k) = dq3dt(i,k) + (gq0(i,k,1)-save_qv(i,k)) * frain
enddo
enddo
endif
endif

t850(1:im) = gt0(1:im,1)

do k = 1, levs-1
Expand All @@ -287,19 +259,21 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
!! and determine explicit rain/snow by snow/ice/graupel coming out directly from MP
!! and convective rainfall from the cumulus scheme if the surface temperature is below
!! \f$0^oC\f$.

if (imp_physics == imp_physics_gfdl .or. imp_physics == imp_physics_thompson) then

! determine convective rain/snow by surface temperature
! determine large-scale rain/snow by rain/snow coming out directly from MP

if (lsm /= lsm_ruc) then
do i = 1, im
!tprcp(i) = max(0.0, rain(i) )! clu: rain -> tprcp ! DH now lines 245-250
srflag(i) = 0. ! clu: default srflag as 'rain' (i.e. 0)
if (tsfc(i) >= 273.15) then
srflag(i) = zero ! clu: default srflag as 'rain' (i.e. 0)
if (tsfc(i) >= 273.15_kind_phys) then
crain = rainc(i)
csnow = 0.0
csnow = zero
else
crain = 0.0
crain = zero
csnow = rainc(i)
endif
! if (snow0(i,1)+ice0(i,1)+graupel0(i,1)+csnow > rain0(i,1)+crain) then
Expand All @@ -319,30 +293,62 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
endif ! lsm==lsm_ruc
elseif( .not. cal_pre) then
if (imp_physics == imp_physics_mg) then ! MG microphysics
tem = con_day / (dtp * con_p001) ! mm / day
do i=1,im
tprcp(i) = max(0.0, rain(i) ) ! clu: rain -> tprcp
if (rain(i)*tem > rainmin) then
srflag(i) = max(zero, min(one, (rain(i)-rainc(i))*sr(i)/rain(i)))
if (rain(i) > rainmin) then
tem1 = max(zero, (rain(i)-rainc(i))) * sr(i)
tem2 = one / rain(i)
if (t850(i) > 273.16_kind_phys) then
srflag(i) = max(zero, min(one, tem1*tem2))
else
srflag(i) = max(zero, min(one, (tem1+rainc(i))*tem2))
endif
else
srflag(i) = 0.0
srflag(i) = zero
rain(i) = zero
rainc(i) = zero
endif
tprcp(i) = max(zero, rain(i))
enddo
else
else ! not GFDL or MG or Thompson microphysics
do i = 1, im
tprcp(i) = max(0.0, rain(i) ) ! clu: rain -> tprcp
srflag(i) = 0.0 ! clu: default srflag as 'rain' (i.e. 0)
if (t850(i) <= 273.16) then
srflag(i) = 1.0 ! clu: set srflag to 'snow' (i.e. 1)
endif
tprcp(i) = max(zero, rain(i))
srflag(i) = sr(i)
enddo
endif
endif

if (lssav) then
! if (Model%me == 0) print *,'in phys drive, kdt=',kdt, &
! 'totprcpb=', totprcpb(1),'totprcp=',totprcp(1), &
! 'rain=',rain(1)
do i=1,im
cnvprcp (i) = cnvprcp (i) + rainc(i)
totprcp (i) = totprcp (i) + rain(i)
totice (i) = totice (i) + ice(i)
totsnw (i) = totsnw (i) + snow(i)
totgrp (i) = totgrp (i) + graupel(i)

cnvprcpb(i) = cnvprcpb(i) + rainc(i)
totprcpb(i) = totprcpb(i) + rain(i)
toticeb (i) = toticeb (i) + ice(i)
totsnwb (i) = totsnwb (i) + snow(i)
totgrpb (i) = totgrpb (i) + graupel(i)
enddo

if (ldiag3d) then
do k=1,levs
do i=1,im
dt3dt(i,k) = dt3dt(i,k) + (gt0(i,k)-save_t(i,k)) * frain
! dq3dt(i,k) = dq3dt(i,k) + (gq0(i,k,1)-save_qv(i,k)) * frain
enddo
enddo
endif
endif

if (cplflx .or. cplchm) then
do i = 1, im
drain_cpl(i) = rain(i) * (one-srflag(i))
dsnow_cpl(i) = rain(i) * srflag(i)
dsnow_cpl(i)= max(zero, rain(i) * srflag(i))
drain_cpl(i)= max(zero, rain(i) - dsnow_cpl(i))
rain_cpl(i) = rain_cpl(i) + drain_cpl(i)
snow_cpl(i) = snow_cpl(i) + dsnow_cpl(i)
enddo
Expand All @@ -354,10 +360,10 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt
enddo
endif

pwat(:) = 0.0
pwat(:) = zero
do k = 1, levs
do i=1, im
work1(i) = 0.0
work1(i) = zero
enddo
if (ncld > 0) then
do ic = ntcw, ntcw+nncl-1
Expand Down
58 changes: 30 additions & 28 deletions physics/GFS_PBL_generic.F90
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac,

implicit none

integer, parameter :: r8 = kind_phys
integer, intent(in) :: im, levs, nvdiff, ntrac
integer, intent(in) :: ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc
integer, intent(in) :: ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef,ntchs, ntchm
Expand All @@ -115,9 +116,11 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac,
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg

real (kind=kind_phys), parameter :: zero = 0.0_r8, one=1.0_r8

! Parameters for canopy heat storage parametrization
real (kind=kind_phys), parameter :: z0min=0.2, z0max=1.0
real (kind=kind_phys), parameter :: u10min=2.5, u10max=7.5
real (kind=kind_phys), parameter :: z0min=0.2_r8, z0max=one
real (kind=kind_phys), parameter :: u10min=2.5_r8, u10max=7.5_r8

! Local variables
integer :: i, k, kk, k1, n
Expand Down Expand Up @@ -283,20 +286,20 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac,
do i=1,im
hflxq(i) = hflx(i)
evapq(i) = evap(i)
hffac(i) = 1.0
hefac(i) = 1.0
hffac(i) = one
hefac(i) = one
enddo
if (lheatstrg) then
do i=1,im
tem = 0.01 * zorl(i) ! change unit from cm to m
tem = 0.01_r8 * zorl(i) ! change unit from cm to m
tem1 = (tem - z0min) / (z0max - z0min)
hffac(i) = z0fac * min(max(tem1, 0.0), 1.0)
tem = sqrt(u10m(i)**2+v10m(i)**2)
hffac(i) = z0fac * min(max(tem1, zero), one)
tem = sqrt(u10m(i)*u10m(i) + v10m(i)*v10m(i))
tem1 = (tem - u10min) / (u10max - u10min)
tem2 = 1.0 - min(max(tem1, 0.0), 1.0)
tem2 = one - min(max(tem1, zero), one)
hffac(i) = tem2 * hffac(i)
hefac(i) = 1. + e0fac * hffac(i)
hffac(i) = 1. + hffac(i)
hefac(i) = one + e0fac * hffac(i)
hffac(i) = one + hffac(i)
hflxq(i) = hflx(i) / hffac(i)
evapq(i) = evap(i) / hefac(i)
enddo
Expand Down Expand Up @@ -339,6 +342,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,

implicit none

integer, parameter :: r8 = kind_phys
integer, intent(in) :: im, levs, nvdiff, ntrac, ntchs, ntchm
integer, intent(in) :: ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef
logical, intent(in) :: trans_aero
Expand All @@ -365,26 +369,26 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
! use assumed-shape arrays. Note that Intel 18 and GNU 6.2.0-8.1.0 tolerate explicit-shape arrays
! as long as these do not get used when not allocated.
real(kind=kind_phys), dimension(:,:), intent(inout) :: dt3dt, du3dt_PBL, du3dt_OGWD, dv3dt_PBL, dv3dt_OGWD, dq3dt, dq3dt_ozone
real(kind=kind_phys), dimension(:), intent(inout) :: dusfc_cpl, dvsfc_cpl, dtsfc_cpl, dqsfc_cpl, dusfci_cpl, dvsfci_cpl, &
real(kind=kind_phys), dimension(:), intent(inout) :: dusfc_cpl, dvsfc_cpl, dtsfc_cpl, dqsfc_cpl, dusfci_cpl, dvsfci_cpl, &
dtsfci_cpl, dqsfci_cpl, dusfc_diag, dvsfc_diag, dtsfc_diag, dqsfc_diag, dusfci_diag, dvsfci_diag, dtsfci_diag, dqsfci_diag

logical, dimension(:),intent(in) :: wet, dry, icy
real(kind=kind_phys), dimension(:), intent(out) :: ushfsfci

real(kind=kind_phys), dimension(:,:), intent(inout) :: dkt_cpl
real(kind=kind_phys), dimension(:,:), intent(in) :: dkt
real(kind=kind_phys), dimension(:,:), intent(in) :: dkt

! From canopy heat storage - reduction factors in latent/sensible heat flux due to surface roughness
real(kind=kind_phys), dimension(im), intent(in) :: hffac, hefac

character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg

real(kind=kind_phys), parameter :: zero = 0.0d0
real(kind=kind_phys), parameter :: one = 1.0d0
real(kind=kind_phys), parameter :: zero = 0.0_r8, one = 1.0_r8
real(kind=kind_phys), parameter :: huge = 9.9692099683868690E36 ! NetCDF float FillValue, same as in GFS_typedefs.F90
real(kind=kind_phys), parameter :: qmin = 1.0e-8_r8
integer :: i, k, kk, k1, n
real(kind=kind_phys) :: tem, tem1, rho
real(kind=kind_phys) :: tem, rho

! Initialize CCPP error handling variables
errmsg = ''
Expand Down Expand Up @@ -437,12 +441,12 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
! Ferrier-Aligo
do k=1,levs
do i=1,im
dqdt(i,k,ntqv) = dvdftra(i,k,1)
dqdt(i,k,ntcw) = dvdftra(i,k,2)
dqdt(i,k,ntiw) = dvdftra(i,k,3)
dqdt(i,k,ntrw) = dvdftra(i,k,4)
dqdt(i,k,ntqv) = dvdftra(i,k,1)
dqdt(i,k,ntcw) = dvdftra(i,k,2)
dqdt(i,k,ntiw) = dvdftra(i,k,3)
dqdt(i,k,ntrw) = dvdftra(i,k,4)
dqdt(i,k,nqrimef) = dvdftra(i,k,5)
dqdt(i,k,ntoz) = dvdftra(i,k,6)
dqdt(i,k,ntoz) = dvdftra(i,k,6)
enddo
enddo

Expand Down Expand Up @@ -536,8 +540,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,

if (cplchm) then
do i = 1, im
tem1 = max(q1(i), 1.e-8)
tem = prsl(i,1) / (rd*t1(i)*(one+fvirt*tem1))
tem = prsl(i,1) / (rd*t1(i)*(one+fvirt*max(q1(i), qmin)))
ushfsfci(i) = -cp * tem * hflx(i) ! upward sensible heat flux
enddo
! dkt_cpl has dimensions (1:im,1:levs), but dkt has (1:im,1:levs-1)
Expand All @@ -563,8 +566,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
dqsfci_cpl(i) = dqsfc1(i)
end if
elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point
tem1 = max(q1(i), 1.e-8)
rho = prsl(i,1) / (rd*t1(i)*(one+fvirt*tem1))
rho = prsl(i,1) / (rd*t1(i)*(one+fvirt*max(q1(i), qmin)))
if (wind(i) > zero) then
tem = - rho * stress_wat(i) / wind(i)
dusfci_cpl(i) = tem * ugrs1(i) ! U-momentum flux
Expand Down Expand Up @@ -600,14 +602,14 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
!-------------------------------------------------------lssav if loop ----------
if (lssav) then
do i=1,im
dusfc_diag (i) = dusfc_diag(i) + dusfc1(i)*dtf
dvsfc_diag (i) = dvsfc_diag(i) + dvsfc1(i)*dtf
dtsfc_diag (i) = dtsfc_diag(i) + dtsfc1(i)*hffac(i)*dtf
dqsfc_diag (i) = dqsfc_diag(i) + dqsfc1(i)*hefac(i)*dtf
dusfc_diag (i) = dusfc_diag(i) + dusfc1(i) * dtf
dvsfc_diag (i) = dvsfc_diag(i) + dvsfc1(i) * dtf
dusfci_diag(i) = dusfc1(i)
dvsfci_diag(i) = dvsfc1(i)
dtsfci_diag(i) = dtsfc1(i)*hffac(i)
dqsfci_diag(i) = dqsfc1(i)*hefac(i)
dtsfc_diag (i) = dtsfc_diag(i) + dtsfci_diag(i) * dtf
dqsfc_diag (i) = dqsfc_diag(i) + dqsfci_diag(i) * dtf
enddo

if (ldiag3d) then
Expand Down
Loading