Update fms_mixedmode branch of NOAA-EMC/GFDL_atmos_cubed_sphere (#73)

driver/fvGFS/atmosphere.F90

@@ -269,7 +269,7 @@ module atmosphere_mod
   logical :: cold_start = .false.     !  used in initial condition
 
   integer, dimension(:), allocatable :: id_tracerdt_dyn
-  integer :: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, cld_amt  ! condensate species tracer indices
+  integer :: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat, cld_amt  ! condensate species tracer indices
 
   integer :: mygrid = 1
   integer :: p_split = 1
@@ -328,10 +328,12 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
    logical :: dycore_only  = .false.
    logical :: debug        = .false.
    logical :: sync         = .false.
+   logical :: ignore_rst_cksum = .false.
    integer, parameter     :: maxhr = 4096
    real, dimension(maxhr) :: fdiag = 0.
    real                   :: fhmax=384.0, fhmaxhf=120.0, fhout=3.0, fhouthf=1.0,avg_max_length=3600.
-   namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, avg_max_length
+   namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, &
+                            & avg_max_length, ignore_rst_cksum
    ! *DH 20210326
 
    !For regional
@@ -361,7 +363,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 #ifdef MOVING_NEST
    call fv_tracker_init(size(Atm))
    if (mygrid .eq. 2) call allocate_tracker(mygrid, Atm(mygrid)%bd%isc, Atm(mygrid)%bd%iec, Atm(mygrid)%bd%jsc, Atm(mygrid)%bd%jec)
-#endif   
+#endif
 
    Atm(mygrid)%Time_init = Time_init
 
@@ -407,9 +409,10 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
    rainwat = get_tracer_index (MODEL_ATMOS, 'rainwat' )
    snowwat = get_tracer_index (MODEL_ATMOS, 'snowwat' )
    graupel = get_tracer_index (MODEL_ATMOS, 'graupel' )
+   hailwat = get_tracer_index (MODEL_ATMOS, 'hailwat' )
    cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt')
 
-   if (max(sphum,liq_wat,ice_wat,rainwat,snowwat,graupel) > Atm(mygrid)%flagstruct%nwat) then
+   if (max(sphum,liq_wat,ice_wat,rainwat,snowwat,graupel,hailwat) > Atm(mygrid)%flagstruct%nwat) then
       call mpp_error (FATAL,' atmosphere_init: condensate species are not first in the list of &
                             &tracers defined in the field_table')
    endif
@@ -453,6 +456,15 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 !--- allocate pref
    allocate(pref(npz+1,2), dum1d(npz+1))
 
+   ! DH* 20210326
+   ! First, read atmos_model_nml namelist section - this is a workaround to avoid
+   ! unnecessary additional changes to the input namelists, in anticipation of the
+   ! implementation of a generic interface for GFDL and CCPP fast physics soon
+   read(input_nml_file, nml=atmos_model_nml, iostat=io)
+   ierr = check_nml_error(io, 'atmos_model_nml')
+   !write(0,'(a)') "It's me, and my physics suite is '" // trim(ccpp_suite) // "'"
+   ! *DH 20210326
+
    call fv_restart(Atm(mygrid)%domain, Atm, seconds, days, cold_start, Atm(mygrid)%gridstruct%grid_type, mygrid)
 
    fv_time = Time
@@ -494,15 +506,6 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 
    ! Do CCPP fast physics initialization before call to adiabatic_init (since this calls fv_dynamics)
 
-   ! DH* 20210326
-   ! First, read atmos_model_nml namelist section - this is a workaround to avoid
-   ! unnecessary additional changes to the input namelists, in anticipation of the
-   ! implementation of a generic interface for GFDL and CCPP fast physics soon
-   read(input_nml_file, nml=atmos_model_nml, iostat=io)
-   ierr = check_nml_error(io, 'atmos_model_nml')
-   !write(0,'(a)') "It's me, and my physics suite is '" // trim(ccpp_suite) // "'"
-   ! *DH 20210326
-
    ! For fast physics running over the entire domain, block
    ! and thread number are not used; set to safe values
    cdata%blk_no = 1
@@ -957,10 +960,10 @@ end subroutine get_nth_domain_info
 !! the "domain2d" variable associated with the coupling grid and the
 !! decomposition for the current cubed-sphere tile.
 !>@detail Coupling is done using the mass/temperature grid with no halos.
- subroutine atmosphere_domain ( fv_domain, layout, regional, nested, &
+ subroutine atmosphere_domain ( fv_domain, rd_domain, layout, regional, nested, &
                                 moving_nest_parent, is_moving_nest, &
                                 ngrids_atmos, mygrid_atmos, pelist )
-   type(domain2d), intent(out) :: fv_domain
+   type(domain2d), intent(out) :: fv_domain, rd_domain
    integer, intent(out) :: layout(2)
    logical, intent(out) :: regional
    logical, intent(out) :: nested
@@ -973,6 +976,7 @@ subroutine atmosphere_domain ( fv_domain, layout, regional, nested, &
    integer :: n
 
    fv_domain = Atm(mygrid)%domain_for_coupler
+   rd_domain = Atm(mygrid)%domain_for_read
    layout(1:2) =  Atm(mygrid)%layout(1:2)
    regional = Atm(mygrid)%flagstruct%regional
    nested = ngrids > 1

model/dyn_core.F90

@@ -1248,7 +1248,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
 
                                                      call timing_on('COMM_TOTAL')
 #ifndef ROT3
-    if ( it/=n_split)   &
+    if ( .not. flagstruct%regional .and. it/=n_split)   &
          call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE)
 #endif
                                                      call timing_off('COMM_TOTAL')
@@ -1351,7 +1351,10 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
                                        isd, ied, jsd, jed,      &
                                        reg_bc_update_time,it )
 
-         call mpp_update_domains(u, v, domain, gridtype=DGRID_NE)
+#ifndef ROT3
+         if (it/=n_split)   &
+            call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE)
+#endif
 
     endif
 

model/fv_arrays.F90

@@ -516,6 +516,7 @@ module fv_arrays_mod
 !-----------------------------------------------------------------------------------------------
 
    logical :: reset_eta = .false.
+   logical :: ignore_rst_cksum = .false. !< enfore (.false.) or override (.true.) data integrity restart checksums
    real    :: p_fac = 0.05  !< Safety factor for minimum nonhydrostatic pressures, which
                             !< will be limited so the full pressure is no less than p_fac
                             !< times the hydrostatic pressure. This is only of concern in mid-top
@@ -1303,6 +1304,7 @@ module fv_arrays_mod
 #if defined(SPMD)
 
      type(domain2D) :: domain_for_coupler !< domain used in coupled model with halo = 1.
+     type(domain2D) :: domain_for_read    !< domain used for reads to increase performance when io_layout=(1,1)
 
      !global tile and tile_of_mosaic only have a meaning for the CURRENT pe
      integer :: num_contact, npes_per_tile, global_tile, tile_of_mosaic, npes_this_grid

model/fv_control.F90

@@ -198,7 +198,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
      real,                intent(in)    :: dt_atmos
      integer,             intent(OUT)   :: this_grid
      logical, allocatable, intent(OUT) :: grids_on_this_pe(:)
-     character(len=32), optional,      intent(in)    :: nml_filename_in ! alternate nml 
+     character(len=32), optional,      intent(in)    :: nml_filename_in ! alternate nml
      logical, optional,                intent(in)    :: skip_nml_read_in ! use previously loaded nml
 
      integer, intent(INOUT) :: p_split
@@ -291,6 +291,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
      real(kind=R_GRID) , pointer :: target_lon
 
      logical , pointer :: reset_eta
+     logical , pointer :: ignore_rst_cksum
      real    , pointer :: p_fac
      real    , pointer :: a_imp
      integer , pointer :: n_split
@@ -676,7 +677,8 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
           Atm(this_grid)%flagstruct%grid_type,Atm(this_grid)%neststruct%nested, &
           Atm(this_grid)%layout,Atm(this_grid)%io_layout,Atm(this_grid)%bd,Atm(this_grid)%tile_of_mosaic, &
           Atm(this_grid)%gridstruct%square_domain,Atm(this_grid)%npes_per_tile,Atm(this_grid)%domain, &
-          Atm(this_grid)%domain_for_coupler,Atm(this_grid)%num_contact,Atm(this_grid)%pelist)
+          Atm(this_grid)%domain_for_coupler,Atm(this_grid)%domain_for_read,Atm(this_grid)%num_contact, &
+          Atm(this_grid)%pelist)
      call broadcast_domains(Atm,Atm(this_grid)%pelist,size(Atm(this_grid)%pelist))
      do n=1,ngrids
         tile_id = mpp_get_tile_id(Atm(n)%domain)
@@ -857,6 +859,7 @@ subroutine set_namelist_pointers(Atm)
        regional_bcs_from_gsi         => Atm%flagstruct%regional_bcs_from_gsi
        write_restart_with_bcs        => Atm%flagstruct%write_restart_with_bcs
        reset_eta                     => Atm%flagstruct%reset_eta
+       ignore_rst_cksum              => Atm%flagstruct%ignore_rst_cksum
        p_fac                         => Atm%flagstruct%p_fac
        a_imp                         => Atm%flagstruct%a_imp
        n_split                       => Atm%flagstruct%n_split
@@ -1075,7 +1078,7 @@ subroutine read_namelist_fv_core_nml(Atm)
             write_coarse_restart_files,&
             write_coarse_diagnostics,&
             write_only_coarse_intermediate_restarts, &
-            write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst
+            write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum
 
 
        ! Read FVCORE namelist

model/fv_regional_bc.F90

@@ -410,6 +410,17 @@ subroutine setup_regional_BC(Atm                                  &
       else
         nrows_blend=nrows_blend_in_data                                    !<-- # of blending rows in the BC files.
       endif
+
+      IF ( north_bc .or. south_bc ) THEN
+        IF ( nrows_blend_user > jed - nhalo_model - (jsd + nhalo_model) + 1 ) THEN
+        call mpp_error(FATAL,'Number of blending rows is greater than the north-south tile size!')
+        ENDIF
+      ENDIF
+      IF ( west_bc .or. east_bc ) THEN
+        IF ( nrows_blend_user > ied - nhalo_model - (isd + nhalo_model) + 1 ) THEN
+        call mpp_error(FATAL,'Number of blending rows is greater than the east-west tile size!')
+        ENDIF
+      ENDIF
 !
       call check(nf90_close(ncid))                                         !<-- Close the BC file for now.
 !
@@ -4352,7 +4363,7 @@ subroutine regional_boundary_update(array                       &
 !
       real,dimension(:,:,:),pointer :: bc_t0,bc_t1                       !<-- Boundary data at the two bracketing times.
 !
-      logical :: blend,call_interp
+      logical :: blend,call_interp,blendtmp
 !
 !---------------------------------------------------------------------
 !*********************************************************************
@@ -4396,13 +4407,21 @@ subroutine regional_boundary_update(array                       &
               i2=ied+1
             endif
 !
-            j1=jsd
-            j2=js-1
+            j1=jsd  ! -2  -- outermost boundary ghost zone
+            j2=js-1 ! 0  -- first boundary ghost zone
 !
+            IF ( east_bc ) THEN
             i1_blend=is
+            ELSE
+            i1_blend=isd !is-nhalo_model
+            ENDIF
+            IF ( west_bc ) THEN
             i2_blend=ie
+            ELSE
+            i2_blend=ied ! ie+nhalo_model
+            ENDIF
             if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then
-              i2_blend=ie+1
+              i2_blend=i2_blend+1 ! ie+1
             endif
             j1_blend=js
             j2_blend=js+nrows_blend_user-1
@@ -4437,8 +4456,19 @@ subroutine regional_boundary_update(array                       &
 !
             i1_blend=is
             i2_blend=ie
-            if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then
-              i2_blend=ie+1
+            IF ( east_bc ) THEN
+            i1_blend=is
+            ELSE
+            i1_blend=isd !is-nhalo_model
+            ENDIF
+            IF ( west_bc ) THEN
+            i2_blend=ie
+            ELSE
+            i2_blend=ied ! ie+nhalo_model
+            ENDIF
+           if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then
+!              i2_blend=ie+1
+              i2_blend=i2_blend+1
             endif
             j2_blend=je
             if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then
@@ -4483,16 +4513,21 @@ subroutine regional_boundary_update(array                       &
               endif
             endif
 !
-            i1_blend=is
-            i2_blend=is+nrows_blend_user-1
+! Note: Original code checked for corner region and avoided overlap, but changed this to blend corners from both boundaries
+             i1_blend=is
+             i2_blend=is+nrows_blend_user-1
+
+            IF ( north_bc ) THEN
             j1_blend=js
+            ELSE
+            j1_blend=jsd !js-nhalo_model
+            ENDIF
+            IF ( south_bc ) THEN
             j2_blend=je
-            if(north_bc)then
-              j1_blend=js+nrows_blend_user     !<-- North BC already handles nrows_blend_user blending rows
-            endif
-            if(south_bc)then
-              j2_blend=je-nrows_blend_user     !<-- South BC already handles nrows_blend_user blending rows
-            endif
+            ELSE
+            j2_blend=jed ! ie+nhalo_model
+            ENDIF
+
             if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then
               j2_blend=j2_blend+1
             endif
@@ -4538,16 +4573,20 @@ subroutine regional_boundary_update(array                       &
               endif
             endif
 !
+! Note: Original code checked for corner region and avoided overlap, but changed this to blend corners from both boundaries
             i1_blend=i1-nrows_blend_user
             i2_blend=i1-1
+
+            IF ( north_bc ) THEN
             j1_blend=js
+            ELSE
+            j1_blend=jsd !is-nhalo_model
+            ENDIF
+            IF ( south_bc ) THEN
             j2_blend=je
-            if(north_bc)then
-              j1_blend=js+nrows_blend_user   !<-- North BC already handled nrows_blend_user blending rows.
-            endif
-            if(south_bc)then
-              j2_blend=je-nrows_blend_user   !<-- South BC already handled nrows_blend_user blending rows.
-            endif
+            ELSE
+            j2_blend=jed ! ie+nhalo_model
+            ENDIF
             if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then
               j2_blend=j2_blend+1
             endif
@@ -4563,6 +4602,8 @@ subroutine regional_boundary_update(array                       &
 !***  then update the boundary points.
 !---------------------------------------------------------------------
 !
+
+
         if(call_interp)then
 !
           call retrieve_bc_variable_data(bc_vbl_name                  &
@@ -4787,7 +4828,7 @@ subroutine bc_time_interpolation(array                               &
 
 !
 !---------------------------------------------------------------------
-!
+! Set values in the boundary points only
       do k=1,ubnd_z
         do j=j1,j2
         do i=i1,i2
@@ -4818,7 +4859,7 @@ subroutine bc_time_interpolation(array                               &
 !-----------
 !
       if(nside==1.and.north_bc)then
-        rdenom=1./real(j2_blend-j_bc-1)
+        rdenom=1./real(Max(1,j2_blend-j_bc-1))
         do k=1,ubnd_z
           do j=j1_blend,j2_blend
             factor_dist=exp(-(blend_exp1+blend_exp2*(j-j_bc-1)*rdenom)) !<-- Exponential falloff of blending weights.
@@ -4837,7 +4878,7 @@ subroutine bc_time_interpolation(array                               &
 !-----------
 !
       if(nside==2.and.south_bc)then
-        rdenom=1./real(j_bc-j1_blend-1)
+        rdenom=1./real(Max(1,j_bc-j1_blend-1))
         do k=1,ubnd_z
           do j=j1_blend,j2_blend
             factor_dist=exp(-(blend_exp1+blend_exp2*(j_bc-j-1)*rdenom)) !<-- Exponential falloff of blending weights.
@@ -4855,7 +4896,7 @@ subroutine bc_time_interpolation(array                               &
 !----------
 !
       if(nside==3.and.east_bc)then
-        rdenom=1./real(i2_blend-i_bc-1)
+        rdenom=1./real(Max(1,i2_blend-i_bc-1))
         do k=1,ubnd_z
           do j=j1_blend,j2_blend
             do i=i1_blend,i2_blend
@@ -4876,7 +4917,7 @@ subroutine bc_time_interpolation(array                               &
 !----------
 !
       if(nside==4.and.west_bc)then
-        rdenom=1./real(i_bc-i1_blend-1)
+        rdenom=1./real(Max(1, i_bc-i1_blend-1))
         do k=1,ubnd_z
           do j=j1_blend,j2_blend
             do i=i1_blend,i2_blend