Merge pull request #13 from NOAA-GFDL/dev/emc

Regional decomposition test fix (when nrows_blend > 0) (NOAA-GFDL#194)
JiliDong-NOAA · Jun 23, 2022 · 949ca85 · 949ca85
2 parents b642410 + 49f30ee
commit 949ca85
Show file tree

Hide file tree

Showing 4 changed files with 88 additions and 31 deletions.
diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90
@@ -265,7 +265,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
@@ -403,9 +403,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

diff --git a/model/dyn_core.F90 b/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
 

diff --git a/model/fv_regional_bc.F90 b/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

diff --git a/model/fv_sg.F90 b/model/fv_sg.F90
@@ -404,12 +404,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt,    &
             qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)
          enddo
       enddo
-   else
+   elseif ( nwat==6 ) then
       do k=1,kbot
          do i=is,ie
             qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)
          enddo
       enddo
+   elseif ( nwat==7 ) then
+      do k=1,kbot
+         do i=is,ie
+            qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)+q0(i,k,hailwat)
+         enddo
+      enddo
    endif
 
       do k=kbot, 2, -1
@@ -985,12 +991,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt,    &
          do i=is,ie
             qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)
          enddo
-   else
+   elseif ( nwat==6 ) then
       do k=1,kbot
          do i=is,ie
             qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)
          enddo
       enddo
+   elseif ( nwat==7 ) then
+      do k=1,kbot
+         do i=is,ie
+            qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)+q0(i,k,hailwat)
+         enddo
+      enddo
    endif
 
       do k=kbot, 2, -1