Merge pull request #8 from natgeo-wong/main

Adding radiative tendencies and fixed wind fluxes

Author: natgeo-wong

forcing.f90

@@ -322,6 +322,26 @@ subroutine forcing
 
 endif
 
+!----------------------------------------------------------------------------
+! Temperature Tendency Forcing:
+! Simple Radiative Tendencies taken from Pauluis & Garner [2006]
+
+if(doradtendency.and.time.gt.timelargescale) then
+
+  do k = 1,nzm
+    do j=1,ny
+      do i=1,nx
+        if (tabs(i,j,k)>207.5) then
+            t(i,j,k) = t(i,j,k) - dtn * troptend / 86400
+        else
+            t(i,j,k) = t(i,j,k) + dtn * (200 - tabs(i,j,k)) / (5*86400)
+        end if
+      end do
+    end do
+  end do
+
+endif
+
 !-------------------------------------------------------------------------------
 
 if(.not.dosfcforcing.and.dodynamicocean) call sst_evolve()

params.f90

@@ -134,7 +134,7 @@ module params
   real:: bubble_dtemp = 0.
   real:: bubble_dq = 0.
 
-! The following are added by Kuang Lab at Harvard
+!!!!!  The following are added by Kuang Lab at Harvard !!!!!
 
 ! Options
 logical:: dompiensemble = .false. ! Subdomains defined in domains.f90 are run separately
@@ -146,5 +146,13 @@ module params
 real:: tperturbA = 1.     ! Default perturbation 1 time positive
 real:: qperturbA = 1.
 
+! Radiative tendencies as per Pauluis & Garner [2006]
+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
+logical :: dobulksfcflx = .false.
+real :: bulksfcflx_u = 0.
 
 end module params

setparm.f90

@@ -47,7 +47,9 @@ subroutine setparm
 
 ! Parameters added by Kuang Lab at Harvard
 NAMELIST /KUANG_PARAMS/ dompiensemble, &
-                dolayerperturb, tperturbi, qperturbi, tperturbA, qperturbA ! linear response perturbation: layer by layer (Song Qiyu, 2022)
+                dolayerperturb, tperturbi, qperturbi, tperturbA, qperturbA, &
+                doradtendency, troptend, &
+                dobulksfcflx, bulksfcflx_u
 
 !bloss: Create dummy namelist, so that we can figure out error code
 !       for a mising namelist.  This lets us differentiate between

surface.f90

@@ -35,21 +35,26 @@ subroutine surface()
 
           if(landmask(i,j).eq.0) then
 
-           u_h = 0.5*(u(i+1,j,1)+u(i,j,1)) + ug
-           v_h = 0.5*(v(i,j+YES3D,1)+v(i,j,1)) + vg
-           ta_h = tabs(i,j,1)
-           q_h = qv(i,j,1)
-           t_h = tabs(i,j,1)*prespot(1)
-           t_s = (sstxy(i,j)+t00)*prespoti(1)
-           q_s = salt_factor*qsatw(sstxy(i,j)+t00,presi(1))
-
-           call oceflx(rho(1), u_h, v_h, ta_h, q_h, t_h, z(1)-zi(1), t_s, q_s, &
-                        fluxt0, fluxq0, taux0, tauy0)
-
-           fluxbu(i,j) = taux0/rho(1)
-           fluxbv(i,j) = tauy0/rho(1)
-           fluxbt(i,j) = fluxt0
-           fluxbq(i,j) = fluxq0
+            u_h = 0.5*(u(i+1,j,1)+u(i,j,1)) + ug
+            v_h = 0.5*(v(i,j+YES3D,1)+v(i,j,1)) + vg
+            ta_h = tabs(i,j,1)
+            q_h = qv(i,j,1)
+            t_h = tabs(i,j,1)*prespot(1)
+            t_s = (sstxy(i,j)+t00)*prespoti(1)
+            q_s = salt_factor*qsatw(sstxy(i,j)+t00,presi(1))
+
+            if (dobulksfcflx) then
+              call oceflx(rho(1), u_h, v_h, ta_h, q_h, t_h, z(1)-zi(1), t_s, q_s, &
+                fluxt0, fluxq0, taux0, tauy0)
+            else
+              call oceflx(rho(1), bulksfcflx_u, 0, ta_h, q_h, t_h, z(1)-zi(1), t_s, q_s, &
+                fluxt0, fluxq0, taux0, tauy0)
+            endif
+
+            fluxbu(i,j) = taux0/rho(1)
+            fluxbv(i,j) = tauy0/rho(1)
+            fluxbt(i,j) = fluxt0
+            fluxbq(i,j) = fluxq0
 
           end if