From 6404481dc9961fe6ac2563b5aa293eb94f1ecc72 Mon Sep 17 00:00:00 2001
From: Fredrik Jansson <fjansson@abo.fi>
Date: Sat, 31 Aug 2024 21:12:48 +0200
Subject: [PATCH] Optional SST spatial SST profile (RCEMIP II)

The profile is a cosinus wave along x.
Add NAMSURFACE option dsst, default 0, for setting the amplitude
of the perturbation in K. Note the perturbation amplitude is given as
a temperature, internally converted to a thls delta.
---
 src/modsurface.f90  | 13 +++++++++----
 src/modsurfdata.f90 |  1 +
 2 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/modsurface.f90 b/src/modsurface.f90
index a3a37cf11..925a34cd0 100644
--- a/src/modsurface.f90
+++ b/src/modsurface.f90
@@ -91,7 +91,7 @@ subroutine initsurface
       ! Jarvis-Steward related variables
       rsminav, rssoilminav, LAIav, gDav, &
       ! Prescribed values for isurf 2, 3, 4
-      z0, thls, ps, ustin, wtsurf, wqsurf, wsvsurf, &
+      z0, thls, dsst, ps, ustin, wtsurf, wqsurf, wsvsurf, &
       ! Heterogeneous variables
       lhetero, xpatches, ypatches, land_use, loldtable, &
       ! AGS variables
@@ -154,6 +154,7 @@ subroutine initsurface
     call D_MPI_BCAST(wsvsurf(1:nsv),nsv,0,comm3d,mpierr)
     call D_MPI_BCAST(ps         ,1,0,comm3d,mpierr)
     call D_MPI_BCAST(thls       ,1,0,comm3d,mpierr)
+    call D_MPI_BCAST(dsst       ,1,0,comm3d,mpierr)
 
     call D_MPI_BCAST(lhetero                    ,            1,  0, comm3d, mpierr)
     call D_MPI_BCAST(loldtable                  ,            1,  0, comm3d, mpierr)
@@ -863,10 +864,12 @@ end subroutine calc_aerodynamic_resistance
 
   !> Prescribes the skin temperature
   subroutine presc_skin_temperature
-    use modglobal, only: i1, j1
+    use modglobal, only: i1, j1, imax, xsize, dx, pi, rd, cp, pref0
+    use modmpi, only: myidx
     implicit none
 
     integer :: i, j
+    real :: x, dthls
 
     if (lhetero) then
       do j = 2, j1
@@ -875,10 +878,12 @@ subroutine presc_skin_temperature
         end do
       end do
     else
+      dthls = dsst / (ps/pref0)**(rd/cp) ! convert dsst (temperature) to a thls amplitude
       !$acc parallel loop collapse(2) default(present)
       do j = 2, j1
-        do i = 2, i1
-          tskin(i,j) = thls
+         do i = 2, i1
+          x = (myidx * imax + i - 2 + 0.5) * dx
+          tskin(i,j) = thls - dthls/2 * cos(2*pi*x / xsize) ! sst + optional cos-shaped perturbation (RCEMIP II)
         end do
       end do
     end if
diff --git a/src/modsurfdata.f90 b/src/modsurfdata.f90
index 64746f3e8..4b6f4a51e 100644
--- a/src/modsurfdata.f90
+++ b/src/modsurfdata.f90
@@ -231,6 +231,7 @@ module modsurfdata
 
   ! Surface properties in case of prescribed conditions (previous isurf 2, 3 and 4)
   real              :: thls  = -1             !<  Surface liquid water potential temperature [K]
+  real              :: dsst  =  0             !<  Amplitude of SST perturbation [K] (note: in temperature, not thl) (RCEMIP II)
   real              :: qts                    !<  Surface specific humidity [kg/kg]
   real              :: thvs                   !<  Surface virtual temperature [K]
   real(field_r), allocatable :: svs   (:)              !<  Surface scalar concentration [-]