Update fms_mixedmode to gfdl_atmos_cubed_sphere/dev/emc

CMakeLists.txt

@@ -59,6 +59,7 @@ endif()
 list(APPEND model_srcs
   model/a2b_edge.F90
   model/multi_gases.F90
+  model/molecular_diffusion.F90
   model/boundary.F90
   model/dyn_core.F90
   model/fv_arrays.F90
@@ -169,6 +170,7 @@ add_library(FV3::fv3 ALIAS fv3)
 
 set_property(SOURCE model/nh_utils.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "${FAST}")
 set_property(SOURCE model/fv_mapz.F90  APPEND_STRING PROPERTY COMPILE_FLAGS "${FAST}")
+set_property(SOURCE tools/fv_diagnostics.F90  APPEND_STRING PROPERTY COMPILE_FLAGS "-O0")
 
 set_target_properties(fv3 PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/include/fv3)
 target_include_directories(fv3 INTERFACE $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/include/fv3>

driver/fvGFS/fv_nggps_diag.F90

@@ -122,9 +122,7 @@ module fv_nggps_diags_mod
  logical :: use_wrtgridcomp_output=.false.
  integer :: sphum, liq_wat, ice_wat       !< GFDL physics
  integer :: rainwat, snowwat, graupel
- real :: vrange(2) = (/ -330.,  330. /)  !< winds
- real :: wrange(2) = (/ -100.,  100. /)  !< vertical wind
- real :: trange(2) = (/  100.,  350. /)  !< temperature
+ real :: vrange(2), wrange(2), trange(2)
  real :: skrange(2) = (/ -10000000.0,  10000000.0 /)  !< dissipation estimate for SKEB
 
 ! file name
@@ -181,6 +179,20 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
     kstt_tracer(:) = 0
     kend_tracer(:) = 0
 
+    if ( Atm(n)%flagstruct%molecular_diffusion ) then
+         vrange = (/ -830.,  830. /)  !< winds for wam
+         wrange = (/ -300.,  300. /)  !< vertical wind for wam
+         trange = (/  100., 3000. /)  !< temperature for wam
+    else
+         vrange = (/ -330.,  330. /)  !< winds
+         wrange = (/ -100.,  100. /)  !< vertical wind
+#ifdef MULTI_GASES
+         trange = (/  100., 3000. /)  !< temperature
+#else
+         trange = (/  100.,  350. /)  !< temperature
+#endif
+    endif
+
     if (Atm(n)%flagstruct%write_3d_diags) then
 !-------------------
 ! A grid winds (lat-lon)

model/dyn_core.F90

@@ -111,7 +111,9 @@ module dyn_core_mod
   use fv_mp_mod,          only: is_master
   use fv_mp_mod,          only: start_group_halo_update, complete_group_halo_update
   use fv_mp_mod,          only: group_halo_update_type
-  use sw_core_mod,        only: c_sw, d_sw
+  use molecular_diffusion_mod,       &
+                          only: md_time, md_layers, md_consv_te, md_tadj_layers
+  use sw_core_mod,        only: c_sw, d_sw, d_md
   use a2b_edge_mod,       only: a2b_ord2, a2b_ord4
   use nh_core_mod,        only: Riem_Solver3, Riem_Solver_C, update_dz_c, update_dz_d, nh_bc
   use tp_core_mod,        only: copy_corners
@@ -139,7 +141,7 @@ module dyn_core_mod
   use test_cases_mod,      only: test_case, case9_forcing1, case9_forcing2
 #endif
 #ifdef MULTI_GASES
-    use multi_gases_mod,  only:  virqd, vicvqd
+    use multi_gases_mod,  only:  virqd, vicpqd, vicvqd, virq, vicvq
 #endif
   use fv_regional_mod,     only: dump_field, exch_uv, H_STAGGER, U_STAGGER, V_STAGGER
   use fv_regional_mod,     only: a_step, p_step, k_step, n_step
@@ -199,7 +201,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
     real, intent(inout), dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed  ,npz):: v  !< D grid meridional wind (m/s)
     real, intent(inout) :: w(   bd%isd:,bd%jsd:,1:)  !< vertical vel. (m/s)
     real, intent(inout) ::  delz(bd%is:,bd%js:,1:)  !< delta-height (m, negative)
-    real, intent(inout) :: cappa(bd%isd:,bd%jsd:,1:) !< moist kappa
+    real, intent(inout) :: cappa(bd%isd:bd%ied,bd%jsd:bd%jed,1:npz)  !< moist kappa
 #ifdef MULTI_GASES
     real, intent(inout) :: kapad(bd%isd:bd%ied,bd%jsd:bd%jed,1:npz) !< multi_gases kappa
 #endif
@@ -1222,6 +1224,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
 ! Prevent accumulation of rounding errors at overlapped domain edges:
        call mpp_get_boundary(u, v, domain, ebuffery=ebuffer,  &
                              nbufferx=nbuffer, gridtype=DGRID_NE )
+                                                     call timing_off('COMM_TOTAL')
 !$OMP parallel do default(none) shared(is,ie,js,je,npz,u,nbuffer,v,ebuffer)
           do k=1,npz
              do i=is,ie
@@ -1234,6 +1237,16 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
 
     endif
 
+    if ( flagstruct%molecular_diffusion .and. md_time ) then
+! -----------------------------------------------------
+! direct explicit molecular diffusion
+! -----------------------------------------------------
+       call molecular_diffusion_run(u, v, w, delp, pt, pkz, cappa, q, bd,   &
+                 gridstruct, flagstruct, domain, i_pack, npx, npy, npz, nq, dt, it, akap, zvir, cv_air)
+    endif
+! -------------------------------------------------
+
+                                                     call timing_on('COMM_TOTAL')
 #ifndef ROT3
     if ( it/=n_split)   &
          call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE)
@@ -1243,9 +1256,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
 #ifdef SW_DYNAMICS
     endif
 #endif
-      if ( gridstruct%nested ) then
+    if ( gridstruct%nested ) then
          neststruct%nest_timestep = neststruct%nest_timestep + 1
-      endif
+    endif
 
 #ifdef SW_DYNAMICS
 #else
@@ -1298,13 +1311,12 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
     if (gridstruct%nested) then
 
 
-
 #ifndef SW_DYNAMICS
-         if (.not. hydrostatic) then
+        if (.not. hydrostatic) then
                call nested_grid_BC_apply_intT(w, &
                0, 0, npx, npy, npz, bd, split_timestep_BC+1, real(n_split*flagstruct%k_split), &
                neststruct%w_BC, bctype=neststruct%nestbctype  )
-       end if
+        end if
 #endif SW_DYNAMICS
             call nested_grid_BC_apply_intT(u, &
             0, 1, npx, npy, npz, bd, split_timestep_BC+1, real(n_split*flagstruct%k_split), &
@@ -1313,9 +1325,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
             1, 0, npx, npy, npz, bd, split_timestep_BC+1, real(n_split*flagstruct%k_split), &
             neststruct%v_BC, bctype=neststruct%nestbctype )
 
-      end if
+    end if
 
-      if (flagstruct%regional) then
+    if (flagstruct%regional) then
 
 #ifndef SW_DYNAMICS
          if (.not. hydrostatic) then
@@ -1340,7 +1352,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
                                        reg_bc_update_time )
 
          call mpp_update_domains(u, v, domain, gridtype=DGRID_NE)
-      end if
+
+    endif
 
       if ( do_diag_debug_dyn ) then
          call debug_column_dyn( pt, delp, delz, u, v, w, q, heat_source, cappa, akap, &
@@ -2784,5 +2797,181 @@ subroutine gz_bc(gz,delzBC,bd,npx,npy,npz,step,split)
 
  end subroutine gz_bc
 
+  subroutine molecular_diffusion_run(u,v,w,delp,pt,pkz,cappa,q,bd,   &
+                            gridstruct,flagstruct,domain,i_pack,npx,npy,npz, &
+                            nq,dt,it,akap,zvir,cv_air)
+  type(fv_grid_bounds_type), intent(IN) :: bd
+  real, intent(inout):: pkz(bd%is:bd%ie,bd%js:bd%je,npz)
+  real, intent(inout):: pt(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
+  real, intent(inout):: delp(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
+  real, intent(inout):: cappa(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
+  real, intent(inout):: w(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
+  real, intent(inout):: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz)
+  real, intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz)
+  real, intent(inout):: q( bd%isd:bd%ied,bd%jsd:bd%jed,npz, nq)
+
+  integer, intent(in):: npx,npy,npz,nq,it
+  real, intent(in)::akap,zvir,dt,cv_air
+  type(fv_grid_type),  intent(INOUT), target :: gridstruct
+  type(fv_flags_type), intent(IN),    target :: flagstruct
+  type(domain2d), intent(inout) :: domain
+  type(group_halo_update_type), intent(inout) :: i_pack(*)
+
+  real :: pkzf(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
+  real, dimension (bd%isd:bd%ied,bd%jsd:bd%jed) :: p, t, e
+  integer :: i, j, k
+  integer :: is,  ie,  js,  je
+  integer :: isd, ied, jsd, jed
+
+      is  = bd%is
+      ie  = bd%ie
+      js  = bd%js
+      je  = bd%je
+      isd = bd%isd
+      ied = bd%ied
+      jsd = bd%jsd
+      jed = bd%jed
+
+! -----------------------------------------------------
+! ------- update halo to prepare for diffusion -------
+    pkzf = 0.0
+    do k=1,npz
+       do j=js,je
+          pkzf(is:ie,j,k) = pkz(is:ie,j,k)
+       enddo
+    enddo
+
+                             call timing_on('COMM_TOTAL')
+    call start_group_halo_update(i_pack(1),delp,  domain, complete=.false.)
+    call start_group_halo_update(i_pack(1), pt,   domain, complete=.true.)
+    call start_group_halo_update(i_pack(2),pkzf,  domain)
+    call start_group_halo_update(i_pack(7), w, domain)
+    call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE)
+    if ( nq > 0 ) then
+                                       call timing_on('COMM_TRACER')
+                    call start_group_halo_update(i_pack(10), q, domain)
+                                       call timing_off('COMM_TRACER')
+    endif
+#ifdef MOIST_CAPPA
+    call start_group_halo_update(i_pack(12), cappa, domain)
+#endif
+
+    call complete_group_halo_update(i_pack(1), domain)  ! delp. pt
+    call complete_group_halo_update(i_pack(2), domain)  ! pkzf
+    call complete_group_halo_update(i_pack(7), domain)  ! w
+    call complete_group_halo_update(i_pack(8), domain)
+    if ( nq>0 ) then
+                                       call timing_on('COMM_TRACER')
+         call complete_group_halo_update(i_pack(10), domain)
+                                       call timing_off('COMM_TRACER')
+    endif
+#ifdef MOIST_CAPPA
+    call complete_group_halo_update(i_pack(12), domain)
+#endif
+                             call timing_off('COMM_TOTAL')
+
+    if( flagstruct%nord>0 .and. (.not. (flagstruct%regional))) then
+        i=mod(it-1,2)+1        ! alternatively to avoid bias
+        do k=1,npz
+        call copy_corners(pt(isd,jsd,k), npx, npy, i, gridstruct%nested, bd, &
+                          gridstruct%sw_corner, gridstruct%se_corner, &
+                          gridstruct%nw_corner, gridstruct%ne_corner)
+        call copy_corners(pkzf(isd,jsd,k), npx, npy, i, gridstruct%nested, bd, &
+                          gridstruct%sw_corner, gridstruct%se_corner, &
+                          gridstruct%nw_corner, gridstruct%ne_corner)
+        call copy_corners(cappa(isd,jsd,k), npx, npy, i, gridstruct%nested,bd, &
+                          gridstruct%sw_corner, gridstruct%se_corner, &
+                          gridstruct%nw_corner, gridstruct%ne_corner)
+        enddo
+    endif
+                                       call timing_on('d_md')
+
+!$OMP parallel do default(none) shared(npz,flagstruct,gridstruct,bd,      &
+!$OMP                                  it,dt,is,ie,js,je,isd,ied,jsd,jed, &
+!$OMP                                  pt,u,v,w,q,pkz,pkzf,cappa,akap,nq, &
+!$OMP                                  zvir,cv_air,md_layers,md_consv_te) &
+!$OMP                          private(k,i,j,p,t,e)
+! ----------------
+    do k=1, md_layers
+! ----------------
+
+! ------- prepare p and t for molecular diffusion coefficients
+
+       do j=jsd,jed
+          do i=isd,ied
+             t(i,j) = pt(i,j,k) * pkzf(i,j,k)
+#ifdef MULTI_GASES
+             t(i,j) = t(i,j) / virq(q(i,j,k,:))
+#else
+             t(i,j) = t(i,j) / (1+zvir*q(i,j,k,1))
+#endif
+#ifdef MOIST_CAPPA
+             p(i,j) = exp( log(pkzf(i,j,k)) / cappa(i,j,k) )
+#else
+#ifdef MULTI_GASES
+             p(i,j) = exp( log(pkzf(i,j,k)) / &
+                      (akap*virqd(q(i,j,k,:))/vicpqd(q(i,j,k,:))) )
+#else
+             p(i,j) = exp( log(pkzf(i,j,k)) / akap )
+#endif
+#endif
+             if( md_consv_te .gt. 0.0 ) &
+               e(i,j) = 0.5*(w(i,j,k)**2 + 0.5*gridstruct%rsin2(i,j)*(         &
+                        u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 -&
+                       (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*            &
+                        gridstruct%cosa_s(i,j)))
+          enddo
+       enddo
+
+! compute molecular diffusion with implicit time and dimensional splits
+
+       call d_md( p(isd,jsd), t(isd,jsd),         &
+                  u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), q, &
+                  it, nq, k, npz, dt,                            &
+                  gridstruct, flagstruct, bd)
+
+       do j=js,je
+          do i=is,ie
+             if( md_consv_te .gt. 0.0 ) then
+               e(i,j) = e(i,j) -  &
+                        0.5*(w(i,j,k)**2 + 0.5*gridstruct%rsin2(i,j)*(         &
+                        u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 -&
+                       (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*            &
+                        gridstruct%cosa_s(i,j)))
+#ifdef MOIST_CAPPA
+               t(i,j) = t(i,j) + e(i,j) / &
+#ifdef MULTI_GASES
+                     (rdgas* virq(q(i,j,k,:))  *(1./cappa(i,j,k)-1.))
+#else
+                     (rdgas*(1+zvir*q(i,j,k,1))*(1./cappa(i,j,k)-1.))
+#endif
+             endif
+             pkz(i,j,k) = exp( log(p(i,j)) * cappa(i,j,k) )
+#else
+#ifdef MULTI_GASES
+               t(i,j) = t(i,j) + e(i,j) / (cv_air*vicvqd(q(i,j,k,:)))
+             endif
+             pkz(i,j,k) = exp( log(p(i,j)) * &
+                          (akap*virqd(q(i,j,k,:))/vicpqd(q(i,j,k,:))) )
+#else
+               t(i,j) = t(i,j) + e(i,j)/cv_air
+             endif
+             pkz(i,j,k) = exp( log(p(i,j)) * akap )
+#endif
+#endif
+#ifdef MULTI_GASES
+             t(i,j) = t(i,j) * virq(q(i,j,k,:))
+#else
+             t(i,j) = t(i,j) * (1+zvir*q(i,j,k,1))
+#endif
+            pt(i,j,k) = t(i,j) / pkz(i,j,k)
+          enddo
+       enddo
+
+! -------------------------------------------------
+    enddo       ! k loop of 2d molecular diffusion
+! -------------------------------------------------
+ return
+ end subroutine molecular_diffusion_run
 
 end module dyn_core_mod