Merge pull request #1 from adcroft/user/adcroft/add_neverland_winds_t…

…emplate

Added a blank template Neverland_surface_forcing

config_src/solo_driver/MOM_surface_forcing.F90

@@ -90,6 +90,8 @@ module MOM_surface_forcing
 use MOM_variables,           only : surface
 use MESO_surface_forcing,    only : MESO_wind_forcing, MESO_buoyancy_forcing
 use MESO_surface_forcing,    only : MESO_surface_forcing_init, MESO_surface_forcing_CS
+use Neverland_surface_forcing, only : Neverland_wind_forcing, Neverland_buoyancy_forcing
+use Neverland_surface_forcing, only : Neverland_surface_forcing_init, Neverland_surface_forcing_CS
 use user_surface_forcing,    only : USER_wind_forcing, USER_buoyancy_forcing
 use user_surface_forcing,    only : USER_surface_forcing_init, user_surface_forcing_CS
 use user_revise_forcing,     only : user_alter_forcing, user_revise_forcing_init
@@ -207,6 +209,7 @@ module MOM_surface_forcing
   type(user_revise_forcing_CS),  pointer :: urf_CS => NULL()
   type(user_surface_forcing_CS), pointer :: user_forcing_CSp => NULL()
   type(MESO_surface_forcing_CS), pointer :: MESO_forcing_CSp => NULL()
+  type(Neverland_surface_forcing_CS), pointer :: Neverland_forcing_CSp => NULL()
   type(SCM_idealized_hurricane_CS), pointer :: SCM_idealized_hurricane_CSp => NULL()
 
 end type surface_forcing_CS
@@ -265,6 +268,8 @@ subroutine set_forcing(state, fluxes, day_start, day_interval, G, CS)
       call wind_forcing_const(state, fluxes, CS%tau_x0, CS%tau_y0, day_center, G, CS)
     elseif (trim(CS%wind_config) == "MESO") then
       call MESO_wind_forcing(state, fluxes, day_center, G, CS%MESO_forcing_CSp)
+    elseif (trim(CS%wind_config) == "Neverland") then
+      call Neverland_wind_forcing(state, fluxes, day_center, G, CS%Neverland_forcing_CSp)
     elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then
       call SCM_idealized_hurricane_wind_forcing(state, fluxes, day_center, G, CS%SCM_idealized_hurricane_CSp)
     elseif (trim(CS%wind_config) == "USER") then
@@ -318,6 +323,8 @@ subroutine set_forcing(state, fluxes, day_start, day_interval, G, CS)
       call buoyancy_forcing_linear(state, fluxes, day_center, dt, G, CS)
     elseif (trim(CS%buoy_config) == "MESO") then
       call MESO_buoyancy_forcing(state, fluxes, day_center, dt, G, CS%MESO_forcing_CSp)
+    elseif (trim(CS%buoy_config) == "Neverland") then
+      call Neverland_buoyancy_forcing(state, fluxes, day_center, dt, G, CS%Neverland_forcing_CSp)
     elseif (trim(CS%buoy_config) == "USER") then
       call USER_buoyancy_forcing(state, fluxes, day_center, dt, G, CS%user_forcing_CSp)
     elseif (trim(CS%buoy_config) == "NONE") then
@@ -1819,6 +1826,8 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, tracer_flow_CSp)
     call USER_surface_forcing_init(Time, G, param_file, diag, CS%user_forcing_CSp)
   elseif (trim(CS%wind_config) == "MESO" .or. trim(CS%buoy_config) == "MESO" ) then
     call MESO_surface_forcing_init(Time, G, param_file, diag, CS%MESO_forcing_CSp)
+  elseif (trim(CS%wind_config) == "Neverland") then
+    call Neverland_surface_forcing_init(Time, G, param_file, diag, CS%Neverland_forcing_CSp)
   elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then
     call SCM_idealized_hurricane_wind_init(Time, G, param_file, CS%SCM_idealized_hurricane_CSp)
   elseif (trim(CS%wind_config) == "const") then

config_src/solo_driver/Neverland_surface_forcing.F90

@@ -0,0 +1,231 @@
+!> Wind and buoyancy forcing for the Neverland configurations
+module Neverland_surface_forcing
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+use MOM_diag_mediator, only : post_data, query_averaging_enabled
+use MOM_diag_mediator, only : register_diag_field, diag_ctrl
+use MOM_domains, only : pass_var, pass_vector, AGRID
+use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
+use MOM_file_parser, only : get_param, log_version, param_file_type
+use MOM_forcing_type, only : forcing, allocate_forcing_type
+use MOM_grid, only : ocean_grid_type
+use MOM_io, only : file_exists, read_data, slasher
+use MOM_time_manager, only : time_type, operator(+), operator(/), get_time
+use MOM_variables, only : surface
+
+implicit none ; private
+
+public Neverland_wind_forcing
+public Neverland_buoyancy_forcing
+public Neverland_surface_forcing_init
+
+!> This control structure should be used to store any run-time variables
+!! associated with the Neverland forcing.  It can be readily modified
+!! for a specific case, and because it is private there will be no changes
+!! needed in other code (although they will have to be recompiled).
+type, public :: Neverland_surface_forcing_CS ; private
+
+  logical :: use_temperature !< If true, use temperature and salinity.
+  logical :: restorebuoy     !< If true, use restoring surface buoyancy forcing.
+  real :: Rho0               !< The density used in the Boussinesq
+                             !! approximation, in kg m-3.
+  real :: G_Earth            !< The gravitational acceleration in m s-2.
+  real :: flux_const         !<  The restoring rate at the surface, in m s-1.
+  real, dimension(:,:), pointer :: &
+    buoy_restore(:,:) => NULL() !< The pattern to restore buoyancy to.
+  character(len=200) :: inputdir !< The directory where NetCDF input files are.
+  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
+                                   !! timing of diagnostic output.
+  logical :: first_call = .true. !< True until Neverland_buoyancy_forcing has been called
+end type Neverland_surface_forcing_CS
+
+contains
+
+!> Sets the surface wind stresses, fluxes%taux and fluxes%tauy for the
+!! Neverland forcing configuration.
+subroutine Neverland_wind_forcing(state, fluxes, day, G, CS)
+  type(surface),                 intent(inout) :: state !< Fields describing surface state of ocean
+  type(forcing),                 intent(inout) :: fluxes !< Forcing fields.
+  type(time_type),               intent(in)    :: day !< Time used for determining the fluxes.
+  type(ocean_grid_type),         intent(inout) :: G !< Grid structure.
+  type(Neverland_surface_forcing_CS), pointer  :: CS !< Control structure for this module.
+  ! Local variable
+  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
+  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
+  Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
+  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
+  IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
+
+  ! Allocate the forcing arrays, if necessary.
+  call allocate_forcing_type(G, fluxes, stress=.true.)
+
+  !  Set the surface wind stresses, in units of Pa.  A positive taux
+  !  accelerates the ocean to the (pseudo-)east.
+
+  !  The i-loop extends to is-1 so that taux can be used later in the
+  ! calculation of ustar - otherwise the lower bound would be Isq.
+  do j=js,je ; do I=is-1,Ieq
+    fluxes%taux(I,j) = G%mask2dCu(I,j) * 0.0  ! Change this to the desired expression.
+  enddo ; enddo
+  do J=js-1,Jeq ; do i=is,ie
+    fluxes%tauy(i,J) = G%mask2dCv(i,J) * 0.0  ! Change this to the desired expression.
+  enddo ; enddo
+
+  !    Set the surface friction velocity, in units of m s-1.  ustar
+  !  is always positive.
+! if (associated(fluxes%ustar)) then ; do j=js,je ; do i=is,ie
+!   !  This expression can be changed if desired, but need not be.
+!   fluxes%ustar(i,j) = G%mask2dT(i,j) * sqrt(CS%gust_const/CS%Rho0 + &
+!      sqrt(0.5*(fluxes%taux(I-1,j)**2 + fluxes%taux(I,j)**2) + &
+!           0.5*(fluxes%tauy(i,J-1)**2 + fluxes%tauy(i,J)**2))/CS%Rho0)
+! enddo ; enddo ; endif
+
+end subroutine Neverland_wind_forcing
+
+!> Surface fluxes of buoyancy for the Neverland configurations.
+subroutine Neverland_buoyancy_forcing(state, fluxes, day, dt, G, CS)
+  type(surface),                 intent(inout) :: state !< Fields describing surface state of ocean
+  type(forcing),                 intent(inout) :: fluxes !< Forcing fields.
+  type(time_type),               intent(in)    :: day !< Time used for determining the fluxes.
+  real,                          intent(in)    :: dt !< Forcing time step (s).
+  type(ocean_grid_type),         intent(inout) :: G !< Grid structure.
+  type(Neverland_surface_forcing_CS), pointer  :: CS !< Control structure for this module.
+  ! Local variables
+  real :: buoy_rest_const  ! A constant relating density anomalies to the
+                           ! restoring buoyancy flux, in m5 s-3 kg-1.
+  real :: density_restore  ! De
+  integer :: i, j, is, ie, js, je
+  integer :: isd, ied, jsd, jed
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
+  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
+
+  !   When modifying the code, comment out this error message.  It is here
+  ! so that the original (unmodified) version is not accidentally used.
+
+  ! Allocate and zero out the forcing arrays, as necessary.  This portion is
+  ! usually not changed.
+  if (CS%use_temperature) then
+    call MOM_error(FATAL, "Neverland_buoyancy_forcing: " // &
+        "Temperature and salinity mode not coded!" )
+  else
+    ! This is the buoyancy only mode.
+    call alloc_if_needed(fluxes%buoy, isd, ied, jsd, jed)
+  endif
+
+
+  ! MODIFY THE CODE IN THE FOLLOWING LOOPS TO SET THE BUOYANCY FORCING TERMS.
+  if (CS%restorebuoy .and. CS%first_call) then
+    call alloc_if_needed(CS%buoy_restore, isd, ied, jsd, jed)
+    CS%first_call = .false.
+    ! Set CS%buoy_restore(i,j) here
+  endif
+
+  if ( CS%use_temperature ) then
+    call MOM_error(FATAL, "Neverland_buoyancy_surface_forcing: " // &
+      "Temperature/salinity restoring not coded!" )
+  else ! This is the buoyancy only mode.
+    do j=js,je ; do i=is,ie
+      !   fluxes%buoy is the buoyancy flux into the ocean in m2 s-3.  A positive
+      ! buoyancy flux is of the same sign as heating the ocean.
+      fluxes%buoy(i,j) = 0.0 * G%mask2dT(i,j)
+    enddo ; enddo
+  endif
+
+  if (CS%restorebuoy) then
+    if (CS%use_temperature) then
+      call MOM_error(FATAL, "Neverland_buoyancy_surface_forcing: " // &
+        "Temperature/salinity restoring not coded!" )
+    else
+      !   When modifying the code, comment out this error message.  It is here
+      ! so that the original (unmodified) version is not accidentally used.
+
+      ! The -1 is because density has the opposite sign to buoyancy.
+      buoy_rest_const = -1.0 * (CS%G_Earth * CS%flux_const) / CS%Rho0
+      do j=js,je ; do i=is,ie
+       !   Set density_restore to an expression for the surface potential
+       ! density in kg m-3 that is being restored toward.
+        density_restore = 1030.0
+
+        fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * &
+                          (density_restore - state%sfc_density(i,j))
+      enddo ; enddo
+    endif
+  endif                                             ! end RESTOREBUOY
+
+end subroutine Neverland_buoyancy_forcing
+
+!> If ptr is not associated, this routine allocates it with the given size
+!! and zeros out its contents.  This is equivalent to safe_alloc_ptr in
+!! MOM_diag_mediator, but is here so as to be completely transparent.
+subroutine alloc_if_needed(ptr, isd, ied, jsd, jed)
+  real, pointer :: ptr(:,:)
+  integer :: isd, ied, jsd, jed
+  if (.not.ASSOCIATED(ptr)) then
+    allocate(ptr(isd:ied,jsd:jed))
+    ptr(:,:) = 0.0
+  endif
+end subroutine alloc_if_needed
+
+!> Initializes the Neverland control structure.
+subroutine Neverland_surface_forcing_init(Time, G, param_file, diag, CS)
+  type(time_type),         intent(in) :: Time       !< The current model time.
+  type(ocean_grid_type),   intent(in) :: G          !< The ocean's grid structure.
+  type(param_file_type),   intent(in) :: param_file !< A structure indicating the open file to parse for
+                                                    !! model parameter values.
+  type(diag_ctrl), target, intent(in) :: diag       !< A structure that is used to regulate diagnostic output.
+  type(Neverland_surface_forcing_CS), pointer :: CS !< A pointer that is set to point to the control structure
+                                                    !! for this module
+  ! This include declares and sets the variable "version".
+#include "version_variable.h"
+  ! Local variables
+  character(len=40)  :: mod = "Neverland_surface_forcing" ! This module's name.
+
+  if (associated(CS)) then
+    call MOM_error(WARNING, "Neverland_surface_forcing_init called with an associated "// &
+                             "control structure.")
+    return
+  endif
+  allocate(CS)
+  CS%diag => diag
+
+  ! Read all relevant parameters and write them to the model log.
+  call log_version(param_file, mod, version, "")
+  call get_param(param_file, mod, "ENABLE_THERMODYNAMICS", CS%use_temperature, &
+                 "If true, Temperature and salinity are used as state \n"//&
+                 "variables.", default=.true.)
+
+  call get_param(param_file, mod, "G_EARTH", CS%G_Earth, &
+                 "The gravitational acceleration of the Earth.", &
+                 units="m s-2", default = 9.80)
+  call get_param(param_file, mod, "RHO_0", CS%Rho0, &
+                 "The mean ocean density used with BOUSSINESQ true to \n"//&
+                 "calculate accelerations and the mass for conservation \n"//&
+                 "properties, or with BOUSSINSEQ false to convert some \n"//&
+                 "parameters from vertical units of m to kg m-2.", &
+                 units="kg m-3", default=1035.0)
+! call get_param(param_file, mod, "GUST_CONST", CS%gust_const, &
+!                "The background gustiness in the winds.", units="Pa", &
+!                default=0.02)
+
+  call get_param(param_file, mod, "RESTOREBUOY", CS%restorebuoy, &
+                 "If true, the buoyancy fluxes drive the model back \n"//&
+                 "toward some specified surface state with a rate \n"//&
+                 "given by FLUXCONST.", default= .false.)
+
+  if (CS%restorebuoy) then
+    call get_param(param_file, mod, "FLUXCONST", CS%flux_const, &
+                 "The constant that relates the restoring surface fluxes \n"//&
+                 "to the relative surface anomalies (akin to a piston \n"//&
+                 "velocity).  Note the non-MKS units.", units="m day-1", &
+                 fail_if_missing=.true.)
+    ! Convert CS%flux_const from m day-1 to m s-1.
+    CS%flux_const = CS%flux_const / 86400.0
+  endif
+
+end subroutine Neverland_surface_forcing_init
+
+end module Neverland_surface_forcing