MOM6
MOM_surface_forcing.F90
Go to the documentation of this file.
1 module mom_surface_forcing
2 !***********************************************************************
3 !* GNU General Public License *
4 !* This file is a part of MOM. *
5 !* *
6 !* MOM is free software; you can redistribute it and/or modify it and *
7 !* are expected to follow the terms of the GNU General Public License *
8 !* as published by the Free Software Foundation; either version 2 of *
9 !* the License, or (at your option) any later version. *
10 !* *
11 !* MOM is distributed in the hope that it will be useful, but WITHOUT *
12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
13 !* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public *
14 !* License for more details. *
15 !* *
16 !* For the full text of the GNU General Public License, *
17 !* write to: Free Software Foundation, Inc., *
18 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
19 !* or see: http://www.gnu.org/licenses/gpl.html *
20 !***********************************************************************
21 !
22 !********+*********+*********+*********+*********+*********+*********+**
23 !* *
24 !* By Robert Hallberg, November 1998 - May 2002 *
25 !* Edited by Stephen Griffies, June 2014 *
26 !* *
27 !* This program contains the subroutines that calculate the *
28 !* surface wind stresses and fluxes of buoyancy or temp/saln and *
29 !* fresh water. These subroutines are called every time step, *
30 !* even if the wind stresses or buoyancy fluxes are constant in time *
31 !* - in that case these routines return quickly without doing *
32 !* anything. In addition, any I/O of forcing fields is controlled *
33 !* by surface_forcing_init, located in this file. *
34 !* *
35 !* set_forcing is a small entry subroutine for the subroutines in *
36 !* this file. It provides the external access to these subroutines. *
37 !* *
38 !* wind_forcing determines the wind stresses and places them into *
39 !* fluxes%taux and fluxes%tauy. Often wind_forcing must be tailored *
40 !* for a particular application - either by specifying file and input *
41 !* variable names or by providing appropriate internal expressions *
42 !* for the stresses within a modified version of USER_wind_forcing. *
43 !* *
44 !* buoyancy_forcing determines the surface fluxes of heat, fresh *
45 !* water and salt, as appropriate. A restoring boundary *
46 !* condition plus a specified flux from a file is implemented here, *
47 !* but a user-provided internal expression can be set by modifying *
48 !* and calling USER_buoyancy_forcing. *
49 !* *
50 !* Macros written all in capital letters are defined in MOM_memory.h. *
51 !* *
52 !* A small fragment of the grid is shown below: *
53 !* *
54 !* j+1 x ^ x ^ x At x: q *
55 !* j+1 > o > o > At ^: v, tauy *
56 !* j x ^ x ^ x At >: u, taux *
57 !* j > o > o > At o: h, fluxes. *
58 !* j-1 x ^ x ^ x *
59 !* i-1 i i+1 At x & ^: *
60 !* i i+1 At > & o: *
61 !* *
62 !* The boundaries always run through q grid points (x). *
63 !* *
64 !********+*********+*********+*********+*********+*********+*********+**
65 !### use MOM_controlled_forcing, only : apply_ctrl_forcing, register_ctrl_forcing_restarts
66 !### use MOM_controlled_forcing, only : controlled_forcing_init, controlled_forcing_end
67 !### use MOM_controlled_forcing, only : ctrl_forcing_CS
68 use mom_constants, only : hlv, hlf
69 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
70 use mom_cpu_clock, only : clock_module
71 use mom_diag_mediator, only : post_data, query_averaging_enabled
72 use mom_diag_mediator, only : diag_ctrl, safe_alloc_ptr
73 use mom_domains, only : pass_var, pass_vector, agrid, to_south, to_west, to_all
74 use mom_domains, only : fill_symmetric_edges, cgrid_ne
75 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
76 use mom_error_handler, only : calltree_enter, calltree_leave
77 use mom_file_parser, only : get_param, log_version, param_file_type
78 use mom_string_functions, only : uppercase
79 use mom_forcing_type, only : forcing, forcing_diags, register_forcing_type_diags, deallocate_forcing_type
80 use mom_forcing_type, only : allocate_forcing_type, deallocate_forcing_type
81 use mom_grid, only : ocean_grid_type
82 use mom_get_input, only : get_mom_input, directories
83 use mom_io, only : file_exists, read_data, slasher, num_timelevels
84 use mom_io, only : east_face, north_face
85 use mom_restart, only : register_restart_field, restart_init, mom_restart_cs
86 use mom_restart, only : restart_init_end, save_restart, restore_state
87 use mom_time_manager, only : time_type, operator(+), operator(/), get_time, set_time
88 use mom_tracer_flow_control, only : call_tracer_set_forcing
89 use mom_tracer_flow_control, only : tracer_flow_control_cs
90 use mom_variables, only : surface
91 use meso_surface_forcing, only : meso_wind_forcing, meso_buoyancy_forcing
92 use meso_surface_forcing, only : meso_surface_forcing_init, meso_surface_forcing_cs
93 use user_surface_forcing, only : user_wind_forcing, user_buoyancy_forcing
94 use user_surface_forcing, only : user_surface_forcing_init, user_surface_forcing_cs
95 use user_revise_forcing, only : user_alter_forcing, user_revise_forcing_init
96 use user_revise_forcing, only : user_revise_forcing_cs
97 use scm_idealized_hurricane, only : scm_idealized_hurricane_wind_init
98 use scm_idealized_hurricane, only : scm_idealized_hurricane_wind_forcing
99 use scm_idealized_hurricane, only : scm_idealized_hurricane_cs
100 use scm_cvmix_tests, only : scm_cvmix_tests_surface_forcing_init
101 use scm_cvmix_tests, only : scm_cvmix_tests_wind_forcing
102 use scm_cvmix_tests, only : scm_cvmix_tests_buoyancy_forcing
103 use scm_cvmix_tests, only : scm_cvmix_tests_cs
104 use bfb_surface_forcing, only : bfb_buoyancy_forcing
105 use bfb_surface_forcing, only : bfb_surface_forcing_init, bfb_surface_forcing_cs
106 
107 use data_override_mod, only : data_override_init, data_override
108 
109 implicit none ; private
110 
111 #include <MOM_memory.h>
112 
113 public set_forcing
114 public surface_forcing_init
115 public forcing_save_restart
116 
117 ! surface_forcing_CS is a structure containing pointers to the forcing fields
118 ! which may be used to drive MOM. All fluxes are positive into the ocean.
119 type, public :: surface_forcing_cs ; private
120 
121  logical :: use_temperature ! if true, temp & salinity used as state variables
122  logical :: restorebuoy ! if true, use restoring surface buoyancy forcing
123  logical :: adiabatic ! if true, no diapycnal mass fluxes or surface buoyancy forcing
124  logical :: variable_winds ! if true, wind stresses vary with time
125  logical :: variable_buoyforce ! if true, buoyancy forcing varies with time.
126  real :: south_lat ! southern latitude of the domain
127  real :: len_lat ! domain length in latitude
128 
129  real :: rho0 ! Boussinesq reference density (kg/m^3)
130  real :: g_earth ! gravitational acceleration (m/s^2)
131  real :: flux_const ! piston velocity for surface restoring (m/s)
132  real :: latent_heat_fusion ! latent heat of fusion (J/kg)
133  real :: latent_heat_vapor ! latent heat of vaporization (J/kg)
134  real :: tau_x0, tau_y0 ! Constant wind stresses used in the WIND_CONFIG="const" forcing
135 
136  real :: gust_const ! constant unresolved background gustiness for ustar (Pa)
137  logical :: read_gust_2d ! if true, use 2-dimensional gustiness supplied from a file
138  real, pointer :: gust(:,:) => null() ! spatially varying unresolved background gustiness (Pa)
139  ! gust is used when read_gust_2d is true.
140 
141  real, pointer :: t_restore(:,:) => null() ! temperature to damp (restore) the SST to (deg C)
142  real, pointer :: s_restore(:,:) => null() ! salinity to damp (restore) the SSS (g/kg)
143  real, pointer :: dens_restore(:,:) => null() ! density to damp (restore) surface density (kg/m^3)
144 
145  integer :: buoy_last_lev_read = -1 ! The last time level read from buoyancy input files
146 
147  real :: gyres_taux_const, gyres_taux_sin_amp, gyres_taux_cos_amp, gyres_taux_n_pis
148  ! if WIND_CONFIG=='gyres' then use
149  ! = A, B, C and n respectively for
150  ! taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L)
151 
152  real :: t_north, t_south ! target temperatures at north and south used in
153  ! buoyancy_forcing_linear
154  real :: s_north, s_south ! target salinity at north and south used in
155  ! buoyancy_forcing_linear
156 
157  logical :: first_call_set_forcing = .true.
158  logical :: archaic_omip_file = .true.
159  logical :: dataoverrideisinitialized = .false.
160 
161  real :: wind_scale ! value by which wind-stresses are scaled, ND.
162  real :: constantheatforcing ! value used for sensible heat flux when buoy_config="const"
163 
164  character(len=8) :: wind_stagger
165  type(tracer_flow_control_cs), pointer :: tracer_flow_csp => null()
166 !### type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL()
167  type(mom_restart_cs), pointer :: restart_csp => null()
168 
169  type(diag_ctrl), pointer :: diag ! structure used to regulate timing of diagnostic output
170 
171  character(len=200) :: inputdir ! directory where NetCDF input files are.
172  character(len=200) :: wind_config ! indicator for wind forcing type (2gyre, USER, FILE..)
173  character(len=200) :: wind_file ! if wind_config is "file", file to use
174  character(len=200) :: buoy_config ! indicator for buoyancy forcing type
175 
176  character(len=200) :: longwave_file = ''
177  character(len=200) :: shortwave_file = ''
178  character(len=200) :: evaporation_file = ''
179  character(len=200) :: sensibleheat_file = ''
180  character(len=200) :: latentheat_file = ''
181 
182  character(len=200) :: rain_file = ''
183  character(len=200) :: snow_file = ''
184  character(len=200) :: runoff_file = ''
185 
186  character(len=200) :: longwaveup_file = ''
187  character(len=200) :: shortwaveup_file = ''
188 
189  character(len=200) :: sstrestore_file = ''
190  character(len=200) :: salinityrestore_file = ''
191 
192  character(len=80) :: & ! Variable names in the input files
193  stress_x_var = '', stress_y_var = '', ustar_var = '', &
194  lw_var = '', sw_var = '', latent_var = '', sens_var = '', evap_var = '', &
195  rain_var = '', snow_var = '', lrunoff_var = '', frunoff_var = '', &
196  sst_restore_var = '', sss_restore_var = ''
197 
198  ! These variables give the number of time levels in the various forcing files.
199  integer :: sw_nlev = -1, lw_nlev = -1, latent_nlev = -1, sens_nlev = -1
200  integer :: wind_nlev = -1, evap_nlev = -1, precip_nlev = -1, runoff_nlev = -1
201  integer :: sst_nlev = -1, sss_nlev = -1
202 
203  ! These variables give the last time level read for the various forcing files.
204  integer :: wind_last_lev = -1
205  integer :: sw_last_lev = -1, lw_last_lev = -1, latent_last_lev = -1
206  integer :: sens_last_lev = -1
207  integer :: evap_last_lev = -1, precip_last_lev = -1, runoff_last_lev = -1
208  integer :: sst_last_lev = -1, sss_last_lev = -1
209 
210  ! Diagnostics handles
211  type(forcing_diags), public :: handles
212 
213  type(user_revise_forcing_cs), pointer :: urf_cs => null()
214  type(user_surface_forcing_cs), pointer :: user_forcing_csp => null()
215  type(bfb_surface_forcing_cs), pointer :: bfb_forcing_csp => null()
216  type(meso_surface_forcing_cs), pointer :: meso_forcing_csp => null()
217  type(scm_idealized_hurricane_cs), pointer :: scm_idealized_hurricane_csp => null()
218  type(scm_cvmix_tests_cs), pointer :: scm_cvmix_tests_csp => null()
219 
220 end type surface_forcing_cs
221 
222 
223 integer :: id_clock_forcing
224 
225 contains
226 
227 subroutine set_forcing(state, fluxes, day_start, day_interval, G, CS)
228  type(surface), intent(inout) :: state
229  type(forcing), intent(inout) :: fluxes
230  type(time_type), intent(in) :: day_start
231  type(time_type), intent(in) :: day_interval
232  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
233  type(surface_forcing_cs), pointer :: CS
234 
235 ! This subroutine calls other subroutines in this file to get surface forcing fields.
236 ! It also allocates and initializes the fields in the flux type.
237 
238 ! Arguments:
239 ! (inout) state = structure describing ocean surface state
240 ! (inout) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
241 ! (in) day_start = Start time of the fluxes
242 ! (in) day_interval = Length of time over which these fluxes applied
243 ! (in) G = ocean grid structure
244 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
245 
246  real :: dt ! length of time in seconds over which fluxes applied
247  type(time_type) :: day_center ! central time of the fluxes.
248  integer :: intdt
249 
250  call cpu_clock_begin(id_clock_forcing)
251  call calltree_enter("set_forcing, MOM_surface_forcing.F90")
252 
253  day_center = day_start + day_interval/2
254  call get_time(day_interval, intdt)
255  dt = real(intdt)
256 
257  ! calls to various wind options
258  if (cs%variable_winds .or. cs%first_call_set_forcing) then
259  if (cs%first_call_set_forcing) call allocate_forcing_type(g, fluxes, stress=.true., ustar=.true.)
260  if (trim(cs%wind_config) == "file") then
261  call wind_forcing_from_file(state, fluxes, day_center, g, cs)
262  elseif (trim(cs%wind_config) == "data_override") then
263  call wind_forcing_by_data_override(state, fluxes, day_center, g, cs)
264  elseif (trim(cs%wind_config) == "2gyre") then
265  call wind_forcing_2gyre(state, fluxes, day_center, g, cs)
266  elseif (trim(cs%wind_config) == "1gyre") then
267  call wind_forcing_1gyre(state, fluxes, day_center, g, cs)
268  elseif (trim(cs%wind_config) == "gyres") then
269  call wind_forcing_gyres(state, fluxes, day_center, g, cs)
270  elseif (trim(cs%wind_config) == "zero") then
271  call wind_forcing_const(state, fluxes, 0., 0., day_center, g, cs)
272  elseif (trim(cs%wind_config) == "const") then
273  call wind_forcing_const(state, fluxes, cs%tau_x0, cs%tau_y0, day_center, g, cs)
274  elseif (trim(cs%wind_config) == "MESO") then
275  call meso_wind_forcing(state, fluxes, day_center, g, cs%MESO_forcing_CSp)
276  elseif (trim(cs%wind_config) == "SCM_ideal_hurr") then
277  call scm_idealized_hurricane_wind_forcing(state, fluxes, day_center, g, cs%SCM_idealized_hurricane_CSp)
278  elseif (trim(cs%wind_config) == "SCM_CVmix_tests") then
279  call scm_cvmix_tests_wind_forcing(state, fluxes, day_center, g, cs%SCM_CVmix_tests_CSp)
280  elseif (trim(cs%wind_config) == "USER") then
281  call user_wind_forcing(state, fluxes, day_center, g, cs%user_forcing_CSp)
282  elseif (cs%variable_winds .and. .not.cs%first_call_set_forcing) then
283  call mom_error(fatal, &
284  "MOM_surface_forcing: Variable winds defined with no wind config")
285  else
286  call mom_error(fatal, &
287  "MOM_surface_forcing:Unrecognized wind config "//trim(cs%wind_config))
288  endif
289  endif
290 
291 ! ! calls to various buoyancy forcing options
292 ! if ((CS%variable_buoyforce .or. CS%first_call_set_forcing) .and. &
293 ! (.not.CS%adiabatic)) then
294 ! if (trim(CS%buoy_config) == "file") then
295 ! if (CS%first_call_set_forcing) call buoyancy_forcing_allocate(fluxes, G, CS)
296 ! elseif (trim(CS%wind_config) == "MESO") then
297 ! call MESO_wind_forcing(state, fluxes, day_center, G, CS%MESO_forcing_CSp)
298 ! elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then
299 ! call SCM_idealized_hurricane_wind_forcing(state, fluxes, day_center, G, CS%SCM_idealized_hurricane_CSp)
300 ! elseif (trim(CS%wind_config) == "USER") then
301 ! call USER_wind_forcing(state, fluxes, day_center, G, CS%user_forcing_CSp)
302 ! elseif (CS%variable_winds .and. .not.CS%first_call_set_forcing) then
303 ! call MOM_error(FATAL, &
304 ! "MOM_surface_forcing: Variable winds defined with no wind config")
305 ! else
306 ! call MOM_error(FATAL, &
307 ! "MOM_surface_forcing:Unrecognized wind config "//trim(CS%wind_config))
308 ! endif
309 ! endif
310 
311  ! calls to various buoyancy forcing options
312  if ((cs%variable_buoyforce .or. cs%first_call_set_forcing) .and. &
313  (.not.cs%adiabatic)) then
314  if (trim(cs%buoy_config) == "file") then
315  if (cs%first_call_set_forcing) call buoyancy_forcing_allocate(fluxes, g, cs)
316  call buoyancy_forcing_from_files(state, fluxes, day_center, dt, g, cs)
317  elseif (trim(cs%buoy_config) == "data_override") then
318  if (cs%first_call_set_forcing) call buoyancy_forcing_allocate(fluxes, g, cs)
319  call buoyancy_forcing_from_data_override(state, fluxes, day_center, dt, g, cs)
320  elseif (trim(cs%buoy_config) == "zero") then
321  if (cs%first_call_set_forcing) call buoyancy_forcing_allocate(fluxes, g, cs)
322  call buoyancy_forcing_zero(state, fluxes, day_center, dt, g, cs)
323  elseif (trim(cs%buoy_config) == "const") then
324  if (cs%first_call_set_forcing) call buoyancy_forcing_allocate(fluxes, g, cs)
325  call buoyancy_forcing_const(state, fluxes, day_center, dt, g, cs)
326  elseif (trim(cs%buoy_config) == "linear") then
327  if (cs%first_call_set_forcing) call buoyancy_forcing_allocate(fluxes, g, cs)
328  call buoyancy_forcing_linear(state, fluxes, day_center, dt, g, cs)
329  elseif (trim(cs%buoy_config) == "MESO") then
330  call meso_buoyancy_forcing(state, fluxes, day_center, dt, g, cs%MESO_forcing_CSp)
331  elseif (trim(cs%buoy_config) == "SCM_CVmix_tests") then
332  call scm_cvmix_tests_buoyancy_forcing(state, fluxes, day_center, g, cs%SCM_CVmix_tests_CSp)
333  elseif (trim(cs%buoy_config) == "USER") then
334  call user_buoyancy_forcing(state, fluxes, day_center, dt, g, cs%user_forcing_CSp)
335  elseif (trim(cs%buoy_config) == "BFB") then
336  call bfb_buoyancy_forcing(state, fluxes, day_center, dt, g, cs%BFB_forcing_CSp)
337  elseif (trim(cs%buoy_config) == "NONE") then
338  call mom_mesg("MOM_surface_forcing: buoyancy forcing has been set to omitted.")
339  elseif (cs%variable_buoyforce .and. .not.cs%first_call_set_forcing) then
340  call mom_error(fatal, &
341  "MOM_surface_forcing: Variable buoy defined with no buoy config.")
342  else
343  call mom_error(fatal, &
344  "MOM_surface_forcing: Unrecognized buoy config "//trim(cs%buoy_config))
345  endif
346  endif
347 
348  if (associated(cs%tracer_flow_CSp)) then
349  if (.not.associated(fluxes%tr_fluxes)) allocate(fluxes%tr_fluxes)
350  call call_tracer_set_forcing(state, fluxes, day_start, day_interval, g, cs%tracer_flow_CSp)
351  endif
352 
353  ! Allow for user-written code to alter the fluxes after all the above
354  call user_alter_forcing(state, fluxes, day_center, g, cs%urf_CS)
355 
356  cs%first_call_set_forcing = .false.
357 
358  call cpu_clock_end(id_clock_forcing)
359  call calltree_leave("set_forcing")
360 
361 end subroutine set_forcing
362 
363 subroutine buoyancy_forcing_allocate(fluxes, G, CS)
364  type(forcing), intent(inout) :: fluxes
365  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
366  type(surface_forcing_cs), pointer :: CS
367 
368 ! This subroutine allocates arrays for buoyancy forcing.
369 
370 ! Arguments:
371 ! (inout) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
372 ! (in) G = ocean grid structure
373 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
374 
375  integer :: isd, ied, jsd, jed
376  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
377 
378 
379  ! this array is zero for all options
380  if (associated(fluxes%p_surf)) then
381  fluxes%p_surf(:,:) = 0.0
382  endif
383 
384  if ( cs%use_temperature ) then
385 
386  ! specify surface freshwater forcing by setting the following (kg/(m^2 * s))
387  ! with convention that positive values for water entering ocean.
388  ! specify surface heat fluxes by setting the following (Watts/m^2)
389  ! with convention that positive values for heat fluxes into the ocean.
390  call allocate_forcing_type(g, fluxes, water=.true., heat=.true.)
391 
392  ! surface restoring fields
393  if (cs%restorebuoy) then
394  call safe_alloc_ptr(cs%T_Restore,isd,ied,jsd,jed)
395  call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
396  call safe_alloc_ptr(cs%S_Restore,isd,ied,jsd,jed)
397  endif
398 
399  else ! CS%use_temperature false.
400 
401  call safe_alloc_ptr(fluxes%buoy,isd,ied,jsd,jed)
402  if (cs%restorebuoy) call safe_alloc_ptr(cs%Dens_Restore,isd,ied,jsd,jed)
403 
404  endif ! endif for CS%use_temperature
405 
406 end subroutine buoyancy_forcing_allocate
407 
408 
409 subroutine wind_forcing_const(state, fluxes, tau_x0, tau_y0, day, G, CS)
410  type(surface), intent(inout) :: state
411  type(forcing), intent(inout) :: fluxes
412  real, intent(in) :: tau_x0
413  real, intent(in) :: tau_y0
414  type(time_type), intent(in) :: day
415  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
416  type(surface_forcing_cs), pointer :: CS
417 
418 ! subroutine sets the surface wind stresses to zero
419 
420 ! Arguments:
421 ! state = structure describing ocean surface state
422 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
423 ! (in) day = time of the fluxes
424 ! (in) G = ocean grid structure
425 ! (in) CS = pointer to control returned by previous surface_forcing_init call
426 
427  real :: mag_tau
428  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
429  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
430 
431  call calltree_enter("wind_forcing_const, MOM_surface_forcing.F90")
432  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
433  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
434  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
435  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
436 
437  !set steady surface wind stresses, in units of Pa.
438  mag_tau = sqrt( tau_x0**2 + tau_y0**2)
439 
440  do j=js,je ; do i=is-1,ieq
441  fluxes%taux(i,j) = tau_x0
442  enddo ; enddo
443 
444  do j=js-1,jeq ; do i=is,ie
445  fluxes%tauy(i,j) = tau_y0
446  enddo ; enddo
447 
448  if (cs%read_gust_2d) then
449  if (associated(fluxes%ustar)) then ; do j=js,je ; do i=is,ie
450  fluxes%ustar(i,j) = sqrt( ( mag_tau + cs%gust(i,j) ) / cs%Rho0 )
451  enddo ; enddo ; endif
452  else
453  if (associated(fluxes%ustar)) then ; do j=js,je ; do i=is,ie
454  fluxes%ustar(i,j) = sqrt( ( mag_tau + cs%gust_const ) / cs%Rho0 )
455  enddo ; enddo ; endif
456  endif
457 
458  call calltree_leave("wind_forcing_const")
459 end subroutine wind_forcing_const
460 
461 
462 subroutine wind_forcing_2gyre(state, fluxes, day, G, CS)
463  type(surface), intent(inout) :: state
464  type(forcing), intent(inout) :: fluxes
465  type(time_type), intent(in) :: day
466  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
467  type(surface_forcing_cs), pointer :: CS
468 
469 ! This subroutine sets the surface wind stresses according to double gyre.
470 
471 ! Arguments:
472 ! state = structure describing ocean surface state
473 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
474 ! (in) day = time of the fluxes
475 ! (in) G = ocean grid structure
476 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
477 
478  real :: PI
479  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
480  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
481 
482  call calltree_enter("wind_forcing_2gyre, MOM_surface_forcing.F90")
483  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
484  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
485  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
486  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
487 
488  !set the steady surface wind stresses, in units of Pa.
489  pi = 4.0*atan(1.0)
490 
491  do j=js,je ; do i=is-1,ieq
492  fluxes%taux(i,j) = 0.1*(1.0 - cos(2.0*pi*(g%geoLatCu(i,j)-cs%South_lat) / &
493  cs%len_lat))
494  enddo ; enddo
495 
496  do j=js-1,jeq ; do i=is,ie
497  fluxes%tauy(i,j) = 0.0
498  enddo ; enddo
499 
500  call calltree_leave("wind_forcing_2gyre")
501 end subroutine wind_forcing_2gyre
502 
503 
504 subroutine wind_forcing_1gyre(state, fluxes, day, G, CS)
505  type(surface), intent(inout) :: state
506  type(forcing), intent(inout) :: fluxes
507  type(time_type), intent(in) :: day
508  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
509  type(surface_forcing_cs), pointer :: CS
510 
511 ! This subroutine sets the surface wind stresses according to single gyre.
512 
513 ! Arguments:
514 ! state = structure describing ocean surface state
515 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
516 ! (in) day = time of the fluxes
517 ! (in) G = ocean grid structure
518 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
519 
520  real :: PI
521  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
522  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
523 
524  call calltree_enter("wind_forcing_1gyre, MOM_surface_forcing.F90")
525  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
526  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
527  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
528  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
529 
530  ! set the steady surface wind stresses, in units of Pa.
531  pi = 4.0*atan(1.0)
532 
533  do j=js,je ; do i=is-1,ieq
534  fluxes%taux(i,j) =-0.2*cos(pi*(g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat)
535  enddo ; enddo
536 
537  do j=js-1,jeq ; do i=is,ie
538  fluxes%tauy(i,j) = 0.0
539  enddo ; enddo
540 
541  call calltree_leave("wind_forcing_1gyre")
542 end subroutine wind_forcing_1gyre
543 
544 
545 subroutine wind_forcing_gyres(state, fluxes, day, G, CS)
546  type(surface), intent(inout) :: state
547  type(forcing), intent(inout) :: fluxes
548  type(time_type), intent(in) :: day
549  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
550  type(surface_forcing_cs), pointer :: CS
551 
552 ! This subroutine sets the surface wind stresses according to gyres.
553 
554 ! Arguments:
555 ! state = structure describing ocean surface state
556 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
557 ! (in) day = time of the fluxes
558 ! (in) G = ocean grid structure
559 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
560 
561  real :: PI, y
562  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
563  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
564 
565  call calltree_enter("wind_forcing_gyres, MOM_surface_forcing.F90")
566  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
567  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
568  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
569  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
570 
571  ! steady surface wind stresses (Pa)
572  pi = 4.0*atan(1.0)
573 
574  do j=jsd,jed ; do i=is-1,iedb
575  y = (g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat
576  fluxes%taux(i,j) = cs%gyres_taux_const + &
577  ( cs%gyres_taux_sin_amp*sin(cs%gyres_taux_n_pis*pi*y) &
578  + cs%gyres_taux_cos_amp*cos(cs%gyres_taux_n_pis*pi*y) )
579  enddo ; enddo
580 
581  do j=js-1,jedb ; do i=isd,ied
582  fluxes%tauy(i,j) = 0.0
583  enddo ; enddo
584 
585  ! set the friction velocity
586  do j=js,je ; do i=is,ie
587  fluxes%ustar(i,j) = sqrt(sqrt(0.5*(fluxes%tauy(i,j-1)*fluxes%tauy(i,j-1) + &
588  fluxes%tauy(i,j)*fluxes%tauy(i,j) + fluxes%taux(i-1,j)*fluxes%taux(i-1,j) + &
589  fluxes%taux(i,j)*fluxes%taux(i,j)))/cs%Rho0 + (cs%gust_const/cs%Rho0))
590  enddo ; enddo
591 
592  call calltree_leave("wind_forcing_gyres")
593 end subroutine wind_forcing_gyres
594 
595 
596 subroutine wind_forcing_from_file(state, fluxes, day, G, CS)
597  type(surface), intent(inout) :: state
598  type(forcing), intent(inout) :: fluxes
599  type(time_type), intent(in) :: day
600  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
601  type(surface_forcing_cs), pointer :: CS
602 
603 ! This subroutine sets the surface wind stresses.
604 
605 ! Arguments:
606 ! state = structure describing ocean surface state
607 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
608 ! (in) day = time of the fluxes
609 ! (in) G = ocean grid structure
610 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
611 
612  character(len=200) :: filename ! The name of the input file.
613  real :: temp_x(szi_(g),szj_(g)) ! Pseudo-zonal and psuedo-meridional
614  real :: temp_y(szi_(g),szj_(g)) ! wind stresses at h-points, in Pa.
615  integer :: time_lev_daily ! The time levels to read for fields with
616  integer :: time_lev_monthly ! daily and montly cycles.
617  integer :: time_lev ! The time level that is used for a field.
618  integer :: days, seconds
619  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
620  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
621  logical :: read_Ustar
622 
623  call calltree_enter("wind_forcing_from_file, MOM_surface_forcing.F90")
624  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
625  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
626  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
627  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
628 
629  call get_time(day,seconds,days)
630  time_lev_daily = days - 365*floor(real(days) / 365.0)
631 
632  if (time_lev_daily < 31) then ; time_lev_monthly = 0
633  else if (time_lev_daily < 59) then ; time_lev_monthly = 1
634  else if (time_lev_daily < 90) then ; time_lev_monthly = 2
635  else if (time_lev_daily < 120) then ; time_lev_monthly = 3
636  else if (time_lev_daily < 151) then ; time_lev_monthly = 4
637  else if (time_lev_daily < 181) then ; time_lev_monthly = 5
638  else if (time_lev_daily < 212) then ; time_lev_monthly = 6
639  else if (time_lev_daily < 243) then ; time_lev_monthly = 7
640  else if (time_lev_daily < 273) then ; time_lev_monthly = 8
641  else if (time_lev_daily < 304) then ; time_lev_monthly = 9
642  else if (time_lev_daily < 334) then ; time_lev_monthly = 10
643  else ; time_lev_monthly = 11
644  endif
645 
646  time_lev_daily = time_lev_daily+1
647  time_lev_monthly = time_lev_monthly+1
648 
649  select case (cs%wind_nlev)
650  case (12) ; time_lev = time_lev_monthly
651  case (365) ; time_lev = time_lev_daily
652  case default ; time_lev = 1
653  end select
654 
655  if (time_lev /= cs%wind_last_lev) then
656  filename = trim(cs%wind_file)
657  read_ustar = (len_trim(cs%ustar_var) > 0)
658 ! if (is_root_pe()) &
659 ! write(*,'("Wind_forcing Reading time level ",I," last was ",I,".")')&
660 ! time_lev-1,CS%wind_last_lev-1
661  select case ( uppercase(cs%wind_stagger(1:1)) )
662  case ("A")
663  temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
664  call read_data(filename,cs%stress_x_var,temp_x(:,:), &
665  domain=g%Domain%mpp_domain,timelevel=time_lev)
666  call read_data(filename,cs%stress_y_var,temp_y(:,:), &
667  domain=g%Domain%mpp_domain,timelevel=time_lev)
668 
669  call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
670  do j=js,je ; do i=is-1,ieq
671  fluxes%taux(i,j) = 0.5 * cs%wind_scale * (temp_x(i,j) + temp_x(i+1,j))
672  enddo ; enddo
673  do j=js-1,jeq ; do i=is,ie
674  fluxes%tauy(i,j) = 0.5 * cs%wind_scale * (temp_y(i,j) + temp_y(i,j+1))
675  enddo ; enddo
676 
677  if (.not.read_ustar) then
678  if (cs%read_gust_2d) then
679  do j=js,je ; do i=is,ie
680  fluxes%ustar(i,j) = sqrt((sqrt(temp_x(i,j)*temp_x(i,j) + &
681  temp_y(i,j)*temp_y(i,j)) + cs%gust(i,j)) / cs%Rho0)
682  enddo ; enddo
683  else
684  do j=js,je ; do i=is,ie
685  fluxes%ustar(i,j) = sqrt(sqrt(temp_x(i,j)*temp_x(i,j) + &
686  temp_y(i,j)*temp_y(i,j))/cs%Rho0 + (cs%gust_const/cs%Rho0))
687  enddo ; enddo
688  endif
689  endif
690  case ("C")
691  if (g%symmetric) then
692  if (.not.associated(g%Domain_aux)) call mom_error(fatal, &
693  " wind_forcing_from_file with C-grid input and symmetric memory "//&
694  " called with a non-associated auxiliary domain in the grid type.")
695  ! Read the data as though symmetric memory were not being used, and
696  ! then translate it appropriately.
697  temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
698  call read_data(filename, cs%stress_x_var, temp_x(:,:), position=east_face, &
699  domain=g%Domain_aux%mpp_domain, timelevel=time_lev)
700  call read_data(filename, cs%stress_y_var, temp_y(:,:), position=north_face, &
701  domain=g%Domain_aux%mpp_domain, timelevel=time_lev)
702 
703  do j=js,je ; do i=is,ie
704  fluxes%taux(i,j) = cs%wind_scale * temp_x(i,j)
705  fluxes%tauy(i,j) = cs%wind_scale * temp_y(i,j)
706  enddo ; enddo
707  call fill_symmetric_edges(fluxes%taux, fluxes%tauy, g%Domain, stagger=cgrid_ne)
708  else
709  call read_data(filename, cs%stress_x_var, fluxes%taux(:,:), &
710  domain=g%Domain%mpp_domain, position=east_face, &
711  timelevel=time_lev)
712  call read_data(filename, cs%stress_y_var, fluxes%tauy(:,:), &
713  domain=g%Domain%mpp_domain, position=north_face, &
714  timelevel=time_lev)
715 
716  if (cs%wind_scale /= 1.0) then
717  do j=js,je ; do i=isq,ieq
718  fluxes%taux(i,j) = cs%wind_scale * fluxes%taux(i,j)
719  enddo ; enddo
720  do j=jsq,jeq ; do i=is,ie
721  fluxes%tauy(i,j) = cs%wind_scale * fluxes%tauy(i,j)
722  enddo ; enddo
723  endif
724  endif
725 
726  call pass_vector(fluxes%taux, fluxes%tauy, g%Domain, to_all)
727  if (.not.read_ustar) then
728  if (cs%read_gust_2d) then
729  do j=js, je ; do i=is, ie
730  fluxes%ustar(i,j) = sqrt((sqrt(0.5*((fluxes%tauy(i,j-1)**2 + &
731  fluxes%tauy(i,j)**2) + (fluxes%taux(i-1,j)**2 + &
732  fluxes%taux(i,j)**2))) + cs%gust(i,j)) / cs%Rho0 )
733  enddo ; enddo
734  else
735  do j=js, je ; do i=is, ie
736  fluxes%ustar(i,j) = sqrt(sqrt(0.5*((fluxes%tauy(i,j-1)**2 + &
737  fluxes%tauy(i,j)**2) + (fluxes%taux(i-1,j)**2 + &
738  fluxes%taux(i,j)**2)))/cs%Rho0 + (cs%gust_const/cs%Rho0))
739  enddo ; enddo
740  endif
741  endif
742  case default
743  call mom_error(fatal, "wind_forcing_from_file: Unrecognized stagger "//&
744  trim(cs%wind_stagger)//" is not 'A' or 'C'.")
745  end select
746 
747  if (read_ustar) then
748  call read_data(filename, cs%Ustar_var, fluxes%ustar(:,:), &
749  domain=g%Domain%mpp_domain, timelevel=time_lev)
750  endif
751 
752  cs%wind_last_lev = time_lev
753 
754  endif ! time_lev /= CS%wind_last_lev
755 
756  call calltree_leave("wind_forcing_from_file")
757 end subroutine wind_forcing_from_file
758 
759 
760 subroutine wind_forcing_by_data_override(state, fluxes, day, G, CS)
761  type(surface), intent(inout) :: state
762  type(forcing), intent(inout) :: fluxes
763  type(time_type), intent(in) :: day
764  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
765  type(surface_forcing_cs), pointer :: CS
766 ! This subroutine sets the surface wind stresses
767 
768 ! Arguments:
769 ! state = structure describing ocean surface state
770 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
771 ! (in) day = time of the fluxes
772 ! (in) G = ocean grid structure
773 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
774 
775  real :: temp_x(szi_(g),szj_(g)) ! Pseudo-zonal and psuedo-meridional
776  real :: temp_y(szi_(g),szj_(g)) ! wind stresses at h-points, in Pa.
777  integer :: i, j, is_in, ie_in, js_in, je_in
778  logical :: read_uStar
779 
780  call calltree_enter("wind_forcing_by_data_override, MOM_surface_forcing.F90")
781 
782  if (.not.cs%dataOverrideIsInitialized) then
783  call allocate_forcing_type(g, fluxes, stress=.true., ustar=.true.)
784  call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
785  cs%dataOverrideIsInitialized = .true.
786  endif
787 
788  is_in = g%isc - g%isd + 1
789  ie_in = g%iec - g%isd + 1
790  js_in = g%jsc - g%jsd + 1
791  je_in = g%jec - g%jsd + 1
792 
793  temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
794  call data_override('OCN', 'taux', temp_x, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
795  call data_override('OCN', 'tauy', temp_y, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
796  call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
797  ! Ignore CS%wind_scale when using data_override ?????
798  do j=g%jsc,g%jec ; do i=g%isc-1,g%IecB
799  fluxes%taux(i,j) = 0.5 * (temp_x(i,j) + temp_x(i+1,j))
800  enddo ; enddo
801  do j=g%jsc-1,g%JecB ; do i=g%isc,g%iec
802  fluxes%tauy(i,j) = 0.5 * (temp_y(i,j) + temp_y(i,j+1))
803  enddo ; enddo
804 
805  read_ustar = (len_trim(cs%ustar_var) > 0) ! Need better control higher up ????
806  if (read_ustar) then
807  call data_override('OCN', 'ustar', fluxes%ustar, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
808  else
809  if (cs%read_gust_2d) then
810  call data_override('OCN', 'gust', cs%gust, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
811  do j=g%jsc,g%jec ; do i=g%isc,g%iec
812  fluxes%ustar(i,j) = sqrt((sqrt(temp_x(i,j)*temp_x(i,j) + &
813  temp_y(i,j)*temp_y(i,j)) + cs%gust(i,j)) / cs%Rho0)
814  enddo ; enddo
815  else
816  do j=g%jsc,g%jec ; do i=g%isc,g%iec
817  fluxes%ustar(i,j) = sqrt(sqrt(temp_x(i,j)*temp_x(i,j) + &
818  temp_y(i,j)*temp_y(i,j))/cs%Rho0 + (cs%gust_const/cs%Rho0))
819  enddo ; enddo
820  endif
821  endif
822 
823  call pass_vector(fluxes%taux, fluxes%tauy, g%Domain, to_all)
824 ! call pass_var(fluxes%ustar, G%Domain, To_All) Not needed ?????
825 
826  call calltree_leave("wind_forcing_by_data_override")
827 end subroutine wind_forcing_by_data_override
828 
829 
830 subroutine buoyancy_forcing_from_files(state, fluxes, day, dt, G, CS)
831  type(surface), intent(inout) :: state
832  type(forcing), intent(inout) :: fluxes
833  type(time_type), intent(in) :: day
834  real, intent(in) :: dt !< The amount of time over which
835  !! the fluxes apply, in s
836  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
837  type(surface_forcing_cs), pointer :: CS
838 
839 ! This subroutine specifies the current surface fluxes of buoyancy
840 ! temperature and fresh water. It may also be modified to add
841 ! surface fluxes of user provided tracers.
842 ! This case has surface buoyancy forcing from input files.
843 
844 ! Arguments:
845 ! state = structure describing ocean surface state
846 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
847 ! (in) day = time of the fluxes
848 ! (in) dt = amount of time over which the fluxes apply
849 ! (in) G = ocean grid structure
850 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
851 
852  real, dimension(SZI_(G),SZJ_(G)) :: &
853  temp, & ! A 2-d temporary work array with various units.
854  SST_anom, & ! Instantaneous sea surface temperature anomalies from a
855  ! target (observed) value, in deg C.
856  sss_anom, & ! Instantaneous sea surface salinity anomalies from a target
857  ! (observed) value, in g kg-1.
858  sss_mean ! A (mean?) salinity about which to normalize local salinity
859  ! anomalies when calculating restorative precipitation
860  ! anomalies, in g kg-1.
861 
862  real :: rhoXcp ! reference density times heat capacity (J/(m^3 * K))
863  real :: Irho0 ! inverse of the Boussinesq reference density (m^3/kg)
864 
865  integer :: time_lev_daily ! time levels to read for fields with daily cycle
866  integer :: time_lev_monthly ! time levels to read for fields with monthly cycle
867  integer :: time_lev ! time level that for a field
868 
869  integer :: days, seconds
870  integer :: i, j, is, ie, js, je
871 
872  call calltree_enter("buoyancy_forcing_from_files, MOM_surface_forcing.F90")
873 
874  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
875 
876  if (cs%use_temperature) rhoxcp = cs%Rho0 * fluxes%C_p
877  irho0 = 1.0/cs%Rho0
878 
879  ! Read the buoyancy forcing file
880  call get_time(day,seconds,days)
881 
882  time_lev_daily = days - 365*floor(real(days) / 365.0)
883 
884  if (time_lev_daily < 31) then ; time_lev_monthly = 0
885  else if (time_lev_daily < 59) then ; time_lev_monthly = 1
886  else if (time_lev_daily < 90) then ; time_lev_monthly = 2
887  else if (time_lev_daily < 120) then ; time_lev_monthly = 3
888  else if (time_lev_daily < 151) then ; time_lev_monthly = 4
889  else if (time_lev_daily < 181) then ; time_lev_monthly = 5
890  else if (time_lev_daily < 212) then ; time_lev_monthly = 6
891  else if (time_lev_daily < 243) then ; time_lev_monthly = 7
892  else if (time_lev_daily < 273) then ; time_lev_monthly = 8
893  else if (time_lev_daily < 304) then ; time_lev_monthly = 9
894  else if (time_lev_daily < 334) then ; time_lev_monthly = 10
895  else ; time_lev_monthly = 11
896  endif
897 
898  time_lev_daily = time_lev_daily +1
899  time_lev_monthly = time_lev_monthly+1
900 
901  if (time_lev_daily /= cs%buoy_last_lev_read) then
902 
903  ! longwave
904  select case (cs%LW_nlev)
905  case (12) ; time_lev = time_lev_monthly
906  case (365) ; time_lev = time_lev_daily
907  case default ; time_lev = 1
908  end select
909  call read_data(cs%longwave_file, cs%LW_var, fluxes%LW(:,:), &
910  domain=g%Domain%mpp_domain, timelevel=time_lev)
911  if (cs%archaic_OMIP_file) then
912  call read_data(cs%longwaveup_file, &
913  "lwup_sfc",temp(:,:), domain=g%Domain%mpp_domain, timelevel=time_lev)
914  do j=js,je ; do i=is,ie ; fluxes%LW(i,j) = fluxes%LW(i,j) - temp(i,j) ; enddo ; enddo
915  endif
916  cs%LW_last_lev = time_lev
917 
918  ! evaporation
919  select case (cs%evap_nlev)
920  case (12) ; time_lev = time_lev_monthly
921  case (365) ; time_lev = time_lev_daily
922  case default ; time_lev = 1
923  end select
924  if (cs%archaic_OMIP_file) then
925  call read_data(cs%evaporation_file, cs%evap_var, temp(:,:), &
926  domain=g%Domain%mpp_domain, timelevel=time_lev)
927  do j=js,je ; do i=is,ie
928  fluxes%latent(i,j) = -cs%latent_heat_vapor*temp(i,j)
929  fluxes%evap(i,j) = -temp(i,j)
930  fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
931  enddo ; enddo
932  else
933  call read_data(cs%evaporation_file, cs%evap_var, fluxes%evap(:,:), &
934  domain=g%Domain%mpp_domain, timelevel=time_lev)
935  endif
936  cs%evap_last_lev = time_lev
937 
938  select case (cs%latent_nlev)
939  case (12) ; time_lev = time_lev_monthly
940  case (365) ; time_lev = time_lev_daily
941  case default ; time_lev = 1
942  end select
943  if (.not.cs%archaic_OMIP_file) then
944  call read_data(cs%latentheat_file, cs%latent_var, fluxes%latent(:,:), &
945  domain=g%Domain%mpp_domain, timelevel=time_lev)
946  do j=js,je ; do i=is,ie
947  fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
948  enddo ; enddo
949  endif
950  cs%latent_last_lev = time_lev
951 
952  select case (cs%sens_nlev)
953  case (12) ; time_lev = time_lev_monthly
954  case (365) ; time_lev = time_lev_daily
955  case default ; time_lev = 1
956  end select
957  if (cs%archaic_OMIP_file) then
958  call read_data(cs%sensibleheat_file, cs%sens_var, temp(:,:), &
959  domain=g%Domain%mpp_domain, timelevel=time_lev)
960  do j=js,je ; do i=is,ie ; fluxes%sens(i,j) = -temp(i,j) ; enddo ; enddo
961  else
962  call read_data(cs%sensibleheat_file, cs%sens_var, fluxes%sens(:,:), &
963  domain=g%Domain%mpp_domain, timelevel=time_lev)
964  endif
965  cs%sens_last_lev = time_lev
966 
967  select case (cs%SW_nlev)
968  case (12) ; time_lev = time_lev_monthly
969  case (365) ; time_lev = time_lev_daily
970  case default ; time_lev = 1
971  end select
972  call read_data(cs%shortwave_file, cs%SW_var, fluxes%sw(:,:), &
973  domain=g%Domain%mpp_domain, timelevel=time_lev)
974  if (cs%archaic_OMIP_file) then
975  call read_data(cs%shortwaveup_file, "swup_sfc", temp(:,:), &
976  domain=g%Domain%mpp_domain, timelevel=time_lev)
977  do j=js,je ; do i=is,ie
978  fluxes%sw(i,j) = fluxes%sw(i,j) - temp(i,j)
979  enddo ; enddo
980  endif
981  cs%SW_last_lev = time_lev
982 
983  select case (cs%precip_nlev)
984  case (12) ; time_lev = time_lev_monthly
985  case (365) ; time_lev = time_lev_daily
986  case default ; time_lev = 1
987  end select
988  call read_data(cs%snow_file, cs%snow_var, &
989  fluxes%fprec(:,:), domain=g%Domain%mpp_domain, timelevel=time_lev)
990  call read_data(cs%rain_file, cs%rain_var, &
991  fluxes%lprec(:,:), domain=g%Domain%mpp_domain, timelevel=time_lev)
992  if (cs%archaic_OMIP_file) then
993  do j=js,je ; do i=is,ie
994  fluxes%lprec(i,j) = fluxes%lprec(i,j) - fluxes%fprec(i,j)
995  enddo ; enddo
996  endif
997  cs%precip_last_lev = time_lev
998 
999  select case (cs%runoff_nlev)
1000  case (12) ; time_lev = time_lev_monthly
1001  case (365) ; time_lev = time_lev_daily
1002  case default ; time_lev = 1
1003  end select
1004  if (cs%archaic_OMIP_file) then
1005  call read_data(cs%runoff_file, cs%lrunoff_var, temp(:,:), &
1006  domain=g%Domain%mpp_domain, timelevel=time_lev)
1007  do j=js,je ; do i=is,ie
1008  fluxes%lrunoff(i,j) = temp(i,j)*g%IareaT(i,j)
1009  enddo ; enddo
1010  call read_data(cs%runoff_file, cs%frunoff_var, temp(:,:), &
1011  domain=g%Domain%mpp_domain, timelevel=time_lev)
1012  do j=js,je ; do i=is,ie
1013  fluxes%frunoff(i,j) = temp(i,j)*g%IareaT(i,j)
1014  enddo ; enddo
1015  else
1016  call read_data(cs%runoff_file, cs%lrunoff_var, fluxes%lrunoff(:,:), &
1017  domain=g%Domain%mpp_domain, timelevel=time_lev)
1018  call read_data(cs%runoff_file, cs%frunoff_var, fluxes%frunoff(:,:), &
1019  domain=g%Domain%mpp_domain, timelevel=time_lev)
1020  endif
1021  cs%runoff_last_lev = time_lev
1022 
1023 ! Read the SST and SSS fields for damping.
1024  if (cs%restorebuoy) then !### .or. associated(CS%ctrl_forcing_CSp)) then
1025  select case (cs%SST_nlev)
1026  case (12) ; time_lev = time_lev_monthly
1027  case (365) ; time_lev = time_lev_daily
1028  case default ; time_lev = 1
1029  end select
1030  call read_data(cs%SSTrestore_file, cs%SST_restore_var, &
1031  cs%T_Restore(:,:), domain=g%Domain%mpp_domain, timelevel=time_lev)
1032  cs%SST_last_lev = time_lev
1033 
1034  select case (cs%SSS_nlev)
1035  case (12) ; time_lev = time_lev_monthly
1036  case (365) ; time_lev = time_lev_daily
1037  case default ; time_lev = 1
1038  end select
1039  call read_data(cs%salinityrestore_file, cs%SSS_restore_var, &
1040  cs%S_Restore(:,:), domain=g%Domain%mpp_domain, timelevel=time_lev)
1041  cs%SSS_last_lev = time_lev
1042  endif
1043  cs%buoy_last_lev_read = time_lev_daily
1044 
1045  ! mask out land points and compute heat content of water fluxes
1046  ! assume liquid precip enters ocean at SST
1047  ! assume frozen precip enters ocean at 0degC
1048  ! assume liquid runoff enters ocean at SST
1049  ! assume solid runoff (calving) enters ocean at 0degC
1050  ! mass leaving the ocean has heat_content determined in MOM_diabatic_driver.F90
1051  do j=js,je ; do i=is,ie
1052  fluxes%evap(i,j) = fluxes%evap(i,j) * g%mask2dT(i,j)
1053  fluxes%lprec(i,j) = fluxes%lprec(i,j) * g%mask2dT(i,j)
1054  fluxes%fprec(i,j) = fluxes%fprec(i,j) * g%mask2dT(i,j)
1055  fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
1056  fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
1057  fluxes%LW(i,j) = fluxes%LW(i,j) * g%mask2dT(i,j)
1058  fluxes%sens(i,j) = fluxes%sens(i,j) * g%mask2dT(i,j)
1059  fluxes%sw(i,j) = fluxes%sw(i,j) * g%mask2dT(i,j)
1060  fluxes%latent(i,j) = fluxes%latent(i,j) * g%mask2dT(i,j)
1061 
1062  fluxes%latent_evap_diag(i,j) = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
1063  fluxes%latent_fprec_diag(i,j) = -fluxes%fprec(i,j)*cs%latent_heat_fusion
1064  fluxes%latent_frunoff_diag(i,j) = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
1065  enddo ; enddo
1066 
1067  endif ! time_lev /= CS%buoy_last_lev_read
1068 
1069 
1070  ! restoring surface boundary fluxes
1071  if (cs%restorebuoy) then
1072 
1073  if (cs%use_temperature) then
1074  do j=js,je ; do i=is,ie
1075  if (g%mask2dT(i,j) > 0) then
1076  fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
1077  ((cs%T_Restore(i,j) - state%SST(i,j)) * rhoxcp * cs%Flux_const)
1078  fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const) * &
1079  (cs%S_Restore(i,j) - state%SSS(i,j)) / &
1080  (0.5*(state%SSS(i,j) + cs%S_Restore(i,j)))
1081  else
1082  fluxes%heat_added(i,j) = 0.0
1083  fluxes%vprec(i,j) = 0.0
1084  endif
1085  enddo ; enddo
1086  else
1087  do j=js,je ; do i=is,ie
1088  if (g%mask2dT(i,j) > 0) then
1089  fluxes%buoy(i,j) = (cs%Dens_Restore(i,j) - state%sfc_density(i,j)) * &
1090  (cs%G_Earth*cs%Flux_const/cs%Rho0)
1091  else
1092  fluxes%buoy(i,j) = 0.0
1093  endif
1094  enddo ; enddo
1095  endif
1096 
1097  else ! not RESTOREBUOY
1098  if (.not.cs%use_temperature) then
1099  call mom_error(fatal, "buoyancy_forcing in MOM_surface_forcing: "// &
1100  "The fluxes need to be defined without RESTOREBUOY.")
1101  endif
1102 
1103  endif ! end RESTOREBUOY
1104 
1105 !### if (associated(CS%ctrl_forcing_CSp)) then
1106 !### do j=js,je ; do i=is,ie
1107 !### SST_anom(i,j) = state%SST(i,j) - CS%T_Restore(i,j)
1108 !### SSS_anom(i,j) = state%SSS(i,j) - CS%S_Restore(i,j)
1109 !### SSS_mean(i,j) = 0.5*(state%SSS(i,j) + CS%S_Restore(i,j))
1110 !### enddo ; enddo
1111 !### call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_added, &
1112 !### fluxes%vprec, day, dt, G, CS%ctrl_forcing_CSp)
1113 !### endif
1114 
1115  call calltree_leave("buoyancy_forcing_from_files")
1116 end subroutine buoyancy_forcing_from_files
1117 
1118 
1119 subroutine buoyancy_forcing_from_data_override(state, fluxes, day, dt, G, CS)
1120  type(surface), intent(inout) :: state
1121  type(forcing), intent(inout) :: fluxes
1122  type(time_type), intent(in) :: day
1123  real, intent(in) :: dt !< The amount of time over which
1124  !! the fluxes apply, in s
1125  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1126  type(surface_forcing_cs), pointer :: CS
1127 
1128 ! This subroutine specifies the current surface fluxes of buoyancy
1129 ! temperature and fresh water. It may also be modified to add
1130 ! surface fluxes of user provided tracers.
1131 ! This case has surface buoyancy forcing from data over-ride.
1132 
1133 ! Arguments:
1134 ! state = structure describing ocean surface state
1135 ! (out) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
1136 ! (in) day = time of the fluxes
1137 ! (in) dt = amount of time over which the fluxes apply
1138 ! (in) G = ocean grid structure
1139 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
1140 
1141  real, dimension(SZI_(G),SZJ_(G)) :: &
1142  temp, & ! A 2-d temporary work array with various units.
1143  SST_anom, & ! Instantaneous sea surface temperature anomalies from a
1144  ! target (observed) value, in deg C.
1145  sss_anom, & ! Instantaneous sea surface salinity anomalies from a target
1146  ! (observed) value, in g kg-1.
1147  sss_mean ! A (mean?) salinity about which to normalize local salinity
1148  ! anomalies when calculating restorative precipitation
1149  ! anomalies, in g kg-1.
1150  real :: rhoXcp ! The mean density times the heat capacity, in J m-3 K-1.
1151  real :: Irho0 ! The inverse of the Boussinesq density, in m3 kg-1.
1152 
1153  integer :: time_lev_daily ! The time levels to read for fields with
1154  integer :: time_lev_monthly ! daily and montly cycles.
1155  integer :: itime_lev ! The time level that is used for a field.
1156 
1157  integer :: days, seconds
1158  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1159  integer :: is_in, ie_in, js_in, je_in
1160 
1161  call calltree_enter("buoyancy_forcing_from_data_override, MOM_surface_forcing.F90")
1162 
1163  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1164  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1165 
1166  if (cs%use_temperature) rhoxcp = cs%Rho0 * fluxes%C_p
1167  irho0 = 1.0/cs%Rho0
1168 
1169  if (.not.cs%dataOverrideIsInitialized) then
1170  call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
1171  cs%dataOverrideIsInitialized = .true.
1172  endif
1173 
1174  is_in = g%isc - g%isd + 1
1175  ie_in = g%iec - g%isd + 1
1176  js_in = g%jsc - g%jsd + 1
1177  je_in = g%jec - g%jsd + 1
1178 
1179  call data_override('OCN', 'lw', fluxes%LW(:,:), day, &
1180  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1181  call data_override('OCN', 'evap', fluxes%evap(:,:), day, &
1182  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1183 
1184  ! note the sign convention
1185  do j=js,je ; do i=is,ie
1186  fluxes%evap(i,j) = -fluxes%evap(i,j) ! Normal convention is positive into the ocean
1187  ! but evap is normally a positive quantity in the files
1188  fluxes%latent(i,j) = cs%latent_heat_vapor*fluxes%evap(i,j)
1189  fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
1190  enddo; enddo
1191 
1192  call data_override('OCN', 'sens', fluxes%sens(:,:), day, &
1193  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1194 
1195  ! note the sign convention
1196  do j=js,je ; do i=is,ie
1197  fluxes%sens(i,j) = -fluxes%sens(i,j) ! Normal convention is positive into the ocean
1198  ! but sensible is normally a positive quantity in the files
1199  enddo; enddo
1200 
1201  call data_override('OCN', 'sw', fluxes%sw(:,:), day, &
1202  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1203 
1204  call data_override('OCN', 'snow', fluxes%fprec(:,:), day, &
1205  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1206 
1207  call data_override('OCN', 'rain', fluxes%lprec(:,:), day, &
1208  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1209 
1210  call data_override('OCN', 'runoff', fluxes%lrunoff(:,:), day, &
1211  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1212 
1213  call data_override('OCN', 'calving', fluxes%frunoff(:,:), day, &
1214  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1215 
1216 ! Read the SST and SSS fields for damping.
1217  if (cs%restorebuoy) then !### .or. associated(CS%ctrl_forcing_CSp)) then
1218  call data_override('OCN', 'SST_restore', cs%T_restore(:,:), day, &
1219  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1220 
1221  call data_override('OCN', 'SSS_restore', cs%S_restore(:,:), day, &
1222  is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1223 
1224  endif
1225 
1226  ! restoring boundary fluxes
1227  if (cs%restorebuoy) then
1228  if (cs%use_temperature) then
1229  do j=js,je ; do i=is,ie
1230  if (g%mask2dT(i,j) > 0) then
1231  fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
1232  ((cs%T_Restore(i,j) - state%SST(i,j)) * rhoxcp * cs%Flux_const)
1233  fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const) * &
1234  (cs%S_Restore(i,j) - state%SSS(i,j)) / &
1235  (0.5*(state%SSS(i,j) + cs%S_Restore(i,j)))
1236  else
1237  fluxes%heat_added(i,j) = 0.0
1238  fluxes%vprec(i,j) = 0.0
1239  endif
1240  enddo ; enddo
1241  else
1242  do j=js,je ; do i=is,ie
1243  if (g%mask2dT(i,j) > 0) then
1244  fluxes%buoy(i,j) = (cs%Dens_Restore(i,j) - state%sfc_density(i,j)) * &
1245  (cs%G_Earth*cs%Flux_const/cs%Rho0)
1246  else
1247  fluxes%buoy(i,j) = 0.0
1248  endif
1249  enddo ; enddo
1250  endif
1251  else ! not RESTOREBUOY
1252  if (.not.cs%use_temperature) then
1253  call mom_error(fatal, "buoyancy_forcing in MOM_surface_forcing: "// &
1254  "The fluxes need to be defined without RESTOREBUOY.")
1255  endif
1256  endif ! end RESTOREBUOY
1257 
1258 
1259  ! mask out land points and compute heat content of water fluxes
1260  ! assume liquid precip enters ocean at SST
1261  ! assume frozen precip enters ocean at 0degC
1262  ! assume liquid runoff enters ocean at SST
1263  ! assume solid runoff (calving) enters ocean at 0degC
1264  ! mass leaving ocean has heat_content determined in MOM_diabatic_driver.F90
1265  do j=js,je ; do i=is,ie
1266  fluxes%evap(i,j) = fluxes%evap(i,j) * g%mask2dT(i,j)
1267  fluxes%lprec(i,j) = fluxes%lprec(i,j) * g%mask2dT(i,j)
1268  fluxes%fprec(i,j) = fluxes%fprec(i,j) * g%mask2dT(i,j)
1269  fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
1270  fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
1271  fluxes%LW(i,j) = fluxes%LW(i,j) * g%mask2dT(i,j)
1272  fluxes%latent(i,j) = fluxes%latent(i,j) * g%mask2dT(i,j)
1273  fluxes%sens(i,j) = fluxes%sens(i,j) * g%mask2dT(i,j)
1274  fluxes%sw(i,j) = fluxes%sw(i,j) * g%mask2dT(i,j)
1275 
1276  fluxes%latent_evap_diag(i,j) = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
1277  fluxes%latent_fprec_diag(i,j) = -fluxes%fprec(i,j)*cs%latent_heat_fusion
1278  fluxes%latent_frunoff_diag(i,j) = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
1279  enddo ; enddo
1280 
1281 
1282 !### if (associated(CS%ctrl_forcing_CSp)) then
1283 !### do j=js,je ; do i=is,ie
1284 !### SST_anom(i,j) = state%SST(i,j) - CS%T_Restore(i,j)
1285 !### SSS_anom(i,j) = state%SSS(i,j) - CS%S_Restore(i,j)
1286 !### SSS_mean(i,j) = 0.5*(state%SSS(i,j) + CS%S_Restore(i,j))
1287 !### enddo ; enddo
1288 !### call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_added, &
1289 !### fluxes%vprec, day, dt, G, CS%ctrl_forcing_CSp)
1290 !### endif
1291 
1292  call calltree_leave("buoyancy_forcing_from_data_override")
1293 end subroutine buoyancy_forcing_from_data_override
1294 
1295 
1296 subroutine buoyancy_forcing_zero(state, fluxes, day, dt, G, CS)
1297  type(surface), intent(inout) :: state
1298  type(forcing), intent(inout) :: fluxes
1299  type(time_type), intent(in) :: day
1300  real, intent(in) :: dt !< The amount of time over which
1301  !! the fluxes apply, in s
1302  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1303  type(surface_forcing_cs), pointer :: CS
1304 
1305 ! This subroutine specifies the current surface fluxes of buoyancy
1306 ! temperature and fresh water. It may also be modified to add
1307 ! surface fluxes of user provided tracers.
1308 ! This case has zero surface buoyancy forcing.
1309 
1310 ! Arguments:
1311 ! (inout) state = structure describing ocean surface state
1312 ! (inout) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
1313 ! (in) day = time of the fluxes
1314 ! (in) dt = amount of time over which the fluxes apply
1315 ! (in) G = ocean grid structure
1316 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
1317 
1318  integer :: i, j, is, ie, js, je
1319 
1320  call calltree_enter("buoyancy_forcing_zero, MOM_surface_forcing.F90")
1321  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1322 
1323  if (cs%use_temperature) then
1324  do j=js,je ; do i=is,ie
1325  fluxes%evap(i,j) = 0.0
1326  fluxes%lprec(i,j) = 0.0
1327  fluxes%fprec(i,j) = 0.0
1328  fluxes%vprec(i,j) = 0.0
1329  fluxes%lrunoff(i,j) = 0.0
1330  fluxes%frunoff(i,j) = 0.0
1331  fluxes%lw(i,j) = 0.0
1332  fluxes%latent(i,j) = 0.0
1333  fluxes%sens(i,j) = 0.0
1334  fluxes%sw(i,j) = 0.0
1335  fluxes%latent_evap_diag(i,j) = 0.0
1336  fluxes%latent_fprec_diag(i,j) = 0.0
1337  fluxes%latent_frunoff_diag(i,j) = 0.0
1338  enddo ; enddo
1339  else
1340  do j=js,je ; do i=is,ie
1341  fluxes%buoy(i,j) = 0.0
1342  enddo ; enddo
1343  endif
1344 
1345  call calltree_leave("buoyancy_forcing_zero")
1346 end subroutine buoyancy_forcing_zero
1347 
1348 
1349 subroutine buoyancy_forcing_const(state, fluxes, day, dt, G, CS)
1350  type(surface), intent(inout) :: state
1351  type(forcing), intent(inout) :: fluxes
1352  type(time_type), intent(in) :: day
1353  real, intent(in) :: dt !< The amount of time over which
1354  !! the fluxes apply, in s
1355  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1356  type(surface_forcing_cs), pointer :: CS
1357 
1358 ! This subroutine specifies the current surface fluxes of buoyancy
1359 ! temperature and fresh water. It may also be modified to add
1360 ! surface fluxes of user provided tracers.
1361 ! We here define a constant surface heat flux.
1362 
1363 ! Arguments:
1364 ! (inout) state = structure describing ocean surface state
1365 ! (inout) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
1366 ! (in) day = time of the fluxes
1367 ! (in) dt = amount of time over which the fluxes apply
1368 ! (in) G = ocean grid structure
1369 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
1370 
1371  integer :: i, j, is, ie, js, je
1372  call calltree_enter("buoyancy_forcing_const, MOM_surface_forcing.F90")
1373  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1374 
1375  if (cs%use_temperature) then
1376  do j=js,je ; do i=is,ie
1377  fluxes%evap(i,j) = 0.0
1378  fluxes%lprec(i,j) = 0.0
1379  fluxes%fprec(i,j) = 0.0
1380  fluxes%vprec(i,j) = 0.0
1381  fluxes%lrunoff(i,j) = 0.0
1382  fluxes%frunoff(i,j) = 0.0
1383  fluxes%lw(i,j) = 0.0
1384  fluxes%latent(i,j) = 0.0
1385  fluxes%sens(i,j) = cs%constantHeatForcing * g%mask2dT(i,j)
1386  fluxes%sw(i,j) = 0.0
1387  fluxes%latent_evap_diag(i,j) = 0.0
1388  fluxes%latent_fprec_diag(i,j) = 0.0
1389  fluxes%latent_frunoff_diag(i,j) = 0.0
1390  enddo ; enddo
1391  else
1392  do j=js,je ; do i=is,ie
1393  fluxes%buoy(i,j) = 0.0
1394  enddo ; enddo
1395  endif
1396 
1397  call calltree_leave("buoyancy_forcing_const")
1398 end subroutine buoyancy_forcing_const
1399 
1400 
1401 subroutine buoyancy_forcing_linear(state, fluxes, day, dt, G, CS)
1402  type(surface), intent(inout) :: state
1403  type(forcing), intent(inout) :: fluxes
1404  type(time_type), intent(in) :: day
1405  real, intent(in) :: dt !< The amount of time over which
1406  !! the fluxes apply, in s
1407  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1408  type(surface_forcing_cs), pointer :: CS
1409 
1410 ! This subroutine specifies the current surface fluxes of buoyancy
1411 ! temperature and fresh water. It may also be modified to add
1412 ! surface fluxes of user provided tracers.
1413 
1414 ! Arguments:
1415 ! (inout) state = structure describing ocean surface state
1416 ! (inout) fluxes = structure with pointers to forcing fields; unused have NULL ptrs
1417 ! (in) day = time of the fluxes
1418 ! (in) dt = amount of time over which the fluxes apply
1419 ! (in) G = ocean grid structure
1420 ! (in) CS = pointer to control struct returned by previous surface_forcing_init call
1421 
1422  real :: y, T_restore, S_restore
1423  integer :: i, j, is, ie, js, je
1424 
1425  call calltree_enter("buoyancy_forcing_linear, MOM_surface_forcing.F90")
1426  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1427 
1428  ! This case has no surface buoyancy forcing.
1429  if (cs%use_temperature) then
1430  do j=js,je ; do i=is,ie
1431  fluxes%evap(i,j) = 0.0
1432  fluxes%lprec(i,j) = 0.0
1433  fluxes%fprec(i,j) = 0.0
1434  fluxes%vprec(i,j) = 0.0
1435  fluxes%lrunoff(i,j) = 0.0
1436  fluxes%frunoff(i,j) = 0.0
1437  fluxes%lw(i,j) = 0.0
1438  fluxes%latent(i,j) = 0.0
1439  fluxes%sens(i,j) = 0.0
1440  fluxes%sw(i,j) = 0.0
1441  fluxes%latent_evap_diag(i,j) = 0.0
1442  fluxes%latent_fprec_diag(i,j) = 0.0
1443  fluxes%latent_frunoff_diag(i,j) = 0.0
1444  enddo ; enddo
1445  else
1446  do j=js,je ; do i=is,ie
1447  fluxes%buoy(i,j) = 0.0
1448  enddo ; enddo
1449  endif
1450 
1451  if (cs%restorebuoy) then
1452  if (cs%use_temperature) then
1453  do j=js,je ; do i=is,ie
1454  y = (g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat
1455  t_restore = cs%T_south + (cs%T_north-cs%T_south)*y
1456  s_restore = cs%S_south + (cs%S_north-cs%S_south)*y
1457  if (g%mask2dT(i,j) > 0) then
1458  fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
1459  ((t_restore - state%SST(i,j)) * ((cs%Rho0 * fluxes%C_p) * cs%Flux_const))
1460  fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const) * &
1461  (s_restore - state%SSS(i,j)) / &
1462  (0.5*(state%SSS(i,j) + s_restore))
1463  else
1464  fluxes%heat_added(i,j) = 0.0
1465  fluxes%vprec(i,j) = 0.0
1466  endif
1467  enddo ; enddo
1468  else
1469  call mom_error(fatal, "buoyancy_forcing_linear in MOM_surface_forcing: "// &
1470  "RESTOREBUOY to linear not written yet.")
1471  !do j=js,je ; do i=is,ie
1472  ! if (G%mask2dT(i,j) > 0) then
1473  ! fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - state%sfc_density(i,j)) * &
1474  ! (CS%G_Earth*CS%Flux_const/CS%Rho0)
1475  ! else
1476  ! fluxes%buoy(i,j) = 0.0
1477  ! endif
1478  !enddo ; enddo
1479  endif
1480  else ! not RESTOREBUOY
1481  if (.not.cs%use_temperature) then
1482  call mom_error(fatal, "buoyancy_forcing_linear in MOM_surface_forcing: "// &
1483  "The fluxes need to be defined without RESTOREBUOY.")
1484  endif
1485  endif ! end RESTOREBUOY
1486 
1487  call calltree_leave("buoyancy_forcing_linear")
1488 end subroutine buoyancy_forcing_linear
1489 
1490 
1491 subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
1492  filename_suffix)
1493  type(surface_forcing_cs), pointer :: CS
1494  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1495  type(time_type), intent(in) :: Time
1496  character(len=*), intent(in) :: directory
1497  logical, optional, intent(in) :: time_stamped
1498  character(len=*), optional, intent(in) :: filename_suffix
1499 
1500 ! Arguments:
1501 ! CS = pointer to control structure from previous surface_forcing_init call
1502 ! (in) G = ocean grid structure
1503 ! (in) Time = model time at this call; needed for mpp_write calls
1504 ! (in, opt) directory = optional directory into which to write these restart files
1505 ! (in, opt) time_stamped = if true, the restart file names include a unique time stamp
1506 ! default is false.
1507 ! (in, opt) filename_suffix = optional suffix (e.g., a time-stamp) to append to the restart fname
1508 
1509  if (.not.associated(cs)) return
1510  if (.not.associated(cs%restart_CSp)) return
1511 
1512  call save_restart(directory, time, g, cs%restart_CSp, time_stamped)
1513 
1514 end subroutine forcing_save_restart
1515 
1516 
1517 subroutine surface_forcing_init(Time, G, param_file, diag, CS, tracer_flow_CSp)
1518  type(time_type), intent(in) :: Time
1519  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1520  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1521  type(diag_ctrl), target, intent(inout) :: diag
1522  type(surface_forcing_cs), pointer :: CS
1523  type(tracer_flow_control_cs), pointer :: tracer_flow_CSp
1524 
1525 ! Arguments:
1526 ! Time = current model time
1527 ! (in) G = ocean grid structure
1528 ! (in) param_file = structure indicating the open file to parse for model parameter values
1529 ! (in) diag = structure used to regulate diagnostic output
1530 ! (in/out) CS = pointer set to point to the control structure for this module
1531 ! (in) tracer_flow_CSp = pointer to the control structure of the tracer flow control module
1532 
1533  type(directories) :: dirs
1534  logical :: new_sim
1535  type(time_type) :: Time_frc
1536 ! This include declares and sets the variable "version".
1537 #include "version_variable.h"
1538  character(len=40) :: mdl = "MOM_surface_forcing" ! This module's name.
1539  character(len=200) :: filename, gust_file ! The name of the gustiness input file.
1540 
1541  if (associated(cs)) then
1542  call mom_error(warning, "surface_forcing_init called with an associated "// &
1543  "control structure.")
1544  return
1545  endif
1546  allocate(cs)
1547 
1548  id_clock_forcing=cpu_clock_id('(Ocean surface forcing)', grain=clock_module)
1549  call cpu_clock_begin(id_clock_forcing)
1550 
1551  cs%diag => diag
1552  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
1553 
1554  ! Read all relevant parameters and write them to the model log.
1555  call log_version(param_file, mdl, version, '')
1556  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
1557  "If true, Temperature and salinity are used as state \n"//&
1558  "variables.", default=.true.)
1559  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, &
1560  "The directory in which all input files are found.", &
1561  default=".")
1562  cs%inputdir = slasher(cs%inputdir)
1563 
1564  call get_param(param_file, mdl, "ADIABATIC", cs%adiabatic, &
1565  "There are no diapycnal mass fluxes if ADIABATIC is \n"//&
1566  "true. This assumes that KD = KDML = 0.0 and that \n"//&
1567  "there is no buoyancy forcing, but makes the model \n"//&
1568  "faster by eliminating subroutine calls.", default=.false.)
1569  call get_param(param_file, mdl, "VARIABLE_WINDS", cs%variable_winds, &
1570  "If true, the winds vary in time after the initialization.", &
1571  default=.true.)
1572  call get_param(param_file, mdl, "VARIABLE_BUOYFORCE", cs%variable_buoyforce, &
1573  "If true, the buoyancy forcing varies in time after the \n"//&
1574  "initialization of the model.", default=.true.)
1575 
1576  call get_param(param_file, mdl, "BUOY_CONFIG", cs%buoy_config, &
1577  "The character string that indicates how buoyancy forcing \n"//&
1578  "is specified. Valid options include (file), (zero), \n"//&
1579  "(linear), (USER), (BFB) and (NONE).", fail_if_missing=.true.)
1580  if (trim(cs%buoy_config) == "file") then
1581  call get_param(param_file, mdl, "ARCHAIC_OMIP_FORCING_FILE", cs%archaic_OMIP_file, &
1582  "If true, use the forcing variable decomposition from \n"//&
1583  "the old German OMIP prescription that predated CORE. If \n"//&
1584  "false, use the variable groupings available from MOM \n"//&
1585  "output diagnostics of forcing variables.", default=.true.)
1586  if (cs%archaic_OMIP_file) then
1587  call get_param(param_file, mdl, "LONGWAVEDOWN_FILE", cs%longwave_file, &
1588  "The file with the downward longwave heat flux, in \n"//&
1589  "variable lwdn_sfc.", fail_if_missing=.true.)
1590  call get_param(param_file, mdl, "LONGWAVEUP_FILE", cs%longwaveup_file, &
1591  "The file with the upward longwave heat flux, in \n"//&
1592  "variable lwup_sfc.", fail_if_missing=.true.)
1593  call get_param(param_file, mdl, "EVAPORATION_FILE", cs%evaporation_file, &
1594  "The file with the evaporative moisture flux, in \n"//&
1595  "variable evap.", fail_if_missing=.true.)
1596  call get_param(param_file, mdl, "SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
1597  "The file with the sensible heat flux, in \n"//&
1598  "variable shflx.", fail_if_missing=.true.)
1599  call get_param(param_file, mdl, "SHORTWAVEUP_FILE", cs%shortwaveup_file, &
1600  "The file with the upward shortwave heat flux.", &
1601  fail_if_missing=.true.)
1602  call get_param(param_file, mdl, "SHORTWAVEDOWN_FILE", cs%shortwave_file, &
1603  "The file with the downward shortwave heat flux.", &
1604  fail_if_missing=.true.)
1605  call get_param(param_file, mdl, "SNOW_FILE", cs%snow_file, &
1606  "The file with the downward frozen precip flux, in \n"//&
1607  "variable snow.", fail_if_missing=.true.)
1608  call get_param(param_file, mdl, "PRECIP_FILE", cs%rain_file, &
1609  "The file with the downward total precip flux, in \n"//&
1610  "variable precip.", fail_if_missing=.true.)
1611  call get_param(param_file, mdl, "FRESHDISCHARGE_FILE", cs%runoff_file, &
1612  "The file with the fresh and frozen runoff/calving fluxes, \n"//&
1613  "invariables disch_w and disch_s.", fail_if_missing=.true.)
1614 
1615  ! These variable names are hard-coded, per the archaic OMIP conventions.
1616  cs%latentheat_file = cs%evaporation_file ; cs%latent_var = "evap"
1617  cs%LW_var = "lwdn_sfc"; cs%SW_var = "swdn_sfc"; cs%sens_var = "shflx"
1618  cs%evap_var = "evap"; cs%rain_var = "precip"; cs%snow_var = "snow"
1619  cs%lrunoff_var = "disch_w"; cs%frunoff_var = "disch_s"
1620 
1621  else
1622  call get_param(param_file, mdl, "LONGWAVE_FILE", cs%longwave_file, &
1623  "The file with the longwave heat flux, in the variable \n"//&
1624  "given by LONGWAVE_FORCING_VAR.", fail_if_missing=.true.)
1625  call get_param(param_file, mdl, "LONGWAVE_FORCING_VAR", cs%LW_var, &
1626  "The variable with the longwave forcing field.", default="LW")
1627 
1628  call get_param(param_file, mdl, "SHORTWAVE_FILE", cs%shortwave_file, &
1629  "The file with the shortwave heat flux, in the variable \n"//&
1630  "given by SHORTWAVE_FORCING_VAR.", fail_if_missing=.true.)
1631  call get_param(param_file, mdl, "SHORTWAVE_FORCING_VAR", cs%SW_var, &
1632  "The variable with the shortwave forcing field.", default="SW")
1633 
1634  call get_param(param_file, mdl, "EVAPORATION_FILE", cs%evaporation_file, &
1635  "The file with the evaporative moisture flux, in the \n"//&
1636  "variable given by EVAP_FORCING_VAR.", fail_if_missing=.true.)
1637  call get_param(param_file, mdl, "EVAP_FORCING_VAR", cs%evap_var, &
1638  "The variable with the evaporative moisture flux.", &
1639  default="evap")
1640 
1641  call get_param(param_file, mdl, "LATENTHEAT_FILE", cs%latentheat_file, &
1642  "The file with the latent heat flux, in the variable \n"//&
1643  "given by LATENT_FORCING_VAR.", fail_if_missing=.true.)
1644  call get_param(param_file, mdl, "LATENT_FORCING_VAR", cs%latent_var, &
1645  "The variable with the latent heat flux.", default="latent")
1646 
1647  call get_param(param_file, mdl, "SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
1648  "The file with the sensible heat flux, in the variable \n"//&
1649  "given by SENSIBLE_FORCING_VAR.", fail_if_missing=.true.)
1650  call get_param(param_file, mdl, "SENSIBLE_FORCING_VAR", cs%sens_var, &
1651  "The variable with the sensible heat flux.", default="sensible")
1652 
1653  call get_param(param_file, mdl, "RAIN_FILE", cs%rain_file, &
1654  "The file with the liquid precipitation flux, in the \n"//&
1655  "variable given by RAIN_FORCING_VAR.", fail_if_missing=.true.)
1656  call get_param(param_file, mdl, "RAIN_FORCING_VAR", cs%rain_var, &
1657  "The variable with the liquid precipitation flux.", &
1658  default="liq_precip")
1659  call get_param(param_file, mdl, "SNOW_FILE", cs%snow_file, &
1660  "The file with the frozen precipitation flux, in the \n"//&
1661  "variable given by SNOW_FORCING_VAR.", fail_if_missing=.true.)
1662  call get_param(param_file, mdl, "SNOW_FORCING_VAR", cs%snow_var, &
1663  "The variable with the frozen precipitation flux.", &
1664  default="froz_precip")
1665 
1666  call get_param(param_file, mdl, "RUNOFF_FILE", cs%runoff_file, &
1667  "The file with the fresh and frozen runoff/calving \n"//&
1668  "fluxes, in variables given by LIQ_RUNOFF_FORCING_VAR \n"//&
1669  "and FROZ_RUNOFF_FORCING_VAR.", fail_if_missing=.true.)
1670  call get_param(param_file, mdl, "LIQ_RUNOFF_FORCING_VAR", cs%lrunoff_var, &
1671  "The variable with the liquid runoff flux.", &
1672  default="liq_runoff")
1673  call get_param(param_file, mdl, "FROZ_RUNOFF_FORCING_VAR", cs%frunoff_var, &
1674  "The variable with the frozen runoff flux.", &
1675  default="froz_runoff")
1676  endif
1677 
1678  call get_param(param_file, mdl, "SSTRESTORE_FILE", cs%SSTrestore_file, &
1679  "The file with the SST toward which to restore in the \n"//&
1680  "variable given by SST_RESTORE_VAR.", fail_if_missing=.true.)
1681  call get_param(param_file, mdl, "SALINITYRESTORE_FILE", cs%salinityrestore_file, &
1682  "The file with the surface salinity toward which to \n"//&
1683  "restore in the variable given by SSS_RESTORE_VAR.", &
1684  fail_if_missing=.true.)
1685 
1686  if (cs%archaic_OMIP_file) then
1687  cs%SST_restore_var = "TEMP" ; cs%SSS_restore_var = "SALT"
1688  else
1689  call get_param(param_file, mdl, "SST_RESTORE_VAR", cs%SST_restore_var, &
1690  "The variable with the SST toward which to restore.", &
1691  default="SST")
1692  call get_param(param_file, mdl, "SSS_RESTORE_VAR", cs%SSS_restore_var, &
1693  "The variable with the SSS toward which to restore.", &
1694  default="SSS")
1695  endif
1696 
1697  ! Add inputdir to the file names.
1698  cs%shortwave_file = trim(cs%inputdir)//trim(cs%shortwave_file)
1699  cs%longwave_file = trim(cs%inputdir)//trim(cs%longwave_file)
1700  cs%sensibleheat_file = trim(cs%inputdir)//trim(cs%sensibleheat_file)
1701  cs%latentheat_file = trim(cs%inputdir)//trim(cs%latentheat_file)
1702  cs%evaporation_file = trim(cs%inputdir)//trim(cs%evaporation_file)
1703  cs%snow_file = trim(cs%inputdir)//trim(cs%snow_file)
1704  cs%rain_file = trim(cs%inputdir)//trim(cs%rain_file)
1705  cs%runoff_file = trim(cs%inputdir)//trim(cs%runoff_file)
1706 
1707  cs%shortwaveup_file = trim(cs%inputdir)//trim(cs%shortwaveup_file)
1708  cs%longwaveup_file = trim(cs%inputdir)//trim(cs%longwaveup_file)
1709 
1710  cs%SSTrestore_file = trim(cs%inputdir)//trim(cs%SSTrestore_file)
1711  cs%salinityrestore_file = trim(cs%inputdir)//trim(cs%salinityrestore_file)
1712  elseif (trim(cs%buoy_config) == "const") then
1713  call get_param(param_file, mdl, "SENSIBLE_HEAT_FLUX", cs%constantHeatForcing, &
1714  "A constant heat forcing (positive into ocean) applied \n"//&
1715  "through the sensible heat flux field. ", &
1716  units='W/m2', fail_if_missing=.true.)
1717  endif
1718  call get_param(param_file, mdl, "WIND_CONFIG", cs%wind_config, &
1719  "The character string that indicates how wind forcing \n"//&
1720  "is specified. Valid options include (file), (2gyre), \n"//&
1721  "(1gyre), (gyres), (zero), and (USER).", fail_if_missing=.true.)
1722  if (trim(cs%wind_config) == "file") then
1723  call get_param(param_file, mdl, "WIND_FILE", cs%wind_file, &
1724  "The file in which the wind stresses are found in \n"//&
1725  "variables STRESS_X and STRESS_Y.", fail_if_missing=.true.)
1726  call get_param(param_file, mdl, "WINDSTRESS_X_VAR",cs%stress_x_var, &
1727  "The name of the x-wind stress variable in WIND_FILE.", &
1728  default="STRESS_X")
1729  call get_param(param_file, mdl, "WINDSTRESS_Y_VAR", cs%stress_y_var, &
1730  "The name of the y-wind stress variable in WIND_FILE.", &
1731  default="STRESS_Y")
1732  call get_param(param_file, mdl, "WINDSTRESS_STAGGER",cs%wind_stagger, &
1733  "A character indicating how the wind stress components \n"//&
1734  "are staggered in WIND_FILE. This may be A or C for now.", &
1735  default="A")
1736  call get_param(param_file, mdl, "WINDSTRESS_SCALE", cs%wind_scale, &
1737  "A value by which the wind stresses in WIND_FILE are rescaled.", &
1738  default=1.0, units="nondim")
1739  call get_param(param_file, mdl, "USTAR_FORCING_VAR", cs%ustar_var, &
1740  "The name of the friction velocity variable in WIND_FILE \n"//&
1741  "or blank to get ustar from the wind stresses plus the \n"//&
1742  "gustiness.", default=" ", units="nondim")
1743  cs%wind_file = trim(cs%inputdir) // trim(cs%wind_file)
1744  endif
1745  if (trim(cs%wind_config) == "gyres") then
1746  call get_param(param_file, mdl, "TAUX_CONST", cs%gyres_taux_const, &
1747  "With the gyres wind_config, the constant offset in the \n"//&
1748  "zonal wind stress profile: \n"//&
1749  " A in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1750  units="Pa", default=0.0)
1751  call get_param(param_file, mdl, "TAUX_SIN_AMP",cs%gyres_taux_sin_amp, &
1752  "With the gyres wind_config, the sine amplitude in the \n"//&
1753  "zonal wind stress profile: \n"//&
1754  " B in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1755  units="Pa", default=0.0)
1756  call get_param(param_file, mdl, "TAUX_COS_AMP",cs%gyres_taux_cos_amp, &
1757  "With the gyres wind_config, the cosine amplitude in \n"//&
1758  "the zonal wind stress profile: \n"//&
1759  " C in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1760  units="Pa", default=0.0)
1761  call get_param(param_file, mdl, "TAUX_N_PIS",cs%gyres_taux_n_pis, &
1762  "With the gyres wind_config, the number of gyres in \n"//&
1763  "the zonal wind stress profile: \n"//&
1764  " n in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1765  units="nondim", default=0.0)
1766  endif
1767  if ((trim(cs%wind_config) == "2gyre") .or. &
1768  (trim(cs%wind_config) == "1gyre") .or. &
1769  (trim(cs%wind_config) == "gyres") .or. &
1770  (trim(cs%buoy_config) == "linear")) then
1771  cs%south_lat = g%south_lat
1772  cs%len_lat = g%len_lat
1773  endif
1774  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
1775  "The mean ocean density used with BOUSSINESQ true to \n"//&
1776  "calculate accelerations and the mass for conservation \n"//&
1777  "properties, or with BOUSSINSEQ false to convert some \n"//&
1778  "parameters from vertical units of m to kg m-2.", &
1779  units="kg m-3", default=1035.0)
1780  call get_param(param_file, mdl, "RESTOREBUOY", cs%restorebuoy, &
1781  "If true, the buoyancy fluxes drive the model back \n"//&
1782  "toward some specified surface state with a rate \n"//&
1783  "given by FLUXCONST.", default= .false.)
1784  call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
1785  "The latent heat of fusion.", units="J/kg", default=hlf)
1786  call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
1787  "The latent heat of fusion.", units="J/kg", default=hlv)
1788  if (cs%restorebuoy) then
1789  call get_param(param_file, mdl, "FLUXCONST", cs%Flux_const, &
1790  "The constant that relates the restoring surface fluxes \n"//&
1791  "to the relative surface anomalies (akin to a piston \n"//&
1792  "velocity). Note the non-MKS units.", units="m day-1", &
1793  fail_if_missing=.true.)
1794  ! Convert CS%Flux_const from m day-1 to m s-1.
1795  cs%Flux_const = cs%Flux_const / 86400.0
1796  if (trim(cs%buoy_config) == "linear") then
1797  call get_param(param_file, mdl, "SST_NORTH", cs%T_north, &
1798  "With buoy_config linear, the sea surface temperature \n"//&
1799  "at the northern end of the domain toward which to \n"//&
1800  "to restore.", units="deg C", default=0.0)
1801  call get_param(param_file, mdl, "SST_SOUTH", cs%T_south, &
1802  "With buoy_config linear, the sea surface temperature \n"//&
1803  "at the southern end of the domain toward which to \n"//&
1804  "to restore.", units="deg C", default=0.0)
1805  call get_param(param_file, mdl, "SSS_NORTH", cs%S_north, &
1806  "With buoy_config linear, the sea surface salinity \n"//&
1807  "at the northern end of the domain toward which to \n"//&
1808  "to restore.", units="PSU", default=35.0)
1809  call get_param(param_file, mdl, "SSS_SOUTH", cs%S_south, &
1810  "With buoy_config linear, the sea surface salinity \n"//&
1811  "at the southern end of the domain toward which to \n"//&
1812  "to restore.", units="PSU", default=35.0)
1813  endif
1814  endif
1815  call get_param(param_file, mdl, "G_EARTH", cs%G_Earth, &
1816  "The gravitational acceleration of the Earth.", &
1817  units="m s-2", default = 9.80)
1818 
1819  call get_param(param_file, mdl, "GUST_CONST", cs%gust_const, &
1820  "The background gustiness in the winds.", units="Pa", &
1821  default=0.02)
1822  call get_param(param_file, mdl, "READ_GUST_2D", cs%read_gust_2d, &
1823  "If true, use a 2-dimensional gustiness supplied from \n"//&
1824  "an input file", default=.false.)
1825  if (cs%read_gust_2d) then
1826  call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
1827  "The file in which the wind gustiness is found in \n"//&
1828  "variable gustiness.", fail_if_missing=.true.)
1829  call safe_alloc_ptr(cs%gust,g%isd,g%ied,g%jsd,g%jed)
1830  filename = trim(cs%inputdir) // trim(gust_file)
1831  call read_data(filename,'gustiness',cs%gust,domain=g%domain%mpp_domain, &
1832  timelevel=1) ! units should be Pa
1833  endif
1834 
1835 ! All parameter settings are now known.
1836 
1837  if (trim(cs%wind_config) == "USER" .or. trim(cs%buoy_config) == "USER" ) then
1838  call user_surface_forcing_init(time, g, param_file, diag, cs%user_forcing_CSp)
1839  elseif (trim(cs%buoy_config) == "BFB" ) then
1840  call bfb_surface_forcing_init(time, g, param_file, diag, cs%BFB_forcing_CSp)
1841  elseif (trim(cs%wind_config) == "MESO" .or. trim(cs%buoy_config) == "MESO" ) then
1842  call meso_surface_forcing_init(time, g, param_file, diag, cs%MESO_forcing_CSp)
1843  elseif (trim(cs%wind_config) == "SCM_ideal_hurr") then
1844  call scm_idealized_hurricane_wind_init(time, g, param_file, cs%SCM_idealized_hurricane_CSp)
1845  elseif (trim(cs%wind_config) == "const") then
1846  call get_param(param_file, mdl, "CONST_WIND_TAUX", cs%tau_x0, &
1847  "With wind_config const, this is the constant zonal\n"//&
1848  "wind-stress", units="Pa", fail_if_missing=.true.)
1849  call get_param(param_file, mdl, "CONST_WIND_TAUY", cs%tau_y0, &
1850  "With wind_config const, this is the constant zonal\n"//&
1851  "wind-stress", units="Pa", fail_if_missing=.true.)
1852  elseif (trim(cs%wind_config) == "SCM_CVmix_tests" .or. &
1853  trim(cs%buoy_config) == "SCM_CVmix_tests") then
1854  call scm_cvmix_tests_surface_forcing_init(time, g, param_file, cs%SCM_CVmix_tests_CSp)
1855  cs%SCM_CVmix_tests_CSp%Rho0 = cs%Rho0 !copy reference density for pass
1856  endif
1857 
1858  call register_forcing_type_diags(time, diag, cs%use_temperature, cs%handles)
1859 
1860  ! Set up any restart fields associated with the forcing.
1861  call restart_init(param_file, cs%restart_CSp, "MOM_forcing.res")
1862 !### call register_ctrl_forcing_restarts(G, param_file, CS%ctrl_forcing_CSp, &
1863 !### CS%restart_CSp)
1864  call restart_init_end(cs%restart_CSp)
1865 
1866  if (associated(cs%restart_CSp)) then
1867  call get_mom_input(dirs=dirs)
1868 
1869  new_sim = .false.
1870  if ((dirs%input_filename(1:1) == 'n') .and. &
1871  (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1872  if (.not.new_sim) then
1873  call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
1874  g, cs%restart_CSp)
1875  endif
1876  endif
1877 
1878  ! Determine how many time levels are in each forcing variable.
1879  if (trim(cs%buoy_config) == "file") then
1880  cs%SW_nlev = num_timelevels(cs%shortwave_file, cs%SW_var, min_dims=3)
1881  cs%LW_nlev = num_timelevels(cs%longwave_file, cs%LW_var, min_dims=3)
1882  cs%latent_nlev = num_timelevels(cs%latentheat_file, cs%latent_var, 3)
1883  cs%sens_nlev = num_timelevels(cs%sensibleheat_file, cs%sens_var, min_dims=3)
1884 
1885  cs%evap_nlev = num_timelevels(cs%evaporation_file, cs%evap_var, min_dims=3)
1886  cs%precip_nlev = num_timelevels(cs%rain_file, cs%rain_var, min_dims=3)
1887  cs%runoff_nlev = num_timelevels(cs%runoff_file, cs%lrunoff_var, 3)
1888 
1889  cs%SST_nlev = num_timelevels(cs%SSTrestore_file, cs%SST_restore_var, 3)
1890  cs%SSS_nlev = num_timelevels(cs%salinityrestore_file, cs%SSS_restore_var, 3)
1891  endif
1892 
1893  if (trim(cs%wind_config) == "file") &
1894  cs%wind_nlev = num_timelevels(cs%wind_file, cs%stress_x_var, min_dims=3)
1895 
1896 !### call controlled_forcing_init(Time, G, param_file, diag, CS%ctrl_forcing_CSp)
1897 
1898  call user_revise_forcing_init(param_file, cs%urf_CS)
1899 
1900  call cpu_clock_end(id_clock_forcing)
1901 end subroutine surface_forcing_init
1902 
1903 
1904 subroutine surface_forcing_end(CS, fluxes)
1905  type(surface_forcing_cs), pointer :: CS
1906  type(forcing), optional, intent(inout) :: fluxes
1907 ! Arguments: CS - A pointer to the control structure returned by a previous
1908 ! call to surface_forcing_init, it will be deallocated here.
1909 ! (inout) fluxes - A structure containing pointers to any possible
1910 ! forcing fields. Unused fields have NULL ptrs.
1911 
1912  if (present(fluxes)) call deallocate_forcing_type(fluxes)
1913 
1914 !### call controlled_forcing_end(CS%ctrl_forcing_CSp)
1915 
1916  if (associated(cs)) deallocate(cs)
1917  cs => null()
1918 
1919 end subroutine surface_forcing_end
1920 
1921 end module mom_surface_forcing
mom_surface_forcing::set_forcing
subroutine, public set_forcing(state, fluxes, day_start, day_interval, G, CS)
Definition: MOM_surface_forcing.F90:228
bfb_surface_forcing::bfb_surface_forcing_cs
Definition: BFB_surface_forcing.F90:52
mom_surface_forcing::wind_forcing_2gyre
subroutine wind_forcing_2gyre(state, fluxes, day, G, CS)
Definition: MOM_surface_forcing.F90:463
mom_variables::surface
The following structure contains pointers to various fields which may be used describe the surface st...
Definition: MOM_variables.F90:49
scm_idealized_hurricane::scm_idealized_hurricane_wind_init
subroutine, public scm_idealized_hurricane_wind_init(Time, G, param_file, CS)
Initializes wind profile for the SCM idealized hurricane example.
Definition: SCM_idealized_hurricane.F90:97
bfb_surface_forcing::bfb_buoyancy_forcing
subroutine, public bfb_buoyancy_forcing(state, fluxes, day, dt, G, CS)
Definition: BFB_surface_forcing.F90:89
mom_time_manager
Definition: MOM_time_manager.F90:1
mom_surface_forcing::wind_forcing_by_data_override
subroutine wind_forcing_by_data_override(state, fluxes, day, G, CS)
Definition: MOM_surface_forcing.F90:761
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
user_revise_forcing::user_revise_forcing_cs
Definition: user_revise_forcing.F90:45
mom_get_input::get_mom_input
subroutine, public get_mom_input(param_file, dirs, check_params)
Definition: MOM_get_input.F90:55
mom_domains::to_all
integer, parameter, public to_all
Definition: MOM_domains.F90:138
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
mom_surface_forcing::surface_forcing_cs
Definition: MOM_surface_forcing.F90:119
meso_surface_forcing::meso_buoyancy_forcing
subroutine, public meso_buoyancy_forcing(state, fluxes, day, dt, G, CS)
Definition: MESO_surface_forcing.F90:178
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:139
mom_get_input
Definition: MOM_get_input.F90:1
scm_cvmix_tests::scm_cvmix_tests_wind_forcing
subroutine, public scm_cvmix_tests_wind_forcing(state, fluxes, day, G, CS)
Definition: SCM_CVmix_tests.F90:189
mom_file_parser::param_file_type
Definition: MOM_file_parser.F90:91
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
scm_idealized_hurricane::scm_idealized_hurricane_cs
Container for parameters describing idealized wind structure.
Definition: SCM_idealized_hurricane.F90:27
mom_cpu_clock
Definition: MOM_cpu_clock.F90:1
mom_restart::register_restart_field
Definition: MOM_restart.F90:135
mom_forcing_type::allocate_forcing_type
subroutine, public allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, press, iceberg)
Conditionally allocate fields within the forcing type.
Definition: MOM_forcing_type.F90:2334
mom_surface_forcing
Definition: MOM_surface_forcing.F90:1
mom_error_handler::calltree_leave
subroutine, public calltree_leave(mesg)
Writes a message about leaving a subroutine if call tree reporting is active.
Definition: MOM_error_handler.F90:147
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_file_parser
Definition: MOM_file_parser.F90:1
mom_domains::fill_symmetric_edges
Definition: MOM_domains.F90:96
mom_surface_forcing::buoyancy_forcing_from_files
subroutine buoyancy_forcing_from_files(state, fluxes, day, dt, G, CS)
Definition: MOM_surface_forcing.F90:831
mom_surface_forcing::wind_forcing_from_file
subroutine wind_forcing_from_file(state, fluxes, day, G, CS)
Definition: MOM_surface_forcing.F90:597
mom_domains::pass_vector
Definition: MOM_domains.F90:69
mom_constants
Definition: MOM_constants.F90:1
mom_tracer_flow_control
Definition: MOM_tracer_flow_control.F90:1
user_revise_forcing::user_alter_forcing
subroutine, public user_alter_forcing(state, fluxes, day, G, CS)
This subroutine sets the surface wind stresses.
Definition: user_revise_forcing.F90:53
mom_surface_forcing::buoyancy_forcing_from_data_override
subroutine buoyancy_forcing_from_data_override(state, fluxes, day, dt, G, CS)
Definition: MOM_surface_forcing.F90:1120
mom_surface_forcing::id_clock_forcing
integer id_clock_forcing
Definition: MOM_surface_forcing.F90:223
mom_surface_forcing::wind_forcing_1gyre
subroutine wind_forcing_1gyre(state, fluxes, day, G, CS)
Definition: MOM_surface_forcing.F90:505
bfb_surface_forcing::bfb_surface_forcing_init
subroutine, public bfb_surface_forcing_init(Time, G, param_file, diag, CS)
Definition: BFB_surface_forcing.F90:252
meso_surface_forcing::meso_surface_forcing_cs
Definition: MESO_surface_forcing.F90:80
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Definition: MOM_diag_mediator.F90:1014
mom_diag_mediator::post_data
Definition: MOM_diag_mediator.F90:80
user_surface_forcing::user_surface_forcing_cs
Definition: user_surface_forcing.F90:80
mom_restart::mom_restart_cs
Definition: MOM_restart.F90:115
mom_string_functions::uppercase
character(len=len(input_string)) function, public uppercase(input_string)
Definition: MOM_string_functions.F90:66
scm_cvmix_tests::scm_cvmix_tests_surface_forcing_init
subroutine, public scm_cvmix_tests_surface_forcing_init(Time, G, param_file, CS)
Initializes surface forcing for the CVmix test case suite.
Definition: SCM_CVmix_tests.F90:124
mom_forcing_type::register_forcing_type_diags
subroutine, public register_forcing_type_diags(Time, diag, use_temperature, handles, use_berg_fluxes)
Register members of the forcing type for diagnostics.
Definition: MOM_forcing_type.F90:1041
mom_domains
Definition: MOM_domains.F90:1
user_revise_forcing::user_revise_forcing_init
subroutine, public user_revise_forcing_init(param_file, CS)
Definition: user_revise_forcing.F90:67
user_revise_forcing
Definition: user_revise_forcing.F90:1
mom_tracer_flow_control::call_tracer_set_forcing
subroutine, public call_tracer_set_forcing(state, fluxes, day_start, day_interval, G, CS)
This subroutine calls the individual tracer modules&#39; subroutines to specify or read quantities relate...
Definition: MOM_tracer_flow_control.F90:374
mom_error_handler::is_root_pe
logical function, public is_root_pe()
Definition: MOM_error_handler.F90:49
mom_diag_mediator
Definition: MOM_diag_mediator.F90:1
mom_error_handler
Definition: MOM_error_handler.F90:1
mom_surface_forcing::forcing_save_restart
subroutine, public forcing_save_restart(CS, G, Time, directory, time_stamped, filename_suffix)
Definition: MOM_surface_forcing.F90:1493
scm_cvmix_tests::scm_cvmix_tests_buoyancy_forcing
subroutine, public scm_cvmix_tests_buoyancy_forcing(state, fluxes, day, G, CS)
Definition: SCM_CVmix_tests.F90:222
mom_restart
Definition: MOM_restart.F90:1
mom_tracer_flow_control::tracer_flow_control_cs
Definition: MOM_tracer_flow_control.F90:93
mom_forcing_type::forcing
Structure that contains pointers to the boundary forcing used to drive the liquid ocean simulated by ...
Definition: MOM_forcing_type.F90:37
user_surface_forcing::user_surface_forcing_init
subroutine, public user_surface_forcing_init(Time, G, param_file, diag, CS)
Definition: user_surface_forcing.F90:320
meso_surface_forcing::meso_surface_forcing_init
subroutine, public meso_surface_forcing_init(Time, G, param_file, diag, CS)
Definition: MESO_surface_forcing.F90:352
mom_surface_forcing::buoyancy_forcing_zero
subroutine buoyancy_forcing_zero(state, fluxes, day, dt, G, CS)
Definition: MOM_surface_forcing.F90:1297
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
Definition: MOM_error_handler.F90:57
meso_surface_forcing::meso_wind_forcing
subroutine, public meso_wind_forcing(state, fluxes, day, G, CS)
Definition: MESO_surface_forcing.F90:117
mom_surface_forcing::surface_forcing_end
subroutine surface_forcing_end(CS, fluxes)
Definition: MOM_surface_forcing.F90:1905
bfb_surface_forcing
Definition: BFB_surface_forcing.F90:1
mom_file_parser::log_version
Definition: MOM_file_parser.F90:140
mom_io::file_exists
Definition: MOM_io.F90:64
user_surface_forcing::user_buoyancy_forcing
subroutine, public user_buoyancy_forcing(state, fluxes, day, dt, G, CS)
Definition: user_surface_forcing.F90:166
mom_restart::save_restart
subroutine, public save_restart(directory, time, G, CS, time_stamped, filename, GV)
Definition: MOM_restart.F90:711
mom_get_input::directories
Definition: MOM_get_input.F90:43
scm_cvmix_tests
Initial conditions and forcing for the single column model (SCM) CVmix test set.
Definition: SCM_CVmix_tests.F90:3
mom_restart::restart_init_end
subroutine, public restart_init_end(CS)
Definition: MOM_restart.F90:1301
mom_surface_forcing::buoyancy_forcing_allocate
subroutine buoyancy_forcing_allocate(fluxes, G, CS)
Definition: MOM_surface_forcing.F90:364
user_surface_forcing::user_wind_forcing
subroutine, public user_wind_forcing(state, fluxes, day, G, CS)
Definition: user_surface_forcing.F90:105
mom_surface_forcing::buoyancy_forcing_linear
subroutine buoyancy_forcing_linear(state, fluxes, day, dt, G, CS)
Definition: MOM_surface_forcing.F90:1402
scm_idealized_hurricane::scm_idealized_hurricane_wind_forcing
subroutine, public scm_idealized_hurricane_wind_forcing(state, fluxes, day, G, CS)
Definition: SCM_idealized_hurricane.F90:159
mom_surface_forcing::buoyancy_forcing_const
subroutine buoyancy_forcing_const(state, fluxes, day, dt, G, CS)
Definition: MOM_surface_forcing.F90:1350
mom_surface_forcing::wind_forcing_gyres
subroutine wind_forcing_gyres(state, fluxes, day, G, CS)
Definition: MOM_surface_forcing.F90:546
mom_string_functions
Definition: MOM_string_functions.F90:1
user_surface_forcing
Definition: user_surface_forcing.F90:1
mom_surface_forcing::wind_forcing_const
subroutine wind_forcing_const(state, fluxes, tau_x0, tau_y0, day, G, CS)
Definition: MOM_surface_forcing.F90:410
mom_surface_forcing::surface_forcing_init
subroutine, public surface_forcing_init(Time, G, param_file, diag, CS, tracer_flow_CSp)
Definition: MOM_surface_forcing.F90:1518
mom_variables
Definition: MOM_variables.F90:1
mom_forcing_type::deallocate_forcing_type
subroutine, public deallocate_forcing_type(fluxes)
Deallocate the forcing type.
Definition: MOM_forcing_type.F90:2425
meso_surface_forcing
Definition: MESO_surface_forcing.F90:1
mom_restart::restart_init
subroutine, public restart_init(param_file, CS, restart_root)
Definition: MOM_restart.F90:1248
mom_forcing_type::forcing_diags
Structure that defines the id handles for the forcing type.
Definition: MOM_forcing_type.F90:171
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
Definition: MOM_error_handler.F90:73
scm_idealized_hurricane
Initial conditions and forcing for the single column model (SCM) idealized hurricane example...
Definition: SCM_idealized_hurricane.F90:3
mom_io::num_timelevels
integer function, public num_timelevels(filename, varname, min_dims)
This function determines how many time levels a variable has.
Definition: MOM_io.F90:463
mom_domains::pass_var
Definition: MOM_domains.F90:65
scm_cvmix_tests::scm_cvmix_tests_cs
Container for surface forcing parameters.
Definition: SCM_CVmix_tests.F90:27
mom_restart::restore_state
subroutine, public restore_state(filename, directory, day, G, CS)
Definition: MOM_restart.F90:889
mom_file_parser::get_param
Definition: MOM_file_parser.F90:135
mom_error_handler::calltree_enter
subroutine, public calltree_enter(mesg, n)
Writes a message about entering a subroutine if call tree reporting is active.
Definition: MOM_error_handler.F90:126