merge of latest dev work from GFDL Weather and Climate Dynamics Division

GFDL_tools/fv_diag_column.F90

@@ -381,7 +381,7 @@ subroutine debug_column(pt, delp, delz, u, v, w, q, npz, ncnst, sphum, nwat, zvi
        write(unit, *) '==================================================================='
        write(unit, *)
        flush(unit)
- 
+
     enddo
 
   end subroutine debug_column

model/boundary.F90

@@ -2543,7 +2543,7 @@ end subroutine fill_var_coarse
 
   subroutine update_coarse_grid_mpp_vector(u_coarse, v_coarse, u_nest, v_nest, nest_domain, dx, dy, area, &
       bd, isd_p, ied_p, jsd_p, jed_p, is_n, ie_n, js_n, je_n, &
-      isu, ieu, jsu, jeu, npx, npy, npz, istag_u, jstag_u, istag_v, jstag_v, &
+      isu, ieu, jsu, jeu, jeu_stag, iev_stag, npx, npy, npz, istag_u, jstag_u, istag_v, jstag_v, &
       r, nestupdate, upoff, nsponge, &
       parent_proc, child_proc, parent_grid, nest_level, flags, gridtype)
 
@@ -2553,7 +2553,7 @@ subroutine update_coarse_grid_mpp_vector(u_coarse, v_coarse, u_nest, v_nest, nes
 
    type(fv_grid_bounds_type), intent(IN) :: bd
    integer, intent(IN) :: isd_p, ied_p, jsd_p, jed_p, is_n, ie_n, js_n, je_n
-   integer, intent(IN) :: isu, ieu, jsu, jeu
+   integer, intent(IN) :: isu, ieu, jsu, jeu, jeu_stag, iev_stag
    integer, intent(IN) :: istag_u, jstag_u, istag_v, jstag_v
    integer, intent(IN) :: npx, npy, npz, r, nestupdate, upoff, nsponge
    real, intent(IN)    :: u_nest(is_n:ie_n+istag_u,js_n:je_n+jstag_u,npz)
@@ -2614,13 +2614,14 @@ subroutine update_coarse_grid_mpp_vector(u_coarse, v_coarse, u_nest, v_nest, nes
    s = r/2 !rounds down (since r > 0)
    qr = r*upoff + nsponge - s
 
-   if (parent_proc .and. .not. (ieu < isu .or. jeu < jsu)) then
+   if (parent_proc .and. .not. (ieu < isu .or. jeu_stag < jsu)) then
       call fill_var_coarse(u_coarse, coarse_dat_recv_u, isd_p, ied_p, jsd_p, jed_p, &
-           isu, ieu, jsu, jeu, npx, npy, npz, istag_u, jstag_u, nestupdate, parent_grid)
+           isu, ieu, jsu, jeu_stag, npx, npy, npz, istag_u, jstag_u, nestupdate, parent_grid)
    endif
-   if (parent_proc .and. .not. (ieu < isu .or. jeu < jsu)) then
+
+   if (parent_proc .and. .not. (iev_stag < isu .or. jeu < jsu)) then
       call fill_var_coarse(v_coarse, coarse_dat_recv_v, isd_p, ied_p, jsd_p, jed_p, &
-           isu, ieu, jsu, jeu, npx, npy, npz, istag_v, jstag_v, nestupdate, parent_grid)
+           isu, iev_stag, jsu, jeu, npx, npy, npz, istag_v, jstag_v, nestupdate, parent_grid)
    endif
 
    if (allocated(coarse_dat_recv_u)) deallocate(coarse_dat_recv_u)

model/fv_arrays.F90

@@ -878,6 +878,11 @@ module fv_arrays_mod
 
    integer :: bc_update_interval = 3   !< Default setting for interval (hours) between external regional BC data files.
 
+  integer :: nrows_blend = 0          !< # of blending rows in the outer integration domain.
+  logical :: write_restart_with_bcs = .false.   !< Default setting for using DA-updated BC files
+  logical :: regional_bcs_from_gsi = .false.    !< Default setting for writing restart files with boundary rows
+
+
   !>Convenience pointers
   integer, pointer :: grid_number
 
@@ -959,8 +964,12 @@ module fv_arrays_mod
      integer :: npx_global
      integer :: upoff = 1 !< currently the same for all variables
      integer :: isu = -999, ieu = -1000, jsu = -999, jeu = -1000 !< limits of update regions on coarse grid
+     integer :: jeu_stag = -1000, iev_stag = -1000 !< limits of update regions on coarse grid for staggered variables in j,i
+     integer :: jeu_stag_boundary = -1000, iev_stag_boundary = -1000 !< BC location
+
      real    :: update_blend = 1. !< option for controlling how much "blending" is done during two-way update
      logical, allocatable :: do_remap_BC(:)
+     logical, allocatable :: do_remap_BC_level(:)
 
      !nest_domain now a global structure defined in fv_mp_mod
      !type(nest_domain_type) :: nest_domain !Structure holding link from this grid to its parent
@@ -1318,7 +1327,7 @@ module fv_arrays_mod
 !>@details It includes an option to define dummy grids that have scalar and
 !! small arrays defined as null 3D arrays.
   subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie_in, js_in, je_in, &
-       npx_in, npy_in, npz_in, ndims_in, ncnst_in, nq_in, dummy, alloc_2d, ngrids_in)
+       npx_in, npy_in, npz_in, ndims_in, ntiles_in, ncnst_in, nq_in, dummy, alloc_2d, ngrids_in)
 
     !WARNING: Before calling this routine, be sure to have set up the
     ! proper domain parameters from the namelists (as is done in
@@ -1327,7 +1336,7 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie
     implicit none
     type(fv_atmos_type), intent(INOUT), target :: Atm
     integer, intent(IN) :: isd_in, ied_in, jsd_in, jed_in, is_in, ie_in, js_in, je_in
-    integer, intent(IN) :: npx_in, npy_in, npz_in, ndims_in, ncnst_in, nq_in
+    integer, intent(IN) :: npx_in, npy_in, npz_in, ndims_in, ntiles_in, ncnst_in, nq_in
     logical, intent(IN) :: dummy, alloc_2d
     integer, intent(IN) :: ngrids_in
     integer:: isd, ied, jsd, jed, is, ie, js, je
@@ -1754,11 +1763,14 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie
     if( ngrids_in > 1 ) then
        if (Atm%flagstruct%grid_type < 4) then
           if (Atm%neststruct%nested) then
-             allocate(Atm%grid_global(1-Atm%ng:npx_2d  +Atm%ng,1-Atm%ng:npy_2d  +Atm%ng,2,1))
+             allocate(Atm%grid_global(1-Atm%ng:npx_2d  +Atm%ng,1-Atm%ng:npy_2d  +Atm%ng,2,ntiles_in))
           else
-             allocate(Atm%grid_global(1-Atm%ng:npx_2d  +Atm%ng,1-Atm%ng:npy_2d  +Atm%ng,2,1:6))
+             allocate(Atm%grid_global(1-Atm%ng:npx_2d  +Atm%ng,1-Atm%ng:npy_2d  +Atm%ng,2,1:ntiles_in))
           endif
        end if
+       if (Atm%flagstruct%grid_type == 4) then
+          allocate(Atm%grid_global(1-Atm%ng:npx_2d  +Atm%ng,1-Atm%ng:npy_2d  +Atm%ng,2,ntiles_in))
+       endif
     endif
 
 

model/fv_cmp.F90

@@ -37,7 +37,7 @@ module fv_cmp_mod
 
     private
 
-    public fv_sat_adj, qs_init
+    public fv_sat_adj, qs_init, c_ice, c_liq
 
     ! real, parameter :: cp_air = cp_air ! 1004.6, heat capacity of dry air at constant pressure, come from constants_mod
     real, parameter :: cp_vap = 4.0 * rvgas ! 1846.0, heat capacity of water vapor at constant pressure

model/fv_control.F90

@@ -563,7 +563,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split)
              Atm(n)%bd%isc, Atm(n)%bd%iec, &
              Atm(n)%bd%jsc, Atm(n)%bd%jec, &
              Atm(n)%flagstruct%npx,    Atm(n)%flagstruct%npy,   Atm(n)%flagstruct%npz, &
-             Atm(n)%flagstruct%ndims,  Atm(n)%flagstruct%ncnst, Atm(n)%flagstruct%ncnst-Atm(n)%flagstruct%pnats, &
+             Atm(n)%flagstruct%ndims, Atm(n)%flagstruct%ntiles,  Atm(n)%flagstruct%ncnst, Atm(n)%flagstruct%ncnst-Atm(n)%flagstruct%pnats, &
              n/=this_grid, n==this_grid, ngrids) !TODO don't need both of the last arguments
      enddo
      if ( (Atm(this_grid)%bd%iec-Atm(this_grid)%bd%isc+1).lt.4 .or. (Atm(this_grid)%bd%jec-Atm(this_grid)%bd%jsc+1).lt.4 ) then
@@ -594,7 +594,9 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split)
      do n=1,ngrids
         Atm(this_grid)%neststruct%child_grids(n) = (grid_coarse(n) == this_grid)
         allocate(Atm(n)%neststruct%do_remap_bc(ngrids))
+        allocate(Atm(n)%neststruct%do_remap_bc_level(Atm(n)%neststruct%num_nest_level))
         Atm(n)%neststruct%do_remap_bc(:) = .false.
+        Atm(n)%neststruct%do_remap_bc_level(:) = .false.
      enddo
      Atm(this_grid)%neststruct%parent_proc = ANY(Atm(this_grid)%neststruct%child_grids) !ANY(tile_coarse == Atm(this_grid)%global_tile)
      Atm(this_grid)%neststruct%child_proc = ASSOCIATED(Atm(this_grid)%parent_grid) !this means a nested grid
@@ -886,7 +888,7 @@ subroutine read_namelist_fv_grid_nml
        ! Read Main namelist
        read (input_nml_file,fv_grid_nml,iostat=ios)
        ierr = check_nml_error(ios,'fv_grid_nml')
- 
+
        call write_version_number ( 'FV_CONTROL_MOD', version )
        unit = stdlog()
        write(unit, nml=fv_grid_nml)
@@ -947,7 +949,7 @@ subroutine read_namelist_fv_core_nml(Atm)
        ierr = check_nml_error(ios,'fv_core_nml')
        ! Reset input_file_nml to default behavior (CHECK do we still need this???)
        !call read_input_nml
- 
+
        call write_version_number ( 'FV_CONTROL_MOD', version )
        unit = stdlog()
        write(unit, nml=fv_core_nml)
@@ -1058,9 +1060,11 @@ end subroutine read_namelist_fv_core_nml
 
      subroutine setup_update_regions
 
-       integer :: isu, ieu, jsu, jeu ! update regions
+       integer :: isu, ieu, jsu, jeu ! update regions for centered variables
        integer :: isc, jsc, iec, jec
        integer :: upoff
+       integer :: isu_stag, jsu_stag, ieu_stag, jeu_stag ! update regions for u
+       integer :: isv_stag, jsv_stag, iev_stag, jev_stag ! update regions for v
 
        isc = Atm(this_grid)%bd%isc
        jsc = Atm(this_grid)%bd%jsc
@@ -1078,35 +1082,108 @@ subroutine setup_update_regions
              jsu = nest_joffsets(n)
              jeu = jsu + jcount_coarse(n) - 1
 
+             isu_stag = isu
+             jsu_stag = jsu
+             ieu_stag = ieu
+             jeu_stag = jeu+1
+
+             isv_stag = isu
+             jsv_stag = jsu
+             iev_stag = ieu+1
+             jev_stag = jeu
+
              !update offset adjustment
              isu = isu + upoff
              ieu = ieu - upoff
              jsu = jsu + upoff
              jeu = jeu - upoff
 
+             isu_stag = isu_stag + upoff
+             ieu_stag = ieu_stag - upoff
+             jsu_stag = jsu_stag + upoff
+             jeu_stag = jeu_stag - upoff
+
+             isv_stag = isv_stag + upoff
+             iev_stag = iev_stag - upoff
+             jsv_stag = jsv_stag + upoff
+             jev_stag = jev_stag - upoff
+
              !restriction to current domain
 !!$             !!! DEBUG CODE
 !!$             if (Atm(this_grid)%flagstruct%fv_debug) then
-!!$                write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS  : ', isu, jsu, ieu, jeu
-!!$                write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 2: ', isc, jsc, iec, jsc
+!!$                write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 1: ', isu, jsu, ieu, jeu
+!!$                write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 11: ', isu_stag, jsu_stag, ieu_stag, jeu_stag
+!!$                write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 111: ', isv_stag, jsv_stag, iev_stag, jev_stag
+!!$                write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 2: ', isc, jsc, iec, jec
 !!$             endif
 !!$             !!! END DEBUG CODE
+
+! Absolute boundary for the staggered point update region on the parent.
+! This is used in remap_uv to control the update of the last staggered point
+! when the the update region coincides with a pe domain to avoid cross-restart repro issues
+
+             Atm(n)%neststruct%jeu_stag_boundary = jeu_stag
+             Atm(n)%neststruct%iev_stag_boundary = iev_stag
+
+             if (isu > iec .or. ieu < isc .or. &
+                 jsu > jec .or. jeu < jsc ) then
+                isu = -999 ; jsu = -999 ; ieu = -1000 ; jeu = -1000
+             else
+                isu = max(isu,isc) ; jsu = max(jsu,jsc)
+                ieu = min(ieu,iec) ; jeu = min(jeu,jec)
+             endif
+
+! Update region for staggered quantity to avoid cross repro issues when the pe domain boundary
+! coincide with the nest. Basically write the staggered update on compute domains
+
+             if (isu_stag > iec .or. ieu_stag < isc .or. &
+                 jsu_stag > jec .or. jeu_stag < jsc ) then
+                isu_stag = -999 ; jsu_stag = -999 ; ieu_stag = -1000 ; jeu_stag = -1000
+             else
+                isu_stag = max(isu_stag,isc) ; jsu_stag = max(jsu_stag,jsc)
+                ieu_stag = min(ieu_stag,iec) ; jeu_stag = min(jeu_stag,jec)
+             endif
+
+             if (isv_stag > iec .or. iev_stag < isc .or. &
+                 jsv_stag > jec .or. jev_stag < jsc ) then
+                isv_stag = -999 ; jsv_stag = -999 ; iev_stag = -1000 ; jev_stag = -1000
+             else
+                isv_stag = max(isv_stag,isc) ; jsv_stag = max(jsv_stag,jsc)
+                iev_stag = min(iev_stag,iec) ; jev_stag = min(jev_stag,jec)
+             endif
              if (isu > iec .or. ieu < isc .or. &
                  jsu > jec .or. jeu < jsc ) then
                 isu = -999 ; jsu = -999 ; ieu = -1000 ; jeu = -1000
              else
                 isu = max(isu,isc) ; jsu = max(jsu,jsc)
                 ieu = min(ieu,iec) ; jeu = min(jeu,jec)
              endif
+
+             ! lump indices
+             isu=max(isu, isu_stag, isv_stag)
+             jsu=max(jsu, jsu_stag, jsv_stag)
+             jeu_stag=max(jeu, jeu_stag)
+             jev_stag=max(jeu, jev_stag)
+             ieu_stag=max(ieu ,ieu_stag)
+             iev_stag=max(ieu ,iev_stag)
+
 !!$             !!! DEBUG CODE
 !!$             if (Atm(this_grid)%flagstruct%fv_debug) &
-!!$                  write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 3: ', isu, jsu, ieu, jeu
+!!$                 write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 3: ', isu, jsu, ieu, jeu
+!!$                 write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 4: ', jeu_stag, iev_stag, ieu_stag, jev_stag
 !!$             !!! END DEBUG CODE
 
              Atm(n)%neststruct%isu = isu
-             Atm(n)%neststruct%ieu = ieu
+             Atm(n)%neststruct%ieu = ieu_stag
              Atm(n)%neststruct%jsu = jsu
-             Atm(n)%neststruct%jeu = jeu
+             Atm(n)%neststruct%jeu = jev_stag
+
+             Atm(n)%neststruct%jeu_stag = jeu_stag
+             Atm(n)%neststruct%iev_stag = iev_stag
+
+!!$                 write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 5: ', Atm(n)%neststruct%isu, Atm(n)%neststruct%jsu, Atm(n)%neststruct%ieu, Atm(n)%neststruct%jeu
+!!$                 write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 6: ', Atm(n)%neststruct%jeu_stag, Atm(n)%neststruct%iev_stag
+
           endif
        enddo
 

model/fv_dynamics.F90

@@ -227,7 +227,14 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot,  ng, bdt, consv_te, fill,
 
         reg_bc_update_time=current_time_in_seconds
         call set_regional_BCs          & !<-- Insert values into the boundary region valid for the start of this large timestep.
-              (delp,delz,w,pt,q_con,cappa,q,u,v,uc,vc, bd, npz, ncnst, reg_bc_update_time )
+              (delp,delz,w,pt                                     &
+#ifdef USE_COND
+               ,q_con                                             &
+#endif
+#ifdef MOIST_CAPPA
+               ,cappa                                             &
+#endif
+               ,q,u,v,uc,vc, bd, npz,  reg_bc_update_time )
 
         call timing_off('Regional_BCs')
       endif