-
Notifications
You must be signed in to change notification settings - Fork 159
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
Changes from 4 commits
9a5e4aa
7faffc0
a8b3501
01ae48c
f292206
5fd3e1c
5d5e441
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
+17 −17 | physics/GFS_rrtmgp_lw_post.F90 | |
+30 −28 | physics/GFS_rrtmgp_sw_post.F90 | |
+7 −46 | physics/GFS_surface_composites.F90 | |
+0 −48 | physics/GFS_surface_composites.meta | |
+56 −55 | physics/cires_ugwp.F90 | |
+111 −111 | physics/micro_mg3_0.F90 | |
+4 −4 | physics/sfc_sice.f |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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)) | ||
Sfcprop(nb)%snodl(ix) = Sfcprop(nb)%snowd(ix) * tem | ||
else | ||
Sfcprop(nb)%snodl(ix) = zero | ||
|
@@ -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 | ||
|
@@ -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.317e-2) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @SMoorthi-emc Moorthi, how is this 0.317e-2 number decided ? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Fanglin, There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Well, I may have messed up the units. It may have to be multiplied by 100. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Jun, is it ok to make this unit correction now or is it too late? |
||
endif | ||
enddo | ||
enddo | ||
endif | ||
|
@@ -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 | ||
|
@@ -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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Jun,
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. 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. There was a problem hiding this comment. Choose a reason for hiding this commentThe 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? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I am not sure I completely understand this. |
||
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) | ||
|
@@ -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 | ||
|
@@ -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) | ||
|
@@ -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 | ||
|
@@ -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", & | ||
|
@@ -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) | ||
|
||
|
There was a problem hiding this comment.
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?