8 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
13 use mom_domains,
only : global_field_sum, bitwise_exact_sum
15 use mom_domains,
only : to_north, to_east, omit_corners
34 use coupler_types_mod,
only : coupler_2d_bc_type, coupler_type_write_chksums
35 use coupler_types_mod,
only : coupler_type_initialized, coupler_type_spawn
36 use coupler_types_mod,
only : coupler_type_copy_data
37 use data_override_mod,
only : data_override_init, data_override
38 use fms_mod,
only : stdout
39 use mpp_mod,
only : mpp_chksum
40 use time_interp_external_mod,
only : init_external_field, time_interp_external
41 use time_interp_external_mod,
only : time_interp_external_init
43 implicit none ;
private
45 #include <MOM_memory.h>
60 integer :: wind_stagger
64 logical :: use_temperature
65 real :: wind_stress_multiplier
68 real :: area_surf = -1.0
69 real :: latent_heat_fusion
70 real :: latent_heat_vapor
78 logical :: use_limited_p_ssh
84 logical :: read_gust_2d
86 real,
pointer,
dimension(:,:) :: &
87 tke_tidal => null(), &
96 logical :: read_tideamp
98 logical :: rigid_sea_ice
103 real :: density_sea_ice
106 real :: rigid_sea_ice_mass
109 logical :: allow_flux_adjustments
110 logical :: liquid_runoff_from_data
113 logical :: salt_restore_as_sflux
114 logical :: adjust_net_srestore_to_zero
115 logical :: adjust_net_srestore_by_scaling
116 logical :: adjust_net_fresh_water_to_zero
117 logical :: use_net_fw_adjustment_sign_bug
118 logical :: adjust_net_fresh_water_by_scaling
119 logical :: mask_srestore_under_ice
121 real :: ice_salt_concentration
122 logical :: mask_srestore_marginal_seas
123 real :: max_delta_srestore
124 real :: max_delta_trestore
125 real,
pointer,
dimension(:,:) :: basin_mask => null()
128 character(len=200) :: inputdir
129 character(len=200) :: salt_restore_file
130 character(len=30) :: salt_restore_var_name
131 logical :: mask_srestore
135 real,
pointer,
dimension(:,:) :: srestore_mask => null()
136 character(len=200) :: temp_restore_file
137 character(len=30) :: temp_restore_var_name
138 logical :: mask_trestore
142 real,
pointer,
dimension(:,:) :: trestore_mask => null()
143 integer :: id_srestore = -1
144 integer :: id_trestore = -1
156 real,
pointer,
dimension(:,:) :: lrunoff =>null()
157 real,
pointer,
dimension(:,:) :: frunoff =>null()
158 real,
pointer,
dimension(:,:) :: u_flux =>null()
159 real,
pointer,
dimension(:,:) :: v_flux =>null()
160 real,
pointer,
dimension(:,:) :: t_flux =>null()
161 real,
pointer,
dimension(:,:) :: q_flux =>null()
162 real,
pointer,
dimension(:,:) :: salt_flux =>null()
163 real,
pointer,
dimension(:,:) :: seaice_melt_heat =>null()
164 real,
pointer,
dimension(:,:) :: seaice_melt =>null()
165 real,
pointer,
dimension(:,:) :: lw_flux =>null()
166 real,
pointer,
dimension(:,:) :: sw_flux_vis_dir =>null()
167 real,
pointer,
dimension(:,:) :: sw_flux_vis_dif =>null()
168 real,
pointer,
dimension(:,:) :: sw_flux_nir_dir =>null()
169 real,
pointer,
dimension(:,:) :: sw_flux_nir_dif =>null()
170 real,
pointer,
dimension(:,:) :: lprec =>null()
171 real,
pointer,
dimension(:,:) :: fprec =>null()
172 real,
pointer,
dimension(:,:) :: ustar_berg =>null()
173 real,
pointer,
dimension(:,:) :: area_berg =>null()
174 real,
pointer,
dimension(:,:) :: mass_berg =>null()
175 real,
pointer,
dimension(:,:) :: lrunoff_hflx =>null()
176 real,
pointer,
dimension(:,:) :: frunoff_hflx =>null()
177 real,
pointer,
dimension(:,:) :: p =>null()
179 real,
pointer,
dimension(:,:) :: mi =>null()
180 real,
pointer,
dimension(:,:) :: ice_rigidity =>null()
185 type(coupler_2d_bc_type) :: fluxes
187 integer :: wind_stagger = -999
202 sfc_state, restore_salt, restore_temp)
204 target,
intent(in) :: iob
206 type(
forcing),
intent(inout) :: fluxes
209 integer,
dimension(4),
intent(in) :: index_bounds
210 type(time_type),
intent(in) :: time
212 real,
intent(in) :: valid_time
218 type(
surface),
intent(in) :: sfc_state
220 logical,
optional,
intent(in) :: restore_salt
221 logical,
optional,
intent(in) :: restore_temp
224 real,
dimension(SZI_(G),SZJ_(G)) :: &
225 data_restore, & !< The surface value toward which to restore [g/kg or degC]
226 sst_anom, & !< Instantaneous sea surface temperature anomalies from a target value [deg C]
227 sss_anom, & !< Instantaneous sea surface salinity anomalies from a target value [g/kg]
228 sss_mean, & !< A (mean?) salinity about which to normalize local salinity
238 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
239 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
240 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
242 logical :: restore_salinity
244 logical :: restore_sst
248 real :: kg_m2_s_conversion
252 real :: sign_for_net_fw_bug
256 isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
257 jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
258 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
259 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
260 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
261 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
262 isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
264 kg_m2_s_conversion = us%kg_m3_to_R*us%m_to_Z*us%T_to_s
266 open_ocn_mask(:,:) = 1.0
268 fluxes%vPrecGlobalAdj = 0.0
269 fluxes%vPrecGlobalScl = 0.0
270 fluxes%saltFluxGlobalAdj = 0.0
271 fluxes%saltFluxGlobalScl = 0.0
272 fluxes%netFWGlobalAdj = 0.0
273 fluxes%netFWGlobalScl = 0.0
275 restore_salinity = .false.
276 if (
present(restore_salt)) restore_salinity = restore_salt
277 restore_sst = .false.
278 if (
present(restore_temp)) restore_sst = restore_temp
282 if (fluxes%dt_buoy_accum < 0)
then
283 call allocate_forcing_type(g, fluxes, water=.true., heat=.true., &
284 ustar=.true., press=.true.)
286 call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
287 call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
288 call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed)
289 call safe_alloc_ptr(fluxes%sw_nir_dif,isd,ied,jsd,jed)
291 call safe_alloc_ptr(fluxes%p_surf ,isd,ied,jsd,jed)
292 call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed)
293 if (cs%use_limited_P_SSH)
then
294 fluxes%p_surf_SSH => fluxes%p_surf
296 fluxes%p_surf_SSH => fluxes%p_surf_full
299 call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed)
300 call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed)
301 call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)
303 call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed)
304 call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed)
306 if (cs%allow_flux_adjustments)
then
307 call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
308 call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)
311 do j=js-2,je+2 ;
do i=is-2,ie+2
312 fluxes%TKE_tidal(i,j) = cs%TKE_tidal(i,j)
313 fluxes%ustar_tidal(i,j) = cs%ustar_tidal(i,j)
316 if (restore_temp)
call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
320 if (((
associated(iob%ustar_berg) .and. (.not.
associated(fluxes%ustar_berg))) &
321 .or. (
associated(iob%area_berg) .and. (.not.
associated(fluxes%area_berg)))) &
322 .or. (
associated(iob%mass_berg) .and. (.not.
associated(fluxes%mass_berg)))) &
323 call allocate_forcing_type(g, fluxes, iceberg=.true.)
325 if ((.not.coupler_type_initialized(fluxes%tr_fluxes)) .and. &
326 coupler_type_initialized(iob%fluxes)) &
327 call coupler_type_spawn(iob%fluxes, fluxes%tr_fluxes, &
328 (/is,is,ie,ie/), (/js,js,je,je/))
334 if (cs%area_surf < 0.0)
then
335 do j=js,je ;
do i=is,ie
336 work_sum(i,j) = us%L_to_m**2*g%areaT(i,j) * g%mask2dT(i,j)
343 fluxes%fluxes_used = .false.
344 fluxes%dt_buoy_accum = us%s_to_T*valid_time
346 if (cs%allow_flux_adjustments)
then
347 fluxes%heat_added(:,:)=0.0
348 fluxes%salt_flux_added(:,:)=0.0
351 do j=js,je ;
do i=is,ie
352 fluxes%salt_flux(i,j) = 0.0
353 fluxes%vprec(i,j) = 0.0
357 if (restore_salinity)
then
358 call time_interp_external(cs%id_srestore,time,data_restore)
360 open_ocn_mask(:,:) = 1.0
361 if (cs%mask_srestore_under_ice)
then
362 do j=js,je ;
do i=is,ie
363 if (sfc_state%SST(i,j) <= -0.0539*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
366 if (cs%salt_restore_as_sflux)
then
367 do j=js,je ;
do i=is,ie
368 delta_sss = data_restore(i,j)- sfc_state%SSS(i,j)
369 delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),cs%max_delta_srestore)
370 fluxes%salt_flux(i,j) = 1.e-3*g%mask2dT(i,j) * (cs%Rho0*us%m_to_Z*us%T_to_s*cs%Flux_const)* &
371 (cs%basin_mask(i,j)*open_ocn_mask(i,j)*cs%srestore_mask(i,j)) *delta_sss
373 if (cs%adjust_net_srestore_to_zero)
then
374 if (cs%adjust_net_srestore_by_scaling)
then
376 unit_scale=us%R_to_kg_m3*us%Z_to_m*us%s_to_T)
377 fluxes%saltFluxGlobalAdj = 0.
379 work_sum(is:ie,js:je) = us%L_to_m**2*us%R_to_kg_m3*us%Z_to_m*us%s_to_T * &
380 g%areaT(is:ie,js:je)*fluxes%salt_flux(is:ie,js:je)
381 fluxes%saltFluxGlobalAdj =
reproducing_sum(work_sum(:,:), isr,ier, jsr,jer)/cs%area_surf
382 fluxes%salt_flux(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) - kg_m2_s_conversion * fluxes%saltFluxGlobalAdj
385 fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je)
387 do j=js,je ;
do i=is,ie
388 if (g%mask2dT(i,j) > 0.5)
then
389 delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
390 delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),cs%max_delta_srestore)
391 fluxes%vprec(i,j) = (cs%basin_mask(i,j)*open_ocn_mask(i,j)*cs%srestore_mask(i,j))* &
392 (us%m_to_Z*us%T_to_s * cs%Rho0*cs%Flux_const) * &
393 delta_sss / (0.5*(sfc_state%SSS(i,j) + data_restore(i,j)))
396 if (cs%adjust_net_srestore_to_zero)
then
397 if (cs%adjust_net_srestore_by_scaling)
then
399 unit_scale=us%R_to_kg_m3*us%Z_to_m*us%s_to_T)
400 fluxes%vPrecGlobalAdj = 0.
402 work_sum(is:ie,js:je) = us%L_to_m**2*g%areaT(is:ie,js:je) * &
403 us%R_to_kg_m3*us%Z_to_m*us%s_to_T*fluxes%vprec(is:ie,js:je)
404 fluxes%vPrecGlobalAdj =
reproducing_sum(work_sum(:,:), isr, ier, jsr, jer) / cs%area_surf
405 do j=js,je ;
do i=is,ie
406 fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion * fluxes%vPrecGlobalAdj ) * g%mask2dT(i,j)
414 if (restore_sst)
then
415 call time_interp_external(cs%id_trestore,time,data_restore)
416 do j=js,je ;
do i=is,ie
417 delta_sst = data_restore(i,j)- sfc_state%SST(i,j)
418 delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),cs%max_delta_trestore)
419 fluxes%heat_added(i,j) = g%mask2dT(i,j) * cs%trestore_mask(i,j) * &
420 (us%R_to_kg_m3*cs%Rho0*fluxes%C_p) * delta_sst * cs%Flux_const
425 if (cs%liquid_runoff_from_data .and. .not.
associated(iob%lrunoff))
then
426 call mom_error(fatal,
"liquid runoff is being added via data_override but "// &
427 "there is no associated runoff in the IOB")
432 i0 = is - isc_bnd ; j0 = js - jsc_bnd
433 do j=js,je ;
do i=is,ie
435 if (
associated(iob%lprec)) &
436 fluxes%lprec(i,j) = kg_m2_s_conversion * iob%lprec(i-i0,j-j0) * g%mask2dT(i,j)
438 if (
associated(iob%fprec)) &
439 fluxes%fprec(i,j) = kg_m2_s_conversion * iob%fprec(i-i0,j-j0) * g%mask2dT(i,j)
441 if (
associated(iob%q_flux)) &
442 fluxes%evap(i,j) = kg_m2_s_conversion * iob%q_flux(i-i0,j-j0) * g%mask2dT(i,j)
445 if (
associated(iob%lrunoff))
then
446 if(cs%liquid_runoff_from_data)
call data_override(
'OCN',
'runoff', iob%lrunoff, time)
447 fluxes%lrunoff(i,j) = kg_m2_s_conversion * iob%lrunoff(i-i0,j-j0) * g%mask2dT(i,j)
451 if (
associated(iob%frunoff))
then
452 fluxes%frunoff(i,j) = kg_m2_s_conversion * iob%frunoff(i-i0,j-j0) * g%mask2dT(i,j)
455 if (
associated(iob%ustar_berg)) &
456 fluxes%ustar_berg(i,j) = us%m_to_Z*us%T_to_s * iob%ustar_berg(i-i0,j-j0) * g%mask2dT(i,j)
458 if (
associated(iob%area_berg)) &
459 fluxes%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
461 if (
associated(iob%mass_berg)) &
462 fluxes%mass_berg(i,j) = iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
464 if (
associated(iob%lrunoff_hflx)) &
465 fluxes%heat_content_lrunoff(i,j) = iob%lrunoff_hflx(i-i0,j-j0) * g%mask2dT(i,j)
467 if (
associated(iob%frunoff_hflx)) &
468 fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion * iob%frunoff_hflx(i-i0,j-j0) * g%mask2dT(i,j)
470 if (
associated(iob%lw_flux)) &
471 fluxes%LW(i,j) = iob%lw_flux(i-i0,j-j0) * g%mask2dT(i,j)
473 if (
associated(iob%t_flux)) &
474 fluxes%sens(i,j) = iob%t_flux(i-i0,j-j0) * g%mask2dT(i,j)
477 if (
associated(iob%seaice_melt_heat)) &
478 fluxes%seaice_melt_heat(i,j) = g%mask2dT(i,j) * iob%seaice_melt_heat(i-i0,j-j0)
481 if (
associated(iob%seaice_melt)) &
482 fluxes%seaice_melt(i,j) = kg_m2_s_conversion * g%mask2dT(i,j) * iob%seaice_melt(i-i0,j-j0)
484 fluxes%latent(i,j) = 0.0
485 if (
associated(iob%fprec))
then
486 fluxes%latent(i,j) = fluxes%latent(i,j) + iob%fprec(i-i0,j-j0)*cs%latent_heat_fusion
487 fluxes%latent_fprec_diag(i,j) = g%mask2dT(i,j) * iob%fprec(i-i0,j-j0)*cs%latent_heat_fusion
489 if (
associated(iob%frunoff))
then
490 fluxes%latent(i,j) = fluxes%latent(i,j) + iob%frunoff(i-i0,j-j0)*cs%latent_heat_fusion
491 fluxes%latent_frunoff_diag(i,j) = g%mask2dT(i,j) * iob%frunoff(i-i0,j-j0)*cs%latent_heat_fusion
493 if (
associated(iob%q_flux))
then
494 fluxes%latent(i,j) = fluxes%latent(i,j) + iob%q_flux(i-i0,j-j0)*cs%latent_heat_vapor
495 fluxes%latent_evap_diag(i,j) = g%mask2dT(i,j) * iob%q_flux(i-i0,j-j0)*cs%latent_heat_vapor
497 fluxes%latent(i,j) = g%mask2dT(i,j) * fluxes%latent(i,j)
499 if (
associated(iob%sw_flux_vis_dir)) &
500 fluxes%sw_vis_dir(i,j) = g%mask2dT(i,j) * iob%sw_flux_vis_dir(i-i0,j-j0)
502 if (
associated(iob%sw_flux_vis_dif)) &
503 fluxes%sw_vis_dif(i,j) = g%mask2dT(i,j) * iob%sw_flux_vis_dif(i-i0,j-j0)
505 if (
associated(iob%sw_flux_nir_dir)) &
506 fluxes%sw_nir_dir(i,j) = g%mask2dT(i,j) * iob%sw_flux_nir_dir(i-i0,j-j0)
508 if (
associated(iob%sw_flux_nir_dif)) &
509 fluxes%sw_nir_dif(i,j) = g%mask2dT(i,j) * iob%sw_flux_nir_dif(i-i0,j-j0)
511 fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + &
512 fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j)
517 if (
associated(iob%p))
then
518 if (cs%max_p_surf >= 0.0)
then
519 do j=js,je ;
do i=is,ie
520 fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
521 fluxes%p_surf(i,j) = min(fluxes%p_surf_full(i,j),cs%max_p_surf)
524 do j=js,je ;
do i=is,ie
525 fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
526 fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
529 fluxes%accumulate_p_surf = .true.
532 if (
associated(iob%salt_flux))
then
533 do j=js,je ;
do i=is,ie
534 fluxes%salt_flux(i,j) = g%mask2dT(i,j)*(fluxes%salt_flux(i,j) + kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0))
535 fluxes%salt_flux_in(i,j) = g%mask2dT(i,j)*( kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0) )
540 if (cs%adjust_net_fresh_water_to_zero)
then
541 sign_for_net_fw_bug = 1.
542 if (cs%use_net_FW_adjustment_sign_bug) sign_for_net_fw_bug = -1.
543 do j=js,je ;
do i=is,ie
544 net_fw(i,j) = us%R_to_kg_m3*us%Z_to_m*us%s_to_T * &
545 (((fluxes%lprec(i,j) + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + &
546 (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
547 (fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * us%L_to_m**2*g%areaT(i,j)
548 net_fw2(i,j) = net_fw(i,j) / (us%L_to_m**2*g%areaT(i,j))
551 if (cs%adjust_net_fresh_water_by_scaling)
then
553 do j=js,je ;
do i=is,ie
554 fluxes%vprec(i,j) = fluxes%vprec(i,j) + us%kg_m3_to_R*us%m_to_Z*us%T_to_s * &
555 (net_fw2(i,j) - net_fw(i,j)/(us%L_to_m**2*g%areaT(i,j))) * g%mask2dT(i,j)
558 fluxes%netFWGlobalAdj =
reproducing_sum(net_fw(:,:), isr, ier, jsr, jer) / cs%area_surf
559 do j=js,je ;
do i=is,ie
560 fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion * fluxes%netFWGlobalAdj ) * g%mask2dT(i,j)
566 if (coupler_type_initialized(fluxes%tr_fluxes) .and. &
567 coupler_type_initialized(iob%fluxes)) &
568 call coupler_type_copy_data(iob%fluxes, fluxes%tr_fluxes)
570 if (cs%allow_flux_adjustments)
then
576 call user_alter_forcing(sfc_state, fluxes, time, g, cs%urf_CS)
587 target,
intent(in) :: iob
590 integer,
dimension(4),
intent(in) :: index_bounds
591 type(time_type),
intent(in) :: time
599 real,
dimension(SZIB_(G),SZJB_(G)) :: &
600 taux_at_q, & !< Zonal wind stresses at q points [Pa]
603 real,
dimension(SZI_(G),SZJ_(G)) :: &
604 rigidity_at_h, & !< Ice rigidity at tracer points (m3 s-1)
605 taux_at_h, & !< Zonal wind stresses at h points [Pa]
612 real :: pa_conversion
613 real :: stress_conversion
619 integer :: wind_stagger
620 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
621 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
622 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
626 isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
627 jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
628 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
629 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
630 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
631 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
632 isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
633 i0 = is - isc_bnd ; j0 = js - jsc_bnd
635 irho0 = us%L_to_Z / cs%Rho0
636 pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
637 stress_conversion = pa_conversion * cs%wind_stress_multiplier
641 if (.not.forces%initialized)
then
645 call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed)
646 call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed)
648 if (cs%rigid_sea_ice)
then
649 call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
650 call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
653 forces%initialized = .true.
656 if ( (
associated(iob%area_berg) .and. (.not.
associated(forces%area_berg))) .or. &
657 (
associated(iob%mass_berg) .and. (.not.
associated(forces%mass_berg))) ) &
659 if (
associated(iob%ice_rigidity))
then
660 rigidity_at_h(:,:) = 0.0
661 call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
662 call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
665 forces%accumulate_rigidity = .true.
666 if (
associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0
667 if (
associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0
670 if (cs%use_limited_P_SSH)
then
671 forces%p_surf_SSH => forces%p_surf
673 forces%p_surf_SSH => forces%p_surf_full
675 if (
associated(iob%p))
then
676 if (cs%max_p_surf >= 0.0)
then
677 do j=js,je ;
do i=is,ie
678 forces%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
679 forces%p_surf(i,j) = min(forces%p_surf_full(i,j),cs%max_p_surf)
682 do j=js,je ;
do i=is,ie
683 forces%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
684 forces%p_surf(i,j) = forces%p_surf_full(i,j)
688 do j=js,je ;
do i=is,ie
689 forces%p_surf_full(i,j) = 0.0
690 forces%p_surf(i,j) = 0.0
693 forces%accumulate_p_surf = .true.
699 wind_stagger = cs%wind_stagger
702 if (wind_stagger == bgrid_ne)
then
704 taux_at_q(:,:) = 0.0 ; tauy_at_q(:,:) = 0.0
706 if (wind_stagger == agrid)
then
708 taux_at_h(:,:) = 0.0 ; tauy_at_h(:,:) = 0.0
712 do j=js,je ;
do i=is,ie
713 if (
associated(iob%area_berg)) &
714 forces%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
716 if (
associated(iob%mass_berg)) &
717 forces%mass_berg(i,j) = iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
719 if (
associated(iob%ice_rigidity)) &
720 rigidity_at_h(i,j) = iob%ice_rigidity(i-i0,j-j0) * g%mask2dT(i,j)
722 if (wind_stagger == bgrid_ne)
then
723 if (
associated(iob%u_flux)) taux_at_q(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
724 if (
associated(iob%v_flux)) tauy_at_q(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
725 elseif (wind_stagger == agrid)
then
726 if (
associated(iob%u_flux)) taux_at_h(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
727 if (
associated(iob%v_flux)) tauy_at_h(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
729 if (
associated(iob%u_flux)) forces%taux(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
730 if (
associated(iob%v_flux)) forces%tauy(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
736 if (wind_stagger == bgrid_ne)
then
739 call pass_vector(taux_at_q, tauy_at_q, g%Domain, stagger=bgrid_ne, halo=1)
741 do j=js,je ;
do i=isq,ieq
742 forces%taux(i,j) = 0.0
743 if ((g%mask2dBu(i,j) + g%mask2dBu(i,j-1)) > 0) &
744 forces%taux(i,j) = (g%mask2dBu(i,j)*taux_at_q(i,j) + &
745 g%mask2dBu(i,j-1)*taux_at_q(i,j-1)) / &
746 (g%mask2dBu(i,j) + g%mask2dBu(i,j-1))
749 do j=jsq,jeq ;
do i=is,ie
750 forces%tauy(i,j) = 0.0
751 if ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j)) > 0) &
752 forces%tauy(i,j) = (g%mask2dBu(i,j)*tauy_at_q(i,j) + &
753 g%mask2dBu(i-1,j)*tauy_at_q(i-1,j)) / &
754 (g%mask2dBu(i,j) + g%mask2dBu(i-1,j))
761 do j=js,je ;
do i=is,ie
762 tau_mag = 0.0 ; gustiness = cs%gust_const
763 if (((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + &
764 (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) > 0)
then
765 tau_mag = sqrt(((g%mask2dBu(i,j)*(taux_at_q(i,j)**2 + tauy_at_q(i,j)**2) + &
766 g%mask2dBu(i-1,j-1)*(taux_at_q(i-1,j-1)**2 + tauy_at_q(i-1,j-1)**2)) + &
767 (g%mask2dBu(i,j-1)*(taux_at_q(i,j-1)**2 + tauy_at_q(i,j-1)**2) + &
768 g%mask2dBu(i-1,j)*(taux_at_q(i-1,j)**2 + tauy_at_q(i-1,j)**2)) ) / &
769 ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) )
770 if (cs%read_gust_2d) gustiness = cs%gust(i,j)
772 forces%ustar(i,j) = sqrt(gustiness*irho0 + irho0*tau_mag)
775 elseif (wind_stagger == agrid)
then
776 call pass_vector(taux_at_h, tauy_at_h, g%Domain, to_all+omit_corners, stagger=agrid, halo=1)
778 do j=js,je ;
do i=isq,ieq
779 forces%taux(i,j) = 0.0
780 if ((g%mask2dT(i,j) + g%mask2dT(i+1,j)) > 0) &
781 forces%taux(i,j) = (g%mask2dT(i,j)*taux_at_h(i,j) + &
782 g%mask2dT(i+1,j)*taux_at_h(i+1,j)) / &
783 (g%mask2dT(i,j) + g%mask2dT(i+1,j))
786 do j=jsq,jeq ;
do i=is,ie
787 forces%tauy(i,j) = 0.0
788 if ((g%mask2dT(i,j) + g%mask2dT(i,j+1)) > 0) &
789 forces%tauy(i,j) = (g%mask2dT(i,j)*tauy_at_h(i,j) + &
790 g%mask2dT(i,j+1)*tauy_at_h(i,j+1)) / &
791 (g%mask2dT(i,j) + g%mask2dT(i,j+1))
794 do j=js,je ;
do i=is,ie
795 gustiness = cs%gust_const
796 if (cs%read_gust_2d .and. (g%mask2dT(i,j) > 0)) gustiness = cs%gust(i,j)
797 forces%ustar(i,j) = sqrt(gustiness*irho0 + irho0 * g%mask2dT(i,j) * &
798 sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
804 call pass_vector(forces%taux, forces%tauy, g%Domain, halo=1)
806 do j=js,je ;
do i=is,ie
808 if ((g%mask2dCu(i-1,j) + g%mask2dCu(i,j)) > 0) &
809 taux2 = (g%mask2dCu(i-1,j)*forces%taux(i-1,j)**2 + &
810 g%mask2dCu(i,j)*forces%taux(i,j)**2) / (g%mask2dCu(i-1,j) + g%mask2dCu(i,j))
813 if ((g%mask2dCv(i,j-1) + g%mask2dCv(i,j)) > 0) &
814 tauy2 = (g%mask2dCv(i,j-1)*forces%tauy(i,j-1)**2 + &
815 g%mask2dCv(i,j)*forces%tauy(i,j)**2) / (g%mask2dCv(i,j-1) + g%mask2dCv(i,j))
817 if (cs%read_gust_2d)
then
818 forces%ustar(i,j) = sqrt(cs%gust(i,j)*irho0 + irho0*sqrt(taux2 + tauy2))
820 forces%ustar(i,j) = sqrt(cs%gust_const*irho0 + irho0*sqrt(taux2 + tauy2))
827 if (
associated(iob%ice_rigidity))
then
828 call pass_var(rigidity_at_h, g%Domain, halo=1)
829 do i=is-1,ie ;
do j=js,je
830 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
831 min(rigidity_at_h(i,j), rigidity_at_h(i+1,j))
833 do i=is,ie ;
do j=js-1,je
834 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
835 min(rigidity_at_h(i,j), rigidity_at_h(i,j+1))
839 if (cs%rigid_sea_ice)
then
840 call pass_var(forces%p_surf_full, g%Domain, halo=1)
841 i_gearth = 1.0 / cs%g_Earth
842 kv_rho_ice = (cs%kv_sea_ice / cs%density_sea_ice)
843 do i=is-1,ie ;
do j=js,je
844 mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * i_gearth
846 if (mass_ice > cs%rigid_sea_ice_mass)
then
847 mass_eff = (mass_ice - cs%rigid_sea_ice_mass) **2 / &
848 (mass_ice + cs%rigid_sea_ice_mass)
850 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * mass_eff
852 do i=is,ie ;
do j=js-1,je
853 mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i,j+1)) * i_gearth
855 if (mass_ice > cs%rigid_sea_ice_mass)
then
856 mass_eff = (mass_ice - cs%rigid_sea_ice_mass) **2 / &
857 (mass_ice + cs%rigid_sea_ice_mass)
859 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * mass_eff
863 if (cs%allow_flux_adjustments)
then
884 type(time_type),
intent(in) :: time
885 type(
forcing),
intent(inout) :: fluxes
888 real,
dimension(SZI_(G),SZJ_(G)) :: temp_at_h
890 integer :: isc, iec, jsc, jec, i, j
891 logical :: overrode_h
893 isc = g%isc; iec = g%iec ; jsc = g%jsc; jec = g%jec
896 call data_override(
'OCN',
'hflx_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
898 if (overrode_h)
then ;
do j=jsc,jec ;
do i=isc,iec
899 fluxes%heat_added(i,j) = fluxes%heat_added(i,j) + temp_at_h(i,j)* g%mask2dT(i,j)
900 enddo ;
enddo ;
endif
904 call data_override(
'OCN',
'sflx_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
906 if (overrode_h)
then ;
do j=jsc,jec ;
do i=isc,iec
907 fluxes%salt_flux_added(i,j) = fluxes%salt_flux_added(i,j) + &
908 us%kg_m3_to_R*us%m_to_Z*us%T_to_s * temp_at_h(i,j)* g%mask2dT(i,j)
909 enddo ;
enddo ;
endif
913 call data_override(
'OCN',
'prcme_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
915 if (overrode_h)
then ;
do j=jsc,jec ;
do i=isc,iec
916 fluxes%vprec(i,j) = fluxes%vprec(i,j) + us%kg_m3_to_R*us%m_to_Z*us%T_to_s * temp_at_h(i,j)* g%mask2dT(i,j)
917 enddo ;
enddo ;
endif
930 type(time_type),
intent(in) :: time
934 real,
dimension(SZI_(G),SZJ_(G)) :: tempx_at_h
935 real,
dimension(SZI_(G),SZJ_(G)) :: tempy_at_h
937 integer :: isc, iec, jsc, jec, i, j
938 real :: dlondx, dlondy, rdlon, cosa, sina, zonal_tau, merid_tau
939 real :: pa_conversion
940 logical :: overrode_x, overrode_y
942 isc = g%isc; iec = g%iec ; jsc = g%jsc; jec = g%jec
943 pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
945 tempx_at_h(:,:) = 0.0 ; tempy_at_h(:,:) = 0.0
947 overrode_x = .false. ; overrode_y = .false.
948 call data_override(
'OCN',
'taux_adj', tempx_at_h(isc:iec,jsc:jec), time, override=overrode_x)
949 call data_override(
'OCN',
'tauy_adj', tempy_at_h(isc:iec,jsc:jec), time, override=overrode_y)
951 if (overrode_x .or. overrode_y)
then
952 if (.not. (overrode_x .and. overrode_y))
call mom_error(fatal,
"apply_flux_adjustments: "//&
953 "Both taux_adj and tauy_adj must be specified, or neither, in data_table")
956 call pass_vector(tempx_at_h, tempy_at_h, g%Domain, to_all, agrid, halo=1)
957 do j=jsc-1,jec+1 ;
do i=isc-1,iec+1
958 dlondx = g%geoLonCu(i,j) - g%geoLonCu(i-1,j)
959 dlondy = g%geoLonCv(i,j) - g%geoLonCv(i,j-1)
960 rdlon = sqrt( dlondx * dlondx + dlondy * dlondy )
961 if (rdlon > 0.) rdlon = 1. / rdlon
962 cosa = dlondx * rdlon
963 sina = dlondy * rdlon
964 zonal_tau = pa_conversion * tempx_at_h(i,j)
965 merid_tau = pa_conversion * tempy_at_h(i,j)
966 tempx_at_h(i,j) = cosa * zonal_tau - sina * merid_tau
967 tempy_at_h(i,j) = sina * zonal_tau + cosa * merid_tau
971 do j=jsc,jec ;
do i=isc-1,iec
972 forces%taux(i,j) = forces%taux(i,j) + 0.5 * ( tempx_at_h(i,j) + tempx_at_h(i+1,j) )
975 do j=jsc-1,jec ;
do i=isc,iec
976 forces%tauy(i,j) = forces%tauy(i,j) + 0.5 * ( tempy_at_h(i,j) + tempy_at_h(i,j+1) )
988 type(time_type),
intent(in) :: time
989 character(len=*),
intent(in) :: directory
991 logical,
optional,
intent(in) :: time_stamped
993 character(len=*),
optional,
intent(in) :: filename_suffix
996 if (.not.
associated(cs))
return
997 if (.not.
associated(cs%restart_CSp))
return
998 call save_restart(directory, time, g, cs%restart_CSp, time_stamped)
1004 type(time_type),
intent(in) :: time
1008 type(
diag_ctrl),
target,
intent(inout) :: diag
1012 logical,
optional,
intent(in) :: restore_salt
1014 logical,
optional,
intent(in) :: restore_temp
1020 logical :: new_sim, iceberg_flux_diags
1021 type(time_type) :: time_frc
1022 character(len=200) :: tideamp_file, gust_file, salt_file, temp_file
1024 #include "version_variable.h"
1025 character(len=40) :: mdl =
"MOM_surface_forcing_nuopc"
1026 character(len=48) :: stagger
1027 character(len=48) :: flnam
1028 character(len=240) :: basin_file
1029 integer :: i, j, isd, ied, jsd, jed
1031 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1033 if (
associated(cs))
then
1034 call mom_error(warning,
"surface_forcing_init called with an associated "// &
1035 "control structure.")
1040 id_clock_forcing=cpu_clock_id(
'Ocean surface forcing', grain=clock_subcomponent)
1045 call write_version_number(version)
1049 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, &
1050 "The directory in which all input files are found.", &
1052 cs%inputdir = slasher(cs%inputdir)
1053 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
1054 "If true, Temperature and salinity are used as state "//&
1055 "variables.", default=.true.)
1056 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
1057 "The mean ocean density used with BOUSSINESQ true to "//&
1058 "calculate accelerations and the mass for conservation "//&
1059 "properties, or with BOUSSINSEQ false to convert some "//&
1060 "parameters from vertical units of m to kg m-2.", &
1061 units=
"kg m-3", default=1035.0, scale=us%kg_m3_to_R)
1062 call get_param(param_file, mdl,
"LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
1063 "The latent heat of fusion.", units=
"J/kg", default=hlf)
1064 call get_param(param_file, mdl,
"LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
1065 "The latent heat of fusion.", units=
"J/kg", default=hlv)
1066 call get_param(param_file, mdl,
"MAX_P_SURF", cs%max_p_surf, &
1067 "The maximum surface pressure that can be exerted by the "//&
1068 "atmosphere and floating sea-ice or ice shelves. This is "//&
1069 "needed because the FMS coupling structure does not "//&
1070 "limit the water that can be frozen out of the ocean and "//&
1071 "the ice-ocean heat fluxes are treated explicitly. No "//&
1072 "limit is applied if a negative value is used.", units=
"Pa", &
1074 call get_param(param_file, mdl,
"ADJUST_NET_SRESTORE_TO_ZERO", &
1075 cs%adjust_net_srestore_to_zero, &
1076 "If true, adjusts the salinity restoring seen to zero "//&
1077 "whether restoring is via a salt flux or virtual precip.",&
1078 default=restore_salt)
1079 call get_param(param_file, mdl,
"ADJUST_NET_SRESTORE_BY_SCALING", &
1080 cs%adjust_net_srestore_by_scaling, &
1081 "If true, adjustments to salt restoring to achieve zero net are "//&
1082 "made by scaling values without moving the zero contour.",&
1084 call get_param(param_file, mdl,
"ADJUST_NET_FRESH_WATER_TO_ZERO", &
1085 cs%adjust_net_fresh_water_to_zero, &
1086 "If true, adjusts the net fresh-water forcing seen "//&
1087 "by the ocean (including restoring) to zero.", default=.false.)
1088 if (cs%adjust_net_fresh_water_to_zero) &
1089 call get_param(param_file, mdl,
"USE_NET_FW_ADJUSTMENT_SIGN_BUG", &
1090 cs%use_net_FW_adjustment_sign_bug, &
1091 "If true, use the wrong sign for the adjustment to "//&
1092 "the net fresh-water.", default=.false.)
1093 call get_param(param_file, mdl,
"ADJUST_NET_FRESH_WATER_BY_SCALING", &
1094 cs%adjust_net_fresh_water_by_scaling, &
1095 "If true, adjustments to net fresh water to achieve zero net are "//&
1096 "made by scaling values without moving the zero contour.",&
1098 call get_param(param_file, mdl,
"ICE_SALT_CONCENTRATION", &
1099 cs%ice_salt_concentration, &
1100 "The assumed sea-ice salinity needed to reverse engineer the "//&
1101 "melt flux (or ice-ocean fresh-water flux).", &
1102 units=
"kg/kg", default=0.005)
1103 call get_param(param_file, mdl,
"USE_LIMITED_PATM_SSH", cs%use_limited_P_SSH, &
1104 "If true, return the sea surface height with the "//&
1105 "correction for the atmospheric (and sea-ice) pressure "//&
1106 "limited by max_p_surf instead of the full atmospheric "//&
1107 "pressure.", default=.true.)
1109 call get_param(param_file, mdl,
"WIND_STAGGER", stagger, &
1110 "A case-insensitive character string to indicate the "//&
1111 "staggering of the input wind stress field. Valid "//&
1112 "values are 'A', 'B', or 'C'.", default=
"C")
1113 if (
uppercase(stagger(1:1)) ==
'A')
then ; cs%wind_stagger = agrid
1114 elseif (
uppercase(stagger(1:1)) ==
'B')
then ; cs%wind_stagger = bgrid_ne
1115 elseif (
uppercase(stagger(1:1)) ==
'C')
then ; cs%wind_stagger = cgrid_ne
1116 else ;
call mom_error(fatal,
"surface_forcing_init: WIND_STAGGER = "// &
1117 trim(stagger)//
" is invalid.") ;
endif
1118 call get_param(param_file, mdl,
"WIND_STRESS_MULTIPLIER", cs%wind_stress_multiplier, &
1119 "A factor multiplying the wind-stress given to the ocean by the "//&
1120 "coupler. This is used for testing and should be =1.0 for any "//&
1121 "production runs.", default=1.0)
1123 if (restore_salt)
then
1124 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
1125 "The constant that relates the restoring surface fluxes "//&
1126 "to the relative surface anomalies (akin to a piston "//&
1127 "velocity). Note the non-MKS units.", units=
"m day-1", &
1128 fail_if_missing=.true.)
1129 call get_param(param_file, mdl,
"SALT_RESTORE_FILE", cs%salt_restore_file, &
1130 "A file in which to find the surface salinity to use for restoring.", &
1131 default=
"salt_restore.nc")
1132 call get_param(param_file, mdl,
"SALT_RESTORE_VARIABLE", cs%salt_restore_var_name, &
1133 "The name of the surface salinity variable to read from "//&
1134 "SALT_RESTORE_FILE for restoring salinity.", &
1137 cs%Flux_const = cs%Flux_const / 86400.0
1139 call get_param(param_file, mdl,
"SRESTORE_AS_SFLUX", cs%salt_restore_as_sflux, &
1140 "If true, the restoring of salinity is applied as a salt "//&
1141 "flux instead of as a freshwater flux.", default=.false.)
1142 call get_param(param_file, mdl,
"MAX_DELTA_SRESTORE", cs%max_delta_srestore, &
1143 "The maximum salinity difference used in restoring terms.", &
1144 units=
"PSU or g kg-1", default=999.0)
1145 call get_param(param_file, mdl,
"MASK_SRESTORE_UNDER_ICE", &
1146 cs%mask_srestore_under_ice, &
1147 "If true, disables SSS restoring under sea-ice based on a frazil "//&
1148 "criteria (SST<=Tf). Only used when RESTORE_SALINITY is True.", &
1150 call get_param(param_file, mdl,
"MASK_SRESTORE_MARGINAL_SEAS", &
1151 cs%mask_srestore_marginal_seas, &
1152 "If true, disable SSS restoring in marginal seas. Only used when "//&
1153 "RESTORE_SALINITY is True.", default=.false.)
1154 call get_param(param_file, mdl,
"BASIN_FILE", basin_file, &
1155 "A file in which to find the basin masks, in variable 'basin'.", &
1157 basin_file = trim(cs%inputdir) // trim(basin_file)
1158 call safe_alloc_ptr(cs%basin_mask,isd,ied,jsd,jed) ; cs%basin_mask(:,:) = 1.0
1159 if (cs%mask_srestore_marginal_seas)
then
1160 call mom_read_data(basin_file,
'basin',cs%basin_mask,g%domain, timelevel=1)
1161 do j=jsd,jed ;
do i=isd,ied
1162 if (cs%basin_mask(i,j) >= 6.0)
then ; cs%basin_mask(i,j) = 0.0
1163 else ; cs%basin_mask(i,j) = 1.0 ;
endif
1166 call get_param(param_file, mdl,
"MASK_SRESTORE", cs%mask_srestore, &
1167 "If true, read a file (salt_restore_mask) containing "//&
1168 "a mask for SSS restoring.", default=.false.)
1171 if (restore_temp)
then
1172 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
1173 "The constant that relates the restoring surface fluxes "//&
1174 "to the relative surface anomalies (akin to a piston "//&
1175 "velocity). Note the non-MKS units.", units=
"m day-1", &
1176 fail_if_missing=.true.)
1177 call get_param(param_file, mdl,
"SST_RESTORE_FILE", cs%temp_restore_file, &
1178 "A file in which to find the surface temperature to use for restoring.", &
1179 default=
"temp_restore.nc")
1180 call get_param(param_file, mdl,
"SST_RESTORE_VARIABLE", cs%temp_restore_var_name, &
1181 "The name of the surface temperature variable to read from "//&
1182 "SST_RESTORE_FILE for restoring sst.", &
1185 cs%Flux_const = cs%Flux_const / 86400.0
1187 call get_param(param_file, mdl,
"MAX_DELTA_TRESTORE", cs%max_delta_trestore, &
1188 "The maximum sst difference used in restoring terms.", &
1189 units=
"degC ", default=999.0)
1190 call get_param(param_file, mdl,
"MASK_TRESTORE", cs%mask_trestore, &
1191 "If true, read a file (temp_restore_mask) containing "//&
1192 "a mask for SST restoring.", default=.false.)
1200 call get_param(param_file, mdl,
"CD_TIDES", cs%cd_tides, &
1201 "The drag coefficient that applies to the tides.", &
1202 units=
"nondim", default=1.0e-4)
1203 call get_param(param_file, mdl,
"READ_TIDEAMP", cs%read_TIDEAMP, &
1204 "If true, read a file (given by TIDEAMP_FILE) containing "//&
1205 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1206 if (cs%read_TIDEAMP)
then
1207 call get_param(param_file, mdl,
"TIDEAMP_FILE", tideamp_file, &
1208 "The path to the file containing the spatially varying "//&
1209 "tidal amplitudes with INT_TIDE_DISSIPATION.", &
1210 default=
"tideamp.nc")
1213 call get_param(param_file, mdl,
"UTIDE", cs%utide, &
1214 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1215 units=
"m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1218 call safe_alloc_ptr(cs%TKE_tidal,isd,ied,jsd,jed)
1219 call safe_alloc_ptr(cs%ustar_tidal,isd,ied,jsd,jed)
1221 if (cs%read_TIDEAMP)
then
1222 tideamp_file = trim(cs%inputdir) // trim(tideamp_file)
1223 call mom_read_data(tideamp_file,
'tideamp',cs%TKE_tidal,g%domain,timelevel=1, scale=us%m_to_Z*us%T_to_s)
1224 do j=jsd, jed;
do i=isd, ied
1225 utide = cs%TKE_tidal(i,j)
1226 cs%TKE_tidal(i,j) = g%mask2dT(i,j)*cs%Rho0*cs%cd_tides*(utide*utide*utide)
1227 cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1230 do j=jsd,jed;
do i=isd,ied
1232 cs%TKE_tidal(i,j) = cs%Rho0*cs%cd_tides*(utide*utide*utide)
1233 cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1237 call time_interp_external_init
1242 call get_param(param_file, mdl,
"READ_GUST_2D", cs%read_gust_2d, &
1243 "If true, use a 2-dimensional gustiness supplied from "//&
1244 "an input file", default=.false.)
1245 call get_param(param_file, mdl,
"GUST_CONST", cs%gust_const, &
1246 "The background gustiness in the winds.", &
1247 units=
"Pa", default=0.02, scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z)
1248 if (cs%read_gust_2d)
then
1249 call get_param(param_file, mdl,
"GUST_2D_FILE", gust_file, &
1250 "The file in which the wind gustiness is found in "//&
1251 "variable gustiness.")
1253 call safe_alloc_ptr(cs%gust,isd,ied,jsd,jed)
1254 gust_file = trim(cs%inputdir) // trim(gust_file)
1255 call mom_read_data(gust_file,
'gustiness',cs%gust,g%domain, timelevel=1, &
1256 scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z)
1260 call get_param(param_file, mdl,
"USE_RIGID_SEA_ICE", cs%rigid_sea_ice, &
1261 "If true, sea-ice is rigid enough to exert a "//&
1262 "nonhydrostatic pressure that resist vertical motion.", &
1264 if (cs%rigid_sea_ice)
then
1265 call get_param(param_file, mdl,
"G_EARTH", cs%g_Earth, &
1266 "The gravitational acceleration of the Earth.", &
1267 units=
"m s-2", default = 9.80)
1268 call get_param(param_file, mdl,
"SEA_ICE_MEAN_DENSITY", cs%density_sea_ice, &
1269 "A typical density of sea ice, used with the kinematic "//&
1270 "viscosity, when USE_RIGID_SEA_ICE is true.", units=
"kg m-3", &
1272 call get_param(param_file, mdl,
"SEA_ICE_VISCOSITY", cs%Kv_sea_ice, &
1273 "The kinematic viscosity of sufficiently thick sea ice "//&
1274 "for use in calculating the rigidity of sea ice.", &
1275 units=
"m2 s-1", default=1.0e9)
1276 call get_param(param_file, mdl,
"SEA_ICE_RIGID_MASS", cs%rigid_sea_ice_mass, &
1277 "The mass of sea-ice per unit area at which the sea-ice "//&
1278 "starts to exhibit rigidity", units=
"kg m-2", default=1000.0)
1281 call get_param(param_file, mdl,
"ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, &
1282 "If true, makes available diagnostics of fluxes from icebergs "//&
1283 "as seen by MOM6.", default=.false.)
1284 call register_forcing_type_diags(time, diag, us, cs%use_temperature, cs%handles, &
1285 use_berg_fluxes=iceberg_flux_diags)
1287 call get_param(param_file, mdl,
"ALLOW_FLUX_ADJUSTMENTS", cs%allow_flux_adjustments, &
1288 "If true, allows flux adjustments to specified via the "//&
1289 "data_table using the component name 'OCN'.", default=.false.)
1291 call get_param(param_file, mdl,
"LIQUID_RUNOFF_FROM_DATA", cs%liquid_runoff_from_data, &
1292 "If true, allows liquid river runoff to be specified via the "//&
1293 "data_table using the component name 'OCN'.", default=.false.)
1295 if (cs%allow_flux_adjustments .or. cs%liquid_runoff_from_data)
then
1296 call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
1299 if (
present(restore_salt))
then ;
if (restore_salt)
then
1300 salt_file = trim(cs%inputdir) // trim(cs%salt_restore_file)
1301 cs%id_srestore = init_external_field(salt_file, cs%salt_restore_var_name, domain=g%Domain%mpp_domain)
1302 call safe_alloc_ptr(cs%srestore_mask,isd,ied,jsd,jed); cs%srestore_mask(:,:) = 1.0
1303 if (cs%mask_srestore)
then
1304 flnam = trim(cs%inputdir) //
'salt_restore_mask.nc'
1305 call mom_read_data(flnam,
'mask', cs%srestore_mask, g%domain, timelevel=1)
1309 if (
present(restore_temp))
then ;
if (restore_temp)
then
1310 temp_file = trim(cs%inputdir) // trim(cs%temp_restore_file)
1311 cs%id_trestore = init_external_field(temp_file, cs%temp_restore_var_name, domain=g%Domain%mpp_domain)
1312 call safe_alloc_ptr(cs%trestore_mask,isd,ied,jsd,jed); cs%trestore_mask(:,:) = 1.0
1313 if (cs%mask_trestore)
then
1314 flnam = trim(cs%inputdir) //
'temp_restore_mask.nc'
1315 call mom_read_data(flnam,
'mask', cs%trestore_mask, g%domain, timelevel=1)
1320 call restart_init(param_file, cs%restart_CSp,
"MOM_forcing.res")
1323 if (
associated(cs%restart_CSp))
then
1327 if ((dirs%input_filename(1:1) ==
'n') .and. &
1328 (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1329 if (.not.new_sim)
then
1330 call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
1335 call user_revise_forcing_init(param_file, cs%urf_CS)
1345 type(
forcing),
optional,
intent(inout) :: fluxes
1349 if (
present(fluxes))
call deallocate_forcing_type(fluxes)
1351 if (
associated(cs))
deallocate(cs)
1359 character(len=*),
intent(in) :: id
1360 integer,
intent(in) :: timestep
1366 integer :: n,m, outunit
1370 write(outunit,*)
"BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep
1371 write(outunit,100)
'iobt%u_flux ' , mpp_chksum( iobt%u_flux )
1372 write(outunit,100)
'iobt%v_flux ' , mpp_chksum( iobt%v_flux )
1373 write(outunit,100)
'iobt%t_flux ' , mpp_chksum( iobt%t_flux )
1374 write(outunit,100)
'iobt%q_flux ' , mpp_chksum( iobt%q_flux )
1375 write(outunit,100)
'iobt%salt_flux ' , mpp_chksum( iobt%salt_flux )
1376 write(outunit,100)
'iobt%seaice_melt_heat' , mpp_chksum( iobt%seaice_melt_heat)
1377 write(outunit,100)
'iobt%seaice_melt ' , mpp_chksum( iobt%seaice_melt )
1378 write(outunit,100)
'iobt%lw_flux ' , mpp_chksum( iobt%lw_flux )
1379 write(outunit,100)
'iobt%sw_flux_vis_dir' , mpp_chksum( iobt%sw_flux_vis_dir)
1380 write(outunit,100)
'iobt%sw_flux_vis_dif' , mpp_chksum( iobt%sw_flux_vis_dif)
1381 write(outunit,100)
'iobt%sw_flux_nir_dir' , mpp_chksum( iobt%sw_flux_nir_dir)
1382 write(outunit,100)
'iobt%sw_flux_nir_dif' , mpp_chksum( iobt%sw_flux_nir_dif)
1383 write(outunit,100)
'iobt%lprec ' , mpp_chksum( iobt%lprec )
1384 write(outunit,100)
'iobt%fprec ' , mpp_chksum( iobt%fprec )
1385 write(outunit,100)
'iobt%lrunoff ' , mpp_chksum( iobt%lrunoff )
1386 write(outunit,100)
'iobt%frunoff ' , mpp_chksum( iobt%frunoff )
1387 write(outunit,100)
'iobt%p ' , mpp_chksum( iobt%p )
1388 if (
associated(iobt%ustar_berg)) &
1389 write(outunit,100)
'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg )
1390 if (
associated(iobt%area_berg)) &
1391 write(outunit,100)
'iobt%area_berg ' , mpp_chksum( iobt%area_berg )
1392 if (
associated(iobt%mass_berg)) &
1393 write(outunit,100)
'iobt%mass_berg ' , mpp_chksum( iobt%mass_berg )
1394 100
FORMAT(
" CHECKSUM::",a20,
" = ",z20)
1396 call coupler_type_write_chksums(iobt%fluxes, outunit,
'iobt%')