Fixes to run 32BIT MULTI_GASES

docs/FV3_citations.bib

@@ -33,6 +33,16 @@ @article{chen2013seasonal
   doi={10.1175/JCLI-D-12-00061.1}
 }
 
+@article{zhou2019toward,
+  title={Toward Convective-Scale Prediction within the Next Generation Global Prediction System},
+  author={Zhou, Linjiong and Lin, Shian-Jiann and Chen, Jan-Huey and Harris, Lucas M. and Chen, Xi and Rees, Shannon L.},
+  journal={Bulletin of the American Meteorological Society},
+  volume={100},
+  issue={7},
+  year={2019},
+  doi={10.1175/bams-d-17-0246.1}
+}
+
 @article {deng2008cirrus,
 author = {Deng, Min and Mace, Gerald G.},
 title = {Cirrus cloud microphysical properties and air motion statistics using cloud radar Doppler moments: Water content, particle size, and sedimentation relationships},

driver/fvGFS/atmosphere.F90

@@ -112,7 +112,7 @@ module atmosphere_mod
 !   <tr>
 !     <td>mpp_mod</td>
 !     <td>mpp_error, stdout, FATAL, NOTE, input_nml_file, mpp_root_pe,
-!                    mpp_npes, mpp_pe, mpp_chksum,mpp_get_current_pelist,     
+!                    mpp_npes, mpp_pe, mpp_chksum,mpp_get_current_pelist,
 !                    mpp_set_current_pelist</td>
 !   </tr>
 !   <tr>
@@ -183,7 +183,7 @@ module atmosphere_mod
 use fv_nggps_diags_mod, only: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
 use fv_restart_mod,     only: fv_restart, fv_write_restart
 use fv_timing_mod,      only: timing_on, timing_off
-use fv_mp_mod,          only: switch_current_Atm
+use fv_mp_mod,          only: switch_current_Atm, is_master
 use fv_sg_mod,          only: fv_subgrid_z
 use fv_update_phys_mod, only: fv_update_phys
 use fv_nwp_nudge_mod,   only: fv_nwp_nudge_init, fv_nwp_nudge_end, do_adiabatic_init
@@ -194,6 +194,7 @@ module atmosphere_mod
                               a_step, p_step, current_time_in_seconds
 
 use mpp_domains_mod,    only:  mpp_get_data_domain, mpp_get_compute_domain
+use fv_grid_utils_mod,  only: g_sum
 
 implicit none
 private
@@ -269,30 +270,13 @@ module atmosphere_mod
 !! and diagnostics.  
  subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 #ifdef CCPP
-#ifdef STATIC
-! For static builds, the ccpp_physics_{init,run,finalize} calls
-! are not pointing to code in the CCPP framework, but to auto-generated
-! ccpp_suite_cap and ccpp_group_*_cap modules behind a ccpp_static_api
-   use ccpp_api,          only: ccpp_init
    use ccpp_static_api,   only: ccpp_physics_init
-#else
-   use iso_c_binding,     only: c_loc
-   use ccpp_api,          only: ccpp_init,           &
-                                ccpp_physics_init,   &
-                                ccpp_field_add,      &
-                                ccpp_error
-#endif
    use CCPP_data,         only: ccpp_suite,          &
                                 cdata => cdata_tile, &
                                 CCPP_interstitial
 #ifdef OPENMP
    use omp_lib
 #endif
-#ifndef STATIC
-! Begin include auto-generated list of modules for ccpp
-#include "ccpp_modules_fast_physics.inc"
-! End include auto-generated list of modules for ccpp
-#endif
 #endif
    type (time_type),    intent(in)    :: Time_init, Time, Time_step
    type(grid_box_type), intent(inout) :: Grid_box
@@ -442,15 +426,8 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 #ifdef CCPP
    ! Do CCPP fast physics initialization before call to adiabatic_init (since this calls fv_dynamics)
 
-   ! Initialize the cdata structure
-   call ccpp_init(trim(ccpp_suite), cdata, ierr)
-   if (ierr/=0) then
-      cdata%errmsg = ' atmosphere_dynamics: error in ccpp_init: ' // trim(cdata%errmsg)
-      call mpp_error (FATAL, cdata%errmsg)
-   end if
-
-   ! For fast physics running over the entire domain, block and thread
-   ! number are not used; set to safe values
+   ! For fast physics running over the entire domain, block
+   ! and thread number are not used; set to safe values
    cdata%blk_no = 1
    cdata%thrd_no = 1
 
@@ -486,18 +463,9 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 #endif
                                  mpirank=mpp_pe(), mpiroot=mpp_root_pe())
 
-#ifndef STATIC
-! Populate cdata structure with fields required to run fast physics (auto-generated).
-#include "ccpp_fields_fast_physics.inc"
-#endif
-
    if (Atm(mytile)%flagstruct%do_sat_adj) then
       ! Initialize fast physics
-#ifdef STATIC
       call ccpp_physics_init(cdata, suite_name=trim(ccpp_suite), group_name="fast_physics", ierr=ierr)
-#else
-      call ccpp_physics_init(cdata, group_name="fast_physics", ierr=ierr)
-#endif
       if (ierr/=0) then
          cdata%errmsg = ' atmosphere_dynamics: error in ccpp_physics_init for group fast_physics: ' // trim(cdata%errmsg)
          call mpp_error (FATAL, cdata%errmsg)
@@ -508,7 +476,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area)
 !  --- initiate the start for a restarted regional forecast
    if ( Atm(mytile)%gridstruct%regional .and. Atm(mytile)%flagstruct%warm_start ) then
 
-     call start_regional_restart(Atm(1),             &
+     call start_regional_restart(Atm(1), dt_atmos,   &
                                  isc, iec, jsc, jec, &
                                  isd, ied, jsd, jed )
    endif
@@ -712,32 +680,22 @@ end subroutine atmosphere_dynamics
 
 !>@brief The subroutine 'atmosphere_end' is an API for the termination of the
 !! FV3 dynamical core responsible for writing out a restart and final diagnostic state.
- subroutine atmosphere_end (Time, Grid_box)
+ subroutine atmosphere_end (Time, Grid_box, restart_endfcst)
 #ifdef CCPP
-#ifdef STATIC
-! For static builds, the ccpp_physics_{init,run,finalize} calls
-! are not pointing to code in the CCPP framework, but to auto-generated
-! ccpp_suite_cap and ccpp_group_*_cap modules behind a ccpp_static_api
    use ccpp_static_api,   only: ccpp_physics_finalize
    use CCPP_data,         only: ccpp_suite
-#else
-   use ccpp_api,          only: ccpp_physics_finalize
-#endif
    use CCPP_data,         only: cdata => cdata_tile
 #endif
    type (time_type),      intent(in)    :: Time
    type(grid_box_type),   intent(inout) :: Grid_box
+   logical,               intent(in)    :: restart_endfcst
 
 #ifdef CCPP
    integer :: ierr
 
    if (Atm(mytile)%flagstruct%do_sat_adj) then
       ! Finalize fast physics
-#ifdef STATIC
       call ccpp_physics_finalize(cdata, suite_name=trim(ccpp_suite), group_name="fast_physics", ierr=ierr)
-#else
-      call ccpp_physics_finalize(cdata, group_name="fast_physics", ierr=ierr)
-#endif
       if (ierr/=0) then
          cdata%errmsg = ' atmosphere_dynamics: error in ccpp_physics_finalize for group fast_physics: ' // trim(cdata%errmsg)
          call mpp_error (FATAL, cdata%errmsg)
@@ -754,7 +712,7 @@ subroutine atmosphere_end (Time, Grid_box)
       call timing_off('FV_DIAG')
    endif
 
-   call fv_end(Atm, grids_on_this_pe)
+   call fv_end(Atm, grids_on_this_pe, restart_endfcst)
    deallocate (Atm)
 
    deallocate( u_dt, v_dt, t_dt, pref, dum1d )
@@ -1267,7 +1225,6 @@ subroutine atmosphere_get_bottom_layer (Atm_block, DYCORE_Data)
    rrg  = rdgas / grav
 
    if (first_time) then
-     print *, 'calculating slp kr value'
      ! determine 0.8 sigma reference level
      sigtop = Atm(mytile)%ak(1)/pstd_mks+Atm(mytile)%bk(1)
      do k = 1, npz
@@ -1394,6 +1351,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    integer :: i, j, ix, k, k1, n, w_diff, nt_dyn, iq
    integer :: nb, blen, nwat, dnats, nq_adv
    real(kind=kind_phys):: rcp, q0, qwat(nq), qt, rdt
+   real psum, qsum, psumb, qsumb, betad
    Time_prev = Time
    Time_next = Time + Time_step_atmos
    rdt = 1.d0 / dt_atmos
@@ -1406,6 +1364,18 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
    if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers')
 
    if (IAU_Data%in_interval) then
+      if (IAU_Data%drymassfixer) then
+         ! global mean total pressure and water before IAU
+         psumb = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
+         qsumb = g_sum(Atm(n)%domain,&
+                 sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+                 isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
+         if (is_master()) then
+           print *,'dry ps before IAU/physics',psumb+Atm(n)%ptop-qsumb
+         endif
+      endif
+
 !     IAU increments are in units of 1/sec
 
 !     add analysis increment to u,v,t tendencies
@@ -1468,11 +1438,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
      call fill_gfs(blen, npz, IPD_Data(nb)%Statein%prsi, IPD_Data(nb)%Stateout%gq0, 1.e-9_kind_phys)
 
      do k = 1, npz
-           if(flip_vc) then
-             k1 = npz+1-k !reverse the k direction 
-           else
-             k1 = k
-           endif
+       if(flip_vc) then
+         k1 = npz+1-k !reverse the k direction 
+       else
+         k1 = k
+       endif
        do ix = 1, blen
          i = Atm_block%index(nb)%ii(ix)
          j = Atm_block%index(nb)%jj(ix)
@@ -1514,11 +1484,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
      !--- See Note in statein...
      do iq = nq+1, ncnst
        do k = 1, npz
-           if(flip_vc) then
-             k1 = npz+1-k !reverse the k direction 
-           else
-             k1 = k
-           endif
+         if(flip_vc) then
+           k1 = npz+1-k !reverse the k direction 
+         else
+           k1 = k
+         endif
          do ix = 1, blen
            i = Atm_block%index(nb)%ii(ix)
            j = Atm_block%index(nb)%jj(ix)
@@ -1529,6 +1499,29 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 
    enddo  ! nb-loop
 
+! dry mass fixer in IAU interval following
+! https://onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0870.2007.00299.x
+   if (IAU_Data%in_interval .and. IAU_data%drymassfixer) then
+      ! global mean total pressure
+      psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
+      ! global mean total water (before adjustment)
+      qsum = g_sum(Atm(n)%domain,&
+             sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+             isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.) 
+      betad = (psum - (psumb - qsumb))/qsum
+      if (is_master()) then
+        print *,'dry ps after IAU/physics',psum+Atm(n)%ptop-qsum
+      endif
+      Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat)
+      !qsum = g_sum(Atm(n)%domain,&
+      !       sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
+      !       isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 
+      !if (is_master()) then
+      !  print *,'dry ps after iau_drymassfixer',psum+Atm(n)%ptop-qsum
+      !endif
+   endif
+
    call timing_off('GFS_TENDENCIES')
 
    w_diff = get_tracer_index (MODEL_ATMOS, 'w_diff' )
@@ -1592,6 +1585,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
 
      call nullify_domain()
      call timing_on('FV_DIAG')
+
      call fv_diag(Atm(mytile:mytile), zvir, fv_time, Atm(mytile)%flagstruct%print_freq)
      first_diag = .false.
      call timing_off('FV_DIAG')
@@ -1905,6 +1899,10 @@ subroutine atmos_phys_driver_statein (IPD_Data, Atm_block,flip_vc)
    real(kind=kind_phys) :: pk0inv, ptop, pktop
    real(kind=kind_phys) :: rTv, dm, qgrs_rad
    integer :: nb, blen, npz, i, j, k, ix, k1, kz, dnats, nq_adv
+#ifdef MULTI_GASES
+   real :: q_grs(nq), q_min
+#endif
+
 
 !!! NOTES: lmh 6nov15
 !!! - "Layer" means "layer mean", ie. the average value in a layer
@@ -1925,7 +1923,13 @@ subroutine atmos_phys_driver_statein (IPD_Data, Atm_block,flip_vc)
 !$OMP             shared  (Atm_block, Atm, IPD_Data, npz, nq, ncnst, sphum, liq_wat, &
 !$OMP                      ice_wat, rainwat, snowwat, graupel, pk0inv, ptop,   &
 !$OMP                      pktop, zvir, mytile, dnats, nq_adv, flip_vc) &
+#ifdef MULTI_GASES
+
+!$OMP             private (dm, nb, blen, i, j, ix, k1, kz, rTv, qgrs_rad, q_min, q_grs)
+
+#else
 !$OMP             private (dm, nb, blen, i, j, ix, k1, kz, rTv, qgrs_rad)
+#endif
 
    do nb = 1,Atm_block%nblks
 ! gas_phase_mass <-- prsl
@@ -2032,7 +2036,9 @@ subroutine atmos_phys_driver_statein (IPD_Data, Atm_block,flip_vc)
         do i=1,blen
 ! Geo-potential at interfaces:
 #ifdef MULTI_GASES
-           rTv = rdgas*IPD_Data(nb)%Statein%tgrs(i,k)*virq_max(IPD_Data(nb)%Statein%qgrs(i,k,:),qmin)
+           q_grs(1:nq_adv) = IPD_Data(nb)%Statein%qgrs(i,k,1:nq_adv)
+           q_min = qmin
+           rTv = rdgas*IPD_Data(nb)%Statein%tgrs(i,k)*virq_max(q_grs(:),q_min)
 #else
            qgrs_rad = max(qmin,IPD_Data(nb)%Statein%qgrs(i,k,sphum))
            rTv = rdgas*IPD_Data(nb)%Statein%tgrs(i,k)*(1.+zvir*qgrs_rad)

driver/fvGFS/fv_nggps_diag.F90

@@ -93,13 +93,13 @@ module fv_nggps_diags_mod
  real, parameter:: stndrd_atmos_lapse = 0.0065
 
  logical master
- integer :: id_ua, id_va, id_pt, id_delp, id_pfhy, id_pfnh 
- integer :: id_w, id_delz, id_diss, id_ps, id_hs, id_dbz
+ integer :: id_ua, id_va, id_pt, id_delp, id_pfhy, id_pfnh
+ integer :: id_w, id_delz, id_diss, id_ps, id_hs, id_dbz, id_omga
  integer :: kstt_ua, kstt_va, kstt_pt, kstt_delp, kstt_pfhy
  integer :: kstt_pfnh, kstt_w, kstt_delz, kstt_diss, kstt_ps,kstt_hs
  integer :: kend_ua, kend_va, kend_pt, kend_delp, kend_pfhy
  integer :: kend_pfnh, kend_w, kend_delz, kend_diss, kend_ps,kend_hs
- integer :: kstt_dbz, kend_dbz
+ integer :: kstt_dbz, kend_dbz, kstt_omga, kend_omga
  integer :: kstt_windvect, kend_windvect
  integer :: id_wmaxup,id_wmaxdn,kstt_wup, kend_wup,kstt_wdn,kend_wdn
  integer :: id_uhmax03,id_uhmin03,id_uhmax25,id_uhmin25,id_maxvort01
@@ -236,6 +236,13 @@ subroutine fv_nggps_diag_init(Atm, axes, Time)
           endif
        endif
 
+       id_omga = register_diag_field ( trim(file_name), 'omga', axes(1:3), Time,        &
+            'Vertical pressure velocity', 'pa/sec', missing_value=missing_value )
+       if (id_omga>0) then
+          kstt_omga = nlevs+1; kend_omga = nlevs+npzo
+          nlevs = nlevs + npzo
+       endif
+
        id_pt   = register_diag_field ( trim(file_name), 'temp', axes(1:3), Time,       &
             'temperature', 'K', missing_value=missing_value, range=trange )
        if (id_pt>0) then
@@ -482,6 +489,11 @@ subroutine fv_nggps_diag(Atm, zvir, Time)
        call store_data(id_w, Atm(n)%w(isco:ieco,jsco:jeco,:), Time, kstt_w, kend_w)
     endif
 
+    !--- OMGA (non-hydrostatic)
+    if ( id_omga>0  ) then
+       call store_data(id_omga, Atm(n)%omga(isco:ieco,jsco:jeco,:), Time, kstt_omga, kend_omga)
+    endif
+
     !--- TEMPERATURE
     if(id_pt   > 0) call store_data(id_pt, Atm(n)%pt(isco:ieco,jsco:jeco,:), Time, kstt_pt, kend_pt)
 
@@ -502,7 +514,7 @@ subroutine fv_nggps_diag(Atm, zvir, Time)
        do k=1,npzo
          do j=jsco,jeco
            do i=isco,ieco
-             wk(i,j,k) = -Atm(n)%delz(i,j,k)
+             wk(i,j,k) = Atm(n)%delz(i,j,k)
            enddo
          enddo
        enddo
@@ -1172,6 +1184,13 @@ subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
        if(rc==0)  num_field_dyn=num_field_dyn+1
      endif
    endif
+!
+   if( id_omga>0  ) then
+     call find_outputname(trim(file_name),'omga',output_name)
+     call add_field_to_bundle(trim(output_name),'Vertical pressure velocity', 'pa/sec', "time: point",   &
+          axes(1:3), fcst_grid, kstt_omga,kend_omga, dyn_bundle, output_file, rcd=rc)
+     if(rc==0)  num_field_dyn=num_field_dyn+1
+   endif
 !
    if(id_pt > 0) then
      call find_outputname(trim(file_name),'temp',output_name)