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

* modify FV3GFS_io.F90 by fixing errors associated with initializing snow depth over ice in the case where both land and water coexist (i.e. fractional grid case) * z0ice is changed to 1.0 cm from 1.1cm in atmos_model.F90
NOAA-EMC · Jan 10, 2022 · c2f06d5 · c2f06d5
1 parent 7eb0ee2
commit c2f06d5
Show file tree

Hide file tree

Showing 4 changed files with 71 additions and 24 deletions.
diff --git a/atmos_model.F90 b/atmos_model.F90
@@ -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

diff --git a/ccpp/physics b/ccpp/physics
diff --git a/cpl/module_block_data.F90 b/cpl/module_block_data.F90
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)
@@ -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)

diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90
@@ -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.317)
+          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)
+          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
+          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)
+17 −17		physics/GFS_rrtmgp_lw_post.F90
+30 −28		physics/GFS_rrtmgp_sw_post.F90
+18 −17		physics/GFS_surface_composites.F90
+56 −55		physics/cires_ugwp.F90
+111 −111		physics/micro_mg3_0.F90
+4 −4		physics/sfc_sice.f