15 use mom_diag_mediator,
only : diag_save_grids, diag_restore_grids, diag_copy_storage_to_diag
19 use mom_eos,
only : gsw_sp_from_sr, gsw_pt_from_ct
32 use coupler_types_mod,
only : coupler_type_send_data
34 implicit none ;
private
36 #include <MOM_memory.h>
38 public calculate_diagnostic_fields, register_time_deriv, write_static_fields
40 public register_surface_diags, post_surface_dyn_diags, post_surface_thermo_diags
41 public register_transport_diags, post_transport_diagnostics
42 public mom_diagnostics_init, mom_diagnostics_end
51 real :: mono_n2_column_fraction = 0.
54 real :: mono_n2_depth = -1.
63 real,
pointer,
dimension(:,:,:) :: &
68 real,
pointer,
dimension(:,:,:) :: &
74 real,
pointer,
dimension(:,:,:) :: h_rlay => null()
76 real,
pointer,
dimension(:,:,:) :: uh_rlay => null()
78 real,
pointer,
dimension(:,:,:) :: vh_rlay => null()
80 real,
pointer,
dimension(:,:,:) :: uhgm_rlay => null()
82 real,
pointer,
dimension(:,:,:) :: vhgm_rlay => null()
86 real,
pointer,
dimension(:,:) :: &
90 cfl_cg1_x => null(), &
94 real,
pointer,
dimension(:,:,:) :: &
98 ke_coradv => null(), &
104 ke_horvisc => null(), &
108 integer :: id_u = -1, id_v = -1, id_h = -1
109 integer :: id_e = -1, id_e_d = -1
110 integer :: id_du_dt = -1, id_dv_dt = -1
111 integer :: id_col_ht = -1, id_dh_dt = -1
112 integer :: id_ke = -1, id_dkedt = -1
113 integer :: id_pe_to_ke = -1, id_ke_coradv = -1
114 integer :: id_ke_adv = -1, id_ke_visc = -1
115 integer :: id_ke_horvisc = -1, id_ke_dia = -1
116 integer :: id_uh_rlay = -1, id_vh_rlay = -1
117 integer :: id_uhgm_rlay = -1, id_vhgm_rlay = -1
118 integer :: id_h_rlay = -1, id_rd1 = -1
119 integer :: id_rml = -1, id_rcv = -1
120 integer :: id_cg1 = -1, id_cfl_cg1 = -1
121 integer :: id_cfl_cg1_x = -1, id_cfl_cg1_y = -1
122 integer :: id_cg_ebt = -1, id_rd_ebt = -1
123 integer :: id_p_ebt = -1
124 integer :: id_temp_int = -1, id_salt_int = -1
125 integer :: id_mass_wt = -1, id_col_mass = -1
126 integer :: id_masscello = -1, id_masso = -1
127 integer :: id_volcello = -1
128 integer :: id_tpot = -1, id_sprac = -1
129 integer :: id_tob = -1, id_sob = -1
130 integer :: id_thetaoga = -1, id_soga = -1
131 integer :: id_sosga = -1, id_tosga = -1
132 integer :: id_temp_layer_ave = -1, id_salt_layer_ave = -1
133 integer :: id_pbo = -1
134 integer :: id_thkcello = -1, id_rhoinsitu = -1
135 integer :: id_rhopot0 = -1, id_rhopot2 = -1
136 integer :: id_h_pre_sync = -1
140 type(
p3d) :: var_ptr(max_fields_)
142 type(
p3d) :: deriv(max_fields_)
143 type(
p3d) :: prev_val(max_fields_)
146 integer :: nlay(max_fields_)
147 integer :: num_time_deriv = 0
149 type(group_pass_type) :: pass_ke_uv
158 integer :: id_zos = -1, id_zossq = -1
159 integer :: id_volo = -1, id_speed = -1
160 integer :: id_ssh = -1, id_ssh_ga = -1
161 integer :: id_sst = -1, id_sst_sq = -1, id_sstcon = -1
162 integer :: id_sss = -1, id_sss_sq = -1, id_sssabs = -1
163 integer :: id_ssu = -1, id_ssv = -1
166 integer :: id_fraz = -1
167 integer :: id_salt_deficit = -1
168 integer :: id_heat_pme = -1
169 integer :: id_intern_heat = -1
177 integer :: id_uhtr = -1, id_umo = -1, id_umo_2d = -1
178 integer :: id_vhtr = -1, id_vmo = -1, id_vmo_2d = -1
179 integer :: id_dynamics_h = -1, id_dynamics_h_tendency = -1
185 subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
186 dt, diag_pre_sync, G, GV, US, CS, eta_bt)
190 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
192 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
194 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
196 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
199 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
208 real,
dimension(:,:),
pointer :: p_surf
211 real,
intent(in) :: dt
216 real,
dimension(SZI_(G),SZJ_(G)), &
217 optional,
intent(in) :: eta_bt
222 integer i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
225 real :: rcv(szi_(g),szj_(g),szk_(g))
227 real :: work_3d(szi_(g),szj_(g),szk_(g))
228 real :: work_2d(szi_(g),szj_(g))
231 real :: surface_field(szi_(g),szj_(g))
232 real :: pressure_1d(szi_(g))
243 real,
parameter :: absurdly_small_freq2 = 1e-34
247 real,
dimension(SZK_(G)) :: temp_layer_ave, salt_layer_ave
248 real :: thetaoga, soga, masso, tosga, sosga
250 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
251 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
252 nz = g%ke ; nkmb = gv%nk_rho_varies
255 call diag_save_grids(cs%diag)
256 call diag_copy_storage_to_diag(cs%diag, diag_pre_sync)
258 if (cs%id_h_pre_sync > 0) &
259 call post_data(cs%id_h_pre_sync, diag_pre_sync%h_state, cs%diag, alt_h = diag_pre_sync%h_state)
261 if (cs%id_du_dt>0)
call post_data(cs%id_du_dt, cs%du_dt, cs%diag, alt_h = diag_pre_sync%h_state)
263 if (cs%id_dv_dt>0)
call post_data(cs%id_dv_dt, cs%dv_dt, cs%diag, alt_h = diag_pre_sync%h_state)
265 if (cs%id_dh_dt>0)
call post_data(cs%id_dh_dt, cs%dh_dt, cs%diag, alt_h = diag_pre_sync%h_state)
267 call diag_restore_grids(cs%diag)
269 call calculate_energy_diagnostics(u, v, h, uh, vh, adp, cdp, g, gv, us, cs)
278 if (nkmb==0 .and. nz > 1) nkmb = nz
280 if (loc(cs)==0)
call mom_error(fatal, &
281 "calculate_diagnostic_fields: Module must be initialized before used.")
283 call calculate_derivs(dt, g, cs)
285 if (cs%id_u > 0)
call post_data(cs%id_u, u, cs%diag)
287 if (cs%id_v > 0)
call post_data(cs%id_v, v, cs%diag)
289 if (cs%id_h > 0)
call post_data(cs%id_h, h, cs%diag)
291 if (
associated(cs%e))
then
292 call find_eta(h, tv, g, gv, us, cs%e, eta_bt)
293 if (cs%id_e > 0)
call post_data(cs%id_e, cs%e, cs%diag)
296 if (
associated(cs%e_D))
then
297 if (
associated(cs%e))
then
298 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
299 cs%e_D(i,j,k) = cs%e(i,j,k) + g%bathyT(i,j)
300 enddo ;
enddo ;
enddo
302 call find_eta(h, tv, g, gv, us, cs%e_D, eta_bt)
303 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
304 cs%e_D(i,j,k) = cs%e_D(i,j,k) + g%bathyT(i,j)
305 enddo ;
enddo ;
enddo
308 if (cs%id_e_D > 0)
call post_data(cs%id_e_D, cs%e_D, cs%diag)
312 if (cs%id_masscello > 0)
then
313 do k=1,nz;
do j=js,je ;
do i=is,ie
314 work_3d(i,j,k) = gv%H_to_kg_m2*h(i,j,k)
315 enddo ;
enddo ;
enddo
316 call post_data(cs%id_masscello, work_3d, cs%diag)
320 if (cs%id_masso > 0)
then
322 do k=1,nz ;
do j=js,je ;
do i=is,ie
323 work_2d(i,j) = work_2d(i,j) + (gv%H_to_kg_m2*h(i,j,k)) * g%areaT(i,j)
324 enddo ;
enddo ;
enddo
326 call post_data(cs%id_masso, masso, cs%diag)
330 if (cs%id_thkcello>0 .or. cs%id_volcello>0)
then
331 if (gv%Boussinesq)
then
332 if (cs%id_thkcello > 0)
then ;
if (gv%H_to_m == 1.0)
then
333 call post_data(cs%id_thkcello, h, cs%diag)
335 do k=1,nz;
do j=js,je ;
do i=is,ie
336 work_3d(i,j,k) = gv%H_to_m*h(i,j,k)
337 enddo ;
enddo ;
enddo
338 call post_data(cs%id_thkcello, work_3d, cs%diag)
340 if (cs%id_volcello > 0)
then
341 do k=1,nz;
do j=js,je ;
do i=is,ie
342 work_3d(i,j,k) = ( gv%H_to_m*h(i,j,k) ) * g%areaT(i,j)
343 enddo ;
enddo ;
enddo
344 call post_data(cs%id_volcello, work_3d, cs%diag)
348 if (
associated(p_surf))
then
350 pressure_1d(i) = p_surf(i,j)
359 pressure_1d(i) = pressure_1d(i) + 0.5*gv%H_to_Pa*h(i,j,k)
363 work_3d(:,j,k), is, ie-is+1, tv%eqn_of_state)
365 work_3d(i,j,k) = (gv%H_to_kg_m2*h(i,j,k)) / work_3d(i,j,k)
368 pressure_1d(i) = pressure_1d(i) + 0.5*gv%H_to_Pa*h(i,j,k)
372 if (cs%id_thkcello > 0)
call post_data(cs%id_thkcello, work_3d, cs%diag)
373 if (cs%id_volcello > 0)
then
374 do k=1,nz;
do j=js,je ;
do i=is,ie
375 work_3d(i,j,k) = g%areaT(i,j) * work_3d(i,j,k)
376 enddo ;
enddo ;
enddo
377 call post_data(cs%id_volcello, work_3d, cs%diag)
383 if (tv%T_is_conT)
then
387 if ((cs%id_Tpot > 0) .or. (cs%id_tob > 0))
then
388 do k=1,nz ;
do j=js,je ;
do i=is,ie
389 work_3d(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k))
390 enddo ;
enddo ;
enddo
391 if (cs%id_Tpot > 0)
call post_data(cs%id_Tpot, work_3d, cs%diag)
392 if (cs%id_tob > 0)
call post_data(cs%id_tob, work_3d(:,:,nz), cs%diag, mask=g%mask2dT)
396 if (cs%id_tob > 0)
call post_data(cs%id_tob, tv%T(:,:,nz), cs%diag, mask=g%mask2dT)
400 if (tv%S_is_absS)
then
404 if ((cs%id_Sprac > 0) .or. (cs%id_sob > 0))
then
405 do k=1,nz ;
do j=js,je ;
do i=is,ie
406 work_3d(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k))
407 enddo ;
enddo ;
enddo
408 if (cs%id_Sprac > 0)
call post_data(cs%id_Sprac, work_3d, cs%diag)
409 if (cs%id_sob > 0)
call post_data(cs%id_sob, work_3d(:,:,nz), cs%diag, mask=g%mask2dT)
413 if (cs%id_sob > 0)
call post_data(cs%id_sob, tv%S(:,:,nz), cs%diag, mask=g%mask2dT)
417 if (cs%id_thetaoga>0)
then
418 thetaoga = global_volume_mean(tv%T, h, g, gv)
419 call post_data(cs%id_thetaoga, thetaoga, cs%diag)
423 if (cs%id_tosga > 0)
then
424 do j=js,je ;
do i=is,ie
425 surface_field(i,j) = tv%T(i,j,1)
427 tosga = global_area_mean(surface_field, g)
428 call post_data(cs%id_tosga, tosga, cs%diag)
432 if (cs%id_soga>0)
then
433 soga = global_volume_mean(tv%S, h, g, gv)
434 call post_data(cs%id_soga, soga, cs%diag)
438 if (cs%id_sosga > 0)
then
439 do j=js,je ;
do i=is,ie
440 surface_field(i,j) = tv%S(i,j,1)
442 sosga = global_area_mean(surface_field, g)
443 call post_data(cs%id_sosga, sosga, cs%diag)
447 if (cs%id_temp_layer_ave>0)
then
448 temp_layer_ave = global_layer_mean(tv%T, h, g, gv)
449 call post_data(cs%id_temp_layer_ave, temp_layer_ave, cs%diag)
453 if (cs%id_salt_layer_ave>0)
then
454 salt_layer_ave = global_layer_mean(tv%S, h, g, gv)
455 call post_data(cs%id_salt_layer_ave, salt_layer_ave, cs%diag)
458 call calculate_vertical_integrals(h, tv, p_surf, g, gv, us, cs)
460 if ((cs%id_Rml > 0) .or. (cs%id_Rcv > 0) .or.
associated(cs%h_Rlay) .or. &
461 associated(cs%uh_Rlay) .or.
associated(cs%vh_Rlay) .or. &
462 associated(cs%uhGM_Rlay) .or.
associated(cs%vhGM_Rlay))
then
464 if (
associated(tv%eqn_of_state))
then
465 pressure_1d(:) = tv%P_Ref
467 do k=1,nz ;
do j=js-1,je+1
469 rcv(:,j,k), is-1, ie-is+3, tv%eqn_of_state)
472 do k=1,nz ;
do j=js-1,je+1 ;
do i=is-1,ie+1
473 rcv(i,j,k) = gv%Rlay(k)
474 enddo ;
enddo ;
enddo
476 if (cs%id_Rml > 0)
call post_data(cs%id_Rml, rcv, cs%diag)
477 if (cs%id_Rcv > 0)
call post_data(cs%id_Rcv, rcv, cs%diag)
479 if (
associated(cs%h_Rlay))
then
484 do k=1,nkmb ;
do i=is,ie
485 cs%h_Rlay(i,j,k) = 0.0
487 do k=nkmb+1,nz ;
do i=is,ie
488 cs%h_Rlay(i,j,k) = h(i,j,k)
490 do k=1,nkmb ;
do i=is,ie
491 call find_weights(gv%Rlay, rcv(i,j,k), k_list, nz, wt, wt_p)
492 cs%h_Rlay(i,j,k_list) = cs%h_Rlay(i,j,k_list) + h(i,j,k)*wt
493 cs%h_Rlay(i,j,k_list+1) = cs%h_Rlay(i,j,k_list+1) + h(i,j,k)*wt_p
497 if (cs%id_h_Rlay > 0)
call post_data(cs%id_h_Rlay, cs%h_Rlay, cs%diag)
500 if (
associated(cs%uh_Rlay))
then
505 do k=1,nkmb ;
do i=isq,ieq
506 cs%uh_Rlay(i,j,k) = 0.0
508 do k=nkmb+1,nz ;
do i=isq,ieq
509 cs%uh_Rlay(i,j,k) = uh(i,j,k)
512 do k=1,nkmb ;
do i=isq,ieq
513 call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
514 cs%uh_Rlay(i,j,k_list) = cs%uh_Rlay(i,j,k_list) + uh(i,j,k)*wt
515 cs%uh_Rlay(i,j,k_list+1) = cs%uh_Rlay(i,j,k_list+1) + uh(i,j,k)*wt_p
519 if (cs%id_uh_Rlay > 0)
call post_data(cs%id_uh_Rlay, cs%uh_Rlay, cs%diag)
522 if (
associated(cs%vh_Rlay))
then
527 do k=1,nkmb ;
do i=is,ie
528 cs%vh_Rlay(i,j,k) = 0.0
530 do k=nkmb+1,nz ;
do i=is,ie
531 cs%vh_Rlay(i,j,k) = vh(i,j,k)
533 do k=1,nkmb ;
do i=is,ie
534 call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
535 cs%vh_Rlay(i,j,k_list) = cs%vh_Rlay(i,j,k_list) + vh(i,j,k)*wt
536 cs%vh_Rlay(i,j,k_list+1) = cs%vh_Rlay(i,j,k_list+1) + vh(i,j,k)*wt_p
540 if (cs%id_vh_Rlay > 0)
call post_data(cs%id_vh_Rlay, cs%vh_Rlay, cs%diag)
543 if (
associated(cs%uhGM_Rlay) .and.
associated(cdp%uhGM))
then
548 do k=1,nkmb ;
do i=isq,ieq
549 cs%uhGM_Rlay(i,j,k) = 0.0
551 do k=nkmb+1,nz ;
do i=isq,ieq
552 cs%uhGM_Rlay(i,j,k) = cdp%uhGM(i,j,k)
554 do k=1,nkmb ;
do i=isq,ieq
555 call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
556 cs%uhGM_Rlay(i,j,k_list) = cs%uhGM_Rlay(i,j,k_list) + cdp%uhGM(i,j,k)*wt
557 cs%uhGM_Rlay(i,j,k_list+1) = cs%uhGM_Rlay(i,j,k_list+1) + cdp%uhGM(i,j,k)*wt_p
561 if (cs%id_uh_Rlay > 0)
call post_data(cs%id_uhGM_Rlay, cs%uhGM_Rlay, cs%diag)
564 if (
associated(cs%vhGM_Rlay) .and.
associated(cdp%vhGM))
then
569 do k=1,nkmb ;
do i=is,ie
570 cs%vhGM_Rlay(i,j,k) = 0.0
572 do k=nkmb+1,nz ;
do i=is,ie
573 cs%vhGM_Rlay(i,j,k) = cdp%vhGM(i,j,k)
575 do k=1,nkmb ;
do i=is,ie
576 call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
577 cs%vhGM_Rlay(i,j,k_list) = cs%vhGM_Rlay(i,j,k_list) + cdp%vhGM(i,j,k)*wt
578 cs%vhGM_Rlay(i,j,k_list+1) = cs%vhGM_Rlay(i,j,k_list+1) + cdp%vhGM(i,j,k)*wt_p
582 if (cs%id_vhGM_Rlay > 0)
call post_data(cs%id_vhGM_Rlay, cs%vhGM_Rlay, cs%diag)
586 if (
associated(tv%eqn_of_state))
then
587 if (cs%id_rhopot0 > 0)
then
590 do k=1,nz ;
do j=js,je
592 rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
594 if (cs%id_rhopot0 > 0)
call post_data(cs%id_rhopot0, rcv, cs%diag)
596 if (cs%id_rhopot2 > 0)
then
597 pressure_1d(:) = 2.e7
599 do k=1,nz ;
do j=js,je
601 rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
603 if (cs%id_rhopot2 > 0)
call post_data(cs%id_rhopot2, rcv, cs%diag)
605 if (cs%id_rhoinsitu > 0)
then
610 pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa
612 rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
613 pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa
616 if (cs%id_rhoinsitu > 0)
call post_data(cs%id_rhoinsitu, rcv, cs%diag)
620 if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
621 (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0))
then
622 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
623 if (cs%id_cg1>0)
call post_data(cs%id_cg1, cs%cg1, cs%diag)
624 if (cs%id_Rd1>0)
then
627 do j=js,je ;
do i=is,ie
629 f2_h = absurdly_small_freq2 + 0.25 * us%s_to_T**2 * &
630 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
631 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
632 mag_beta = sqrt(0.5 * us%s_to_T**2 * ( &
633 (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
634 ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
635 (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
636 ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
637 cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)
640 call post_data(cs%id_Rd1, cs%Rd1, cs%diag)
642 if (cs%id_cfl_cg1>0)
then
643 do j=js,je ;
do i=is,ie
644 cs%cfl_cg1(i,j) = (dt*cs%cg1(i,j)) * (g%IdxT(i,j) + g%IdyT(i,j))
646 call post_data(cs%id_cfl_cg1, cs%cfl_cg1, cs%diag)
648 if (cs%id_cfl_cg1_x>0)
then
649 do j=js,je ;
do i=is,ie
650 cs%cfl_cg1_x(i,j) = (dt*cs%cg1(i,j)) * g%IdxT(i,j)
652 call post_data(cs%id_cfl_cg1_x, cs%cfl_cg1_x, cs%diag)
654 if (cs%id_cfl_cg1_y>0)
then
655 do j=js,je ;
do i=is,ie
656 cs%cfl_cg1_y(i,j) = (dt*cs%cg1(i,j)) * g%IdyT(i,j)
658 call post_data(cs%id_cfl_cg1_y, cs%cfl_cg1_y, cs%diag)
661 if ((cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0))
then
662 if (cs%id_p_ebt>0)
then
663 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
664 mono_n2_column_fraction=cs%mono_N2_column_fraction, &
665 mono_n2_depth=cs%mono_N2_depth, modal_structure=cs%p_ebt)
666 call post_data(cs%id_p_ebt, cs%p_ebt, cs%diag)
668 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
669 mono_n2_column_fraction=cs%mono_N2_column_fraction, &
670 mono_n2_depth=cs%mono_N2_depth)
672 if (cs%id_cg_ebt>0)
call post_data(cs%id_cg_ebt, cs%cg1, cs%diag)
673 if (cs%id_Rd_ebt>0)
then
676 do j=js,je ;
do i=is,ie
678 f2_h = absurdly_small_freq2 + 0.25 * us%s_to_T**2 * &
679 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
680 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
681 mag_beta = sqrt(0.5 * us%s_to_T**2 * ( &
682 (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
683 ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
684 (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
685 ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
686 cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)
689 call post_data(cs%id_Rd_ebt, cs%Rd1, cs%diag)
693 end subroutine calculate_diagnostic_fields
699 subroutine find_weights(Rlist, R_in, k, nz, wt, wt_p)
700 real,
dimension(:), &
702 real,
intent(in) :: R_in
703 integer,
intent(inout) :: k
705 integer,
intent(in) :: nz
706 real,
intent(out) :: wt
707 real,
intent(out) :: wt_p
714 integer :: k_upper, k_lower, k_new, inc
717 if ((k < 1) .or. (k > nz)) k = nz/2
719 k_upper = k ; k_lower = k ; inc = 1
720 if (r_in < rlist(k))
then
722 k_lower = max(k_lower-inc, 1)
723 if ((k_lower == 1) .or. (r_in >= rlist(k_lower)))
exit
729 k_upper = min(k_upper+inc, nz)
730 if ((k_upper == nz) .or. (r_in < rlist(k_upper)))
exit
736 if ((k_lower == 1) .and. (r_in <= rlist(k_lower)))
then
737 k = 1 ; wt = 1.0 ; wt_p = 0.0
738 elseif ((k_upper == nz) .and. (r_in >= rlist(k_upper)))
then
739 k = nz-1 ; wt = 0.0 ; wt_p = 1.0
742 if (k_upper <= k_lower+1)
exit
743 k_new = (k_upper + k_lower) / 2
744 if (r_in < rlist(k_new))
then
756 wt = (rlist(k_upper) - r_in) / (rlist(k_upper) - rlist(k_lower))
761 end subroutine find_weights
766 subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
770 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
774 real,
dimension(:,:),
pointer :: p_surf
780 real,
dimension(SZI_(G), SZJ_(G)) :: &
796 integer :: i, j, k, is, ie, js, je, nz
797 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
799 if (cs%id_mass_wt > 0)
then
800 do j=js,je ;
do i=is,ie ; mass(i,j) = 0.0 ;
enddo ;
enddo
801 do k=1,nz ;
do j=js,je ;
do i=is,ie
802 mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
803 enddo ;
enddo ;
enddo
804 call post_data(cs%id_mass_wt, mass, cs%diag)
807 if (cs%id_temp_int > 0)
then
808 do j=js,je ;
do i=is,ie ; tr_int(i,j) = 0.0 ;
enddo ;
enddo
809 do k=1,nz ;
do j=js,je ;
do i=is,ie
810 tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%T(i,j,k)
811 enddo ;
enddo ;
enddo
812 call post_data(cs%id_temp_int, tr_int, cs%diag)
815 if (cs%id_salt_int > 0)
then
816 do j=js,je ;
do i=is,ie ; tr_int(i,j) = 0.0 ;
enddo ;
enddo
817 do k=1,nz ;
do j=js,je ;
do i=is,ie
818 tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%S(i,j,k)
819 enddo ;
enddo ;
enddo
820 call post_data(cs%id_salt_int, tr_int, cs%diag)
823 if (cs%id_col_ht > 0)
then
824 call find_eta(h, tv, g, gv, us, z_top)
825 do j=js,je ;
do i=is,ie
826 z_bot(i,j) = z_top(i,j) + g%bathyT(i,j)
828 call post_data(cs%id_col_ht, z_bot, cs%diag)
831 if (cs%id_col_mass > 0 .or. cs%id_pbo > 0)
then
832 do j=js,je ;
do i=is,ie ; mass(i,j) = 0.0 ;
enddo ;
enddo
833 if (gv%Boussinesq)
then
834 if (
associated(tv%eqn_of_state))
then
835 ig_earth = 1.0 / gv%mks_g_Earth
837 do j=js,je ;
do i=is,ie ; z_bot(i,j) = 0.0 ;
enddo ;
enddo
839 do j=js,je ;
do i=is,ie
840 z_top(i,j) = z_bot(i,j)
841 z_bot(i,j) = z_top(i,j) - gv%H_to_Z*h(i,j,k)
843 call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), &
844 z_top, z_bot, 0.0, gv%Rho0, gv%mks_g_Earth*us%Z_to_m, &
845 g%HI, g%HI, tv%eqn_of_state, dpress)
846 do j=js,je ;
do i=is,ie
847 mass(i,j) = mass(i,j) + dpress(i,j) * ig_earth
851 do k=1,nz ;
do j=js,je ;
do i=is,ie
852 mass(i,j) = mass(i,j) + (gv%H_to_m*gv%Rlay(k))*h(i,j,k)
853 enddo ;
enddo ;
enddo
856 do k=1,nz ;
do j=js,je ;
do i=is,ie
857 mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
858 enddo ;
enddo ;
enddo
860 if (cs%id_col_mass > 0)
then
861 call post_data(cs%id_col_mass, mass, cs%diag)
863 if (cs%id_pbo > 0)
then
864 do j=js,je ;
do i=is,ie ; btm_pres(i,j) = 0.0 ;
enddo ;
enddo
868 do j=js,je ;
do i=is,ie
869 btm_pres(i,j) = mass(i,j) * gv%mks_g_Earth
870 if (
associated(p_surf))
then
871 btm_pres(i,j) = btm_pres(i,j) + p_surf(i,j)
874 call post_data(cs%id_pbo, btm_pres, cs%diag)
878 end subroutine calculate_vertical_integrals
881 subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS)
884 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
886 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
888 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
890 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
893 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
905 real :: KE_u(SZIB_(G),SZJ_(G))
906 real :: KE_v(SZI_(G),SZJB_(G))
907 real :: KE_h(SZI_(G),SZJ_(G))
909 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
910 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
911 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
913 do j=js-1,je ;
do i=is-1,ie
914 ke_u(i,j) = 0.0 ; ke_v(i,j) = 0.0
917 if (
associated(cs%KE))
then
918 do k=1,nz ;
do j=js,je ;
do i=is,ie
919 cs%KE(i,j,k) = ((u(i,j,k)*u(i,j,k) + u(i-1,j,k)*u(i-1,j,k)) + &
920 (v(i,j,k)*v(i,j,k) + v(i,j-1,k)*v(i,j-1,k)))*0.25
924 enddo ;
enddo ;
enddo
925 if (cs%id_KE > 0)
call post_data(cs%id_KE, cs%KE, cs%diag)
928 if (.not.g%symmetric)
then
929 if (
associated(cs%dKE_dt) .OR.
associated(cs%PE_to_KE) .OR.
associated(cs%KE_CorAdv) .OR. &
930 associated(cs%KE_adv) .OR.
associated(cs%KE_visc) .OR.
associated(cs%KE_horvisc).OR. &
931 associated(cs%KE_dia) )
then
936 if (
associated(cs%dKE_dt))
then
938 do j=js,je ;
do i=isq,ieq
939 ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*cs%du_dt(i,j,k)
941 do j=jsq,jeq ;
do i=is,ie
942 ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*cs%dv_dt(i,j,k)
944 do j=js,je ;
do i=is,ie
945 ke_h(i,j) = cs%KE(i,j,k)*cs%dh_dt(i,j,k)
947 if (.not.g%symmetric) &
948 call do_group_pass(cs%pass_KE_uv, g%domain)
949 do j=js,je ;
do i=is,ie
950 cs%dKE_dt(i,j,k) = gv%H_to_m * (ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
951 (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1)))
954 if (cs%id_dKEdt > 0)
call post_data(cs%id_dKEdt, cs%dKE_dt, cs%diag)
957 if (
associated(cs%PE_to_KE))
then
959 do j=js,je ;
do i=isq,ieq
960 ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%PFu(i,j,k)
962 do j=jsq,jeq ;
do i=is,ie
963 ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%PFv(i,j,k)
965 if (.not.g%symmetric) &
966 call do_group_pass(cs%pass_KE_uv, g%domain)
967 do j=js,je ;
do i=is,ie
968 cs%PE_to_KE(i,j,k) = gv%H_to_m * 0.5 * g%IareaT(i,j) * &
969 (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
972 if (cs%id_PE_to_KE > 0)
call post_data(cs%id_PE_to_KE, cs%PE_to_KE, cs%diag)
975 if (
associated(cs%KE_CorAdv))
then
977 do j=js,je ;
do i=isq,ieq
978 ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%CAu(i,j,k)
980 do j=jsq,jeq ;
do i=is,ie
981 ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%CAv(i,j,k)
983 do j=js,je ;
do i=is,ie
984 ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) * &
985 (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
987 if (.not.g%symmetric) &
988 call do_group_pass(cs%pass_KE_uv, g%domain)
989 do j=js,je ;
do i=is,ie
990 cs%KE_CorAdv(i,j,k) = gv%H_to_m * (ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
991 (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1)))
994 if (cs%id_KE_Coradv > 0)
call post_data(cs%id_KE_Coradv, cs%KE_Coradv, cs%diag)
997 if (
associated(cs%KE_adv))
then
1001 ke_u(:,:) = 0. ; ke_v(:,:) = 0.
1003 do j=js,je ;
do i=isq,ieq
1004 if (g%mask2dCu(i,j) /= 0.) &
1005 ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%gradKEu(i,j,k)
1007 do j=jsq,jeq ;
do i=is,ie
1008 if (g%mask2dCv(i,j) /= 0.) &
1009 ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%gradKEv(i,j,k)
1011 do j=js,je ;
do i=is,ie
1012 ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) * &
1013 (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
1015 if (.not.g%symmetric) &
1016 call do_group_pass(cs%pass_KE_uv, g%domain)
1017 do j=js,je ;
do i=is,ie
1018 cs%KE_adv(i,j,k) = gv%H_to_m * (ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
1019 (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1)))
1022 if (cs%id_KE_adv > 0)
call post_data(cs%id_KE_adv, cs%KE_adv, cs%diag)
1025 if (
associated(cs%KE_visc))
then
1027 do j=js,je ;
do i=isq,ieq
1028 ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%du_dt_visc(i,j,k)
1030 do j=jsq,jeq ;
do i=is,ie
1031 ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%dv_dt_visc(i,j,k)
1033 if (.not.g%symmetric) &
1034 call do_group_pass(cs%pass_KE_uv, g%domain)
1035 do j=js,je ;
do i=is,ie
1036 cs%KE_visc(i,j,k) = gv%H_to_m * (0.5 * g%IareaT(i,j) * &
1037 (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1)))
1040 if (cs%id_KE_visc > 0)
call post_data(cs%id_KE_visc, cs%KE_visc, cs%diag)
1043 if (
associated(cs%KE_horvisc))
then
1045 do j=js,je ;
do i=isq,ieq
1046 ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*us%s_to_T*adp%diffu(i,j,k)
1048 do j=jsq,jeq ;
do i=is,ie
1049 ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*us%s_to_T*adp%diffv(i,j,k)
1051 if (.not.g%symmetric) &
1052 call do_group_pass(cs%pass_KE_uv, g%domain)
1053 do j=js,je ;
do i=is,ie
1054 cs%KE_horvisc(i,j,k) = gv%H_to_m * 0.5 * g%IareaT(i,j) * &
1055 (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
1058 if (cs%id_KE_horvisc > 0)
call post_data(cs%id_KE_horvisc, cs%KE_horvisc, cs%diag)
1061 if (
associated(cs%KE_dia))
then
1063 do j=js,je ;
do i=isq,ieq
1064 ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%du_dt_dia(i,j,k)
1066 do j=jsq,jeq ;
do i=is,ie
1067 ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%dv_dt_dia(i,j,k)
1069 do j=js,je ;
do i=is,ie
1070 ke_h(i,j) = cs%KE(i,j,k) * &
1071 (cdp%diapyc_vel(i,j,k) - cdp%diapyc_vel(i,j,k+1))
1073 if (.not.g%symmetric) &
1074 call do_group_pass(cs%pass_KE_uv, g%domain)
1075 do j=js,je ;
do i=is,ie
1076 cs%KE_dia(i,j,k) = ke_h(i,j) + gv%H_to_m * 0.5 * g%IareaT(i,j) * &
1077 (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
1080 if (cs%id_KE_dia > 0)
call post_data(cs%id_KE_dia, cs%KE_dia, cs%diag)
1083 end subroutine calculate_energy_diagnostics
1086 subroutine register_time_deriv(lb, f_ptr, deriv_ptr, CS)
1087 integer,
intent(in),
dimension(3) :: lb
1088 real,
dimension(lb(1):,lb(2):,:),
target :: f_ptr
1090 real,
dimension(lb(1):,lb(2):,:),
target :: deriv_ptr
1102 if (.not.
associated(cs))
call mom_error(fatal, &
1103 "register_time_deriv: Module must be initialized before it is used.")
1105 if (cs%num_time_deriv >= max_fields_)
then
1106 call mom_error(warning,
"MOM_diagnostics: Attempted to register more than " // &
1107 "MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.")
1111 m = cs%num_time_deriv+1 ; cs%num_time_deriv = m
1113 ub(:) = lb(:) + shape(f_ptr) - 1
1115 cs%nlay(m) =
size(f_ptr, 3)
1116 cs%deriv(m)%p => deriv_ptr
1117 allocate(cs%prev_val(m)%p(lb(1):ub(1), lb(2):ub(2), cs%nlay(m)))
1119 cs%var_ptr(m)%p => f_ptr
1120 cs%prev_val(m)%p(:,:,:) = f_ptr(:,:,:)
1122 end subroutine register_time_deriv
1125 subroutine calculate_derivs(dt, G, CS)
1126 real,
intent(in) :: dt
1135 if (dt > 0.0)
then ; idt = 1.0/dt
1136 else ;
return ;
endif
1147 do m=1,cs%num_time_deriv
1148 do k=1,cs%nlay(m) ;
do j=g%jsc-1,g%jec ;
do i=g%isc-1,g%iec
1149 cs%deriv(m)%p(i,j,k) = (cs%var_ptr(m)%p(i,j,k) - cs%prev_val(m)%p(i,j,k)) * idt
1150 cs%prev_val(m)%p(i,j,k) = cs%var_ptr(m)%p(i,j,k)
1151 enddo ;
enddo ;
enddo
1154 end subroutine calculate_derivs
1159 subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh)
1163 type(
surface),
intent(in) :: sfc_state
1164 real,
dimension(SZI_(G),SZJ_(G)), &
1168 real,
dimension(SZI_(G),SZJ_(G)) :: work_2d
1169 integer :: i, j, is, ie, js, je
1171 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1173 if (ids%id_ssh > 0) &
1174 call post_data(ids%id_ssh, ssh, diag, mask=g%mask2dT)
1176 if (ids%id_ssu > 0) &
1177 call post_data(ids%id_ssu, sfc_state%u, diag, mask=g%mask2dCu)
1179 if (ids%id_ssv > 0) &
1180 call post_data(ids%id_ssv, sfc_state%v, diag, mask=g%mask2dCv)
1182 if (ids%id_speed > 0)
then
1183 do j=js,je ;
do i=is,ie
1184 work_2d(i,j) = sqrt(0.5*(sfc_state%u(i-1,j)**2 + sfc_state%u(i,j)**2) + &
1185 0.5*(sfc_state%v(i,j-1)**2 + sfc_state%v(i,j)**2))
1187 call post_data(ids%id_speed, work_2d, diag, mask=g%mask2dT)
1190 end subroutine post_surface_dyn_diags
1195 subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv, &
1202 real,
intent(in) :: dt_int
1203 type(
surface),
intent(in) :: sfc_state
1205 real,
dimension(SZI_(G),SZJ_(G)), &
1207 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: ssh_ibc
1210 real,
dimension(SZI_(G),SZJ_(G)) :: work_2d
1211 real,
dimension(SZI_(G),SZJ_(G)) :: &
1214 real :: zos_area_mean, volo, ssh_ga
1215 integer :: i, j, is, ie, js, je
1217 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1220 if (ids%id_ssh_ga > 0)
then
1221 ssh_ga = global_area_mean(ssh, g)
1222 call post_data(ids%id_ssh_ga, ssh_ga, diag)
1228 if (ids%id_zos > 0 .or. ids%id_zossq > 0)
then
1230 do j=js,je ;
do i=is,ie
1231 zos(i,j) = ssh_ibc(i,j)
1233 zos_area_mean = global_area_mean(zos, g)
1234 do j=js,je ;
do i=is,ie
1235 zos(i,j) = zos(i,j) - g%mask2dT(i,j)*zos_area_mean
1237 if (ids%id_zos > 0)
call post_data(ids%id_zos, zos, diag, mask=g%mask2dT)
1238 if (ids%id_zossq > 0)
then
1239 do j=js,je ;
do i=is,ie
1240 work_2d(i,j) = zos(i,j)*zos(i,j)
1242 call post_data(ids%id_zossq, work_2d, diag, mask=g%mask2dT)
1247 if (ids%id_volo > 0)
then
1248 do j=js,je ;
do i=is,ie
1249 work_2d(i,j) = g%mask2dT(i,j)*(ssh(i,j) + us%Z_to_m*g%bathyT(i,j))
1251 volo = global_area_integral(work_2d, g)
1256 i_time_int = 0.0 ;
if (dt_int > 0.0) i_time_int = 1.0 / dt_int
1259 if (
associated(tv%frazil) .and. (ids%id_fraz > 0))
then
1260 do j=js,je ;
do i=is,ie
1261 work_2d(i,j) = tv%frazil(i,j) * i_time_int
1263 call post_data(ids%id_fraz, work_2d, diag, mask=g%mask2dT)
1267 if (
associated(tv%salt_deficit) .and. (ids%id_salt_deficit > 0))
then
1268 do j=js,je ;
do i=is,ie
1269 work_2d(i,j) = tv%salt_deficit(i,j) * i_time_int
1271 call post_data(ids%id_salt_deficit, work_2d, diag, mask=g%mask2dT)
1275 if (
associated(tv%TempxPmE) .and. (ids%id_Heat_PmE > 0))
then
1276 do j=js,je ;
do i=is,ie
1277 work_2d(i,j) = tv%TempxPmE(i,j) * (tv%C_p * i_time_int)
1279 call post_data(ids%id_Heat_PmE, work_2d, diag, mask=g%mask2dT)
1283 if (
associated(tv%internal_heat) .and. (ids%id_intern_heat > 0))
then
1284 do j=js,je ;
do i=is,ie
1285 work_2d(i,j) = tv%internal_heat(i,j) * (tv%C_p * i_time_int)
1287 call post_data(ids%id_intern_heat, work_2d, diag, mask=g%mask2dT)
1290 if (tv%T_is_conT)
then
1292 if (ids%id_sstcon > 0)
call post_data(ids%id_sstcon, sfc_state%SST, diag, mask=g%mask2dT)
1295 do j=js,je ;
do i=is,ie
1296 work_2d(i,j) = gsw_pt_from_ct(sfc_state%SSS(i,j),sfc_state%SST(i,j))
1298 if (ids%id_sst > 0)
call post_data(ids%id_sst, work_2d, diag, mask=g%mask2dT)
1301 if (ids%id_sst > 0)
call post_data(ids%id_sst, sfc_state%SST, diag, mask=g%mask2dT)
1304 if (tv%S_is_absS)
then
1306 if (ids%id_sssabs > 0)
call post_data(ids%id_sssabs, sfc_state%SSS, diag, mask=g%mask2dT)
1309 do j=js,je ;
do i=is,ie
1310 work_2d(i,j) = gsw_sp_from_sr(sfc_state%SSS(i,j))
1312 if (ids%id_sss > 0)
call post_data(ids%id_sss, work_2d, diag, mask=g%mask2dT)
1315 if (ids%id_sss > 0)
call post_data(ids%id_sss, sfc_state%SSS, diag, mask=g%mask2dT)
1318 if (ids%id_sst_sq > 0)
then
1319 do j=js,je ;
do i=is,ie
1320 work_2d(i,j) = sfc_state%SST(i,j)*sfc_state%SST(i,j)
1322 call post_data(ids%id_sst_sq, work_2d, diag, mask=g%mask2dT)
1324 if (ids%id_sss_sq > 0)
then
1325 do j=js,je ;
do i=is,ie
1326 work_2d(i,j) = sfc_state%SSS(i,j)*sfc_state%SSS(i,j)
1328 call post_data(ids%id_sss_sq, work_2d, diag, mask=g%mask2dT)
1331 call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag))
1333 end subroutine post_surface_thermo_diags
1338 subroutine post_transport_diagnostics(G, GV, uhtr, vhtr, h, IDs, diag_pre_dyn, &
1339 diag, dt_trans, Reg)
1342 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: uhtr
1344 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: vhtr
1346 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1351 real,
intent(in) :: dt_trans
1355 real,
dimension(SZIB_(G), SZJ_(G)) :: umo2d
1356 real,
dimension(SZI_(G), SZJB_(G)) :: vmo2d
1357 real,
dimension(SZIB_(G), SZJ_(G), SZK_(G)) :: umo
1358 real,
dimension(SZI_(G), SZJB_(G), SZK_(G)) :: vmo
1359 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tend
1362 real :: h_to_kg_m2_dt
1364 integer :: i, j, k, is, ie, js, je, nz
1365 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1368 h_to_kg_m2_dt = gv%H_to_kg_m2 * idt
1370 call diag_save_grids(diag)
1371 call diag_copy_storage_to_diag(diag, diag_pre_dyn)
1373 if (ids%id_umo_2d > 0)
then
1375 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
1376 umo2d(i,j) = umo2d(i,j) + uhtr(i,j,k) * h_to_kg_m2_dt
1377 enddo ;
enddo ;
enddo
1378 call post_data(ids%id_umo_2d, umo2d, diag)
1380 if (ids%id_umo > 0)
then
1382 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
1383 umo(i,j,k) = uhtr(i,j,k) * h_to_kg_m2_dt
1384 enddo ;
enddo ;
enddo
1385 call post_data(ids%id_umo, umo, diag, alt_h = diag_pre_dyn%h_state)
1387 if (ids%id_vmo_2d > 0)
then
1389 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
1390 vmo2d(i,j) = vmo2d(i,j) + vhtr(i,j,k) * h_to_kg_m2_dt
1391 enddo ;
enddo ;
enddo
1392 call post_data(ids%id_vmo_2d, vmo2d, diag)
1394 if (ids%id_vmo > 0)
then
1396 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
1397 vmo(i,j,k) = vhtr(i,j,k) * h_to_kg_m2_dt
1398 enddo ;
enddo ;
enddo
1399 call post_data(ids%id_vmo, vmo, diag, alt_h = diag_pre_dyn%h_state)
1402 if (ids%id_uhtr > 0)
call post_data(ids%id_uhtr, uhtr, diag, alt_h = diag_pre_dyn%h_state)
1403 if (ids%id_vhtr > 0)
call post_data(ids%id_vhtr, vhtr, diag, alt_h = diag_pre_dyn%h_state)
1404 if (ids%id_dynamics_h > 0)
call post_data(ids%id_dynamics_h, diag_pre_dyn%h_state, diag, &
1405 alt_h = diag_pre_dyn%h_state)
1407 if (ids%id_dynamics_h_tendency > 0)
then
1409 do k=1,nz ;
do j=js,je ;
do i=is,ie
1410 h_tend(i,j,k) = (h(i,j,k) - diag_pre_dyn%h_state(i,j,k))*idt
1411 enddo ;
enddo ;
enddo
1412 call post_data(ids%id_dynamics_h_tendency, h_tend, diag, alt_h = diag_pre_dyn%h_state)
1415 call post_tracer_transport_diagnostics(g, gv, reg, diag_pre_dyn%h_state, diag)
1417 call diag_restore_grids(diag)
1419 end subroutine post_transport_diagnostics
1423 subroutine mom_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag, CS, tv)
1431 type(time_type),
intent(in) :: time
1437 type(
diag_ctrl),
target,
intent(inout) :: diag
1444 real :: omega, f2_min, convert_h
1446 # include "version_variable.h"
1447 character(len=40) :: mdl =
"MOM_diagnostics"
1448 character(len=48) :: thickness_units, flux_units
1449 logical :: use_temperature, adiabatic
1450 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nkml, nkbl
1451 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j
1453 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1454 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1455 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1456 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1458 if (
associated(cs))
then
1459 call mom_error(warning,
"MOM_diagnostics_init called with an associated "// &
1460 "control structure.")
1466 use_temperature =
associated(tv%T)
1467 call get_param(param_file, mdl,
"ADIABATIC", adiabatic, default=.false., &
1472 call get_param(param_file, mdl,
"DIAG_EBT_MONO_N2_COLUMN_FRACTION", cs%mono_N2_column_fraction, &
1473 "The lower fraction of water column over which N2 is limited as monotonic "// &
1474 "for the purposes of calculating the equivalent barotropic wave speed.", &
1475 units=
'nondim', default=0.)
1476 call get_param(param_file, mdl,
"DIAG_EBT_MONO_N2_DEPTH", cs%mono_N2_depth, &
1477 "The depth below which N2 is limited as monotonic for the "// &
1478 "purposes of calculating the equivalent barotropic wave speed.", &
1479 units=
'm', default=-1.)
1481 if (gv%Boussinesq)
then
1482 thickness_units =
"m" ; flux_units =
"m3 s-1" ; convert_h = gv%H_to_m
1484 thickness_units =
"kg m-2" ; flux_units =
"kg s-1" ; convert_h = gv%H_to_kg_m2
1487 cs%id_masscello = register_diag_field(
'ocean_model',
'masscello', diag%axesTL,&
1488 time,
'Mass per unit area of liquid ocean grid cell',
'kg m-2', &
1489 standard_name=
'sea_water_mass_per_unit_area', v_extensive=.true.)
1491 cs%id_masso = register_scalar_field(
'ocean_model',
'masso', time, &
1492 diag,
'Mass of liquid ocean',
'kg', standard_name=
'sea_water_mass')
1494 cs%id_thkcello = register_diag_field(
'ocean_model',
'thkcello', diag%axesTL, time, &
1495 long_name =
'Cell Thickness', standard_name=
'cell_thickness', units=
'm', v_extensive=.true.)
1496 cs%id_h_pre_sync = register_diag_field(
'ocean_model',
'h_pre_sync', diag%axesTL, time, &
1497 long_name =
'Cell thickness from the previous timestep', units=
'm', v_extensive=.true.)
1502 cs%id_volcello = diag_get_volume_cell_measure_dm_id(diag)
1504 if (use_temperature)
then
1505 if (tv%T_is_conT)
then
1506 cs%id_Tpot = register_diag_field(
'ocean_model',
'temp', diag%axesTL, &
1507 time,
'Potential Temperature',
'degC')
1509 if (tv%S_is_absS)
then
1510 cs%id_Sprac = register_diag_field(
'ocean_model',
'salt', diag%axesTL, &
1511 time,
'Salinity',
'psu')
1514 cs%id_tob = register_diag_field(
'ocean_model',
'tob', diag%axesT1, time, &
1515 long_name=
'Sea Water Potential Temperature at Sea Floor', &
1516 standard_name=
'sea_water_potential_temperature_at_sea_floor', units=
'degC')
1517 cs%id_sob = register_diag_field(
'ocean_model',
'sob',diag%axesT1, time, &
1518 long_name=
'Sea Water Salinity at Sea Floor', &
1519 standard_name=
'sea_water_salinity_at_sea_floor', units=
'psu')
1521 cs%id_temp_layer_ave = register_diag_field(
'ocean_model',
'temp_layer_ave', &
1522 diag%axesZL, time,
'Layer Average Ocean Temperature',
'degC')
1523 cs%id_salt_layer_ave = register_diag_field(
'ocean_model',
'salt_layer_ave', &
1524 diag%axesZL, time,
'Layer Average Ocean Salinity',
'psu')
1526 cs%id_thetaoga = register_scalar_field(
'ocean_model',
'thetaoga', &
1527 time, diag,
'Global Mean Ocean Potential Temperature',
'degC',&
1528 standard_name=
'sea_water_potential_temperature')
1529 cs%id_soga = register_scalar_field(
'ocean_model',
'soga', &
1530 time, diag,
'Global Mean Ocean Salinity',
'psu', &
1531 standard_name=
'sea_water_salinity')
1533 cs%id_tosga = register_scalar_field(
'ocean_model',
'sst_global', time, diag,&
1534 long_name=
'Global Area Average Sea Surface Temperature', &
1535 units=
'degC', standard_name=
'sea_surface_temperature', &
1536 cmor_field_name=
'tosga', cmor_standard_name=
'sea_surface_temperature', &
1537 cmor_long_name=
'Sea Surface Temperature')
1538 cs%id_sosga = register_scalar_field(
'ocean_model',
'sss_global', time, diag,&
1539 long_name=
'Global Area Average Sea Surface Salinity', &
1540 units=
'psu', standard_name=
'sea_surface_salinity', &
1541 cmor_field_name=
'sosga', cmor_standard_name=
'sea_surface_salinity', &
1542 cmor_long_name=
'Sea Surface Salinity')
1545 cs%id_u = register_diag_field(
'ocean_model',
'u', diag%axesCuL, time, &
1546 'Zonal velocity',
'm s-1', cmor_field_name=
'uo', &
1547 cmor_standard_name=
'sea_water_x_velocity', cmor_long_name=
'Sea Water X Velocity')
1548 cs%id_v = register_diag_field(
'ocean_model',
'v', diag%axesCvL, time, &
1549 'Meridional velocity',
'm s-1', cmor_field_name=
'vo', &
1550 cmor_standard_name=
'sea_water_y_velocity', cmor_long_name=
'Sea Water Y Velocity')
1551 cs%id_h = register_diag_field(
'ocean_model',
'h', diag%axesTL, time, &
1552 'Layer Thickness', thickness_units, v_extensive=.true., conversion=convert_h)
1554 cs%id_e = register_diag_field(
'ocean_model',
'e', diag%axesTi, time, &
1555 'Interface Height Relative to Mean Sea Level',
'm', conversion=us%Z_to_m)
1556 if (cs%id_e>0)
call safe_alloc_ptr(cs%e,isd,ied,jsd,jed,nz+1)
1558 cs%id_e_D = register_diag_field(
'ocean_model',
'e_D', diag%axesTi, time, &
1559 'Interface Height above the Seafloor',
'm', conversion=us%Z_to_m)
1560 if (cs%id_e_D>0)
call safe_alloc_ptr(cs%e_D,isd,ied,jsd,jed,nz+1)
1562 cs%id_Rml = register_diag_field(
'ocean_model',
'Rml', diag%axesTL, time, &
1563 'Mixed Layer Coordinate Potential Density',
'kg m-3')
1565 cs%id_Rcv = register_diag_field(
'ocean_model',
'Rho_cv', diag%axesTL, time, &
1566 'Coordinate Potential Density',
'kg m-3')
1568 cs%id_rhopot0 = register_diag_field(
'ocean_model',
'rhopot0', diag%axesTL, time, &
1569 'Potential density referenced to surface',
'kg m-3')
1570 cs%id_rhopot2 = register_diag_field(
'ocean_model',
'rhopot2', diag%axesTL, time, &
1571 'Potential density referenced to 2000 dbar',
'kg m-3')
1572 cs%id_rhoinsitu = register_diag_field(
'ocean_model',
'rhoinsitu', diag%axesTL, time, &
1573 'In situ density',
'kg m-3')
1575 cs%id_du_dt = register_diag_field(
'ocean_model',
'dudt', diag%axesCuL, time, &
1576 'Zonal Acceleration',
'm s-2')
1577 if ((cs%id_du_dt>0) .and. .not.
associated(cs%du_dt))
then
1578 call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
1579 call register_time_deriv(lbound(mis%u), mis%u, cs%du_dt, cs)
1582 cs%id_dv_dt = register_diag_field(
'ocean_model',
'dvdt', diag%axesCvL, time, &
1583 'Meridional Acceleration',
'm s-2')
1584 if ((cs%id_dv_dt>0) .and. .not.
associated(cs%dv_dt))
then
1585 call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
1586 call register_time_deriv(lbound(mis%v), mis%v, cs%dv_dt, cs)
1589 cs%id_dh_dt = register_diag_field(
'ocean_model',
'dhdt', diag%axesTL, time, &
1590 'Thickness tendency', trim(thickness_units)//
" s-1", v_extensive = .true.)
1591 if ((cs%id_dh_dt>0) .and. .not.
associated(cs%dh_dt))
then
1592 call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
1593 call register_time_deriv(lbound(mis%h), mis%h, cs%dh_dt, cs)
1598 cs%id_h_Rlay = register_diag_field(
'ocean_model',
'h_rho', diag%axesTL, time, &
1599 'Layer thicknesses in pure potential density coordinates', thickness_units, &
1600 conversion=convert_h)
1601 if (cs%id_h_Rlay>0)
call safe_alloc_ptr(cs%h_Rlay,isd,ied,jsd,jed,nz)
1603 cs%id_uh_Rlay = register_diag_field(
'ocean_model',
'uh_rho', diag%axesCuL, time, &
1604 'Zonal volume transport in pure potential density coordinates', flux_units, &
1605 conversion=convert_h)
1606 if (cs%id_uh_Rlay>0)
call safe_alloc_ptr(cs%uh_Rlay,isdb,iedb,jsd,jed,nz)
1608 cs%id_vh_Rlay = register_diag_field(
'ocean_model',
'vh_rho', diag%axesCvL, time, &
1609 'Meridional volume transport in pure potential density coordinates', flux_units, &
1610 conversion=convert_h)
1611 if (cs%id_vh_Rlay>0)
call safe_alloc_ptr(cs%vh_Rlay,isd,ied,jsdb,jedb,nz)
1613 cs%id_uhGM_Rlay = register_diag_field(
'ocean_model',
'uhGM_rho', diag%axesCuL, time, &
1614 'Zonal volume transport due to interface height diffusion in pure potential &
1615 &density coordinates', flux_units, conversion=convert_h)
1616 if (cs%id_uhGM_Rlay>0)
call safe_alloc_ptr(cs%uhGM_Rlay,isdb,iedb,jsd,jed,nz)
1618 cs%id_vhGM_Rlay = register_diag_field(
'ocean_model',
'vhGM_rho', diag%axesCvL, time, &
1619 'Meridional volume transport due to interface height diffusion in pure &
1620 &potential density coordinates', flux_units, conversion=convert_h)
1621 if (cs%id_vhGM_Rlay>0)
call safe_alloc_ptr(cs%vhGM_Rlay,isd,ied,jsdb,jedb,nz)
1626 cs%id_KE = register_diag_field(
'ocean_model',
'KE', diag%axesTL, time, &
1627 'Layer kinetic energy per unit mass',
'm2 s-2')
1628 if (cs%id_KE>0)
call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)
1630 cs%id_dKEdt = register_diag_field(
'ocean_model',
'dKE_dt', diag%axesTL, time, &
1631 'Kinetic Energy Tendency of Layer',
'm3 s-3')
1632 if (cs%id_dKEdt>0)
call safe_alloc_ptr(cs%dKE_dt,isd,ied,jsd,jed,nz)
1634 cs%id_PE_to_KE = register_diag_field(
'ocean_model',
'PE_to_KE', diag%axesTL, time, &
1635 'Potential to Kinetic Energy Conversion of Layer',
'm3 s-3')
1636 if (cs%id_PE_to_KE>0)
call safe_alloc_ptr(cs%PE_to_KE,isd,ied,jsd,jed,nz)
1638 cs%id_KE_Coradv = register_diag_field(
'ocean_model',
'KE_Coradv', diag%axesTL, time, &
1639 'Kinetic Energy Source from Coriolis and Advection',
'm3 s-3')
1640 if (cs%id_KE_Coradv>0)
call safe_alloc_ptr(cs%KE_Coradv,isd,ied,jsd,jed,nz)
1642 cs%id_KE_adv = register_diag_field(
'ocean_model',
'KE_adv', diag%axesTL, time, &
1643 'Kinetic Energy Source from Advection',
'm3 s-3')
1644 if (cs%id_KE_adv>0)
call safe_alloc_ptr(cs%KE_adv,isd,ied,jsd,jed,nz)
1646 cs%id_KE_visc = register_diag_field(
'ocean_model',
'KE_visc', diag%axesTL, time, &
1647 'Kinetic Energy Source from Vertical Viscosity and Stresses',
'm3 s-3')
1648 if (cs%id_KE_visc>0)
call safe_alloc_ptr(cs%KE_visc,isd,ied,jsd,jed,nz)
1650 cs%id_KE_horvisc = register_diag_field(
'ocean_model',
'KE_horvisc', diag%axesTL, time, &
1651 'Kinetic Energy Source from Horizontal Viscosity',
'm3 s-3')
1652 if (cs%id_KE_horvisc>0)
call safe_alloc_ptr(cs%KE_horvisc,isd,ied,jsd,jed,nz)
1654 if (.not. adiabatic)
then
1655 cs%id_KE_dia = register_diag_field(
'ocean_model',
'KE_dia', diag%axesTL, time, &
1656 'Kinetic Energy Source from Diapycnal Diffusion',
'm3 s-3')
1657 if (cs%id_KE_dia>0)
call safe_alloc_ptr(cs%KE_dia,isd,ied,jsd,jed,nz)
1661 cs%id_cg1 = register_diag_field(
'ocean_model',
'cg1', diag%axesT1, time, &
1662 'First baroclinic gravity wave speed',
'm s-1')
1663 cs%id_Rd1 = register_diag_field(
'ocean_model',
'Rd1', diag%axesT1, time, &
1664 'First baroclinic deformation radius',
'm')
1665 cs%id_cfl_cg1 = register_diag_field(
'ocean_model',
'CFL_cg1', diag%axesT1, time, &
1666 'CFL of first baroclinic gravity wave = dt*cg1*(1/dx+1/dy)',
'nondim')
1667 cs%id_cfl_cg1_x = register_diag_field(
'ocean_model',
'CFL_cg1_x', diag%axesT1, time, &
1668 'i-component of CFL of first baroclinic gravity wave = dt*cg1*/dx',
'nondim')
1669 cs%id_cfl_cg1_y = register_diag_field(
'ocean_model',
'CFL_cg1_y', diag%axesT1, time, &
1670 'j-component of CFL of first baroclinic gravity wave = dt*cg1*/dy',
'nondim')
1671 cs%id_cg_ebt = register_diag_field(
'ocean_model',
'cg_ebt', diag%axesT1, time, &
1672 'Equivalent barotropic gravity wave speed',
'm s-1')
1673 cs%id_Rd_ebt = register_diag_field(
'ocean_model',
'Rd_ebt', diag%axesT1, time, &
1674 'Equivalent barotropic deformation radius',
'm')
1675 cs%id_p_ebt = register_diag_field(
'ocean_model',
'p_ebt', diag%axesTL, time, &
1676 'Equivalent barotropic modal strcuture',
'nondim')
1678 if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
1679 (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0) .or. &
1680 (cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0))
then
1681 call wave_speed_init(cs%wave_speed_CSp)
1682 call safe_alloc_ptr(cs%cg1,isd,ied,jsd,jed)
1683 if (cs%id_Rd1>0)
call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
1684 if (cs%id_Rd_ebt>0)
call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
1685 if (cs%id_cfl_cg1>0)
call safe_alloc_ptr(cs%cfl_cg1,isd,ied,jsd,jed)
1686 if (cs%id_cfl_cg1_x>0)
call safe_alloc_ptr(cs%cfl_cg1_x,isd,ied,jsd,jed)
1687 if (cs%id_cfl_cg1_y>0)
call safe_alloc_ptr(cs%cfl_cg1_y,isd,ied,jsd,jed)
1688 if (cs%id_p_ebt>0)
call safe_alloc_ptr(cs%p_ebt,isd,ied,jsd,jed,nz)
1691 cs%id_mass_wt = register_diag_field(
'ocean_model',
'mass_wt', diag%axesT1, time, &
1692 'The column mass for calculating mass-weighted average properties',
'kg m-2')
1694 if (use_temperature)
then
1695 cs%id_temp_int = register_diag_field(
'ocean_model',
'temp_int', diag%axesT1, time, &
1696 'Density weighted column integrated potential temperature',
'degC kg m-2', &
1697 cmor_field_name=
'opottempmint', &
1698 cmor_long_name=
'integral_wrt_depth_of_product_of_sea_water_density_and_potential_temperature',&
1699 cmor_standard_name=
'Depth integrated density times potential temperature')
1701 cs%id_salt_int = register_diag_field(
'ocean_model',
'salt_int', diag%axesT1, time, &
1702 'Density weighted column integrated salinity',
'psu kg m-2', &
1703 cmor_field_name=
'somint', &
1704 cmor_long_name=
'integral_wrt_depth_of_product_of_sea_water_density_and_salinity',&
1705 cmor_standard_name=
'Depth integrated density times salinity')
1708 cs%id_col_mass = register_diag_field(
'ocean_model',
'col_mass', diag%axesT1, time, &
1709 'The column integrated in situ density',
'kg m-2')
1711 cs%id_col_ht = register_diag_field(
'ocean_model',
'col_height', diag%axesT1, time, &
1712 'The height of the water column',
'm', conversion=us%Z_to_m)
1713 cs%id_pbo = register_diag_field(
'ocean_model',
'pbo', diag%axesT1, time, &
1714 long_name=
'Sea Water Pressure at Sea Floor', standard_name=
'sea_water_pressure_at_sea_floor', &
1717 call set_dependent_diagnostics(mis, adp, cdp, g, cs)
1719 end subroutine mom_diagnostics_init
1723 subroutine register_surface_diags(Time, G, IDs, diag, tv)
1724 type(time_type),
intent(in) :: time
1731 ids%id_volo = register_scalar_field(
'ocean_model',
'volo', time, diag,&
1732 long_name=
'Total volume of liquid ocean', units=
'm3', &
1733 standard_name=
'sea_water_volume')
1734 ids%id_zos = register_diag_field(
'ocean_model',
'zos', diag%axesT1, time,&
1735 standard_name =
'sea_surface_height_above_geoid', &
1736 long_name=
'Sea surface height above geoid', units=
'm')
1737 ids%id_zossq = register_diag_field(
'ocean_model',
'zossq', diag%axesT1, time,&
1738 standard_name=
'square_of_sea_surface_height_above_geoid', &
1739 long_name=
'Square of sea surface height above geoid', units=
'm2')
1740 ids%id_ssh = register_diag_field(
'ocean_model',
'SSH', diag%axesT1, time, &
1741 'Sea Surface Height',
'm')
1742 ids%id_ssh_ga = register_scalar_field(
'ocean_model',
'ssh_ga', time, diag,&
1743 long_name=
'Area averaged sea surface height', units=
'm', &
1744 standard_name=
'area_averaged_sea_surface_height')
1745 ids%id_ssu = register_diag_field(
'ocean_model',
'SSU', diag%axesCu1, time, &
1746 'Sea Surface Zonal Velocity',
'm s-1')
1747 ids%id_ssv = register_diag_field(
'ocean_model',
'SSV', diag%axesCv1, time, &
1748 'Sea Surface Meridional Velocity',
'm s-1')
1749 ids%id_speed = register_diag_field(
'ocean_model',
'speed', diag%axesT1, time, &
1750 'Sea Surface Speed',
'm s-1')
1752 if (
associated(tv%T))
then
1753 ids%id_sst = register_diag_field(
'ocean_model',
'SST', diag%axesT1, time, &
1754 'Sea Surface Temperature',
'degC', cmor_field_name=
'tos', &
1755 cmor_long_name=
'Sea Surface Temperature', &
1756 cmor_standard_name=
'sea_surface_temperature')
1757 ids%id_sst_sq = register_diag_field(
'ocean_model',
'SST_sq', diag%axesT1, time, &
1758 'Sea Surface Temperature Squared',
'degC2', cmor_field_name=
'tossq', &
1759 cmor_long_name=
'Square of Sea Surface Temperature ', &
1760 cmor_standard_name=
'square_of_sea_surface_temperature')
1761 ids%id_sss = register_diag_field(
'ocean_model',
'SSS', diag%axesT1, time, &
1762 'Sea Surface Salinity',
'psu', cmor_field_name=
'sos', &
1763 cmor_long_name=
'Sea Surface Salinity', &
1764 cmor_standard_name=
'sea_surface_salinity')
1765 ids%id_sss_sq = register_diag_field(
'ocean_model',
'SSS_sq', diag%axesT1, time, &
1766 'Sea Surface Salinity Squared',
'psu', cmor_field_name=
'sossq', &
1767 cmor_long_name=
'Square of Sea Surface Salinity ', &
1768 cmor_standard_name=
'square_of_sea_surface_salinity')
1769 if (tv%T_is_conT)
then
1770 ids%id_sstcon = register_diag_field(
'ocean_model',
'conSST', diag%axesT1, time, &
1771 'Sea Surface Conservative Temperature',
'Celsius')
1773 if (tv%S_is_absS)
then
1774 ids%id_sssabs = register_diag_field(
'ocean_model',
'absSSS', diag%axesT1, time, &
1775 'Sea Surface Absolute Salinity',
'g kg-1')
1777 if (
associated(tv%frazil))
then
1778 ids%id_fraz = register_diag_field(
'ocean_model',
'frazil', diag%axesT1, time, &
1779 'Heat from frazil formation',
'W m-2', cmor_field_name=
'hfsifrazil', &
1780 cmor_standard_name=
'heat_flux_into_sea_water_due_to_frazil_ice_formation', &
1781 cmor_long_name=
'Heat Flux into Sea Water due to Frazil Ice Formation')
1785 ids%id_salt_deficit = register_diag_field(
'ocean_model',
'salt_deficit', diag%axesT1, time, &
1786 'Salt sink in ocean due to ice flux',
'psu m-2 s-1')
1787 ids%id_Heat_PmE = register_diag_field(
'ocean_model',
'Heat_PmE', diag%axesT1, time, &
1788 'Heat flux into ocean from mass flux into ocean',
'W m-2')
1789 ids%id_intern_heat = register_diag_field(
'ocean_model',
'internal_heat', diag%axesT1, time,&
1790 'Heat flux into ocean from geothermal or other internal sources',
'W m-2')
1792 end subroutine register_surface_diags
1795 subroutine register_transport_diags(Time, G, GV, IDs, diag)
1796 type(time_type),
intent(in) :: time
1803 character(len=48) :: thickness_units
1805 thickness_units = get_thickness_units(gv)
1806 if (gv%Boussinesq)
then
1807 h_convert = gv%H_to_m
1809 h_convert = gv%H_to_kg_m2
1813 ids%id_uhtr = register_diag_field(
'ocean_model',
'uhtr', diag%axesCuL, time, &
1814 'Accumulated zonal thickness fluxes to advect tracers',
'kg', &
1815 y_cell_method=
'sum', v_extensive=.true., conversion=h_convert)
1816 ids%id_vhtr = register_diag_field(
'ocean_model',
'vhtr', diag%axesCvL, time, &
1817 'Accumulated meridional thickness fluxes to advect tracers',
'kg', &
1818 x_cell_method=
'sum', v_extensive=.true., conversion=h_convert)
1819 ids%id_umo = register_diag_field(
'ocean_model',
'umo', &
1820 diag%axesCuL, time,
'Ocean Mass X Transport',
'kg s-1', &
1821 standard_name=
'ocean_mass_x_transport', y_cell_method=
'sum', v_extensive=.true.)
1822 ids%id_vmo = register_diag_field(
'ocean_model',
'vmo', &
1823 diag%axesCvL, time,
'Ocean Mass Y Transport',
'kg s-1', &
1824 standard_name=
'ocean_mass_y_transport', x_cell_method=
'sum', v_extensive=.true.)
1825 ids%id_umo_2d = register_diag_field(
'ocean_model',
'umo_2d', &
1826 diag%axesCu1, time,
'Ocean Mass X Transport Vertical Sum',
'kg s-1', &
1827 standard_name=
'ocean_mass_x_transport_vertical_sum', y_cell_method=
'sum')
1828 ids%id_vmo_2d = register_diag_field(
'ocean_model',
'vmo_2d', &
1829 diag%axesCv1, time,
'Ocean Mass Y Transport Vertical Sum',
'kg s-1', &
1830 standard_name=
'ocean_mass_y_transport_vertical_sum', x_cell_method=
'sum')
1831 ids%id_dynamics_h = register_diag_field(
'ocean_model',
'dynamics_h', &
1832 diag%axesTl, time,
'Change in layer thicknesses due to horizontal dynamics', &
1833 'm s-1', v_extensive = .true.)
1834 ids%id_dynamics_h_tendency = register_diag_field(
'ocean_model',
'dynamics_h_tendency', &
1835 diag%axesTl, time,
'Change in layer thicknesses due to horizontal dynamics', &
1836 'm s-1', v_extensive = .true.)
1838 end subroutine register_transport_diags
1841 subroutine write_static_fields(G, GV, US, tv, diag)
1846 type(
diag_ctrl),
target,
intent(inout) :: diag
1851 id = register_static_field(
'ocean_model',
'geolat', diag%axesT1, &
1852 'Latitude of tracer (T) points',
'degrees_north')
1853 if (id > 0)
call post_data(id, g%geoLatT, diag, .true.)
1855 id = register_static_field(
'ocean_model',
'geolon', diag%axesT1, &
1856 'Longitude of tracer (T) points',
'degrees_east')
1857 if (id > 0)
call post_data(id, g%geoLonT, diag, .true.)
1859 id = register_static_field(
'ocean_model',
'geolat_c', diag%axesB1, &
1860 'Latitude of corner (Bu) points',
'degrees_north', interp_method=
'none')
1861 if (id > 0)
call post_data(id, g%geoLatBu, diag, .true.)
1863 id = register_static_field(
'ocean_model',
'geolon_c', diag%axesB1, &
1864 'Longitude of corner (Bu) points',
'degrees_east', interp_method=
'none')
1865 if (id > 0)
call post_data(id, g%geoLonBu, diag, .true.)
1867 id = register_static_field(
'ocean_model',
'geolat_v', diag%axesCv1, &
1868 'Latitude of meridional velocity (Cv) points',
'degrees_north', interp_method=
'none')
1869 if (id > 0)
call post_data(id, g%geoLatCv, diag, .true.)
1871 id = register_static_field(
'ocean_model',
'geolon_v', diag%axesCv1, &
1872 'Longitude of meridional velocity (Cv) points',
'degrees_east', interp_method=
'none')
1873 if (id > 0)
call post_data(id, g%geoLonCv, diag, .true.)
1875 id = register_static_field(
'ocean_model',
'geolat_u', diag%axesCu1, &
1876 'Latitude of zonal velocity (Cu) points',
'degrees_north', interp_method=
'none')
1877 if (id > 0)
call post_data(id, g%geoLatCu, diag, .true.)
1879 id = register_static_field(
'ocean_model',
'geolon_u', diag%axesCu1, &
1880 'Longitude of zonal velocity (Cu) points',
'degrees_east', interp_method=
'none')
1881 if (id > 0)
call post_data(id, g%geoLonCu, diag, .true.)
1883 id = register_static_field(
'ocean_model',
'area_t', diag%axesT1, &
1884 'Surface area of tracer (T) cells',
'm2', &
1885 cmor_field_name=
'areacello', cmor_standard_name=
'cell_area', &
1886 cmor_long_name=
'Ocean Grid-Cell Area', &
1887 x_cell_method=
'sum', y_cell_method=
'sum', area_cell_method=
'sum')
1889 call post_data(id, g%areaT, diag, .true.)
1890 call diag_register_area_ids(diag, id_area_t=id)
1893 id = register_static_field(
'ocean_model',
'area_u', diag%axesCu1, &
1894 'Surface area of x-direction flow (U) cells',
'm2', &
1895 cmor_field_name=
'areacello_cu', cmor_standard_name=
'cell_area', &
1896 cmor_long_name=
'Ocean Grid-Cell Area', &
1897 x_cell_method=
'sum', y_cell_method=
'sum', area_cell_method=
'sum')
1898 if (id > 0)
call post_data(id, g%areaCu, diag, .true.)
1900 id = register_static_field(
'ocean_model',
'area_v', diag%axesCv1, &
1901 'Surface area of y-direction flow (V) cells',
'm2', &
1902 cmor_field_name=
'areacello_cv', cmor_standard_name=
'cell_area', &
1903 cmor_long_name=
'Ocean Grid-Cell Area', &
1904 x_cell_method=
'sum', y_cell_method=
'sum', area_cell_method=
'sum')
1905 if (id > 0)
call post_data(id, g%areaCv, diag, .true.)
1907 id = register_static_field(
'ocean_model',
'area_q', diag%axesB1, &
1908 'Surface area of B-grid flow (Q) cells',
'm2', &
1909 cmor_field_name=
'areacello_bu', cmor_standard_name=
'cell_area', &
1910 cmor_long_name=
'Ocean Grid-Cell Area', &
1911 x_cell_method=
'sum', y_cell_method=
'sum', area_cell_method=
'sum')
1912 if (id > 0)
call post_data(id, g%areaBu, diag, .true.)
1914 id = register_static_field(
'ocean_model',
'depth_ocean', diag%axesT1, &
1915 'Depth of the ocean at tracer points',
'm', &
1916 standard_name=
'sea_floor_depth_below_geoid', &
1917 cmor_field_name=
'deptho', cmor_long_name=
'Sea Floor Depth', &
1918 cmor_standard_name=
'sea_floor_depth_below_geoid', area=diag%axesT1%id_area, &
1919 x_cell_method=
'mean', y_cell_method=
'mean', area_cell_method=
'mean', &
1920 conversion=us%Z_to_m)
1921 if (id > 0)
call post_data(id, g%bathyT, diag, .true., mask=g%mask2dT)
1923 id = register_static_field(
'ocean_model',
'wet', diag%axesT1, &
1924 '0 if land, 1 if ocean at tracer points',
'none', area=diag%axesT1%id_area)
1925 if (id > 0)
call post_data(id, g%mask2dT, diag, .true.)
1927 id = register_static_field(
'ocean_model',
'wet_c', diag%axesB1, &
1928 '0 if land, 1 if ocean at corner (Bu) points',
'none', interp_method=
'none')
1929 if (id > 0)
call post_data(id, g%mask2dBu, diag, .true.)
1931 id = register_static_field(
'ocean_model',
'wet_u', diag%axesCu1, &
1932 '0 if land, 1 if ocean at zonal velocity (Cu) points',
'none', interp_method=
'none')
1933 if (id > 0)
call post_data(id, g%mask2dCu, diag, .true.)
1935 id = register_static_field(
'ocean_model',
'wet_v', diag%axesCv1, &
1936 '0 if land, 1 if ocean at meridional velocity (Cv) points',
'none', interp_method=
'none')
1937 if (id > 0)
call post_data(id, g%mask2dCv, diag, .true.)
1939 id = register_static_field(
'ocean_model',
'Coriolis', diag%axesB1, &
1940 'Coriolis parameter at corner (Bu) points',
's-1', interp_method=
'none', conversion=us%s_to_T)
1941 if (id > 0)
call post_data(id, g%CoriolisBu, diag, .true.)
1943 id = register_static_field(
'ocean_model',
'dxt', diag%axesT1, &
1944 'Delta(x) at thickness/tracer points (meter)',
'm', interp_method=
'none')
1945 if (id > 0)
call post_data(id, g%dxt, diag, .true.)
1947 id = register_static_field(
'ocean_model',
'dyt', diag%axesT1, &
1948 'Delta(y) at thickness/tracer points (meter)',
'm', interp_method=
'none')
1949 if (id > 0)
call post_data(id, g%dyt, diag, .true.)
1951 id = register_static_field(
'ocean_model',
'dxCu', diag%axesCu1, &
1952 'Delta(x) at u points (meter)',
'm', interp_method=
'none')
1953 if (id > 0)
call post_data(id, g%dxCu, diag, .true.)
1955 id = register_static_field(
'ocean_model',
'dyCu', diag%axesCu1, &
1956 'Delta(y) at u points (meter)',
'm', interp_method=
'none')
1957 if (id > 0)
call post_data(id, g%dyCu, diag, .true.)
1959 id = register_static_field(
'ocean_model',
'dxCv', diag%axesCv1, &
1960 'Delta(x) at v points (meter)',
'm', interp_method=
'none')
1961 if (id > 0)
call post_data(id, g%dxCv, diag, .true.)
1963 id = register_static_field(
'ocean_model',
'dyCv', diag%axesCv1, &
1964 'Delta(y) at v points (meter)',
'm', interp_method=
'none')
1965 if (id > 0)
call post_data(id, g%dyCv, diag, .true.)
1967 id = register_static_field(
'ocean_model',
'dyCuo', diag%axesCu1, &
1968 'Open meridional grid spacing at u points (meter)',
'm', interp_method=
'none')
1969 if (id > 0)
call post_data(id, g%dy_Cu, diag, .true.)
1971 id = register_static_field(
'ocean_model',
'dxCvo', diag%axesCv1, &
1972 'Open zonal grid spacing at v points (meter)',
'm', interp_method=
'none')
1973 if (id > 0)
call post_data(id, g%dx_Cv, diag, .true.)
1975 id = register_static_field(
'ocean_model',
'sin_rot', diag%axesT1, &
1976 'sine of the clockwise angle of the ocean grid north to true north',
'none')
1977 if (id > 0)
call post_data(id, g%sin_rot, diag, .true.)
1979 id = register_static_field(
'ocean_model',
'cos_rot', diag%axesT1, &
1980 'cosine of the clockwise angle of the ocean grid north to true north',
'none')
1981 if (id > 0)
call post_data(id, g%cos_rot, diag, .true.)
1986 id = register_static_field(
'ocean_model',
'area_t_percent', diag%axesT1, &
1987 'Percentage of cell area covered by ocean',
'%', &
1988 cmor_field_name=
'sftof', cmor_standard_name=
'SeaAreaFraction', &
1989 cmor_long_name=
'Sea Area Fraction', &
1990 x_cell_method=
'mean', y_cell_method=
'mean', area_cell_method=
'mean', &
1992 if (id > 0)
call post_data(id, g%mask2dT, diag, .true.)
1995 id = register_static_field(
'ocean_model',
'Rho_0', diag%axesNull, &
1996 'mean ocean density used with the Boussinesq approximation', &
1997 'kg m-3', cmor_field_name=
'rhozero', &
1998 cmor_standard_name=
'reference_sea_water_density_for_boussinesq_approximation', &
1999 cmor_long_name=
'reference sea water density for boussinesq approximation')
2000 if (id > 0)
call post_data(id, gv%Rho0, diag, .true.)
2002 id = register_static_field(
'ocean_model',
'C_p', diag%axesNull, &
2003 'heat capacity of sea water',
'J kg-1 K-1', cmor_field_name=
'cpocean', &
2004 cmor_standard_name=
'specific_heat_capacity_of_sea_water', &
2005 cmor_long_name=
'specific_heat_capacity_of_sea_water')
2006 if (id > 0)
call post_data(id, tv%C_p, diag, .true.)
2008 end subroutine write_static_fields
2012 subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, CS)
2025 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
2026 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
2027 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2029 if (
associated(cs%dKE_dt) .or.
associated(cs%PE_to_KE) .or. &
2030 associated(cs%KE_CorAdv) .or.
associated(cs%KE_adv) .or. &
2031 associated(cs%KE_visc) .or.
associated(cs%KE_horvisc) .or. &
2032 associated(cs%KE_dia)) &
2033 call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)
2035 if (
associated(cs%dKE_dt))
then
2036 if (.not.
associated(cs%du_dt))
then
2037 call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
2038 call register_time_deriv(lbound(mis%u), mis%u, cs%du_dt, cs)
2040 if (.not.
associated(cs%dv_dt))
then
2041 call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
2042 call register_time_deriv(lbound(mis%v), mis%v, cs%dv_dt, cs)
2044 if (.not.
associated(cs%dh_dt))
then
2045 call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
2046 call register_time_deriv(lbound(mis%h), mis%h, cs%dh_dt, cs)
2050 if (
associated(cs%KE_adv))
then
2051 call safe_alloc_ptr(adp%gradKEu,isdb,iedb,jsd,jed,nz)
2052 call safe_alloc_ptr(adp%gradKEv,isd,ied,jsdb,jedb,nz)
2055 if (
associated(cs%KE_visc))
then
2056 call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
2057 call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
2060 if (
associated(cs%KE_dia))
then
2061 call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
2062 call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
2065 if (
associated(cs%uhGM_Rlay))
call safe_alloc_ptr(cdp%uhGM,isdb,iedb,jsd,jed,nz)
2066 if (
associated(cs%vhGM_Rlay))
call safe_alloc_ptr(cdp%vhGM,isd,ied,jsdb,jedb,nz)
2068 end subroutine set_dependent_diagnostics
2071 subroutine mom_diagnostics_end(CS, ADp)
2078 if (
associated(cs%e))
deallocate(cs%e)
2079 if (
associated(cs%e_D))
deallocate(cs%e_D)
2080 if (
associated(cs%KE))
deallocate(cs%KE)
2081 if (
associated(cs%dKE_dt))
deallocate(cs%dKE_dt)
2082 if (
associated(cs%PE_to_KE))
deallocate(cs%PE_to_KE)
2083 if (
associated(cs%KE_Coradv))
deallocate(cs%KE_Coradv)
2084 if (
associated(cs%KE_adv))
deallocate(cs%KE_adv)
2085 if (
associated(cs%KE_visc))
deallocate(cs%KE_visc)
2086 if (
associated(cs%KE_horvisc))
deallocate(cs%KE_horvisc)
2087 if (
associated(cs%KE_dia))
deallocate(cs%KE_dia)
2088 if (
associated(cs%dv_dt))
deallocate(cs%dv_dt)
2089 if (
associated(cs%dh_dt))
deallocate(cs%dh_dt)
2090 if (
associated(cs%du_dt))
deallocate(cs%du_dt)
2091 if (
associated(cs%h_Rlay))
deallocate(cs%h_Rlay)
2092 if (
associated(cs%uh_Rlay))
deallocate(cs%uh_Rlay)
2093 if (
associated(cs%vh_Rlay))
deallocate(cs%vh_Rlay)
2094 if (
associated(cs%uhGM_Rlay))
deallocate(cs%uhGM_Rlay)
2095 if (
associated(cs%vhGM_Rlay))
deallocate(cs%vhGM_Rlay)
2097 if (
associated(adp%gradKEu))
deallocate(adp%gradKEu)
2098 if (
associated(adp%gradKEu))
deallocate(adp%gradKEu)
2099 if (
associated(adp%du_dt_visc))
deallocate(adp%du_dt_visc)
2100 if (
associated(adp%dv_dt_visc))
deallocate(adp%dv_dt_visc)
2101 if (
associated(adp%du_dt_dia))
deallocate(adp%du_dt_dia)
2102 if (
associated(adp%dv_dt_dia))
deallocate(adp%dv_dt_dia)
2103 if (
associated(adp%du_other))
deallocate(adp%du_other)
2104 if (
associated(adp%dv_other))
deallocate(adp%dv_other)
2106 do m=1,cs%num_time_deriv ;
deallocate(cs%prev_val(m)%p) ;
enddo
2110 end subroutine mom_diagnostics_end