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

Bugfixes #8

Merged
merged 12 commits into from
Nov 15, 2019
7 changes: 4 additions & 3 deletions .gitmodules
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
[submodule "atmos_cubed_sphere"]
path = atmos_cubed_sphere
url = https://github.com/NOAA-EMC/GFDL_atmos_cubed_sphere
branch = dev/emc
url = https://github.com/junwang-noaa/GFDL_atmos_cubed_sphere
branch = regfv3_rst
[submodule "ccpp/framework"]
path = ccpp/framework
url = https://github.com/NCAR/ccpp-framework
[submodule "ccpp/physics"]
path = ccpp/physics
url = https://github.com/NCAR/ccpp-physics
url = https://github.com/junwang-noaa/ccpp-physics
branch = gfdlmp_update
2 changes: 1 addition & 1 deletion atmos_cubed_sphere
2 changes: 1 addition & 1 deletion ccpp/physics
2 changes: 1 addition & 1 deletion gfsphysics/GFS_layer/GFS_physics_driver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1149,7 +1149,7 @@ subroutine GFS_physics_driver &
if (fice(i) < one) then
wet(i) = .true.
! Sfcprop%tsfco(i) = tgice
Sfcprop%tsfco(i) = max(Sfcprop%tisfc(i), tgice)
if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tisfc(i), tgice)
! Sfcprop%tsfco(i) = max((Sfcprop%tsfc(i) - fice(i)*sfcprop%tisfc(i)) &
! / (one - fice(i)), tgice)
endif
Expand Down
11 changes: 10 additions & 1 deletion gfsphysics/GFS_layer/GFS_restart.F90
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,9 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
#endif

Restart%num3d = Model%ntot3d
if(Model%lrefres) then
Restart%num3d = Model%ntot3d+1
endif
#ifdef CCPP
! GF
if (Model%imfdeepcnv == 3) then
Expand Down Expand Up @@ -252,7 +255,13 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
Restart%data(nb,num)%var3p => Tbd(nb)%phy_f3d(:,:,num)
enddo
enddo

if (Model%lrefres) then
num = Model%ntot3d+1
restart%name3d(num) = 'ref_f3d'
do nb = 1,nblks
Restart%data(nb,num)%var3p => IntDiag(nb)%refl_10cm(:,:)
enddo
endif
#ifdef CCPP
!--- RAP/HRRR-specific variables, 3D
num = Model%ntot3d
Expand Down
12 changes: 11 additions & 1 deletion gfsphysics/GFS_layer/GFS_typedefs.F90
Original file line number Diff line number Diff line change
Expand Up @@ -701,6 +701,9 @@ module GFS_typedefs
!--- GFDL microphysical paramters
logical :: lgfdlmprad !< flag for GFDL mp scheme and radiation consistency

!--- Thompson,GFDL mp parameter
logical :: lrefres !< flag for radar reflectivity in restart file

!--- land/surface model parameters
integer :: lsm !< flag for land surface model lsm=1 for noah lsm
integer :: lsm_noah=1 !< flag for NOAH land surface model
Expand Down Expand Up @@ -2740,6 +2743,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
!--- GFDL microphysical parameters
logical :: lgfdlmprad = .false. !< flag for GFDLMP radiation interaction

!--- Thompson,GFDL microphysical parameter
logical :: lrefres = .false. !< flag for radar reflectivity in restart file

!--- land/surface model parameters
integer :: lsm = 1 !< flag for land surface model to use =0 for osu lsm; =1 for noah lsm; =2 for noah mp lsm; =3 for RUC lsm
integer :: lsoil = 4 !< number of soil layers
Expand Down Expand Up @@ -3023,7 +3029,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
mg_do_graupel, mg_do_hail, mg_nccons, mg_nicons, mg_ngcons, &
mg_ncnst, mg_ninst, mg_ngnst, sed_supersat, do_sb_physics, &
mg_alf, mg_qcmin, mg_do_ice_gmao, mg_do_liq_liu, &
ltaerosol, lradar, ttendlim, lgfdlmprad, &
ltaerosol, lradar, lrefres, ttendlim, lgfdlmprad, &
!--- max hourly
avg_max_length, &
!--- land/surface model control
Expand Down Expand Up @@ -3312,6 +3318,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%ttendlim = ttendlim
!--- gfdl MP parameters
Model%lgfdlmprad = lgfdlmprad
!--- Thompson,GFDL MP parameter
Model%lrefres = lrefres

!--- land/surface model parameters
Model%lsm = lsm
Expand Down Expand Up @@ -4310,6 +4318,7 @@ subroutine control_print(Model)
print *, ' Thompson microphysical parameters'
print *, ' ltaerosol : ', Model%ltaerosol
print *, ' lradar : ', Model%lradar
print *, ' lrefres : ', Model%lrefres
print *, ' ttendlim : ', Model%ttendlim
print *, ' '
endif
Expand All @@ -4327,6 +4336,7 @@ subroutine control_print(Model)
if (Model%imp_physics == Model%imp_physics_gfdl) then
print *, ' GFDL microphysical parameters'
print *, ' GFDL MP radiation inter: ', Model%lgfdlmprad
print *, ' lrefres : ', Model%lrefres
print *, ' '
endif

Expand Down
8 changes: 4 additions & 4 deletions gfsphysics/physics/gfdl_cloud_microphys.F90
Original file line number Diff line number Diff line change
Expand Up @@ -4688,7 +4688,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
real :: qmin = 1.0e-12, beta = 1.22
real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6

do k = ks, ke
do i = is, ie
Expand Down Expand Up @@ -4718,7 +4718,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
! cloud ice (Heymsfield and Mcfarquhar, 1996)
! -----------------------------------------------------------------------

if (qmi (i, k) .gt. qmin) then
if (qmi (i, k) .gt. qmin1) then
qci (i, k) = dpg * qmi (i, k) * 1.0e3
rei_fac = log (1.0e3 * qmi (i, k) * den (i, k))
if (t (i, k) - tice .lt. - 50) then
Expand All @@ -4744,7 +4744,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
! cloud ice (Wyser, 1998)
! -----------------------------------------------------------------------

if (qmi (i, k) .gt. qmin) then
if (qmi (i, k) .gt. qmin1) then
qci (i, k) = dpg * qmi (i, k) * 1.0e3
bw = - 2. + 1.e-3 * log10 (den (i, k) * qmi (i, k) / rho_0) * max (0.0, tice - t (i, k)) ** 1.5
rei (i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
Expand Down Expand Up @@ -4774,7 +4774,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms,
! snow (Lin et al., 1983)
! -----------------------------------------------------------------------

if (qms (i, k) .gt. qmin) then
if (qms (i, k) .gt. qmin1) then
qcs (i, k) = dpg * qms (i, k) * 1.0e3
lambdas = exp (0.25 * log (pi * rhos * n0s / qms (i, k) / den (i, k)))
res (i, k) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6
Expand Down
53 changes: 27 additions & 26 deletions gfsphysics/physics/satmedmfvdifq.f
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,7 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
& rlmn, rlmn1, rlmx, elmx,
& ttend, utend, vtend, qtend,
& zfac, zfmin, vk, spdk2,
& tkmin, xkzinv, xkgdx,
& tkmin, tkminx, xkzinv, xkgdx,
& zlup, zldn, bsum,
& tem, tem1, tem2,
& ptem, ptem0, ptem1, ptem2
Expand All @@ -176,11 +176,11 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
parameter(prmin=0.25,prmax=4.0)
parameter(pr0=1.0,prtke=1.0,prscu=0.67)
parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
parameter(tkmin=1.e-9,dspmax=10.0)
parameter(tkmin=1.e-9,tkminx=0.2,dspmax=10.0)
parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8)
parameter(aphi5=5.,aphi16=16.)
parameter(elmfac=1.0,elefac=1.0,cql=100.)
parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=25000.)
parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=5000.)
parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.1)
parameter(h1=0.33333333)
parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15)
Expand Down Expand Up @@ -273,20 +273,20 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
xkzo(i,k) = 0.0
xkzmo(i,k) = 0.0
if (k < kinver(i)) then
! vertical background diffusivity
ptem = prsi(i,k+1) * tx1(i)
tem1 = 1.0 - ptem
tem2 = tem1 * tem1 * 10.0
tem2 = min(1.0, exp(-tem2))
xkzo(i,k) = xkzm_hx(i) * tem2
!
! minimum turbulent mixing length
ptem = prsl(i,k) * tx1(i)
tem1 = 1.0 - ptem
tem2 = tem1 * tem1 * 2.5
tem2 = min(1.0, exp(-tem2))
rlmnz(i,k)= rlmn * tem2
rlmnz(i,k)= max(rlmnz(i,k), rlmn1)
! vertical background diffusivity for momentum
! vertical background diffusivity
ptem = prsi(i,k+1) * tx1(i)
tem1 = 1.0 - ptem
tem2 = tem1 * tem1 * 10.0
tem2 = min(1.0, exp(-tem2))
xkzo(i,k) = xkzm_hx(i) * tem2
! vertical background diffusivity for momentum
if (ptem >= xkzm_s) then
xkzmo(i,k) = xkzm_mx(i)
kx1(i) = k + 1
Expand Down Expand Up @@ -674,20 +674,20 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
!
! background diffusivity decreasing with increasing surface layer stability
!
do i = 1, im
if(.not.sfcflg(i)) then
tem = (1. + 5. * rbsoil(i))**2.
! tem = (1. + 5. * zol(i))**2.
frik(i) = 0.1 + 0.9 / tem
endif
enddo
!
do k = 1,km1
do i=1,im
xkzo(i,k) = frik(i) * xkzo(i,k)
xkzmo(i,k)= frik(i) * xkzmo(i,k)
enddo
enddo
! do i = 1, im
! if(.not.sfcflg(i)) then
! tem = (1. + 5. * rbsoil(i))**2.
!! tem = (1. + 5. * zol(i))**2.
! frik(i) = 0.1 + 0.9 / tem
! endif
! enddo
!
! do k = 1,km1
! do i=1,im
! xkzo(i,k) = frik(i) * xkzo(i,k)
! xkzmo(i,k)= frik(i) * xkzmo(i,k)
! enddo
! enddo
!
! The background vertical diffusivities in the inversion layers are limited
! to be less than or equal to xkzminv
Expand Down Expand Up @@ -867,13 +867,14 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
do i = 1, im
if(k == 1) then
tem = ckz(i,1)
tem1 = xkzmo(i,1)
tem1 = 0.5 * xkzmo(i,1)
else
tem = 0.5 * (ckz(i,k-1) + ckz(i,k))
tem1 = 0.5 * (xkzmo(i,k-1) + xkzmo(i,k))
endif
ptem = tem1 / (tem * elm(i,k))
tkmnz(i,k) = ptem * ptem
tkmnz(i,k) = min(tkmnz(i,k), tkminx)
tkmnz(i,k) = max(tkmnz(i,k), tkmin)
enddo
enddo
Expand Down
9 changes: 2 additions & 7 deletions io/FV3GFS_io.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1119,15 +1119,10 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain)
Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix)
Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorlo(ix)
Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfco(ix)
if (Sfcprop(nb)%slmsk(ix) < 0.1 .or. Sfcprop(nb)%slmsk(ix) > 1.9) then
if (Sfcprop(nb)%slmsk(ix) > 1.9) then
Sfcprop(nb)%landfrac(ix) = 0.0
if (Sfcprop(nb)%oro_uf(ix) > 0.01) then
Sfcprop(nb)%lakefrac(ix) = 1.0 ! lake
else
Sfcprop(nb)%lakefrac(ix) = 0.0 ! ocean
endif
else
Sfcprop(nb)%landfrac(ix) = 1.0 ! land
Sfcprop(nb)%landfrac(ix) = Sfcprop(nb)%slmsk(ix)
endif
enddo
enddo
Expand Down
67 changes: 40 additions & 27 deletions io/module_wrt_grid_comp.F90
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ module module_wrt_grid_comp
logical,save :: first_init=.false.
logical,save :: first_run=.false.
logical,save :: first_getlatlon=.true.
logical,save :: change_wrtidate=.false.
!
!-----------------------------------------------------------------------
!
Expand Down Expand Up @@ -173,7 +174,8 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
real(ESMF_KIND_R4) :: valueR4
real(ESMF_KIND_R8) :: valueR8

integer :: attCount, axeslen, jidx, noutfile
integer :: attCount, axeslen, jidx, idx, noutfile
character(19) :: newdate
character(128) :: FBlist_outfilename(100), outfile_name
character(128),dimension(:,:), allocatable :: outfilename
real(8), dimension(:), allocatable :: slat
Expand All @@ -183,7 +185,6 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
real(ESMF_KIND_R8) :: geo_lon, geo_lat
real(ESMF_KIND_R8) :: lon1_r8, lat1_r8
real(ESMF_KIND_R8) :: x1, y1, x, y
type(ESMF_Time) :: IO_BASETIME_IAU
type(ESMF_TimeInterval) :: IAU_offsetTI
type(ESMF_DataCopy_Flag) :: copyflag=ESMF_DATACOPY_REFERENCE
! real(8),parameter :: PI=3.14159265358979d0
Expand Down Expand Up @@ -469,8 +470,32 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
call ESMF_Finalize(endflag=ESMF_END_ABORT)

endif
!
!-----------------------------------------------------------------------
!*** get write grid component initial time from clock
!-----------------------------------------------------------------------
!
call ESMF_ClockGet(clock =CLOCK & !<-- The ESMF Clock
,startTime=wrt_int_state%IO_BASETIME & !<-- The Clock's starting time
,rc =RC)


call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), &
m=idate(5),s=idate(6),rc=rc)
! if (lprnt) write(0,*) 'in wrt initial, io_baseline time=',idate,'rc=',rc
idate(7) = 1
wrt_int_state%idate = idate
wrt_int_state%fdate = idate
Copy link
Collaborator

@jswhit2 jswhit2 Nov 12, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add a comment here since it's not obvious what's going on. Something like "update idate when IAU is enabled to reference center of IAU window"

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added comments "! update IO-BASETIME and idate on write grid comp when IAU is enabled"

! update IO-BASETIME and idate on write grid comp when IAU is enabled
if(iau_offset > 0 ) then
call ESMF_TimeIntervalSet(IAU_offsetTI, h=iau_offset, rc=rc)
wrt_int_state%IO_BASETIME = wrt_int_state%IO_BASETIME + IAU_offsetTI
call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), &
m=idate(5),s=idate(6),rc=rc)
wrt_int_state%idate = idate
change_wrtidate = .true.
if (lprnt) print *,'in wrt initial, with iau, io_baseline time=',idate,'rc=',rc
endif
!
! Create field bundle
!-------------------------------------------------------------------
!
Expand Down Expand Up @@ -867,6 +892,17 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)

if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return

! update the time:units when idate on write grid component is changed
if ( index(trim(attNameList(i)),'time:units')>0) then
if ( change_wrtidate ) then
idx = index(trim(valueS),' since ')
if(lprnt) print *,'in write grid comp, time:unit=',trim(valueS)
write(newdate,'(I4.4,a,I2.2,a,I2.2,a,I2.2,a,I2.2,a,I2.2)') idate(1),'-', &
idate(2),'-',idate(3),' ',idate(4),':',idate(5),':',idate(6)
valueS = valueS(1:idx+6)//newdate
if(lprnt) print *,'in write grid comp, new time:unit=',trim(valueS)
endif
endif
call ESMF_AttributeSet(wrtgrid, convention="NetCDF", purpose="FV3", &
name=trim(attNameList(i)), value=valueS, rc=rc)

Expand Down Expand Up @@ -1036,28 +1072,6 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
!
!-----------------------------------------------------------------------
!*** SET THE IO_BaseTime TO THE INITIAL CLOCK TIME.
!-----------------------------------------------------------------------
!
call ESMF_ClockGet(clock =CLOCK & !<-- The ESMF Clock
,startTime=wrt_int_state%IO_BASETIME & !<-- The Clock's starting time
,rc =RC)

call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), &
m=idate(5),s=idate(6),rc=rc)
! if (lprnt) write(0,*) 'in wrt initial, io_baseline time=',idate,'rc=',rc
idate(7) = 1
wrt_int_state%idate = idate
wrt_int_state%fdate = idate
if(iau_offset > 0 ) then
call ESMF_TimeIntervalSet(IAU_offsetTI, h=iau_offset, rc=rc)
IO_BASETIME_IAU = wrt_int_state%IO_BASETIME + IAU_offsetTI
call ESMF_TimeGet(time=IO_BASETIME_IAU,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), &
m=idate(5),s=idate(6),rc=rc)
! if (lprnt) write(0,*) 'in wrt initial, with iau, io_baseline time=',idate,'rc=',rc
endif
!
!-----------------------------------------------------------------------
!*** SET THE FIRST HISTORY FILE'S TIME INDEX.
!-----------------------------------------------------------------------
!
Expand Down Expand Up @@ -1250,10 +1264,9 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc)
! 'nseconds_num=',nseconds_num,nseconds_den,'mype=',mype
!
nf_seconds = nf_hours*3600+nf_minuteS*60+nseconds+real(nseconds_num)/real(nseconds_den)
! shift forecast hour by iau_offset if iau is on.
nf_seconds = nf_seconds - iau_offset*3600
wrt_int_state%nfhour = nf_seconds/3600.
nf_hours = int(nf_seconds/3600.)
if(mype == lead_write_task) print *,'in write grid comp, nf_hours=',nf_hours
! if iau_offset > nf_hours, don't write out anything
if (nf_hours < 0) return

Expand Down
Loading