13 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
35 use mom_time_manager,
only : time_type,
operator(+),
operator(/), get_time, time_type_to_real
59 use data_override_mod,
only : data_override_init, data_override
61 implicit none ;
private
63 #include <MOM_memory.h>
73 logical :: use_temperature
74 logical :: restorebuoy
76 logical :: variable_winds
77 logical :: variable_buoyforce
86 real :: latent_heat_fusion
87 real :: latent_heat_vapor
92 logical :: read_gust_2d
93 real,
pointer :: gust(:,:) => null()
96 real,
pointer :: t_restore(:,:) => null()
97 real,
pointer :: s_restore(:,:) => null()
98 real,
pointer :: dens_restore(:,:) => null()
100 integer :: buoy_last_lev_read = -1
104 real :: gyres_taux_const
105 real :: gyres_taux_sin_amp
106 real :: gyres_taux_cos_amp
107 real :: gyres_taux_n_pis
108 logical :: answers_2018
118 logical :: first_call_set_forcing = .true.
119 logical :: archaic_omip_file = .true.
121 logical :: dataoverrideisinitialized = .false.
124 real :: constantheatforcing
126 character(len=8) :: wind_stagger
135 character(len=200) :: inputdir
136 character(len=200) :: wind_config
137 character(len=200) :: wind_file
138 character(len=200) :: buoy_config
140 character(len=200) :: longwave_file =
''
141 character(len=200) :: shortwave_file =
''
142 character(len=200) :: evaporation_file =
''
143 character(len=200) :: sensibleheat_file =
''
144 character(len=200) :: latentheat_file =
''
146 character(len=200) :: rain_file =
''
147 character(len=200) :: snow_file =
''
148 character(len=200) :: runoff_file =
''
150 character(len=200) :: longwaveup_file =
''
151 character(len=200) :: shortwaveup_file =
''
153 character(len=200) :: sstrestore_file =
''
155 character(len=200) :: salinityrestore_file =
''
158 character(len=80) :: stress_x_var =
''
159 character(len=80) :: stress_y_var =
''
160 character(len=80) :: ustar_var =
''
161 character(len=80) :: lw_var =
''
162 character(len=80) :: sw_var =
''
163 character(len=80) :: latent_var =
''
164 character(len=80) :: sens_var =
''
165 character(len=80) :: evap_var =
''
166 character(len=80) :: rain_var =
''
167 character(len=80) :: snow_var =
''
168 character(len=80) :: lrunoff_var =
''
169 character(len=80) :: frunoff_var =
''
170 character(len=80) :: sst_restore_var =
''
171 character(len=80) :: sss_restore_var =
''
174 integer :: wind_nlev = -1
175 integer :: sw_nlev = -1
176 integer :: lw_nlev = -1
177 integer :: latent_nlev = -1
178 integer :: sens_nlev = -1
179 integer :: evap_nlev = -1
180 integer :: precip_nlev = -1
181 integer :: runoff_nlev = -1
182 integer :: sst_nlev = -1
183 integer :: sss_nlev = -1
186 integer :: wind_last_lev = -1
187 integer :: sw_last_lev = -1
188 integer :: lw_last_lev = -1
189 integer :: latent_last_lev = -1
190 integer :: sens_last_lev = -1
191 integer :: evap_last_lev = -1
192 integer :: precip_last_lev = -1
193 integer :: runoff_last_lev = -1
194 integer :: sst_last_lev = -1
195 integer :: sss_last_lev = -1
220 subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US, CS)
224 type(
forcing),
intent(inout) :: fluxes
225 type(time_type),
intent(in) :: day_start
226 type(time_type),
intent(in) :: day_interval
233 type(time_type) :: day_center
234 integer :: isd, ied, jsd, jed
235 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
240 day_center = day_start + day_interval/2
241 dt = time_type_to_real(day_interval)
243 if (cs%first_call_set_forcing)
then
247 call allocate_forcing_type(g, fluxes, ustar=.true.)
248 if (trim(cs%buoy_config) /=
"NONE")
then
249 if ( cs%use_temperature )
then
250 call allocate_forcing_type(g, fluxes, water=.true., heat=.true., press=.true.)
251 if (cs%restorebuoy)
then
252 call safe_alloc_ptr(cs%T_Restore,isd, ied, jsd, jed)
253 call safe_alloc_ptr(fluxes%heat_added, isd, ied, jsd, jed)
254 call safe_alloc_ptr(cs%S_Restore, isd, ied, jsd, jed)
257 call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
259 if (cs%restorebuoy)
call safe_alloc_ptr(cs%Dens_Restore, isd, ied, jsd, jed)
265 if (cs%variable_winds .or. cs%first_call_set_forcing)
then
267 if (trim(cs%wind_config) ==
"file")
then
269 elseif (trim(cs%wind_config) ==
"data_override")
then
271 elseif (trim(cs%wind_config) ==
"2gyre")
then
273 elseif (trim(cs%wind_config) ==
"1gyre")
then
275 elseif (trim(cs%wind_config) ==
"gyres")
then
277 elseif (trim(cs%wind_config) ==
"zero")
then
279 elseif (trim(cs%wind_config) ==
"const")
then
280 call wind_forcing_const(sfc_state, forces, cs%tau_x0, cs%tau_y0, day_center, g, us, cs)
281 elseif (trim(cs%wind_config) ==
"Neverland")
then
282 call neverland_wind_forcing(sfc_state, forces, day_center, g, us, cs%Neverland_forcing_CSp)
283 elseif (trim(cs%wind_config) ==
"ideal_hurr")
then
284 call idealized_hurricane_wind_forcing(sfc_state, forces, day_center, g, us, cs%idealized_hurricane_CSp)
285 elseif (trim(cs%wind_config) ==
"SCM_ideal_hurr")
then
286 call scm_idealized_hurricane_wind_forcing(sfc_state, forces, day_center, g, us, cs%idealized_hurricane_CSp)
287 elseif (trim(cs%wind_config) ==
"SCM_CVmix_tests")
then
288 call scm_cvmix_tests_wind_forcing(sfc_state, forces, day_center, g, us, cs%SCM_CVmix_tests_CSp)
289 elseif (trim(cs%wind_config) ==
"USER")
then
290 call user_wind_forcing(sfc_state, forces, day_center, g, us, cs%user_forcing_CSp)
291 elseif (cs%variable_winds .and. .not.cs%first_call_set_forcing)
then
292 call mom_error(fatal, &
293 "MOM_surface_forcing: Variable winds defined with no wind config")
295 call mom_error(fatal, &
296 "MOM_surface_forcing:Unrecognized wind config "//trim(cs%wind_config))
301 if ((cs%variable_buoyforce .or. cs%first_call_set_forcing) .and. &
302 (.not.cs%adiabatic))
then
303 if (trim(cs%buoy_config) ==
"file")
then
305 elseif (trim(cs%buoy_config) ==
"data_override")
then
307 elseif (trim(cs%buoy_config) ==
"zero")
then
309 elseif (trim(cs%buoy_config) ==
"const")
then
311 elseif (trim(cs%buoy_config) ==
"linear")
then
313 elseif (trim(cs%buoy_config) ==
"MESO")
then
314 call meso_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%MESO_forcing_CSp)
315 elseif (trim(cs%buoy_config) ==
"Neverland")
then
316 call neverland_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%Neverland_forcing_CSp)
317 elseif (trim(cs%buoy_config) ==
"SCM_CVmix_tests")
then
318 call scm_cvmix_tests_buoyancy_forcing(sfc_state, fluxes, day_center, g, us, cs%SCM_CVmix_tests_CSp)
319 elseif (trim(cs%buoy_config) ==
"USER")
then
320 call user_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%user_forcing_CSp)
321 elseif (trim(cs%buoy_config) ==
"BFB")
then
322 call bfb_buoyancy_forcing(sfc_state, fluxes, day_center, dt, g, us, cs%BFB_forcing_CSp)
323 elseif (trim(cs%buoy_config) ==
"dumbbell")
then
325 elseif (trim(cs%buoy_config) ==
"NONE")
then
326 call mom_mesg(
"MOM_surface_forcing: buoyancy forcing has been set to omitted.")
327 elseif (cs%variable_buoyforce .and. .not.cs%first_call_set_forcing)
then
328 call mom_error(fatal, &
329 "MOM_surface_forcing: Variable buoy defined with no buoy config.")
331 call mom_error(fatal, &
332 "MOM_surface_forcing: Unrecognized buoy config "//trim(cs%buoy_config))
336 if (
associated(cs%tracer_flow_CSp))
then
337 call call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, g, cs%tracer_flow_CSp)
341 call user_alter_forcing(sfc_state, fluxes, day_center, g, cs%urf_CS)
344 if (cs%variable_winds .or. cs%first_call_set_forcing)
then
345 call copy_common_forcing_fields(forces, fluxes, g)
346 call set_derived_forcing_fields(forces, fluxes, g, us, cs%Rho0)
349 if ((cs%variable_buoyforce .or. cs%first_call_set_forcing) .and. &
350 (.not.cs%adiabatic))
then
351 call set_net_mass_forcing(fluxes, forces, g, us)
354 cs%first_call_set_forcing = .false.
366 real,
intent(in) :: tau_x0
367 real,
intent(in) :: tau_y0
368 type(time_type),
intent(in) :: day
374 real :: Pa_conversion
376 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
378 call calltree_enter(
"wind_forcing_const, MOM_surface_forcing.F90")
379 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
380 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
381 pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
385 mag_tau = pa_conversion * sqrt( tau_x0**2 + tau_y0**2)
387 do j=js,je ;
do i=is-1,ieq
388 forces%taux(i,j) = tau_x0 * pa_conversion
391 do j=js-1,jeq ;
do i=is,ie
392 forces%tauy(i,j) = tau_y0 * pa_conversion
395 if (cs%read_gust_2d)
then
396 if (
associated(forces%ustar))
then ;
do j=js,je ;
do i=is,ie
397 forces%ustar(i,j) = sqrt( us%L_to_Z * ( mag_tau + cs%gust(i,j) ) / cs%Rho0 )
398 enddo ;
enddo ;
endif
400 if (
associated(forces%ustar))
then ;
do j=js,je ;
do i=is,ie
401 forces%ustar(i,j) = sqrt( us%L_to_Z * ( mag_tau + cs%gust_const ) / cs%Rho0 )
402 enddo ;
enddo ;
endif
414 type(time_type),
intent(in) :: day
421 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
423 call calltree_enter(
"wind_forcing_2gyre, MOM_surface_forcing.F90")
424 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
425 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
430 do j=js,je ;
do i=is-1,ieq
431 forces%taux(i,j) = 0.1*us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
432 (1.0 - cos(2.0*pi*(g%geoLatCu(i,j)-cs%South_lat) / cs%len_lat))
435 do j=js-1,jeq ;
do i=is,ie
436 forces%tauy(i,j) = 0.0
448 type(time_type),
intent(in) :: day
455 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
457 call calltree_enter(
"wind_forcing_1gyre, MOM_surface_forcing.F90")
458 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
459 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
464 do j=js,je ;
do i=is-1,ieq
465 forces%taux(i,j) = -0.2*us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
466 cos(pi*(g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat)
469 do j=js-1,jeq ;
do i=is,ie
470 forces%tauy(i,j) = 0.0
481 type(time_type),
intent(in) :: day
488 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
490 call calltree_enter(
"wind_forcing_gyres, MOM_surface_forcing.F90")
491 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
492 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
497 do j=js-1,je+1 ;
do i=is-1,ieq
498 y = (g%geoLatCu(i,j)-cs%South_lat) / cs%len_lat
499 forces%taux(i,j) = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
500 (cs%gyres_taux_const + &
501 ( cs%gyres_taux_sin_amp*sin(cs%gyres_taux_n_pis*pi*y) &
502 + cs%gyres_taux_cos_amp*cos(cs%gyres_taux_n_pis*pi*y) ))
505 do j=js-1,jeq ;
do i=is-1,ie+1
506 forces%tauy(i,j) = 0.0
510 if (cs%answers_2018)
then
511 do j=js,je ;
do i=is,ie
512 forces%ustar(i,j) = sqrt(us%L_to_Z * ((cs%gust_const/cs%Rho0) + &
513 sqrt(0.5*(forces%tauy(i,j-1)*forces%tauy(i,j-1) + forces%tauy(i,j)*forces%tauy(i,j) + &
514 forces%taux(i-1,j)*forces%taux(i-1,j) + forces%taux(i,j)*forces%taux(i,j)))/cs%Rho0) )
517 i_rho = us%L_to_Z / cs%Rho0
518 do j=js,je ;
do i=is,ie
519 forces%ustar(i,j) = sqrt( (cs%gust_const + &
520 sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
521 (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * i_rho )
534 type(time_type),
intent(in) :: day
540 character(len=200) :: filename
541 real :: temp_x(SZI_(G),SZJ_(G))
542 real :: temp_y(SZI_(G),SZJ_(G))
543 real :: Pa_conversion
545 integer :: time_lev_daily
546 integer :: time_lev_monthly
548 integer :: days, seconds
549 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
550 logical :: read_Ustar
552 call calltree_enter(
"wind_forcing_from_file, MOM_surface_forcing.F90")
553 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
554 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
555 pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
557 call get_time(day, seconds, days)
558 time_lev_daily = days - 365*floor(real(days) / 365.0)
560 if (time_lev_daily < 31)
then ; time_lev_monthly = 0
561 elseif (time_lev_daily < 59)
then ; time_lev_monthly = 1
562 elseif (time_lev_daily < 90)
then ; time_lev_monthly = 2
563 elseif (time_lev_daily < 120)
then ; time_lev_monthly = 3
564 elseif (time_lev_daily < 151)
then ; time_lev_monthly = 4
565 elseif (time_lev_daily < 181)
then ; time_lev_monthly = 5
566 elseif (time_lev_daily < 212)
then ; time_lev_monthly = 6
567 elseif (time_lev_daily < 243)
then ; time_lev_monthly = 7
568 elseif (time_lev_daily < 273)
then ; time_lev_monthly = 8
569 elseif (time_lev_daily < 304)
then ; time_lev_monthly = 9
570 elseif (time_lev_daily < 334)
then ; time_lev_monthly = 10
571 else ; time_lev_monthly = 11
574 time_lev_daily = time_lev_daily+1
575 time_lev_monthly = time_lev_monthly+1
577 select case (cs%wind_nlev)
578 case (12) ; time_lev = time_lev_monthly
579 case (365) ; time_lev = time_lev_daily
580 case default ; time_lev = 1
583 if (time_lev /= cs%wind_last_lev)
then
584 filename = trim(cs%wind_file)
585 read_ustar = (len_trim(cs%ustar_var) > 0)
589 select case (
uppercase(cs%wind_stagger(1:1)) )
591 temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
593 temp_x(:,:), temp_y(:,:), g%Domain, stagger=agrid, &
594 timelevel=time_lev, scale=pa_conversion)
596 call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
597 do j=js,je ;
do i=is-1,ieq
598 forces%taux(i,j) = 0.5 * cs%wind_scale * (temp_x(i,j) + temp_x(i+1,j))
600 do j=js-1,jeq ;
do i=is,ie
601 forces%tauy(i,j) = 0.5 * cs%wind_scale * (temp_y(i,j) + temp_y(i,j+1))
604 if (.not.read_ustar)
then
605 if (cs%read_gust_2d)
then
606 do j=js,je ;
do i=is,ie
607 forces%ustar(i,j) = sqrt((cs%gust(i,j) + &
608 sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j))) * us%L_to_Z / cs%Rho0)
611 do j=js,je ;
do i=is,ie
612 forces%ustar(i,j) = sqrt(us%L_to_Z * (cs%gust_const/cs%Rho0 + &
613 sqrt(temp_x(i,j)*temp_x(i,j) + temp_y(i,j)*temp_y(i,j)) / cs%Rho0) )
618 if (g%symmetric)
then
619 if (.not.
associated(g%Domain_aux))
call mom_error(fatal, &
620 " wind_forcing_from_file with C-grid input and symmetric memory "//&
621 " called with a non-associated auxiliary domain in the grid type.")
624 temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
626 temp_x(:,:), temp_y(:,:), &
627 g%Domain_aux, stagger=cgrid_ne, timelevel=time_lev, &
629 do j=js,je ;
do i=is,ie
630 forces%taux(i,j) = cs%wind_scale * temp_x(i,j)
631 forces%tauy(i,j) = cs%wind_scale * temp_y(i,j)
636 forces%taux(:,:), forces%tauy(:,:), &
637 g%Domain, stagger=cgrid_ne, timelevel=time_lev, &
640 if (cs%wind_scale /= 1.0)
then
641 do j=js,je ;
do i=isq,ieq
642 forces%taux(i,j) = cs%wind_scale * forces%taux(i,j)
644 do j=jsq,jeq ;
do i=is,ie
645 forces%tauy(i,j) = cs%wind_scale * forces%tauy(i,j)
650 call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)
651 if (.not.read_ustar)
then
652 if (cs%read_gust_2d)
then
653 do j=js, je ;
do i=is, ie
654 forces%ustar(i,j) = sqrt((cs%gust(i,j) + &
655 sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
656 (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * us%L_to_Z / cs%Rho0 )
659 do j=js, je ;
do i=is, ie
660 forces%ustar(i,j) = sqrt(us%L_to_Z * ( (cs%gust_const/cs%Rho0) + &
661 sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
662 (forces%taux(i-1,j)**2 + forces%taux(i,j)**2)))/cs%Rho0))
667 call mom_error(fatal,
"wind_forcing_from_file: Unrecognized stagger "//&
668 trim(cs%wind_stagger)//
" is not 'A' or 'C'.")
672 call mom_read_data(filename, cs%Ustar_var, forces%ustar(:,:), &
673 g%Domain, timelevel=time_lev, scale=us%m_to_Z*us%T_to_s)
676 cs%wind_last_lev = time_lev
689 type(time_type),
intent(in) :: day
695 real :: temp_x(SZI_(G),SZJ_(G))
696 real :: temp_y(SZI_(G),SZJ_(G))
697 real :: temp_ustar(SZI_(G),SZJ_(G))
698 real :: Pa_conversion
699 integer :: i, j, is_in, ie_in, js_in, je_in
700 logical :: read_uStar
702 call calltree_enter(
"wind_forcing_by_data_override, MOM_surface_forcing.F90")
704 if (.not.cs%dataOverrideIsInitialized)
then
706 call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
707 cs%dataOverrideIsInitialized = .true.
710 is_in = g%isc - g%isd + 1 ; ie_in = g%iec - g%isd + 1
711 js_in = g%jsc - g%jsd + 1 ; je_in = g%jec - g%jsd + 1
712 pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
714 temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
715 call data_override(
'OCN',
'taux', temp_x, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
716 call data_override(
'OCN',
'tauy', temp_y, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
717 call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
719 do j=g%jsc,g%jec ;
do i=g%isc-1,g%IecB
720 forces%taux(i,j) = pa_conversion * 0.5 * (temp_x(i,j) + temp_x(i+1,j))
722 do j=g%jsc-1,g%JecB ;
do i=g%isc,g%iec
723 forces%tauy(i,j) = pa_conversion * 0.5 * (temp_y(i,j) + temp_y(i,j+1))
726 read_ustar = (len_trim(cs%ustar_var) > 0)
728 do j=g%jsc,g%jec ;
do i=g%isc,g%iec ; temp_ustar(i,j) = us%Z_to_m*us%s_to_T*forces%ustar(i,j) ;
enddo ;
enddo
729 call data_override(
'OCN',
'ustar', temp_ustar, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
730 do j=g%jsc,g%jec ;
do i=g%isc,g%iec ; forces%ustar(i,j) = us%m_to_Z*us%T_to_s*temp_ustar(i,j) ;
enddo ;
enddo
732 if (cs%read_gust_2d)
then
733 call data_override(
'OCN',
'gust', cs%gust, day, is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
734 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
735 forces%ustar(i,j) = sqrt((pa_conversion * sqrt(temp_x(i,j)*temp_x(i,j) + &
736 temp_y(i,j)*temp_y(i,j)) + cs%gust(i,j)) * us%L_to_Z / cs%Rho0)
739 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
740 forces%ustar(i,j) = sqrt(us%L_to_Z * (pa_conversion*sqrt(temp_x(i,j)*temp_x(i,j) + &
741 temp_y(i,j)*temp_y(i,j))/cs%Rho0 + cs%gust_const/cs%Rho0 ))
746 call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)
757 type(
forcing),
intent(inout) :: fluxes
758 type(time_type),
intent(in) :: day
759 real,
intent(in) :: dt
766 real,
dimension(SZI_(G),SZJ_(G)) :: &
776 real :: kg_m2_s_conversion
780 integer :: time_lev_daily
781 integer :: time_lev_monthly
784 integer :: days, seconds
785 integer :: i, j, is, ie, js, je
787 call calltree_enter(
"buoyancy_forcing_from_files, MOM_surface_forcing.F90")
789 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
790 kg_m2_s_conversion = us%kg_m3_to_R*us%m_to_Z*us%T_to_s
792 if (cs%use_temperature) rhoxcp = us%R_to_kg_m3*cs%Rho0 * fluxes%C_p
795 call get_time(day, seconds, days)
797 time_lev_daily = days - 365*floor(real(days) / 365.0)
799 if (time_lev_daily < 31)
then ; time_lev_monthly = 0
800 elseif (time_lev_daily < 59)
then ; time_lev_monthly = 1
801 elseif (time_lev_daily < 90)
then ; time_lev_monthly = 2
802 elseif (time_lev_daily < 120)
then ; time_lev_monthly = 3
803 elseif (time_lev_daily < 151)
then ; time_lev_monthly = 4
804 elseif (time_lev_daily < 181)
then ; time_lev_monthly = 5
805 elseif (time_lev_daily < 212)
then ; time_lev_monthly = 6
806 elseif (time_lev_daily < 243)
then ; time_lev_monthly = 7
807 elseif (time_lev_daily < 273)
then ; time_lev_monthly = 8
808 elseif (time_lev_daily < 304)
then ; time_lev_monthly = 9
809 elseif (time_lev_daily < 334)
then ; time_lev_monthly = 10
810 else ; time_lev_monthly = 11
813 time_lev_daily = time_lev_daily +1
814 time_lev_monthly = time_lev_monthly+1
816 if (time_lev_daily /= cs%buoy_last_lev_read)
then
819 select case (cs%LW_nlev)
820 case (12) ; time_lev = time_lev_monthly
821 case (365) ; time_lev = time_lev_daily
822 case default ; time_lev = 1
824 call mom_read_data(cs%longwave_file, cs%LW_var, fluxes%LW(:,:), &
825 g%Domain, timelevel=time_lev)
826 if (cs%archaic_OMIP_file)
then
827 call mom_read_data(cs%longwaveup_file,
"lwup_sfc", temp(:,:), g%Domain, &
829 do j=js,je ;
do i=is,ie ; fluxes%LW(i,j) = fluxes%LW(i,j) - temp(i,j) ;
enddo ;
enddo
831 cs%LW_last_lev = time_lev
834 select case (cs%evap_nlev)
835 case (12) ; time_lev = time_lev_monthly
836 case (365) ; time_lev = time_lev_daily
837 case default ; time_lev = 1
839 if (cs%archaic_OMIP_file)
then
840 call mom_read_data(cs%evaporation_file, cs%evap_var, temp(:,:), &
841 g%Domain, timelevel=time_lev)
842 do j=js,je ;
do i=is,ie
843 fluxes%latent(i,j) = -cs%latent_heat_vapor*temp(i,j)
844 fluxes%evap(i,j) = -kg_m2_s_conversion*temp(i,j)
845 fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
848 call mom_read_data(cs%evaporation_file, cs%evap_var, fluxes%evap(:,:), &
849 g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
851 cs%evap_last_lev = time_lev
853 select case (cs%latent_nlev)
854 case (12) ; time_lev = time_lev_monthly
855 case (365) ; time_lev = time_lev_daily
856 case default ; time_lev = 1
858 if (.not.cs%archaic_OMIP_file)
then
859 call mom_read_data(cs%latentheat_file, cs%latent_var, fluxes%latent(:,:), &
860 g%Domain, timelevel=time_lev)
861 do j=js,je ;
do i=is,ie
862 fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
865 cs%latent_last_lev = time_lev
867 select case (cs%sens_nlev)
868 case (12) ; time_lev = time_lev_monthly
869 case (365) ; time_lev = time_lev_daily
870 case default ; time_lev = 1
872 if (cs%archaic_OMIP_file)
then
873 call mom_read_data(cs%sensibleheat_file, cs%sens_var, temp(:,:), &
874 g%Domain, timelevel=time_lev)
875 do j=js,je ;
do i=is,ie ; fluxes%sens(i,j) = -temp(i,j) ;
enddo ;
enddo
877 call mom_read_data(cs%sensibleheat_file, cs%sens_var, fluxes%sens(:,:), &
878 g%Domain, timelevel=time_lev)
880 cs%sens_last_lev = time_lev
882 select case (cs%SW_nlev)
883 case (12) ; time_lev = time_lev_monthly
884 case (365) ; time_lev = time_lev_daily
885 case default ; time_lev = 1
887 call mom_read_data(cs%shortwave_file, cs%SW_var, fluxes%sw(:,:), &
888 g%Domain, timelevel=time_lev)
889 if (cs%archaic_OMIP_file)
then
890 call mom_read_data(cs%shortwaveup_file,
"swup_sfc", temp(:,:), &
891 g%Domain, timelevel=time_lev)
892 do j=js,je ;
do i=is,ie
893 fluxes%sw(i,j) = fluxes%sw(i,j) - temp(i,j)
896 cs%SW_last_lev = time_lev
898 select case (cs%precip_nlev)
899 case (12) ; time_lev = time_lev_monthly
900 case (365) ; time_lev = time_lev_daily
901 case default ; time_lev = 1
904 fluxes%fprec(:,:), g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
906 fluxes%lprec(:,:), g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
907 if (cs%archaic_OMIP_file)
then
908 do j=js,je ;
do i=is,ie
909 fluxes%lprec(i,j) = fluxes%lprec(i,j) - fluxes%fprec(i,j)
912 cs%precip_last_lev = time_lev
914 select case (cs%runoff_nlev)
915 case (12) ; time_lev = time_lev_monthly
916 case (365) ; time_lev = time_lev_daily
917 case default ; time_lev = 1
919 if (cs%archaic_OMIP_file)
then
920 call mom_read_data(cs%runoff_file, cs%lrunoff_var, temp(:,:), &
921 g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
922 do j=js,je ;
do i=is,ie
923 fluxes%lrunoff(i,j) = temp(i,j)*us%m_to_L**2*g%IareaT(i,j)
925 call mom_read_data(cs%runoff_file, cs%frunoff_var, temp(:,:), &
926 g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
927 do j=js,je ;
do i=is,ie
928 fluxes%frunoff(i,j) = temp(i,j)*us%m_to_L**2*g%IareaT(i,j)
931 call mom_read_data(cs%runoff_file, cs%lrunoff_var, fluxes%lrunoff(:,:), &
932 g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
933 call mom_read_data(cs%runoff_file, cs%frunoff_var, fluxes%frunoff(:,:), &
934 g%Domain, timelevel=time_lev, scale=kg_m2_s_conversion)
936 cs%runoff_last_lev = time_lev
939 if (cs%restorebuoy)
then
940 select case (cs%SST_nlev)
941 case (12) ; time_lev = time_lev_monthly
942 case (365) ; time_lev = time_lev_daily
943 case default ; time_lev = 1
946 cs%T_Restore(:,:), g%Domain, timelevel=time_lev)
947 cs%SST_last_lev = time_lev
949 select case (cs%SSS_nlev)
950 case (12) ; time_lev = time_lev_monthly
951 case (365) ; time_lev = time_lev_daily
952 case default ; time_lev = 1
954 call mom_read_data(cs%salinityrestore_file, cs%SSS_restore_var, &
955 cs%S_Restore(:,:), g%Domain, timelevel=time_lev)
956 cs%SSS_last_lev = time_lev
958 cs%buoy_last_lev_read = time_lev_daily
966 do j=js,je ;
do i=is,ie
967 fluxes%evap(i,j) = fluxes%evap(i,j) * g%mask2dT(i,j)
968 fluxes%lprec(i,j) = fluxes%lprec(i,j) * g%mask2dT(i,j)
969 fluxes%fprec(i,j) = fluxes%fprec(i,j) * g%mask2dT(i,j)
970 fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
971 fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
972 fluxes%LW(i,j) = fluxes%LW(i,j) * g%mask2dT(i,j)
973 fluxes%sens(i,j) = fluxes%sens(i,j) * g%mask2dT(i,j)
974 fluxes%sw(i,j) = fluxes%sw(i,j) * g%mask2dT(i,j)
975 fluxes%latent(i,j) = fluxes%latent(i,j) * g%mask2dT(i,j)
977 fluxes%latent_evap_diag(i,j) = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
978 fluxes%latent_fprec_diag(i,j) = -fluxes%fprec(i,j)*cs%latent_heat_fusion
979 fluxes%latent_frunoff_diag(i,j) = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
986 if (cs%restorebuoy)
then
988 if (cs%use_temperature)
then
989 do j=js,je ;
do i=is,ie
990 if (g%mask2dT(i,j) > 0)
then
991 fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
992 ((cs%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoxcp * cs%Flux_const_T)
993 fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const_S) * &
994 (cs%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
995 (0.5*(sfc_state%SSS(i,j) + cs%S_Restore(i,j)))
997 fluxes%heat_added(i,j) = 0.0
998 fluxes%vprec(i,j) = 0.0
1002 do j=js,je ;
do i=is,ie
1003 if (g%mask2dT(i,j) > 0)
then
1004 fluxes%buoy(i,j) = us%kg_m3_to_R * (cs%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
1005 (cs%G_Earth * cs%Flux_const / cs%Rho0)
1007 fluxes%buoy(i,j) = 0.0
1013 if (.not.cs%use_temperature)
then
1014 call mom_error(fatal,
"buoyancy_forcing in MOM_surface_forcing: "// &
1015 "The fluxes need to be defined without RESTOREBUOY.")
1037 type(
forcing),
intent(inout) :: fluxes
1038 type(time_type),
intent(in) :: day
1039 real,
intent(in) :: dt
1046 real,
dimension(SZI_(G),SZJ_(G)) :: &
1055 real :: kg_m2_s_conversion
1059 integer :: time_lev_daily
1060 integer :: time_lev_monthly
1061 integer :: itime_lev
1063 integer :: days, seconds
1064 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1065 integer :: is_in, ie_in, js_in, je_in
1067 call calltree_enter(
"buoyancy_forcing_from_data_override, MOM_surface_forcing.F90")
1069 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1070 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1071 kg_m2_s_conversion = us%kg_m3_to_R*us%m_to_Z*us%T_to_s
1073 if (cs%use_temperature) rhoxcp = cs%Rho0 * fluxes%C_p
1075 if (.not.cs%dataOverrideIsInitialized)
then
1076 call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
1077 cs%dataOverrideIsInitialized = .true.
1080 is_in = g%isc - g%isd + 1
1081 ie_in = g%iec - g%isd + 1
1082 js_in = g%jsc - g%jsd + 1
1083 je_in = g%jec - g%jsd + 1
1085 call data_override(
'OCN',
'lw', fluxes%LW(:,:), day, &
1086 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1087 call data_override(
'OCN',
'evap', fluxes%evap(:,:), day, &
1088 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1091 do j=js,je ;
do i=is,ie
1093 fluxes%evap(i,j) = -fluxes%evap(i,j)
1095 fluxes%latent(i,j) = cs%latent_heat_vapor*fluxes%evap(i,j)
1096 fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)
1097 fluxes%evap(i,j) = kg_m2_s_conversion*fluxes%evap(i,j)
1100 call data_override(
'OCN',
'sens', fluxes%sens(:,:), day, &
1101 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1104 do j=js,je ;
do i=is,ie
1105 fluxes%sens(i,j) = -fluxes%sens(i,j)
1109 call data_override(
'OCN',
'sw', fluxes%sw(:,:), day, &
1110 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1112 call data_override(
'OCN',
'snow', fluxes%fprec(:,:), day, &
1113 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1115 call data_override(
'OCN',
'rain', fluxes%lprec(:,:), day, &
1116 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1118 call data_override(
'OCN',
'runoff', fluxes%lrunoff(:,:), day, &
1119 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1121 call data_override(
'OCN',
'calving', fluxes%frunoff(:,:), day, &
1122 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1124 if (kg_m2_s_conversion /= 1.0)
then ;
do j=js,je ;
do i=is,ie
1125 fluxes%lprec(i,j) = fluxes%lprec(i,j) * kg_m2_s_conversion
1126 fluxes%fprec(i,j) = fluxes%fprec(i,j) * kg_m2_s_conversion
1127 fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * kg_m2_s_conversion
1128 fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * kg_m2_s_conversion
1129 enddo ;
enddo ;
endif
1132 if (cs%restorebuoy)
then
1133 call data_override(
'OCN',
'SST_restore', cs%T_restore(:,:), day, &
1134 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1136 call data_override(
'OCN',
'SSS_restore', cs%S_restore(:,:), day, &
1137 is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in)
1142 if (cs%restorebuoy)
then
1143 if (cs%use_temperature)
then
1144 do j=js,je ;
do i=is,ie
1145 if (g%mask2dT(i,j) > 0)
then
1146 fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
1147 ((cs%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoxcp * cs%Flux_const_T)
1148 fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const_S) * &
1149 (cs%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
1150 (0.5*(sfc_state%SSS(i,j) + cs%S_Restore(i,j)))
1152 fluxes%heat_added(i,j) = 0.0
1153 fluxes%vprec(i,j) = 0.0
1157 do j=js,je ;
do i=is,ie
1158 if (g%mask2dT(i,j) > 0)
then
1159 fluxes%buoy(i,j) = us%kg_m3_to_R * (cs%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
1160 (cs%G_Earth * cs%Flux_const / cs%Rho0)
1162 fluxes%buoy(i,j) = 0.0
1167 if (.not.cs%use_temperature)
then
1168 call mom_error(fatal,
"buoyancy_forcing in MOM_surface_forcing: "// &
1169 "The fluxes need to be defined without RESTOREBUOY.")
1180 do j=js,je ;
do i=is,ie
1181 fluxes%evap(i,j) = fluxes%evap(i,j) * g%mask2dT(i,j)
1182 fluxes%lprec(i,j) = fluxes%lprec(i,j) * g%mask2dT(i,j)
1183 fluxes%fprec(i,j) = fluxes%fprec(i,j) * g%mask2dT(i,j)
1184 fluxes%lrunoff(i,j) = fluxes%lrunoff(i,j) * g%mask2dT(i,j)
1185 fluxes%frunoff(i,j) = fluxes%frunoff(i,j) * g%mask2dT(i,j)
1186 fluxes%LW(i,j) = fluxes%LW(i,j) * g%mask2dT(i,j)
1187 fluxes%latent(i,j) = fluxes%latent(i,j) * g%mask2dT(i,j)
1188 fluxes%sens(i,j) = fluxes%sens(i,j) * g%mask2dT(i,j)
1189 fluxes%sw(i,j) = fluxes%sw(i,j) * g%mask2dT(i,j)
1191 fluxes%latent_evap_diag(i,j) = fluxes%latent_evap_diag(i,j) * g%mask2dT(i,j)
1192 fluxes%latent_fprec_diag(i,j) = -fluxes%fprec(i,j)*cs%latent_heat_fusion
1193 fluxes%latent_frunoff_diag(i,j) = -fluxes%frunoff(i,j)*cs%latent_heat_fusion
1214 type(
forcing),
intent(inout) :: fluxes
1215 type(time_type),
intent(in) :: day
1216 real,
intent(in) :: dt
1222 integer :: i, j, is, ie, js, je
1224 call calltree_enter(
"buoyancy_forcing_zero, MOM_surface_forcing.F90")
1225 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1227 if (cs%use_temperature)
then
1228 do j=js,je ;
do i=is,ie
1229 fluxes%evap(i,j) = 0.0
1230 fluxes%lprec(i,j) = 0.0
1231 fluxes%fprec(i,j) = 0.0
1232 fluxes%vprec(i,j) = 0.0
1233 fluxes%lrunoff(i,j) = 0.0
1234 fluxes%frunoff(i,j) = 0.0
1235 fluxes%lw(i,j) = 0.0
1236 fluxes%latent(i,j) = 0.0
1237 fluxes%sens(i,j) = 0.0
1238 fluxes%sw(i,j) = 0.0
1239 fluxes%latent_evap_diag(i,j) = 0.0
1240 fluxes%latent_fprec_diag(i,j) = 0.0
1241 fluxes%latent_frunoff_diag(i,j) = 0.0
1244 do j=js,je ;
do i=is,ie
1245 fluxes%buoy(i,j) = 0.0
1257 type(
forcing),
intent(inout) :: fluxes
1258 type(time_type),
intent(in) :: day
1259 real,
intent(in) :: dt
1265 integer :: i, j, is, ie, js, je
1266 call calltree_enter(
"buoyancy_forcing_const, MOM_surface_forcing.F90")
1267 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1269 if (cs%use_temperature)
then
1270 do j=js,je ;
do i=is,ie
1271 fluxes%evap(i,j) = 0.0
1272 fluxes%lprec(i,j) = 0.0
1273 fluxes%fprec(i,j) = 0.0
1274 fluxes%vprec(i,j) = 0.0
1275 fluxes%lrunoff(i,j) = 0.0
1276 fluxes%frunoff(i,j) = 0.0
1277 fluxes%lw(i,j) = 0.0
1278 fluxes%latent(i,j) = 0.0
1279 fluxes%sens(i,j) = cs%constantHeatForcing * g%mask2dT(i,j)
1280 fluxes%sw(i,j) = 0.0
1281 fluxes%latent_evap_diag(i,j) = 0.0
1282 fluxes%latent_fprec_diag(i,j) = 0.0
1283 fluxes%latent_frunoff_diag(i,j) = 0.0
1286 do j=js,je ;
do i=is,ie
1287 fluxes%buoy(i,j) = 0.0
1299 type(
forcing),
intent(inout) :: fluxes
1300 type(time_type),
intent(in) :: day
1301 real,
intent(in) :: dt
1308 real :: y, T_restore, S_restore
1309 integer :: i, j, is, ie, js, je
1311 call calltree_enter(
"buoyancy_forcing_linear, MOM_surface_forcing.F90")
1312 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1315 if (cs%use_temperature)
then
1316 do j=js,je ;
do i=is,ie
1317 fluxes%evap(i,j) = 0.0
1318 fluxes%lprec(i,j) = 0.0
1319 fluxes%fprec(i,j) = 0.0
1320 fluxes%vprec(i,j) = 0.0
1321 fluxes%lrunoff(i,j) = 0.0
1322 fluxes%frunoff(i,j) = 0.0
1323 fluxes%lw(i,j) = 0.0
1324 fluxes%latent(i,j) = 0.0
1325 fluxes%sens(i,j) = 0.0
1326 fluxes%sw(i,j) = 0.0
1327 fluxes%latent_evap_diag(i,j) = 0.0
1328 fluxes%latent_fprec_diag(i,j) = 0.0
1329 fluxes%latent_frunoff_diag(i,j) = 0.0
1332 do j=js,je ;
do i=is,ie
1333 fluxes%buoy(i,j) = 0.0
1337 if (cs%restorebuoy)
then
1338 if (cs%use_temperature)
then
1339 do j=js,je ;
do i=is,ie
1340 y = (g%geoLatCu(i,j)-cs%South_lat)/cs%len_lat
1341 t_restore = cs%T_south + (cs%T_north-cs%T_south)*y
1342 s_restore = cs%S_south + (cs%S_north-cs%S_south)*y
1343 if (g%mask2dT(i,j) > 0)
then
1344 fluxes%heat_added(i,j) = g%mask2dT(i,j) * (us%R_to_kg_m3*us%Z_to_m*us%s_to_T) * &
1345 ((t_restore - sfc_state%SST(i,j)) * ((cs%Rho0 * fluxes%C_p) * cs%Flux_const))
1346 fluxes%vprec(i,j) = - (cs%Rho0*cs%Flux_const) * &
1347 (s_restore - sfc_state%SSS(i,j)) / &
1348 (0.5*(sfc_state%SSS(i,j) + s_restore))
1350 fluxes%heat_added(i,j) = 0.0
1351 fluxes%vprec(i,j) = 0.0
1355 call mom_error(fatal,
"buoyancy_forcing_linear in MOM_surface_forcing: "// &
1356 "RESTOREBUOY to linear not written yet.")
1367 if (.not.cs%use_temperature)
then
1368 call mom_error(fatal,
"buoyancy_forcing_linear in MOM_surface_forcing: "// &
1369 "The fluxes need to be defined without RESTOREBUOY.")
1382 type(time_type),
intent(in) :: time
1383 character(len=*),
intent(in) :: directory
1384 logical,
optional,
intent(in) :: time_stamped
1386 character(len=*),
optional,
intent(in) :: filename_suffix
1389 if (.not.
associated(cs))
return
1390 if (.not.
associated(cs%restart_CSp))
return
1392 call save_restart(directory, time, g, cs%restart_CSp, time_stamped)
1398 type(time_type),
intent(in) :: time
1402 type(
diag_ctrl),
target,
intent(inout) :: diag
1410 type(time_type) :: time_frc
1412 # include "version_variable.h"
1413 real :: flux_const_default
1414 logical :: default_2018_answers
1415 character(len=40) :: mdl =
"MOM_surface_forcing"
1416 character(len=200) :: filename, gust_file
1418 if (
associated(cs))
then
1419 call mom_error(warning,
"surface_forcing_init called with an associated "// &
1420 "control structure.")
1425 id_clock_forcing=cpu_clock_id(
'(Ocean surface forcing)', grain=clock_module)
1429 if (
associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
1433 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
1434 "If true, Temperature and salinity are used as state "//&
1435 "variables.", default=.true.)
1436 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, &
1437 "The directory in which all input files are found.", &
1439 cs%inputdir = slasher(cs%inputdir)
1441 call get_param(param_file, mdl,
"ADIABATIC", cs%adiabatic, &
1442 "There are no diapycnal mass fluxes if ADIABATIC is "//&
1443 "true. This assumes that KD = KDML = 0.0 and that "//&
1444 "there is no buoyancy forcing, but makes the model "//&
1445 "faster by eliminating subroutine calls.", default=.false.)
1446 call get_param(param_file, mdl,
"VARIABLE_WINDS", cs%variable_winds, &
1447 "If true, the winds vary in time after the initialization.", &
1449 call get_param(param_file, mdl,
"VARIABLE_BUOYFORCE", cs%variable_buoyforce, &
1450 "If true, the buoyancy forcing varies in time after the "//&
1451 "initialization of the model.", default=.true.)
1453 call get_param(param_file, mdl,
"BUOY_CONFIG", cs%buoy_config, &
1454 "The character string that indicates how buoyancy forcing "//&
1455 "is specified. Valid options include (file), (zero), "//&
1456 "(linear), (USER), (BFB) and (NONE).", fail_if_missing=.true.)
1457 if (trim(cs%buoy_config) ==
"file")
then
1458 call get_param(param_file, mdl,
"ARCHAIC_OMIP_FORCING_FILE", cs%archaic_OMIP_file, &
1459 "If true, use the forcing variable decomposition from "//&
1460 "the old German OMIP prescription that predated CORE. If "//&
1461 "false, use the variable groupings available from MOM "//&
1462 "output diagnostics of forcing variables.", default=.true.)
1463 if (cs%archaic_OMIP_file)
then
1464 call get_param(param_file, mdl,
"LONGWAVEDOWN_FILE", cs%longwave_file, &
1465 "The file with the downward longwave heat flux, in "//&
1466 "variable lwdn_sfc.", fail_if_missing=.true.)
1467 call get_param(param_file, mdl,
"LONGWAVEUP_FILE", cs%longwaveup_file, &
1468 "The file with the upward longwave heat flux, in "//&
1469 "variable lwup_sfc.", fail_if_missing=.true.)
1470 call get_param(param_file, mdl,
"EVAPORATION_FILE", cs%evaporation_file, &
1471 "The file with the evaporative moisture flux, in "//&
1472 "variable evap.", fail_if_missing=.true.)
1473 call get_param(param_file, mdl,
"SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
1474 "The file with the sensible heat flux, in "//&
1475 "variable shflx.", fail_if_missing=.true.)
1476 call get_param(param_file, mdl,
"SHORTWAVEUP_FILE", cs%shortwaveup_file, &
1477 "The file with the upward shortwave heat flux.", &
1478 fail_if_missing=.true.)
1479 call get_param(param_file, mdl,
"SHORTWAVEDOWN_FILE", cs%shortwave_file, &
1480 "The file with the downward shortwave heat flux.", &
1481 fail_if_missing=.true.)
1482 call get_param(param_file, mdl,
"SNOW_FILE", cs%snow_file, &
1483 "The file with the downward frozen precip flux, in "//&
1484 "variable snow.", fail_if_missing=.true.)
1485 call get_param(param_file, mdl,
"PRECIP_FILE", cs%rain_file, &
1486 "The file with the downward total precip flux, in "//&
1487 "variable precip.", fail_if_missing=.true.)
1488 call get_param(param_file, mdl,
"FRESHDISCHARGE_FILE", cs%runoff_file, &
1489 "The file with the fresh and frozen runoff/calving fluxes, "//&
1490 "invariables disch_w and disch_s.", fail_if_missing=.true.)
1493 cs%latentheat_file = cs%evaporation_file ; cs%latent_var =
"evap"
1494 cs%LW_var =
"lwdn_sfc"; cs%SW_var =
"swdn_sfc"; cs%sens_var =
"shflx"
1495 cs%evap_var =
"evap"; cs%rain_var =
"precip"; cs%snow_var =
"snow"
1496 cs%lrunoff_var =
"disch_w"; cs%frunoff_var =
"disch_s"
1499 call get_param(param_file, mdl,
"LONGWAVE_FILE", cs%longwave_file, &
1500 "The file with the longwave heat flux, in the variable "//&
1501 "given by LONGWAVE_FORCING_VAR.", fail_if_missing=.true.)
1502 call get_param(param_file, mdl,
"LONGWAVE_FORCING_VAR", cs%LW_var, &
1503 "The variable with the longwave forcing field.", default=
"LW")
1505 call get_param(param_file, mdl,
"SHORTWAVE_FILE", cs%shortwave_file, &
1506 "The file with the shortwave heat flux, in the variable "//&
1507 "given by SHORTWAVE_FORCING_VAR.", fail_if_missing=.true.)
1508 call get_param(param_file, mdl,
"SHORTWAVE_FORCING_VAR", cs%SW_var, &
1509 "The variable with the shortwave forcing field.", default=
"SW")
1511 call get_param(param_file, mdl,
"EVAPORATION_FILE", cs%evaporation_file, &
1512 "The file with the evaporative moisture flux, in the "//&
1513 "variable given by EVAP_FORCING_VAR.", fail_if_missing=.true.)
1514 call get_param(param_file, mdl,
"EVAP_FORCING_VAR", cs%evap_var, &
1515 "The variable with the evaporative moisture flux.", &
1518 call get_param(param_file, mdl,
"LATENTHEAT_FILE", cs%latentheat_file, &
1519 "The file with the latent heat flux, in the variable "//&
1520 "given by LATENT_FORCING_VAR.", fail_if_missing=.true.)
1521 call get_param(param_file, mdl,
"LATENT_FORCING_VAR", cs%latent_var, &
1522 "The variable with the latent heat flux.", default=
"latent")
1524 call get_param(param_file, mdl,
"SENSIBLEHEAT_FILE", cs%sensibleheat_file, &
1525 "The file with the sensible heat flux, in the variable "//&
1526 "given by SENSIBLE_FORCING_VAR.", fail_if_missing=.true.)
1527 call get_param(param_file, mdl,
"SENSIBLE_FORCING_VAR", cs%sens_var, &
1528 "The variable with the sensible heat flux.", default=
"sensible")
1530 call get_param(param_file, mdl,
"RAIN_FILE", cs%rain_file, &
1531 "The file with the liquid precipitation flux, in the "//&
1532 "variable given by RAIN_FORCING_VAR.", fail_if_missing=.true.)
1533 call get_param(param_file, mdl,
"RAIN_FORCING_VAR", cs%rain_var, &
1534 "The variable with the liquid precipitation flux.", &
1535 default=
"liq_precip")
1536 call get_param(param_file, mdl,
"SNOW_FILE", cs%snow_file, &
1537 "The file with the frozen precipitation flux, in the "//&
1538 "variable given by SNOW_FORCING_VAR.", fail_if_missing=.true.)
1539 call get_param(param_file, mdl,
"SNOW_FORCING_VAR", cs%snow_var, &
1540 "The variable with the frozen precipitation flux.", &
1541 default=
"froz_precip")
1543 call get_param(param_file, mdl,
"RUNOFF_FILE", cs%runoff_file, &
1544 "The file with the fresh and frozen runoff/calving "//&
1545 "fluxes, in variables given by LIQ_RUNOFF_FORCING_VAR "//&
1546 "and FROZ_RUNOFF_FORCING_VAR.", fail_if_missing=.true.)
1547 call get_param(param_file, mdl,
"LIQ_RUNOFF_FORCING_VAR", cs%lrunoff_var, &
1548 "The variable with the liquid runoff flux.", &
1549 default=
"liq_runoff")
1550 call get_param(param_file, mdl,
"FROZ_RUNOFF_FORCING_VAR", cs%frunoff_var, &
1551 "The variable with the frozen runoff flux.", &
1552 default=
"froz_runoff")
1555 call get_param(param_file, mdl,
"SSTRESTORE_FILE", cs%SSTrestore_file, &
1556 "The file with the SST toward which to restore in the "//&
1557 "variable given by SST_RESTORE_VAR.", fail_if_missing=.true.)
1558 call get_param(param_file, mdl,
"SALINITYRESTORE_FILE", cs%salinityrestore_file, &
1559 "The file with the surface salinity toward which to "//&
1560 "restore in the variable given by SSS_RESTORE_VAR.", &
1561 fail_if_missing=.true.)
1563 if (cs%archaic_OMIP_file)
then
1564 cs%SST_restore_var =
"TEMP" ; cs%SSS_restore_var =
"SALT"
1566 call get_param(param_file, mdl,
"SST_RESTORE_VAR", cs%SST_restore_var, &
1567 "The variable with the SST toward which to restore.", &
1569 call get_param(param_file, mdl,
"SSS_RESTORE_VAR", cs%SSS_restore_var, &
1570 "The variable with the SSS toward which to restore.", &
1575 cs%shortwave_file = trim(cs%inputdir)//trim(cs%shortwave_file)
1576 cs%longwave_file = trim(cs%inputdir)//trim(cs%longwave_file)
1577 cs%sensibleheat_file = trim(cs%inputdir)//trim(cs%sensibleheat_file)
1578 cs%latentheat_file = trim(cs%inputdir)//trim(cs%latentheat_file)
1579 cs%evaporation_file = trim(cs%inputdir)//trim(cs%evaporation_file)
1580 cs%snow_file = trim(cs%inputdir)//trim(cs%snow_file)
1581 cs%rain_file = trim(cs%inputdir)//trim(cs%rain_file)
1582 cs%runoff_file = trim(cs%inputdir)//trim(cs%runoff_file)
1584 cs%shortwaveup_file = trim(cs%inputdir)//trim(cs%shortwaveup_file)
1585 cs%longwaveup_file = trim(cs%inputdir)//trim(cs%longwaveup_file)
1587 cs%SSTrestore_file = trim(cs%inputdir)//trim(cs%SSTrestore_file)
1588 cs%salinityrestore_file = trim(cs%inputdir)//trim(cs%salinityrestore_file)
1589 elseif (trim(cs%buoy_config) ==
"const")
then
1590 call get_param(param_file, mdl,
"SENSIBLE_HEAT_FLUX", cs%constantHeatForcing, &
1591 "A constant heat forcing (positive into ocean) applied "//&
1592 "through the sensible heat flux field. ", &
1593 units=
'W/m2', fail_if_missing=.true.)
1595 call get_param(param_file, mdl,
"WIND_CONFIG", cs%wind_config, &
1596 "The character string that indicates how wind forcing "//&
1597 "is specified. Valid options include (file), (2gyre), "//&
1598 "(1gyre), (gyres), (zero), and (USER).", fail_if_missing=.true.)
1599 if (trim(cs%wind_config) ==
"file")
then
1600 call get_param(param_file, mdl,
"WIND_FILE", cs%wind_file, &
1601 "The file in which the wind stresses are found in "//&
1602 "variables STRESS_X and STRESS_Y.", fail_if_missing=.true.)
1603 call get_param(param_file, mdl,
"WINDSTRESS_X_VAR",cs%stress_x_var, &
1604 "The name of the x-wind stress variable in WIND_FILE.", &
1606 call get_param(param_file, mdl,
"WINDSTRESS_Y_VAR", cs%stress_y_var, &
1607 "The name of the y-wind stress variable in WIND_FILE.", &
1609 call get_param(param_file, mdl,
"WINDSTRESS_STAGGER",cs%wind_stagger, &
1610 "A character indicating how the wind stress components "//&
1611 "are staggered in WIND_FILE. This may be A or C for now.", &
1613 call get_param(param_file, mdl,
"WINDSTRESS_SCALE", cs%wind_scale, &
1614 "A value by which the wind stresses in WIND_FILE are rescaled.", &
1615 default=1.0, units=
"nondim")
1616 call get_param(param_file, mdl,
"USTAR_FORCING_VAR", cs%ustar_var, &
1617 "The name of the friction velocity variable in WIND_FILE "//&
1618 "or blank to get ustar from the wind stresses plus the "//&
1619 "gustiness.", default=
" ", units=
"nondim")
1620 cs%wind_file = trim(cs%inputdir) // trim(cs%wind_file)
1622 if (trim(cs%wind_config) ==
"gyres")
then
1623 call get_param(param_file, mdl,
"TAUX_CONST", cs%gyres_taux_const, &
1624 "With the gyres wind_config, the constant offset in the "//&
1625 "zonal wind stress profile: "//&
1626 " A in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1627 units=
"Pa", default=0.0)
1628 call get_param(param_file, mdl,
"TAUX_SIN_AMP",cs%gyres_taux_sin_amp, &
1629 "With the gyres wind_config, the sine amplitude in the "//&
1630 "zonal wind stress profile: "//&
1631 " B in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1632 units=
"Pa", default=0.0)
1633 call get_param(param_file, mdl,
"TAUX_COS_AMP",cs%gyres_taux_cos_amp, &
1634 "With the gyres wind_config, the cosine amplitude in "//&
1635 "the zonal wind stress profile: "//&
1636 " C in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1637 units=
"Pa", default=0.0)
1638 call get_param(param_file, mdl,
"TAUX_N_PIS",cs%gyres_taux_n_pis, &
1639 "With the gyres wind_config, the number of gyres in "//&
1640 "the zonal wind stress profile: "//&
1641 " n in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
1642 units=
"nondim", default=0.0)
1643 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1644 "This sets the default value for the various _2018_ANSWERS parameters.", &
1646 call get_param(param_file, mdl,
"WIND_GYRES_2018_ANSWERS", cs%answers_2018, &
1647 "If true, use the order of arithmetic and expressions that recover the answers "//&
1648 "from the end of 2018. Otherwise, use expressions for the gyre friction velocities "//&
1649 "that are rotationally invariant and more likely to be the same between compilers.", &
1650 default=default_2018_answers)
1652 cs%answers_2018 = .false.
1654 if ((trim(cs%wind_config) ==
"2gyre") .or. &
1655 (trim(cs%wind_config) ==
"1gyre") .or. &
1656 (trim(cs%wind_config) ==
"gyres") .or. &
1657 (trim(cs%buoy_config) ==
"linear"))
then
1658 cs%south_lat = g%south_lat
1659 cs%len_lat = g%len_lat
1661 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
1662 "The mean ocean density used with BOUSSINESQ true to "//&
1663 "calculate accelerations and the mass for conservation "//&
1664 "properties, or with BOUSSINSEQ false to convert some "//&
1665 "parameters from vertical units of m to kg m-2.", &
1666 units=
"kg m-3", default=1035.0, scale=us%kg_m3_to_R)
1667 call get_param(param_file, mdl,
"RESTOREBUOY", cs%restorebuoy, &
1668 "If true, the buoyancy fluxes drive the model back "//&
1669 "toward some specified surface state with a rate "//&
1670 "given by FLUXCONST.", default= .false.)
1671 call get_param(param_file, mdl,
"LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
1672 "The latent heat of fusion.", default=hlf, &
1673 units=
"J/kg", scale=us%R_to_kg_m3*us%Z_to_m*us%s_to_T)
1674 call get_param(param_file, mdl,
"LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
1675 "The latent heat of fusion.", units=
"J/kg", default=hlv)
1676 if (cs%restorebuoy)
then
1677 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
1678 "The constant that relates the restoring surface fluxes "//&
1679 "to the relative surface anomalies (akin to a piston "//&
1680 "velocity). Note the non-MKS units.", &
1681 units=
"m day-1", scale=us%m_to_Z*us%T_to_s, &
1682 fail_if_missing=.true., unscaled=flux_const_default)
1684 if (cs%use_temperature)
then
1685 call get_param(param_file, mdl,
"FLUXCONST_T", cs%Flux_const_T, &
1686 "The constant that relates the restoring surface temperature "//&
1687 "flux to the relative surface anomaly (akin to a piston "//&
1688 "velocity). Note the non-MKS units.", &
1689 units=
"m day-1", scale=1.0, &
1690 default=flux_const_default)
1691 call get_param(param_file, mdl,
"FLUXCONST_S", cs%Flux_const_S, &
1692 "The constant that relates the restoring surface salinity "//&
1693 "flux to the relative surface anomaly (akin to a piston "//&
1694 "velocity). Note the non-MKS units.", &
1695 units=
"m day-1", scale=us%m_to_Z*us%T_to_s, &
1696 default=flux_const_default)
1702 cs%Flux_const = cs%Flux_const / 86400.0
1703 cs%Flux_const_T = cs%Flux_const_T / 86400.0
1704 cs%Flux_const_S = cs%Flux_const_S / 86400.0
1706 if (trim(cs%buoy_config) ==
"linear")
then
1707 call get_param(param_file, mdl,
"SST_NORTH", cs%T_north, &
1708 "With buoy_config linear, the sea surface temperature "//&
1709 "at the northern end of the domain toward which to "//&
1710 "to restore.", units=
"deg C", default=0.0)
1711 call get_param(param_file, mdl,
"SST_SOUTH", cs%T_south, &
1712 "With buoy_config linear, the sea surface temperature "//&
1713 "at the southern end of the domain toward which to "//&
1714 "to restore.", units=
"deg C", default=0.0)
1715 call get_param(param_file, mdl,
"SSS_NORTH", cs%S_north, &
1716 "With buoy_config linear, the sea surface salinity "//&
1717 "at the northern end of the domain toward which to "//&
1718 "to restore.", units=
"PSU", default=35.0)
1719 call get_param(param_file, mdl,
"SSS_SOUTH", cs%S_south, &
1720 "With buoy_config linear, the sea surface salinity "//&
1721 "at the southern end of the domain toward which to "//&
1722 "to restore.", units=
"PSU", default=35.0)
1725 call get_param(param_file, mdl,
"G_EARTH", cs%G_Earth, &
1726 "The gravitational acceleration of the Earth.", &
1727 units=
"m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
1729 call get_param(param_file, mdl,
"GUST_CONST", cs%gust_const, &
1730 "The background gustiness in the winds.", &
1731 units=
"Pa", default=0.02, scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z)
1732 call get_param(param_file, mdl,
"READ_GUST_2D", cs%read_gust_2d, &
1733 "If true, use a 2-dimensional gustiness supplied from "//&
1734 "an input file", default=.false.)
1735 if (cs%read_gust_2d)
then
1736 call get_param(param_file, mdl,
"GUST_2D_FILE", gust_file, &
1737 "The file in which the wind gustiness is found in "//&
1738 "variable gustiness.", fail_if_missing=.true.)
1739 call safe_alloc_ptr(cs%gust,g%isd,g%ied,g%jsd,g%jed)
1740 filename = trim(cs%inputdir) // trim(gust_file)
1741 call mom_read_data(filename,
'gustiness',cs%gust,g%domain, timelevel=1, &
1742 scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z)
1747 if (trim(cs%wind_config) ==
"USER" .or. trim(cs%buoy_config) ==
"USER" )
then
1749 elseif (trim(cs%buoy_config) ==
"BFB" )
then
1751 elseif (trim(cs%buoy_config) ==
"dumbbell" )
then
1752 call dumbbell_surface_forcing_init(time, g, us, param_file, diag, cs%dumbbell_forcing_CSp)
1753 elseif (trim(cs%wind_config) ==
"MESO" .or. trim(cs%buoy_config) ==
"MESO" )
then
1755 elseif (trim(cs%wind_config) ==
"Neverland")
then
1757 elseif (trim(cs%wind_config) ==
"ideal_hurr" .or.&
1758 trim(cs%wind_config) ==
"SCM_ideal_hurr")
then
1759 call idealized_hurricane_wind_init(time, g, us, param_file, cs%idealized_hurricane_CSp)
1760 elseif (trim(cs%wind_config) ==
"const")
then
1761 call get_param(param_file, mdl,
"CONST_WIND_TAUX", cs%tau_x0, &
1762 "With wind_config const, this is the constant zonal "//&
1763 "wind-stress", units=
"Pa", fail_if_missing=.true.)
1764 call get_param(param_file, mdl,
"CONST_WIND_TAUY", cs%tau_y0, &
1765 "With wind_config const, this is the constant meridional "//&
1766 "wind-stress", units=
"Pa", fail_if_missing=.true.)
1767 elseif (trim(cs%wind_config) ==
"SCM_CVmix_tests" .or. &
1768 trim(cs%buoy_config) ==
"SCM_CVmix_tests")
then
1769 call scm_cvmix_tests_surface_forcing_init(time, g, param_file, cs%SCM_CVmix_tests_CSp)
1770 cs%SCM_CVmix_tests_CSp%Rho0 = us%R_to_kg_m3*cs%Rho0
1773 call register_forcing_type_diags(time, diag, us, cs%use_temperature, cs%handles)
1776 call restart_init(param_file, cs%restart_CSp,
"MOM_forcing.res")
1781 if (
associated(cs%restart_CSp))
then
1785 if ((dirs%input_filename(1:1) ==
'n') .and. &
1786 (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1787 if (.not.new_sim)
then
1788 call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
1794 if (trim(cs%buoy_config) ==
"file")
then
1795 cs%SW_nlev =
num_timelevels(cs%shortwave_file, cs%SW_var, min_dims=3)
1796 cs%LW_nlev =
num_timelevels(cs%longwave_file, cs%LW_var, min_dims=3)
1797 cs%latent_nlev =
num_timelevels(cs%latentheat_file, cs%latent_var, 3)
1798 cs%sens_nlev =
num_timelevels(cs%sensibleheat_file, cs%sens_var, min_dims=3)
1800 cs%evap_nlev =
num_timelevels(cs%evaporation_file, cs%evap_var, min_dims=3)
1801 cs%precip_nlev =
num_timelevels(cs%rain_file, cs%rain_var, min_dims=3)
1802 cs%runoff_nlev =
num_timelevels(cs%runoff_file, cs%lrunoff_var, 3)
1804 cs%SST_nlev =
num_timelevels(cs%SSTrestore_file, cs%SST_restore_var, 3)
1805 cs%SSS_nlev =
num_timelevels(cs%salinityrestore_file, cs%SSS_restore_var, 3)
1808 if (trim(cs%wind_config) ==
"file") &
1809 cs%wind_nlev =
num_timelevels(cs%wind_file, cs%stress_x_var, min_dims=3)
1813 call user_revise_forcing_init(param_file, cs%urf_CS)
1823 type(
forcing),
optional,
intent(inout) :: fluxes
1829 if (
present(fluxes))
call deallocate_forcing_type(fluxes)
1833 if (
associated(cs))
deallocate(cs)