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

Fixes on initializing snow depth over ice and changes z0ice #461

Merged
merged 7 commits into from
Jan 10, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 1 addition & 1 deletion atmos_model.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1572,7 +1572,7 @@ subroutine assign_importdata(jdat, rc)
type(ESMF_Grid) :: grid
type(ESMF_Field) :: dbgField
character(19) :: currtimestring
real (kind=GFS_kind_phys), parameter :: z0ice=1.1 ! (in cm)
real (kind=GFS_kind_phys), parameter :: z0ice=1.0 ! (in cm)

!
! real(kind=GFS_kind_phys), parameter :: himax = 8.0 !< maximum ice thickness allowed
Expand Down
13 changes: 13 additions & 0 deletions cpl/module_block_data.F90
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ subroutine block_copy_1d_i4_to_2d_r8(destin_ptr, source_ptr, block, block_index,
if (associated(destin_ptr) .and. associated(source_ptr)) then
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -105,6 +106,7 @@ subroutine block_copy_1d_to_2d_r8(destin_ptr, source_ptr, block, block_index, sc
if (associated(destin_ptr) .and. associated(source_ptr)) then
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -144,6 +146,7 @@ subroutine block_copy_1dslice_to_2d_r8(destin_ptr, source_ptr, slice, block, blo
if (slice > 0 .and. slice <= size(source_ptr, dim=2)) then
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -182,6 +185,7 @@ subroutine block_copy_2d_to_3d_r8(destin_ptr, source_ptr, block, block_index, sc
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
do k = 1, size(source_ptr, dim=2)
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -219,6 +223,7 @@ subroutine block_copy_2d_to_2d_r8(destin_ptr, source_ptr, block, block_index, sc
if (associated(destin_ptr) .and. associated(source_ptr)) then
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -253,6 +258,7 @@ subroutine block_array_copy_2d_to_2d_r8(destin_ptr, source_arr, block, block_ind
if (associated(destin_ptr)) then
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -290,6 +296,7 @@ subroutine block_copy_3d_to_3d_r8(destin_ptr, source_ptr, block, block_index, sc
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
do k = 1, size(source_ptr, dim=3)
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -326,6 +333,7 @@ subroutine block_array_copy_3d_to_3d_r8(destin_ptr, source_arr, block, block_ind
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
do k = 1, size(source_arr, dim=3)
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -367,6 +375,7 @@ subroutine block_copy_3dslice_to_3d_r8(destin_ptr, source_ptr, slice, block, blo
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
do k = 1, size(source_ptr, dim=3)
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -407,6 +416,7 @@ subroutine block_array_copy_3dslice_to_3d_r8(destin_ptr, source_arr, slice, bloc
factor = 1._kind_phys
if (present(scale_factor)) factor = scale_factor
do k = 1, size(source_arr, dim=3)
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -441,6 +451,7 @@ subroutine block_fill_2d_r8(destin_ptr, fill_value, block, block_index, rc)
! -- begin
localrc = ESMF_RC_PTR_NOTALLOC
if (associated(destin_ptr)) then
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -474,6 +485,7 @@ subroutine block_fill_3d_r8(destin_ptr, fill_value, block, block_index, rc)
localrc = ESMF_RC_PTR_NOTALLOC
if (associated(destin_ptr)) then
do k = 1, size(destin_ptr, dim=3)
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down Expand Up @@ -586,6 +598,7 @@ subroutine block_combine_frac_1d_to_2d_r8(destin_ptr, fract1_ptr, fract2_ptr, bl
localrc = ESMF_RC_PTR_NOTALLOC
if (associated(destin_ptr) .and. &
associated(fract1_ptr) .and. associated(fract2_ptr)) then
!$omp parallel do private(ix,ib,jb,i,j)
do ix = 1, block%blksz(block_index)
ib = block%index(block_index)%ii(ix)
jb = block%index(block_index)%jj(ix)
Expand Down
78 changes: 56 additions & 22 deletions io/FV3GFS_io.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1462,7 +1462,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta
do nb = 1, Atm_block%nblks
do ix = 1, Atm_block%blksz(nb)
if (Sfcprop(nb)%landfrac(ix) > zero) then
tem = one / Sfcprop(nb)%landfrac(ix)
tem = one / (Sfcprop(nb)%fice(ix)*(one-Sfcprop(nb)%landfrac(ix))+Sfcprop(nb)%landfrac(ix))
Copy link
Collaborator

Choose a reason for hiding this comment

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

Moorthi, why do we include the ice section when computing the snod on land?

Sfcprop(nb)%snodl(ix) = Sfcprop(nb)%snowd(ix) * tem
else
Sfcprop(nb)%snodl(ix) = zero
Expand All @@ -1477,7 +1477,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta
do nb = 1, Atm_block%nblks
do ix = 1, Atm_block%blksz(nb)
if (Sfcprop(nb)%landfrac(ix) > zero) then
tem = one / Sfcprop(nb)%landfrac(ix)
tem = one / (Sfcprop(nb)%fice(ix)*(one-Sfcprop(nb)%landfrac(ix))+Sfcprop(nb)%landfrac(ix))
Sfcprop(nb)%weasdl(ix) = Sfcprop(nb)%weasd(ix) * tem
else
Sfcprop(nb)%weasdl(ix) = zero
Expand All @@ -1501,7 +1501,9 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta
!$omp parallel do default(shared) private(nb, ix)
do nb = 1, Atm_block%nblks
do ix = 1, Atm_block%blksz(nb)
Sfcprop(nb)%zorlw(ix) = Sfcprop(nb)%zorl(ix) !--- compute zorlw from existing variables
if (Sfcprop(nb)%landfrac(ix) < one .and. Sfcprop(nb)%fice(ix) < one) then
Sfcprop(nb)%zorlw(ix) = min(Sfcprop(nb)%zorl(ix), 0.317)
endif
enddo
enddo
endif
Expand All @@ -1521,7 +1523,9 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta
!$omp parallel do default(shared) private(nb, ix)
do nb = 1, Atm_block%nblks
do ix = 1, Atm_block%blksz(nb)
Sfcprop(nb)%zorli(ix) = Sfcprop(nb)%zorl(ix) !--- compute zorli from existing variables
if (Sfcprop(nb)%fice(ix)*(one-Sfcprop(nb)%landfrac(ix)) > zero) then
Sfcprop(nb)%zorli(ix) = one
endif
enddo
enddo
endif
Expand All @@ -1547,6 +1551,36 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta
enddo
endif

if (sfc_var2(i,j,47) < -9990.0_r8) then
if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing snodi')
!$omp parallel do default(shared) private(nb, ix, tem)
do nb = 1, Atm_block%nblks
do ix = 1, Atm_block%blksz(nb)
if (Sfcprop(nb)%fice(ix) > zero) then
tem = one / (Sfcprop(nb)%fice(ix)*(one-Sfcprop(nb)%landfrac(ix))+Sfcprop(nb)%landfrac(ix))
Sfcprop(nb)%snodi(ix) = min(Sfcprop(nb)%snowd(ix) * tem, 3.0)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Also here, why is landfrac included when computing snod on ice?

else
Sfcprop(nb)%snodi(ix) = zero
endif
enddo
enddo
endif

if (sfc_var2(i,j,48) < -9990.0_r8) then
if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing weasdi')
!$omp parallel do default(shared) private(nb, ix, tem)
do nb = 1, Atm_block%nblks
do ix = 1, Atm_block%blksz(nb)
if (Sfcprop(nb)%fice(ix) > zero) then
tem = one / (Sfcprop(nb)%fice(ix)*(one-Sfcprop(nb)%landfrac(ix))+Sfcprop(nb)%landfrac(ix))
Sfcprop(nb)%weasdi(ix) = Sfcprop(nb)%weasd(ix)*tem
Copy link
Collaborator

@junwang-noaa junwang-noaa Jan 10, 2022

Choose a reason for hiding this comment

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

I saw in GFS_typedefs.meta, the weasd is already defined as: "water equiv of acc snow depth over land and sea ice", why do we need to multiply "tem" and not directly set weasdi=weasd?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Jun,
In "GFS_surface_composites.F90" weasd is computed in the loop
" if (frac_grid) then

    do i=1, im

      ! Three-way composites (fields from sfc_diff)
      txl   = landfrac(i)            ! land fraction
      wfrac = one - txl              ! ocean fraction
      txi   = cice(i) * wfrac        ! txi = ice fraction wrt whole cell
      txo   = max(zero, wfrac-txi)   ! txo = open water fraction

     !gflx(i)   = txl*gflx_lnd(i)  + txi*gflx_ice(i)  + txo*gflx_wat(i)
      ep1d(i)   = txl*ep1d_lnd(i)  + txi*ep1d_ice(i)  + txo*ep1d_wat(i)
      weasd(i)  = txl*weasd_lnd(i) + txi*weasd_ice(i)
      snowd(i)  = txl*snowd_lnd(i) + txi*snowd_ice(i)
     !tprcp(i)  = txl*tprcp_lnd(i) + txi*tprcp_ice(i) + txo*tprcp_wat(i)"

So, if you set weasd_lnd=weasd/(txi+txl) and wead_ice=weasd/(txi+txl) and then use the formula above you get the right answer. i.el.
txlweasd_lnd(i) + txiweasd_ice(i) = weasd(i) * (txl+txi) / (txl+txi) = weasd(i)

Please note that this has any impact when "txl+txi) is < 1. If txl+txi=1, then there is no impact - so this affects only the fractional grid.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Moorthi, I see the changes only have impact on fraction grid when txl+txi< 1. Just for discussion. I am wondering since in the code weasd/snowd are computed as composite fields GFS_surface_composites.F90 fields over land and ice, could we just use txl/(txl+txi) and txi/(txl+txi) in weasd/snowd computation?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I am not sure I completely understand this.
For completeness,
weasd = txl * weasd_lnd + txi * weasd_ice + txo * weasd_ocn
and txl+txi+txo=1.0
However, weasd_ocn=0.0; so it is left out from the computation.
What is done in initial condition is going to be changed within one timestep by the snow imported from CICE6 when coupled.
Moorthi

else
Sfcprop(nb)%weasdi(ix) = zero
endif
enddo
enddo
endif

if (Model%use_cice_alb) then
if (sfc_var2(i,j,49) < -9990.0_r8) then
!$omp parallel do default(shared) private(nb, ix)
Expand Down Expand Up @@ -3058,7 +3092,7 @@ subroutine fv_phys_bundle_setup(Diag, axes, phys_bundle, fcst_grid, quilting, nb
!
implicit none
!
type(GFS_externaldiag_type),intent(in) :: Diag(:)
type(GFS_externaldiag_type),intent(in) :: Diag(:)
integer, intent(in) :: axes(:)
type(ESMF_FieldBundle),intent(inout) :: phys_bundle(:)
type(ESMF_Grid),intent(inout) :: fcst_grid
Expand Down Expand Up @@ -3099,7 +3133,7 @@ subroutine fv_phys_bundle_setup(Diag, axes, phys_bundle, fcst_grid, quilting, nb
!------------------------------------------------------------
!
allocate(bdl_intplmethod(nbdlphys), outputfile(nbdlphys))
if(mpp_pe()==mpp_root_pe())print *,'in fv_phys bundle,nbdl=',nbdlphys
if(mpp_pe()==mpp_root_pe()) print *,'in fv_phys bundle,nbdl=',nbdlphys
do ibdl = 1, nbdlphys
loutputfile = .false.
call ESMF_FieldBundleGet(phys_bundle(ibdl), name=physbdl_name,rc=rc)
Expand Down Expand Up @@ -3178,14 +3212,14 @@ subroutine fv_phys_bundle_setup(Diag, axes, phys_bundle, fcst_grid, quilting, nb
allocate(udimList(udimCount))
call ESMF_AttributeGet(fcst_grid, convention="NetCDF", purpose="FV3", &
name="vertical_dim_labels", valueList=udimList, rc=rc)
! if(mpp_pe()==mpp_root_pe())print *,'in fv3gfsio, vertical
! if(mpp_pe()==mpp_root_pe()) print *,'in fv3gfsio, vertical
! list=',udimList(1:udimCount),'rc=',rc

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

else

if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,axis_name_vert=',axis_name_vert
if(mpp_pe()==mpp_root_pe()) print *,'in fv_dyn bundle,axis_name_vert=',axis_name_vert
call ESMF_AttributeAdd(fcst_grid, convention="NetCDF", purpose="FV3", &
attrList=(/"vertical_dim_labels"/), rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
Expand All @@ -3207,13 +3241,13 @@ subroutine fv_phys_bundle_setup(Diag, axes, phys_bundle, fcst_grid, quilting, nb
direction, edges, Domain, DomainU, axis_data, &
num_attributes=num_attributes, attributes=attributes)
!
edgesS=''
edgesS = ''
do i = 1,num_axes_phys
if(axes(i) == edges) edgesS=axis_name(i)
enddo
! Add vertical dimension Attributes to Grid
if( id>2 ) then
! if(mpp_pe()==mpp_root_pe())print *,' in dyn add grid, axis_name=', &
! if(mpp_pe()==mpp_root_pe()) print *,' in dyn add grid, axis_name=', &
! trim(axis_name(id)),'axis_data=',axis_data
if(trim(edgesS)/='') then
call ESMF_AttributeAdd(fcst_grid, convention="NetCDF", purpose="FV3", &
Expand Down Expand Up @@ -3415,62 +3449,62 @@ subroutine add_field_to_phybundle(var_name,long_name,units,cell_methods, axes,ph
!
!*** add field attributes
call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", &
attrList=(/"long_name"/), rc=rc)
attrList=(/"long_name"/), rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
name='long_name',value=trim(long_name),rc=rc)
name='long_name',value=trim(long_name),rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", &
attrList=(/"units"/), rc=rc)
attrList=(/"units"/), rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
name='units',value=trim(units),rc=rc)
name='units',value=trim(units),rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", &
attrList=(/"missing_value"/), rc=rc)
attrList=(/"missing_value"/), rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
name='missing_value',value=missing_value,rc=rc)
name='missing_value',value=missing_value,rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", &
attrList=(/"_FillValue"/), rc=rc)
attrList=(/"_FillValue"/), rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
name='_FillValue',value=missing_value,rc=rc)
name='_FillValue',value=missing_value,rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", &
attrList=(/"cell_methods"/), rc=rc)
attrList=(/"cell_methods"/), rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
name='cell_methods',value=trim(cell_methods),rc=rc)
name='cell_methods',value=trim(cell_methods),rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
!
call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", &
attrList=(/"output_file"/), rc=rc)
attrList=(/"output_file"/), rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
name='output_file',value=trim(output_file),rc=rc)
name='output_file',value=trim(output_file),rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)

Expand Down