Merge pull request #11 from natgeo-wong/main

Merging natgeo-wong changes to KuangLab

forcing.f90

@@ -11,13 +11,17 @@ subroutine forcing
 
 integer i,j,k,n,nn,m,iz,iday0,iday
 real coef, radtend, dayy
-real tt(nzm,2),qq(nzm,2),uu(nzm,2),vv(nzm,2),ww(nzm,2),pp
+real tt(nzm,2),qq(nzm,2),uu(nzm,2),vv(nzm,2),ww(nzm,2),tp(nzm,2),pp(nzm,2)
+real tpm(nzm)
 real ratio1, ratio2, ratio_t1, ratio_t2
 logical zgrid, pgrid
 
 ! linear response perturbation (Song Qiyu, 2022)
 real, save :: delt_t, delt_q    ! Layer by layer perturbation
 
+! ktrop index for tropopause
+integer :: ktrop
+
 call t_startf ('forcing')
 
 
@@ -63,11 +67,12 @@ subroutine forcing
             coef = (z(iz)-zsnd(i-1,m))/(zsnd(i,m)-zsnd(i-1,m)) 
             tt(iz,n)=tsnd(i-1,m)+(tsnd(i,m)-tsnd(i-1,m))*coef
             if(pgrid) then
-               pp=psnd(i-1,m)+(psnd(i,m)-psnd(i-1,m))*coef
-               tt(iz,n)=tt(iz,n)/((1000./pp)**(rgas/cp))
+               pp(iz,n)=psnd(i-1,m)+(psnd(i,m)-psnd(i-1,m))*coef
+               tt(iz,n)=tt(iz,n)/((1000./pp(iz,n))**(rgas/cp))
             else
                tt(iz,n)=tt(iz,n)/prespotb(iz)
             endif
+            tp(iz,n)=tsnd(i-1,m)+(tsnd(i,m)-tsnd(i-1,m))*coef
             qq(iz,n)=qsnd(i-1,m)+(qsnd(i,m)-qsnd(i-1,m))*coef
             uu(iz,n)=usnd(i-1,m)+(usnd(i,m)-usnd(i-1,m))*coef
             vv(iz,n)=vsnd(i-1,m)+(vsnd(i,m)-vsnd(i-1,m))*coef
@@ -78,11 +83,12 @@ subroutine forcing
       do i = 2,nzsnd
          if(pres(iz).ge.psnd(i,m)) then
             coef = (pres(iz)-psnd(i-1,m))/(psnd(i,m)-psnd(i-1,m))
-            tt(iz,n)=tsnd(i-1,m)+(tsnd(i,m)-tsnd(i-1,m))*coef
-            tt(iz,n)=tt(iz,n)/prespotb(iz)
+            tt(iz,n)=tsnd(i-1,m)+(tsnd(i,m)-tsnd(i-1,m))*coef/prespotb(iz)
+            tp(iz,n)=tsnd(i-1,m)+(tsnd(i,m)-tsnd(i-1,m))*coef
             qq(iz,n)=qsnd(i-1,m)+(qsnd(i,m)-qsnd(i-1,m))*coef
             uu(iz,n)=usnd(i-1,m)+(usnd(i,m)-usnd(i-1,m))*coef
             vv(iz,n)=vsnd(i-1,m)+(vsnd(i,m)-vsnd(i-1,m))*coef
+            pp(iz,n)=psnd(i-1,m)+(psnd(i,m)-psnd(i-1,m))*coef
             goto 11
          endif
       end do
@@ -107,11 +113,13 @@ subroutine forcing
 
 do k=1,nzm
    tg0(k)=tt(k,1)+(tt(k,2)-tt(k,1))*coef
+   tp0(k)=tp(k,1)+(tp(k,2)-tp(k,1))*coef
    qg0(k)=qq(k,1)+(qq(k,2)-qq(k,1))*coef
    qg0(k)=qg0(k)*1.e-3
    ! Note that ug0 and vg0 maybe reset if dolargescale is true)
    ug0(k)=uu(k,1)+(uu(k,2)-uu(k,1))*coef - ug
    vg0(k)=vv(k,1)+(vv(k,2)-vv(k,1))*coef - vg
+   pg0(k)=pp(k,1)+(pp(k,2)-pp(k,1))*coef - vg
 end do
 
 ! ---------------------------------------------------------------
@@ -226,6 +234,88 @@ subroutine forcing
       end do
    end if
 
+   !-------------------------------------------------------------------------------
+   ! Kuang Lab Addition
+   ! Save reference copy of large-scale vertical velocity before modification
+   ! by WTG or scaling techniques, similar to Blossey's version of SAM
+   wsub_ref(1:nzm) = wsub(1:nzm)
+
+   if(dodgw) then
+
+      if(wtgscale_time.gt.0) then
+         twtgmax = (nstop * dt - timelargescale) * wtgscale_time
+         twtg = time-timelargescale
+         if(twtg.gt.twtgmax) then
+         am_wtg_time = am_wtg
+         else
+         am_wtg_time = am_wtg * twtgmax / twtg
+         endif
+      else
+         am_wtg_time = am_wtg
+      endif
+
+      if (dowtg_blossey_etal_JAMES2009) call wtg_james2009(nzm, &
+            100.*pres, tg0, qg0, tabs0, qv0, qn0+qp0, &
+            fcor, lambda_wtg, am_wtg_time, am_wtg_exp, o_wtg, ktrop)
+      if (dowtg_decompdgw) then
+         call wtg_james2009(nzm, &
+            100.*pres, tg0, qg0, tabs0, qv0, qn0+qp0, &
+            fcor, lambda_wtg, am_wtg_time, am_wtg_exp, owtgr, ktrop)
+         call wtg_decompdgw(masterproc, &
+            nzm, nz, z, 100.*pg0, tg0, qg0, tabs0, qv0, qn0+qp0, &
+            lambda_wtg, am_wtg_time, wtgscale_vertmodenum, wtgscale_vertmodescl, &
+            o_wtg, wwtgc, ktrop)
+      end if
+
+      ! convert from omega in Pa/s to wsub in m/s
+      w_wtg(1:nzm) = -o_wtg(1:nzm)/rho(1:nzm)/ggr
+      if (dowtg_decompdgw) wwtgr(1:nzm) = -owtgr(1:nzm)/rho(1:nzm)/ggr
+
+   end if
+
+   if (dotgr) then
+
+      if(wtgscale_time.gt.0) then
+         twtgmax = (nstop * dt - timelargescale) * wtgscale_time
+         twtg = time-timelargescale
+         if(twtg.gt.twtgmax) then
+         tau_wtg_time = tau_wtg
+         else
+         tau_wtg_time = tau_wtg * twtg / twtgmax
+         endif
+      else
+         tau_wtg_time = tau_wtg
+      endif
+
+      do k = 1,nzm
+         tpm(k) = tabs0(k) * prespot(k)
+      end do
+
+      if (dowtg_raymondzeng_QJRMS2005)   call wtg_qjrms2005(masterproc, nzm, nz, z, &
+                              tp0, tpm, tabs0, tau_wtg_time, dowtgLBL, boundstatic, &
+                              dthetadz_min, w_wtg, wwtgr)
+      if (dowtg_hermanraymond_JAMES2014) call wtg_james2014(masterproc, nzm, nz, z, &
+                              tp0, tpm, tabs0, tau_wtg_time, dowtgLBL, boundstatic, &
+                              dthetadz_min, wtgscale_vertmodepwr, w_wtg, wwtgr, wwtgc)
+      if (dowtg_decomptgr)               call wtg_decomptgr(masterproc, nzm, nz, z, &
+                              tp0, tpm, tabs0, tau_wtg_time, &
+                              wtgscale_vertmodenum, wtgscale_vertmodescl, &
+                              dowtgLBL, boundstatic, dthetadz_min, w_wtg, wwtgr, wwtgc)
+
+      ! convert from omega in Pa/s to wsub in m/s
+      o_wtg(1:nzm) = -w_wtg(1:nzm)*rho(1:nzm)*ggr
+      owtgr(1:nzm) = -wwtgr(1:nzm)*rho(1:nzm)*ggr
+
+   end if
+
+   if (dotgr.OR.dodgw) then
+
+      ! add to reference large-scale vertical velocity.
+      wsub(1:nzm) = wsub(1:nzm) + w_wtg(1:nzm)
+      dosubsidence = .true.
+
+   end if
+
    if(dosubsidence) call subsidence()
 
 end if 

hbuf_conditionals_init.f90

@@ -1,18 +1,50 @@
 subroutine hbuf_conditionals_init(count,trcount)
   use vars, only: ncondavg, condavgname, condavglongname
   use rad, only: do_output_clearsky_heating_profiles
+  use params, only: dodgw, dotgr, dowtg_decomp, dowtg_decompdgw, dowtg_decomptgr, &
+       dowtg_raymondzeng_QJRMS2005, dowtg_hermanraymond_JAMES2014
   implicit none
 
   ! Initialize the list of UW statistics variables written in statistics.f90
   integer count,trcount, n
 
+  call add_to_namelist(count,trcount,'TBIAS', &
+         'Absolute temperature bias (model-OBS)','K',0)
+  call add_to_namelist(count,trcount,'QBIAS', &
+         'Water vapor mass mixing ratio bias (model-OBS)','g/kg',0)
+
   if(do_output_clearsky_heating_profiles) then
     call add_to_namelist(count,trcount,'RADQRCLW', &
          'Clearsky longwave heating rate','K/d',0)
     call add_to_namelist(count,trcount,'RADQRCSW', &
          'Clearsky shortwave heating rate','K/d',0)
   end if
 
+  if(dodgw.OR.dotgr) then
+    call add_to_namelist(count,trcount,'WWTG', &
+         'Large-scale W induced by weak temperature gradient approx','m/s',0)
+    call add_to_namelist(count,trcount,'OWTG', &
+         'Large-scale Omega induced by weak temperature gradient approx','Pa/s',0)
+    call add_to_namelist(count,trcount,'WOBSREF', &
+         'Reference Large-scale W Before Modifications by WTG/Scaling','m/s',0)
+  end if
+
+  if(dowtg_raymondzeng_QJRMS2005) then
+    call add_to_namelist(count,trcount,'WWTGRAW', &
+         'Raw (Non-Adjusted) Component of the WTG Vertical Velocity','m/s',0)
+    call add_to_namelist(count,trcount,'OWTGRAW', &
+         'Raw (Non-Adjusted) Component of the WTG Pressure Velocity','Pa/s',0)
+  end if
+
+  if(dowtg_hermanraymond_JAMES2014.OR.dowtg_decomp) then
+    call add_to_namelist(count,trcount,'WTGCOEF', &
+         'Coefficients of Vertical Modes for Decomposed WTG Velocities',' ',0)
+    call add_to_namelist(count,trcount,'WWTGRAW', &
+         'Raw (Non-Adjusted) Component of the WTG Vertical Velocity','m/s',0)
+    call add_to_namelist(count,trcount,'OWTGRAW', &
+         'Raw (Non-Adjusted) Component of the WTG Pressure Velocity','Pa/s',0)
+  end if
+
   !bloss: setup to add an arbitrary number of conditional statistics
   do n = 1,ncondavg
 

main.f90

@@ -9,6 +9,7 @@ program crm
 use tracers
 use movies, only: init_movies
 use params, only: dompiensemble
+use simple_ocean, only: set_sst
 implicit none
 
 integer k, icyc, nn, nstatsteps
@@ -60,6 +61,7 @@ program crm
    call micro_init()  !initialize microphysics
    nstep = 0
    day0 = day
+   if((.not.SLM.and..not.dosfcforcing).AND.nrestart_resetsst) call set_sst()
 else
    print *,'Error: confused by value of NRESTART'
    call task_abort() 

params.f90

@@ -134,7 +134,8 @@ module params
   real:: bubble_dtemp = 0.
   real:: bubble_dq = 0.
 
-!!!!!  The following are added by Kuang Lab at Harvard !!!!!
+!=====================================================
+! Kuang-Lab Additions Begin Here
 
 ! Options
 logical:: dompiensemble = .false. ! Subdomains defined in domains.f90 are run separately
@@ -147,12 +148,66 @@ module params
 real:: qperturbA = 1.
 
 ! Radiative tendencies as per Pauluis & Garner [2006]
+! Added by Nathanael Wong on 2023/07/05
 logical :: doradtendency = .false. 
 real :: troptend = 1.5 ! Convective tendency in Pauluis & Garner [2006]
 
 ! Option to fix wind speed used in calculation of bulk surface fluxes
 ! Taken from Peter Blossey's version of SAM
+! Added by Nathanael Wong on 2023/07/05
 logical :: dobulksfcflx = .false.
 real :: bulksfcflx_u = 0.
 
+! Damped Gravity Wave and Temperature Gradient Relaxation Implementations
+! Added by Nathanael Wong on 2023/07/05
+logical :: dodgw = .false.
+logical :: dotgr = .false.
+logical :: dowtg_decomp = .false.
+real :: wtgscale_time = 0. ! period over which theta relaxation timescale scales from infinity to ttheta_wtg.  Express as fraction of time over which WTG large-scale forcing is implemented.  So if WTG/Large-scale is turned on for 100 days, twtg_scale = 1/4 means that the scaling up to WTG occurs over 25 days.
+
+logical :: dowtg_blossey_etal_JAMES2009  = .false.
+logical :: dowtg_raymondzeng_QJRMS2005   = .false. 
+logical :: dowtg_hermanraymond_JAMES2014 = .false.
+logical :: dowtg_decompdgw = .false.
+logical :: dowtg_decomptgr = .false.
+
+real :: am_wtg = 1. ! momentum damping rate in 1/d -- note must be non-zero.
+real :: am_wtg_exp = 0. ! exponent of p/p0 in momentum damping rate.
+real :: lambda_wtg = 650.e3 ! quarter wavelength in m. default = 650.e3 (=650 km).
+
+real :: tau_wtg = 1. ! Relaxation timescale (in hours) for WTG Approximation of Raymond and Zeng [2005]
+logical :: dowtgLBL = .false.
+logical :: boundstatic = .true. ! Restrict the static stability lower bound to prevent unrealistically large values of w_wtg
+real :: dthetadz_min = 1.e-3 ! if boundstatic = .true., what is the minimum bound? Default from Raymond & Zeng [2005] is 1.e-3 K/km
+real :: wtgscale_vertmodepwr = 1. ! Spectral decomposition power, default is 1 as per Herman and Raymond [2014]
+
+integer :: wtgscale_vertmodenum = 2! number of vertical modes
+real, dimension(2) :: wtgscale_vertmodescl = (/1., 1./) ! strength scaling for vertical modes (number of items = wtgscale_vertmodenum)
+
+! Specify a "island" within which SST is allowed to vary
+! If dosstisland = .false. and  dodynamicocean = .true. the entire domain SST varies
+logical :: dosstislands = .false. ! specify an island within which SST is allowed to vary
+real    :: sstislands_oceanmld = 0. ! "ocean" slab depth, if 0, ocean SST is constant
+real    :: sstislands_landmld  = 0. ! "island" slab depth, set to depth_slab_ocean if 0
+
+! Specify round islands using formula. If readlsm = true, then lsm file will override this
+real    :: sstislands_radius   = 0. ! "island" radii in meters
+integer :: sstislands_nrow = 1.  ! number of island rows
+integer :: sstislands_ncol = 1.  ! number of island columns
+real    :: sstislands_sep  = 0.  ! spacing between island centers, should be at least 2*sstisland_radius
+
+! Alternatively, specify a file to read the land-sea mask data
+! File is a binary file, with variables in this order:
+! (1) nx_lsm, which is an integer specifying number of x points
+! (2) ny_lsm, which is an integer specifying number of y points
+! (3) lsm, which is an array of 1s and 0s, with 0s denoting ocean and 1s denoting land
+logical       :: readlsm = .false. ! read land-sea mask from file
+character(80) :: lsmfile = ""
+
+! If nrestart = 2 and dodynamicocean = false, if nrestart_resetsst = true, set all sst back to tabs_s
+logical :: nrestart_resetsst = .false.
+
+! Kuang-Lab Additions End Here
+!=====================================================
+
 end module params