Diagnostics not more naturally calculated elsewhere are computed here.
187 type(ocean_grid_type),
intent(inout) :: G
188 type(verticalGrid_type),
intent(in) :: GV
189 type(unit_scale_type),
intent(in) :: US
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)), &
202 type(thermo_var_ptrs),
intent(in) :: tv
204 type(accel_diag_ptrs),
intent(in) :: ADp
206 type(cont_diag_ptrs),
intent(in) :: CDp
208 real,
dimension(:,:),
pointer :: p_surf
211 real,
intent(in) :: dt
213 type(diag_grid_storage),
intent(in) :: diag_pre_sync
214 type(diagnostics_CS),
intent(inout) :: CS
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
325 masso = reproducing_sum(work_2d)
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)
362 call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, &
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
468 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, &
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
591 call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
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
600 call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
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
611 call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
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)