mpas_state_mod.F90

Go to the documentation of this file.

! (C) Copyright 2017 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 
module mpas_state_mod
 
use fckit_configuration_module, only: fckit_configuration
use fckit_mpi_module, only: fckit_mpi_comm, fckit_mpi_sum
use fckit_log_module, only: fckit_log
 
!oops
use datetime_mod
use kinds, only: kind_real
use oops_variables_mod, only: oops_variables
 
!dcmip initialization
use dcmip_initial_conditions_test_1_2_3, only : test1_advection_deformation, &
       test1_advection_hadley, test3_gravity_wave
use dcmip_initial_conditions_test_4, only : test4_baroclinic_wave
 
!ufo
use ufo_geovals_mod
use ufo_vars_mod
 
!MPAS-Model
use mpas_constants
use mpas_derived_types
use mpas_field_routines
use mpas_kind_types, only: strkind
use mpas_pool_routines
use mpas_dmpar, only: mpas_dmpar_exch_halo_field
 
 
!mpas-jedi
use mpas_constants_mod
use mpas_geom_mod
use mpas_fields_mod
use mpas2ufo_vars_mod
use mpas4da_mod
 
implicit none
 
private
 
public :: add_incr, &
        & analytic_ic, invent_state
 
! ------------------------------------------------------------------------------
 
contains
 
 
! ------------------------------------------------------------------------------
!> add increment to state
!!
!! \details **add_incr()** adds "increment" to "state", such as
!!          state (containing analysis) = state (containing guess) + increment
!!          Here, we also update "theta", "rho", and "u" (edge-normal wind), which are
!!          close to MPAS prognostic variable.
!!          Intermediate 3D pressure is diagnosed with hydrostatic balance.
!!          While conversion to "theta" and "rho" uses full state variables,
!!          conversion to "u" from cell center winds uses their increment to reduce
!!          the smoothing effect.
!!
subroutine add_incr(self, increment)
 
   implicit none
   class(mpas_fields), intent(inout) :: self !< state
   class(mpas_fields), intent(in)    :: increment
   character(len=StrKIND) :: kind_op
 
   integer :: ngrid
   type (mpas_pool_type), pointer :: state, diag, mesh
   type (field2dreal), pointer :: fld2d_pb, fld2d_u, fld2d_u_inc, fld2d_urm, fld2d_urz
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_qv, ptrr2_sh
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_p, ptrr2_rho, ptrr2_t, ptrr2_th, ptrr2_pp
   real(kind=kind_real), dimension(:), pointer :: ptrr1_ps
 
   ! Difference with self_add other is that self%subFields can contain extra fields
   ! beyond increment%subFields and the resolution of increment can be different.
 
   if (self%geom%nCells==increment%geom%nCells .and. self%geom%nVertLevels==increment%geom%nVertLevels) then
      ! First, update subFields that are common between self and increment
      kind_op = 'add'
      call da_operator(trim(kind_op), self%subFields, increment%subFields, fld_select = increment%fldnames_ci)
 
      ! Impose positive-definite limits on hydrometeors and moistureFields
      ! note: nonlinear COV (change of variables) from increment to state
      call da_posdef( self%subFields, mpas_hydrometeor_fields)
      call da_posdef( self%subFields, moisturefields)
 
      !NOTE: second, also update variables which are closely related to MPAS prognostic vars.
      ngrid = self%geom%nCellsSolve
 
      ! Update qv (water vapor mixing ratio) from spechum (specific humidity) [ w = q / (1 - q) ]
      ! note: nonlinear COV
      if ( all(self%has(moisturefields)) .and. &
           increment%has('spechum') .and. .not.increment%has('qv')) then
         call self%get(     'qv', ptrr2_qv)
         call self%get('spechum', ptrr2_sh)
         call q_to_w( ptrr2_sh(:,1:ngrid), ptrr2_qv(:,1:ngrid) )
      endif
 
      ! Enforce a hydrostatic balance to diagnose 3D pressure,
      ! plus update full state theta and rho
      ! note: nonlinear COV
      if ( all(self%has(analysisthermofields)) .and. &
           all(self%has(modelthermofields)) .and. &
           all(increment%has(analysisthermofields)) .and. &
           .not. all(increment%has(modelthermofields)) ) then
         call self%get(              'qv', ptrr2_qv)
         call self%get(        'pressure', ptrr2_p)
         call self%get(             'rho', ptrr2_rho)
         call self%get('surface_pressure', ptrr1_ps)
         call self%get(     'temperature', ptrr2_t)
         call self%get(           'theta', ptrr2_th)
 
         call hydrostatic_balance( ngrid, self%geom%nVertLevels, self%geom%zgrid(:,1:ngrid), &
                   ptrr2_t(:,1:ngrid), ptrr2_qv(:,1:ngrid), ptrr1_ps(1:ngrid), &
                   ptrr2_p(:,1:ngrid), ptrr2_rho(:,1:ngrid), ptrr2_th(:,1:ngrid) )
      endif
 
      ! Update pressure_p (pressure perturbation) , which is a diagnostic variable
      if ( self%has('pressure_p') .and. self%has('pressure') .and. &
           .not.increment%has('pressure_p') ) then
         call self%get(  'pressure', ptrr2_p)
         call self%get('pressure_p', ptrr2_pp)
         call mpas_pool_get_field(self%geom%domain%blocklist%allFields, 'pressure_base', fld2d_pb)
         ptrr2_pp(:,1:ngrid) = ptrr2_p(:,1:ngrid) - fld2d_pb%array(:,1:ngrid)
      endif
 
      ! Update edge normal wind u from uReconstructZonal and uReconstructMeridional "incrementally"
      ! note: linear COV
      if ( self%has('u') .and. &
           all(increment%has(cellcenteredwindfields)) .and. &
           .not.increment%has('u') ) then
         call mpas_pool_get_field(self%subFields, 'u', fld2d_u)
         call mpas_pool_get_field(increment%subFields, 'uReconstructMeridional', fld2d_urm)
         call mpas_pool_get_field(increment%subFields, 'uReconstructZonal', fld2d_urz)
 
         call mpas_duplicate_field(fld2d_u, fld2d_u_inc)
 
         call mpas_dmpar_exch_halo_field(fld2d_urz)
         call mpas_dmpar_exch_halo_field(fld2d_urm)
         call uv_cell_to_edges(self%geom%domain, fld2d_urz, fld2d_urm, fld2d_u_inc, &
                    self%geom%lonCell, self%geom%latCell, self%geom%nCells, &
                    self%geom%edgeNormalVectors, self%geom%nEdgesOnCell, self%geom%edgesOnCell, &
                    self%geom%nVertLevels)
         ngrid = self%geom%nEdgesSolve
         fld2d_u%array(:,1:ngrid) = fld2d_u%array(:,1:ngrid) + fld2d_u_inc%array(:,1:ngrid)
 
         ! TODO: DO we need HALO exchange here or in ModelMPAS::initialize for model integration?
         call mpas_deallocate_field( fld2d_u_inc )
      endif
   else
      call abor1_ftn("mpas_state:add_incr: dimension mismatch")
   endif
 
   return
 
end subroutine add_incr
 
! ------------------------------------------------------------------------------
!> Analytic Initialization for the MPAS Model
!!
!! \details **analytic_IC()** initializes the MPAS Field and State objects using one of
!! several alternative idealized analytic models.  This is intended to facilitate testing by
!! eliminating the need to read in the initial state from a file and by providing exact expressions
!! to test interpolations.  This function is activated by setting the "analytic_init" field in the
!! "initial" or "StateFile" section of the configuration file.
!!
!! Initialization options that begin with "dcmip" refer to tests defined by the multi-institutional
!! 2012 [Dynamical Core Intercomparison Project](https://earthsystealcmcog.org/projects/dcmip-2012)
!! and the associated Summer School, sponsored by NOAA, NSF, DOE, NCAR, and the University of Michigan.
!!
!! Currently implemented options for analytic_init include:
!! * invent-state: Backward compatibility with original analytic init option
!! * dcmip-test-1-1: 3D deformational flow
!! * dcmip-test-1-2: 3D Hadley-like meridional circulation
!! * dcmip-test-3-1: Non-hydrostatic gravity wave
!! * dcmip-test-4-0: Baroclinic instability
!!
!! \author J. Guerrette (adapted from fv3jedi code by M. Miesch)
!! \date July, 2018: Created
!!
subroutine analytic_ic(self, geom, f_conf, vdate)
 
!  !MPAS Test Cases
!  !JJG: This initialization requires the init_atmospher_core core_type
!  !      in the MPAS library for OOPS, but currently it is not included
!  use init_atm_core, only: init_atm_core_run!, init_atm_core_finalize (could be used for cleanup...)
 
  implicit none
 
  class(mpas_fields),        intent(inout) :: self   !< State
  type(mpas_geom), target,   intent(in)    :: geom   !< Geometry
  type(fckit_configuration), intent(in)    :: f_conf !< Configuration
  type(datetime),            intent(inout) :: vdate  !< DateTime
 
  character(len=:), allocatable :: str
  character(len=30) :: ic
  character(len=20) :: sdate
  character(len=1024) :: buf
  Integer :: jlev,ii
  integer :: ierr = 0
  real(kind=kind_real) :: rlat, rlon, z
  real(kind=kind_real) :: pk,pe1,pe2,ps
  real(kind=kind_real) :: u0,v0,w0,t0,phis0,ps0,rho0,hum0,q1,q2,q3,q4
 
  real(kind=kind_real)             :: dtdummy = 900.0
  logical, allocatable             :: grids_on_this_pe(:)
  integer                          :: p_split = 1
  real (kind=kind_real), dimension(:,:), pointer :: &
              u_ptr, v_ptr, temperature_ptr, p_ptr, &
              qv_ptr, qc_ptr, qr_ptr, qi_ptr, qs_ptr, &
              ln_p_ptr
  real(kind=kind_real), dimension(:),pointer :: ps_ptr
  integer, pointer :: index_qv, index_qc, index_qr, index_qi, index_qs
  type (field3dreal), pointer :: field3d
  type (mpas_pool_type), pointer :: pool_a, pool_b, state
  real(kind=kind_real) :: zhalf
 
  ! Pointer to geometry component of field object
  self%geom => geom
 
  If (f_conf%has("analytic_init")) Then
     call f_conf%get_or_die("analytic_init",str)
     ic = str
  Else
     ! This default value is for backward compatibility
     ic = "invent-state"
  EndIf
  call fckit_log%info ("mpas_state:analytic_init: "//ic)
 
! Conflicts with natural log below
!  call log%warning("mpas_state:analytic_init: "//IC)
  call f_conf%get_or_die("date",str)
  sdate = str
  call fckit_log%info ("validity date is: "//sdate)
 
  call datetime_set(sdate, vdate)
 
   ! Need to initialize variables that are used in interpolation/getVals
   ! In "create" and "read" subroutines, subFields are
   ! initialized from geom % domain % blocklist % allFields, zeroed,
   ! reread from file into allFields, then values copied to subFields
   ! -> must initialize allFields here and copy to subFields
 
   call mpas_pool_get_subpool(geom % domain % blocklist % structs, &
                              'state', state)
 
   pool_a => geom % domain % blocklist % allFields
 
   !Diagnostic vars (diag pool)
   call mpas_pool_get_array(pool_a, "pressure", p_ptr)
   call mpas_pool_get_array(pool_a, "uReconstructZonal", u_ptr)
   call mpas_pool_get_array(pool_a, "uReconstructMeridional", v_ptr)
   call mpas_pool_get_array(pool_a, "temperature", temperature_ptr)
   call mpas_pool_get_array(pool_a, "surface_pressure", ps_ptr)
 
 
   !Scalars (state pool)
   call mpas_pool_get_field(pool_a, "scalars", field3d)
   call mpas_pool_get_dimension(state, "index_qv", index_qv)
   if ( index_qv .gt. 0 ) &
      qv_ptr => field3d % array(index_qv,:,:)
 
   call mpas_pool_get_dimension(state, "index_qc", index_qc)
   if ( index_qc .gt. 0 ) &
      qc_ptr => field3d % array(index_qc,:,:)
 
   call mpas_pool_get_dimension(state, "index_qr", index_qr)
   if ( index_qr .gt. 0 ) &
      qr_ptr => field3d % array(index_qr,:,:)
 
   call mpas_pool_get_dimension(state, "index_qi", index_qi)
   if ( index_qi .gt. 0 ) &
      qi_ptr => field3d % array(index_qi,:,:)
 
   call mpas_pool_get_dimension(state, "index_qs", index_qs)
   if ( index_qs .gt. 0 ) &
      qs_ptr => field3d % array(index_qs,:,:)
 
  !===================================================================
  int_option: Select Case (ic)
 
     Case("invent-state")
 
        call invent_state(self,f_conf)
 
 
!     !TODO: This case requires the init_atmospher_core core_type to be
!     !      built as part of the MPAS library.
!     Case("mpas_init_case")
!
!!Would use init_atm_setup_case in MPAS-Release/src/core_init_atmosphere/mpas_init_atm_cases.F
!!mpas_init has already been called at this point from geo_setup
!
!!init_atms_setup_case is normally called from the following set of subroutines:
!!mpas_run => core_run [init_atm_core_run] => init_atm_setup_case => [select from preset cases]
!!Can we bypass the first two somehow?
!!Would use "config_init_case" in the yaml file, then check for matching with one of the ideal cases below... (not 7 or 8)
!
!!if ((config_init_case == 1) .or. (config_init_case == 2) .or. (config_init_case == 3)) then
!!   write(0,*) ' Jablonowski and Williamson baroclinic wave test case '
!!   if (config_init_case == 1) write(0,*) ' no initial perturbation '
!!   if (config_init_case == 2) write(0,*) ' initial perturbation included '
!!   if (config_init_case == 3) write(0,*) ' normal-mode perturbation included '
!!else if ((config_init_case == 4) .or. (config_init_case == 5)) then
!!   write(0,*) ' squall line - super cell test case '
!!   if (config_init_case == 4) write(0,*) ' squall line test case'
!!   if (config_init_case == 5) write(0,*) ' supercell test case'
!!      else if (config_init_case == 6 ) then
!!   write(0,*) ' mountain wave test case '
!!else if (config_init_case == 7 ) then
!!   write(0,*) ' real-data GFS test case '
!!else if (config_init_case == 8 ) then
!!   write(0,*) 'real-data surface (SST) update test case '
!
!       ierr = init_atm_core_run(geom % domain)
!       if ( ierr .ne. 0  ) then
!          call abor1_ftn("mpas_state: init_atm_core_run failed")
!       end if
 
     Case ("dcmip-test-1-1")
        do ii = 1, geom%nCellsSolve
           rlat = geom%latCell(ii)
           rlon = geom%lonCell(ii)
 
           ! Now loop over all levels
           do jlev = 1, geom%nVertLevels
 
              zhalf = mpas_jedi_half_kr * (geom%zgrid(jlev,ii) + geom%zgrid(jlev+1,ii))
              Call test1_advection_deformation(rlon,rlat,pk,zhalf,1,u0,v0,w0,t0,&
                                               phis0,ps0,rho0,hum0,q1,q2,q3,q4)
              p_ptr(jlev,ii) = pk
 
              u_ptr(jlev,ii) = u0 !MMiesch: ATTN Not going to necessary keep a-grid winds, u can be either a-
              v_ptr(jlev,ii) = v0 ! or staggered-grid so this needs to be generic. You cannot drive the model
                                  ! with A grid winds
              if (index_qv.gt.0) qv_ptr(jlev,ii) = hum0 !set to zero for this test
              if (index_qc.gt.0) qc_ptr(jlev,ii) = q1
              if (index_qi.gt.0) qi_ptr(jlev,ii) = q2
              if (index_qr.gt.0) qr_ptr(jlev,ii) = q3
              if (index_qs.gt.0) qs_ptr(jlev,ii) = q4
 
              temperature_ptr(jlev,ii) = t0
           enddo
           ps_ptr(ii) = ps0
        enddo
 
     Case ("dcmip-test-1-2")
 
        do ii = 1, geom%nCellsSolve
           rlat = geom%latCell(ii)
           rlon = geom%lonCell(ii)
 
           ! Now loop over all levels
           do jlev = 1, geom%nVertLevels
 
              zhalf = mpas_jedi_half_kr * (geom%zgrid(jlev,ii) + geom%zgrid(jlev+1,ii))
              Call test1_advection_hadley(rlon,rlat,pk,zhalf,1,u0,v0,w0,&
                                          t0,phis0,ps0,rho0,hum0,q1)
              p_ptr(jlev,ii) = pk
 
              u_ptr(jlev,ii) = u0 !MMiesch: ATTN Not going to necessary keep a-grid winds, u can be either a-
              v_ptr(jlev,ii) = v0 ! or staggered-grid so this needs to be generic. You cannot drive the model
                                  ! with A grid winds
              if (index_qv.gt.0) qv_ptr(jlev,ii) = hum0 !set to zero for this test
              if (index_qc.gt.0) qc_ptr(jlev,ii) = q1
 
              if (index_qi.gt.0) qi_ptr(jlev,ii) = 0._kind_real
              if (index_qr.gt.0) qr_ptr(jlev,ii) = 0._kind_real
              if (index_qs.gt.0) qs_ptr(jlev,ii) = 0._kind_real
 
              temperature_ptr(jlev,ii) = t0
           enddo
           ps_ptr(ii) = ps0
        enddo
 
     Case ("dcmip-test-3-1")
 
        do ii = 1, geom%nCellsSolve
           rlat = geom%latCell(ii)
           rlon = geom%lonCell(ii)
 
           ! Now loop over all levels
           do jlev = 1, geom%nVertLevels
 
              zhalf = mpas_jedi_half_kr * (geom%zgrid(jlev,ii) + geom%zgrid(jlev+1,ii))
              Call test3_gravity_wave(rlon,rlat,pk,zhalf,1,u0,v0,w0,&
                                      t0,phis0,ps0,rho0,hum0)
 
              p_ptr(jlev,ii) = pk
              u_ptr(jlev,ii) = u0 !MMiesch: ATTN Not going to necessary keep a-grid winds, u can be either a-
              v_ptr(jlev,ii) = v0 ! or staggered-grid so this needs to be generic. You cannot drive the model
                                  ! with A grid winds
 
              if (index_qv.gt.0) qv_ptr(jlev,ii) = hum0 !set to zero for this test
              if (index_qc.gt.0) qc_ptr(jlev,ii) = 0._kind_real
              if (index_qi.gt.0) qi_ptr(jlev,ii) = 0._kind_real
              if (index_qr.gt.0) qr_ptr(jlev,ii) = 0._kind_real
              if (index_qs.gt.0) qs_ptr(jlev,ii) = 0._kind_real
 
              temperature_ptr(jlev,ii) = t0
           enddo
           ps_ptr(ii) = ps0
        enddo
 
     Case ("dcmip-test-4-0")
 
        do ii = 1, geom%nCellsSolve
           rlat = geom%latCell(ii)
           rlon = geom%lonCell(ii)
 
           ! Now loop over all levels
           do jlev = 1, geom%nVertLevels
 
              zhalf = mpas_jedi_half_kr * (geom%zgrid(jlev,ii) + geom%zgrid(jlev+1,ii))
              Call test4_baroclinic_wave(0,mpas_jedi_one_kr,rlon,rlat,pk,zhalf,1,u0,v0,w0,&
                                      t0,phis0,ps0,rho0,hum0,q1,q2)
 
              p_ptr(jlev,ii) = pk
 
              u_ptr(jlev,ii) = u0 !MMiesch: ATTN Not going to necessary keep a-grid winds, u can be either a-
 
              v_ptr(jlev,ii) = v0 ! or staggered-grid so this needs to be generic. You cannot drive the model
                                  ! with A grid winds
 
              if (index_qv.gt.0) qv_ptr(jlev,ii) = hum0 !set to zero for this test
              if (index_qc.gt.0) qc_ptr(jlev,ii) = 0._kind_real
              if (index_qi.gt.0) qi_ptr(jlev,ii) = 0._kind_real
              if (index_qr.gt.0) qr_ptr(jlev,ii) = 0._kind_real
              if (index_qs.gt.0) qs_ptr(jlev,ii) = 0._kind_real
 
              temperature_ptr(jlev,ii) = t0
           enddo
           ps_ptr(ii) = ps0
        enddo
 
     Case Default
 
        call invent_state(self,f_conf)
 
     End Select int_option
 
     call da_copy_all2sub_fields(self % geom % domain, self % subFields)
 
     call fckit_log%debug ('==> end mpas_state:analytic_init')
 
end subroutine analytic_ic
 
! ------------------------------------------------------------------------------
 
subroutine invent_state(self,f_conf)
 
   implicit none
 
   class(mpas_fields),        intent(inout) :: self    !< Model fields
   type(fckit_configuration), intent(in)    :: f_conf  !< Configuration structure
   real (kind=kind_real), dimension(:,:), pointer :: r2d_ptr_a, r2d_ptr_b
   real (kind=kind_real), dimension(:), pointer :: r1d_ptr_a, r1d_ptr_b
   integer :: jlev,ii
   type (mpas_pool_type), pointer :: pool_a, state
   type (field3dreal), pointer :: field3d
   integer, pointer :: index_qv
 
   !- read/interp.
 
   call mpas_pool_get_subpool(self % geom % domain % blocklist % structs, &
                              'state', state)
   pool_a => self % geom % domain % blocklist % allFields
 
   !Diagnostic vars (diag pool)
   !u
   call mpas_pool_get_array(pool_a, "uReconstructZonal", r2d_ptr_a)
   do jlev = 1,self % geom % nVertLevels
      do ii = 1, self % geom % nCellsSolve
         r2d_ptr_a(jlev,ii) = cos(0.25*self % geom % lonEdge(ii)) + cos(0.25*self % geom % latEdge(ii))
      enddo
   enddo
 
   !v
   call mpas_pool_get_array(pool_a, "uReconstructMeridional", r2d_ptr_a)
   do jlev = 1,self % geom % nVertLevels
      do ii = 1, self % geom % nCellsSolve
         r2d_ptr_a(jlev,ii) = mpas_jedi_one_kr
      enddo
   enddo
 
   !temperature
   call mpas_pool_get_array(pool_a, "temperature", r2d_ptr_a)
   do jlev = 1,self % geom % nVertLevels
      do ii = 1, self % geom % nCellsSolve
         r2d_ptr_a(jlev,ii) = cos(0.25*self % geom % lonCell(ii)) + cos(0.25*self % geom % latCell(ii))
      enddo
   enddo
 
   !pressure
   call mpas_pool_get_array(pool_a, "pressure", r2d_ptr_a)
   do jlev = 1,self % geom % nVertLevels
      do ii = 1, self % geom % nCellsSolve
         r2d_ptr_a(jlev,ii) = real(jlev,kind_real)
      enddo
   enddo
 
   !surface_pressure
   call mpas_pool_get_array(pool_a, "surface_pressure", r1d_ptr_a)
   do ii = 1, self % geom % nCellsSolve
      r1d_ptr_a(ii) = mpas_jedi_one_kr
   enddo
   !Scalars (state pool)
   !qv
   call mpas_pool_get_field(pool_a, 'scalars', field3d)
   call mpas_pool_get_dimension(state, 'index_qv', index_qv)
   if ( index_qv .gt. 0 ) then
      r2d_ptr_a => field3d % array(index_qv,:,:)
      do jlev = 1,self % geom % nVertLevels
         do ii = 1, self % geom % nCellsSolve
            r2d_ptr_a(jlev,ii) = mpas_jedi_zero_kr
         enddo
      enddo
   end if
 
   return
 
end subroutine invent_state
 
! ------------------------------------------------------------------------------
 
end module mpas_state_mod