Add TCWA2 Microphysics for NCAR MPAS

src/core_atmosphere/dynamics/mpas_atm_time_integration.F

@@ -1920,6 +1920,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
       type (mpas_pool_type), pointer :: tend
       type (mpas_pool_type), pointer :: tend_physics => null()
       type (mpas_pool_type), pointer :: lbc  ! regional_MPAS addition
+      type (mpas_pool_type), pointer :: sfc_input ! for tcwa2
 
       real (kind=RKIND), dimension(:,:), pointer :: w
       real (kind=RKIND), dimension(:,:), pointer :: u, uReconstructZonal, uReconstructMeridional, uReconstructX, uReconstructY, uReconstructZ
@@ -1966,6 +1967,8 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
       call mpas_pool_get_subpool(block % structs, 'diag', diag)
       call mpas_pool_get_subpool(block % structs, 'tend', tend)
       call mpas_pool_get_subpool(block % structs, 'tend_physics', tend_physics)
+      call mpas_pool_get_subpool(block % structs, 'sfc_input', sfc_input)
+
 #ifdef DO_PHYSICS
       call mpas_pool_get_subpool(block % structs, 'diag_physics', diag_physics)
 #endif
@@ -2808,7 +2811,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
         call mpas_timer_start('microphysics')
 !$OMP PARALLEL DO
         do thread=1,nThreads
-           call driver_microphysics ( block % configs, mesh, state, 2, diag, diag_physics, tend_physics, tend, itimestep, &
+           call driver_microphysics ( block % configs, mesh, state, 2, diag, diag_physics, tend_physics, tend, sfc_input, itimestep, &
                                       cellSolveThreadStart(thread), cellSolveThreadEnd(thread))
         end do
 !$OMP END PARALLEL DO

src/core_atmosphere/physics/mpas_atmphys_control.F

@@ -177,6 +177,7 @@ subroutine physics_namelist_check(configs)
            config_microp_scheme .eq. 'mp_kessler'           .or. &
            config_microp_scheme .eq. 'mp_thompson'          .or. &
            config_microp_scheme .eq. 'mp_thompson_aerosols' .or. &
+           config_microp_scheme .eq. 'mp_tcwa2'             .or. &
            config_microp_scheme .eq. 'mp_wsm6')) then
 
     write(mpas_err_message,'(A,A20)') 'illegal value for config_microp_scheme:', &

src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F

@@ -19,6 +19,7 @@ module mpas_atmphys_driver_microphysics
 !wrf physics:
  use module_mp_kessler
  use module_mp_thompson
+ use module_mp_tcwa2
  use module_mp_wsm6,only: wsm6
  use mp_wsm6,only: mp_wsm6_init,refl10cm_wsm6
 
@@ -127,7 +128,7 @@ subroutine allocate_microphysics(configs)
  if(.not.allocated(rainncv_p)) allocate(rainncv_p(ims:ime,jms:jme))
 
  microp_select: select case(trim(microp_scheme))
-    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6")
+    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6","mp_tcwa2")
        !mass mixing ratios:
        if(.not.allocated(qi_p)) allocate(qi_p(ims:ime,kms:kme,jms:jme))
        if(.not.allocated(qs_p)) allocate(qs_p(ims:ime,kms:kme,jms:jme))
@@ -167,6 +168,20 @@ subroutine allocate_microphysics(configs)
              case default
           end select microp3_select
 
+       case("mp_tcwa2")
+          if(.not.allocated(nc_p) ) allocate(nc_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(nr_p) ) allocate(nr_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(ni_p) ) allocate(ni_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(ns_p) ) allocate(ns_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(ng_p) ) allocate(ng_p(ims:ime,kms:kme,jms:jme))
+          !height above sea level (terrain height) and land-ocean mask
+          if(.not.allocated(xland_p)) allocate(xland_p(ims:ime,jms:jme))
+          if(.not.allocated(refl10cm_p)) allocate(refl10cm_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(zdr_p)) allocate(zdr_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(kdp_p)) allocate(kdp_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(smlf_p)) allocate(smlf_p(ims:ime,kms:kme,jms:jme))
+          if(.not.allocated(gmlf_p)) allocate(gmlf_p(ims:ime,kms:kme,jms:jme))
+
        case default
     end select microp2_select
 
@@ -208,7 +223,7 @@ subroutine deallocate_microphysics(configs)
  if(allocated(rainncv_p)) deallocate(rainncv_p)
 
  microp_select: select case(trim(microp_scheme))
-    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6")
+    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6","mp_tcwa2")
        !mass mixing ratios:
        if(allocated(qi_p)) deallocate(qi_p)
        if(allocated(qs_p)) deallocate(qs_p)
@@ -248,6 +263,20 @@ subroutine deallocate_microphysics(configs)
              case default
           end select microp3_select
 
+       case("mp_tcwa2")
+          if(allocated(nc_p) ) deallocate(nc_p )
+          if(allocated(nr_p) ) deallocate(nr_p )
+          if(allocated(ni_p) ) deallocate(ni_p )
+          if(allocated(ns_p) ) deallocate(ns_p )
+          if(allocated(ng_p) ) deallocate(ng_p )
+          !height above sea level (terrain height) and land-ocean mask
+          if(allocated(xland_p)) deallocate(xland_p)
+          if(allocated(refl10cm_p)) deallocate(refl10cm_p)
+          if(allocated(zdr_p)) deallocate(zdr_p)
+          if(allocated(kdp_p)) deallocate(kdp_p)
+          if(allocated(smlf_p)) deallocate(smlf_p)
+          if(allocated(gmlf_p)) deallocate(gmlf_p)
+
        case default
     end select microp2_select
 
@@ -314,14 +343,15 @@ subroutine init_microphysics(dminfo,configs,mesh,state,time_lev,sfc_input,diag_p
  end subroutine init_microphysics
 
 !=================================================================================================================
- subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,tend_physics,tend,itimestep,its,ite)
+ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,tend_physics,tend,sfc_input,itimestep,its,ite)
 !=================================================================================================================
 
  use mpas_constants, only : rvord
 
 !input arguments:
  type(mpas_pool_type),intent(in):: configs
  type(mpas_pool_type),intent(in):: mesh
+ type(mpas_pool_type),intent(in):: sfc_input
 
  integer,intent(in):: time_lev
  integer,intent(in):: itimestep
@@ -364,7 +394,7 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
  call precip_from_MPAS(configs,diag_physics,its,ite)
 
 !... initialization of soundings for non-hydrostatic dynamical cores.
- call microphysics_from_MPAS(configs,mesh,state,time_lev,diag,diag_physics,tend_physics,its,ite)
+ call microphysics_from_MPAS(configs,mesh,state,time_lev,diag,diag_physics,tend_physics,sfc_input,its,ite)
 
 !... call to different cloud microphysics schemes:
  microp_select: select case(trim(microp_scheme))
@@ -460,6 +490,21 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
                )
        call mpas_timer_stop('mp_wsm6')
 
+    case ("mp_tcwa2")
+       call mpas_timer_start('mp_tcwa2')
+       call mp_tcwa2(TH = th_p, P = pres_p, DZ = dz_p, W = w_p, PII = pi_p, DT_MP = dt_microp,  &
+                     XLAND = xland_p, QV = qv_p, QC = qc_p, QR = qr_p, QI = qi_p, QS = qs_p,    &
+                     QG = qg_p, NC = nc_p, NR = nr_p, NI = ni_p, NS = ns_p, NG = ng_p,          &
+                     DBZM = refl10cm_p, ZDR = zdr_p, KDP = kdp_p, SMLF = smlf_p, GMLF = gmlf_p, &
+                     re_cloud = recloud_p, re_ice = reice_p, re_snow = resnow_p,                &
+                     has_reqc = has_reqc, has_reqi = has_reqi, has_reqs = has_reqs, SR = sr_p,  &
+                     RAINNC = rainnc_p, RAINNCV = rainncv_p, SNOWNC = snownc_p,                 &
+                     SNOWNCV = snowncv_p, GRAPNC = graupelnc_p, GRAPNCV = graupelncv_p,         &
+                     IDS = ids, IDE = ide, JDS = jds, JDE = jde, KDS = kds, KDE = kde,          &
+                     IMS = ims, IME = ime, JMS = jms, JME = jme, KMS = kms, KME = kme,          &
+                     ITS = its, ITE = ite, JTS = jts, JTE = jte, KTS = kts, KTE = kte)
+       call mpas_timer_stop('mp_tcwa2')
+
     case default
  end select microp_select
 
@@ -468,6 +513,7 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten
     !ensure that we only call compute_radar_reflectivity() if we are using an MPS that supports
     !the computation of simulated radar reflectivity:
     if(trim(microp_scheme) == "mp_wsm6"     .or. &
+       trim(microp_scheme) == "mp_tcwa2"    .or. &
        trim(microp_scheme) == "mp_thompson" .or. &
        trim(microp_scheme) == "mp_thompson_aerosols") then
        call compute_radar_reflectivity(configs,diag_physics,its,ite)
@@ -541,7 +587,7 @@ subroutine precip_from_MPAS(configs,diag_physics,its,ite)
 
 !variables specific to different cloud microphysics schemes:
  microp_select: select case(trim(microp_scheme))
-    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6")
+    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6","mp_tcwa2")
        do j = jts, jte
        do i = its, ite
           snowncv_p(i,j)    = 0._RKIND
@@ -633,7 +679,7 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite)
 
 !variables specific to different cloud microphysics schemes:
  microp_select: select case(trim(microp_scheme))
-    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6")
+    case ("mp_thompson","mp_thompson_aerosols","mp_wsm6","mp_tcwa2")
        do j = jts,jte
        do i = its,ite
           !time-step precipitation:
@@ -687,6 +733,8 @@ subroutine compute_radar_reflectivity(configs,diag_physics,its,ite)
        call physics_error_fatal('--- calculation of radar reflectivity is not available' // &
                                  'with kessler cloud microphysics')
 
+    case ("mp_tcwa2")
+
     case ("mp_wsm6")
        if(.not.allocated(p1d)  ) allocate(p1d(kts:kte)  )
        if(.not.allocated(t1d)  ) allocate(t1d(kts:kte)  )

src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F

@@ -413,7 +413,7 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi
  radiation_lw_select: select case (trim(radt_lw_scheme))
     case("rrtmg_lw")
        microp_select: select case(microp_scheme)
-          case("mp_thompson","mp_thompson_aerosols","mp_wsm6")
+          case("mp_thompson","mp_thompson_aerosols","mp_wsm6","mp_tcwa2")
              if(config_microp_re) then
                 call mpas_pool_get_array(diag_physics,'re_cloud',re_cloud)
                 call mpas_pool_get_array(diag_physics,'re_ice'  ,re_ice  )
@@ -684,7 +684,7 @@ subroutine radiation_lw_to_MPAS(configs,diag_physics,tend_physics,its,ite)
     case("rrtmg_lw")
 
        microp_select: select case(microp_scheme)
-          case("mp_thompson","mp_thompson_aerosols","mp_wsm6")
+          case("mp_thompson","mp_thompson_aerosols","mp_wsm6","mp_tcwa2")
              call mpas_pool_get_array(diag_physics,'rre_cloud',rre_cloud)
              call mpas_pool_get_array(diag_physics,'rre_ice'  ,rre_ice  )
              call mpas_pool_get_array(diag_physics,'rre_snow' ,rre_snow )

src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F

@@ -494,7 +494,7 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i
        enddo
 
        microp_select: select case(microp_scheme)
-          case("mp_thompson","mp_thompson_aerosols","mp_wsm6")
+          case("mp_thompson","mp_thompson_aerosols","mp_wsm6","mp_tcwa2")
              if(config_microp_re) then
                 call mpas_pool_get_array(diag_physics,'re_cloud',re_cloud)
                 call mpas_pool_get_array(diag_physics,'re_ice'  ,re_ice  )