This subroutine uses a substantially refactored tridiagonal equation for diapycnal mixing of temperature and salinity to estimate the potential energy change due to diapycnal mixing in a column of water. It does this estimate 4 different ways, all of which should be equivalent, but reports only one. The various estimates are taken because they will later be used as templates for other bits of code.
122 type(ocean_grid_type),
intent(in) :: G
123 type(verticalGrid_type),
intent(in) :: GV
124 type(unit_scale_type),
intent(in) :: US
125 real,
dimension(GV%ke),
intent(in) :: h_in
127 real,
dimension(GV%ke),
intent(in) :: T_in
128 real,
dimension(GV%ke),
intent(in) :: S_in
129 real,
dimension(GV%ke+1),
intent(in) :: Kd
131 real,
intent(in) :: dt
132 real,
intent(out) :: energy_Kd
134 type(thermo_var_ptrs),
intent(inout) :: tv
137 logical,
optional,
intent(in) :: may_print
139 type(diapyc_energy_req_CS), &
140 optional,
pointer :: CS
149 real,
dimension(GV%ke) :: &
197 real,
dimension(GV%ke+1) :: &
211 real,
dimension(GV%ke+1,4) :: &
242 real :: dTe_t2, dT_km1_t2, dT_k_t2
243 real :: dSe_t2, dS_km1_t2, dS_k_t2
247 real,
dimension(GV%ke) :: &
249 real,
dimension(GV%ke+1) :: &
250 dPEa_dKd, dPEa_dKd_est, dPEa_dKd_err, dPEa_dKd_trunc, dPEa_dKd_err_norm, &
251 dPEb_dKd, dPEb_dKd_est, dPEb_dKd_err, dPEb_dKd_trunc, dPEb_dKd_err_norm
252 real :: PE_chg_tot1A, PE_chg_tot2A, T_chg_totA
253 real :: PE_chg_tot1B, PE_chg_tot2B, T_chg_totB
254 real :: PE_chg_tot1C, PE_chg_tot2C, T_chg_totC
255 real :: PE_chg_tot1D, PE_chg_tot2D, T_chg_totD
256 real,
dimension(GV%ke+1) :: dPEchg_dKd
258 real,
dimension(6) :: dT_k_itt, dS_k_itt, dT_km1_itt, dS_km1_itt
260 integer :: k, nz, itt, max_itt, k_cent
261 logical :: surface_BL, bottom_BL, central, halves, debug
262 logical :: old_PE_calc
264 h_neglect = gv%H_subroundoff
268 surface_bl = .true. ; bottom_bl = .true. ; halves = .true.
269 central = .true. ; k_cent = nz/2
271 do_print = .false. ;
if (
present(may_print) .and.
present(cs)) do_print = may_print
273 dpea_dkd(:) = 0.0 ; dpea_dkd_est(:) = 0.0 ; dpea_dkd_err(:) = 0.0
274 dpea_dkd_err_norm(:) = 0.0 ; dpea_dkd_trunc(:) = 0.0
275 dpeb_dkd(:) = 0.0 ; dpeb_dkd_est(:) = 0.0 ; dpeb_dkd_err(:) = 0.0
276 dpeb_dkd_err_norm(:) = 0.0 ; dpeb_dkd_trunc(:) = 0.0
278 htot = 0.0 ; pres(1) = 0.0 ; pres_z(1) = 0.0 ; z_int(1) = 0.0
280 t0(k) = t_in(k) ; s0(k) = s_in(k)
282 htot = htot + h_tr(k)
283 pres(k+1) = pres(k) + gv%H_to_Pa * h_tr(k)
284 pres_z(k+1) = us%Z_to_m * pres(k+1)
285 p_lay(k) = 0.5*(pres(k) + pres(k+1))
286 z_int(k+1) = z_int(k) - h_tr(k)
289 h_tr(k) = max(h_tr(k), 1e-15*htot)
294 kddt_h(1) = 0.0 ; kddt_h(nz+1) = 0.0
296 kddt_h(k) = min((gv%Z_to_H**2*dt)*kd(k) / (0.5*(h_tr(k-1) + h_tr(k))), 1e3*htot)
301 call calculate_specific_vol_derivs(t0, s0, p_lay, dsv_dt, dsv_ds, 1, nz, tv%eqn_of_state)
304 dmass = gv%H_to_kg_m2 * h_tr(k)
305 dpres = gv%H_to_Pa * h_tr(k)
306 dt_to_dpe(k) = (dmass * (pres(k) + 0.5*dpres)) * dsv_dt(k)
307 ds_to_dpe(k) = (dmass * (pres(k) + 0.5*dpres)) * dsv_ds(k)
308 dt_to_dcolht(k) = dmass * us%m_to_Z * dsv_dt(k) * cs%ColHt_scaling
309 ds_to_dcolht(k) = dmass * us%m_to_Z * dsv_ds(k) * cs%ColHt_scaling
314 pe_chg_k(:,:) = 0.0 ; colht_cor_k(:,:) = 0.0
318 old_pe_calc = .false.
322 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
323 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
324 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
325 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
326 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
334 if (old_pe_calc)
then
336 dt_km1_t2 = (t0(k)-t0(k-1))
337 ds_km1_t2 = (s0(k)-s0(k-1))
338 dte_t2 = 0.0 ; dse_t2 = 0.0
340 dte_t2 = kddt_h(k-1) * ((t0(k-2) - t0(k-1)) + dte(k-2))
341 dse_t2 = kddt_h(k-1) * ((s0(k-2) - s0(k-1)) + dse(k-2))
342 dt_km1_t2 = (t0(k)-t0(k-1)) - &
343 (kddt_h(k-1) / hp_a(k-1)) * ((t0(k-2) - t0(k-1)) + dte(k-2))
344 ds_km1_t2 = (s0(k)-s0(k-1)) - &
345 (kddt_h(k-1) / hp_a(k-1)) * ((s0(k-2) - s0(k-1)) + dse(k-2))
347 dte_term = dte_t2 + hp_a(k-1) * (t0(k-1)-t0(k))
348 dse_term = dse_t2 + hp_a(k-1) * (s0(k-1)-s0(k))
351 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
353 th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
354 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
356 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
362 kddt_h_guess = kddt_h(k)
363 if (old_pe_calc)
then
364 call find_pe_chg_orig(kddt_h_guess, h_tr(k), hp_a(k-1), &
365 dte_term, dse_term, dt_km1_t2, ds_km1_t2, &
366 dt_to_dpe(k), ds_to_dpe(k), dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
367 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
368 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
369 pe_chg_k(k,1), dpea_dkd(k))
371 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
372 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
373 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
374 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
375 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
376 pe_chg=pe_chg_k(k,1), dpec_dkd=dpea_dkd(k), &
377 colht_cor=colht_cor_k(k,1))
382 kddt_h_guess = (1.0+0.01*(itt-3))*kddt_h(k)
384 if (old_pe_calc)
then
385 call find_pe_chg_orig(kddt_h_guess, h_tr(k), hp_a(k-1), &
386 dte_term, dse_term, dt_km1_t2, ds_km1_t2, &
387 dt_to_dpe(k), ds_to_dpe(k), dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
388 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
389 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
392 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
393 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
394 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
395 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
396 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
401 dpea_dkd_est(k) = (4.0*(pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
402 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))) / 3.0
403 dpea_dkd_trunc(k) = (pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
404 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))
405 dpea_dkd_err(k) = (dpea_dkd_est(k) - dpea_dkd(k))
406 dpea_dkd_err_norm(k) = (dpea_dkd_est(k) - dpea_dkd(k)) / &
407 (abs(dpea_dkd_est(k)) + abs(dpea_dkd(k)) + 1e-100)
413 b1 = 1.0 / (hp_a(k-1) + kddt_h(k))
414 c1_a(k) = kddt_h(k) * b1
416 te(1) = b1*(h_tr(1)*t0(1))
417 se(1) = b1*(h_tr(1)*s0(1))
419 te(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
420 se(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
422 if (old_pe_calc)
then
423 dte(k-1) = b1 * ( kddt_h(k)*(t0(k)-t0(k-1)) + dte_t2 )
424 dse(k-1) = b1 * ( kddt_h(k)*(s0(k)-s0(k-1)) + dse_t2 )
427 hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kddt_h(k)
428 dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
429 ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
430 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
431 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
435 b1 = 1.0 / (hp_a(nz))
436 tf(nz) = b1 * (h_tr(nz) * t0(nz) + kddt_h(nz) * te(nz-1))
437 sf(nz) = b1 * (h_tr(nz) * s0(nz) + kddt_h(nz) * se(nz-1))
438 if (old_pe_calc)
then
439 dte(nz) = b1 * kddt_h(nz) * ((t0(nz-1)-t0(nz)) + dte(nz-1))
440 dse(nz) = b1 * kddt_h(nz) * ((s0(nz-1)-s0(nz)) + dse(nz-1))
444 tf(k) = te(k) + c1_a(k+1)*tf(k+1)
445 sf(k) = se(k) + c1_a(k+1)*sf(k+1)
449 do k=1,nz ; ta(k) = tf(k) ; sa(k) = sf(k) ;
enddo
450 pe_chg_tot1a = 0.0 ; pe_chg_tot2a = 0.0 ; t_chg_tota = 0.0
452 pe_chg_tot1a = pe_chg_tot1a + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
453 ds_to_dpe(k) * (sf(k) - s0(k)))
454 t_chg_tota = t_chg_tota + h_tr(k) * (tf(k) - t0(k))
457 pe_chg_tot2a = pe_chg_tot2a + (pe_chg_k(k,1) - colht_cor_k(k,1))
464 old_pe_calc = .false.
468 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
469 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
470 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
471 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
472 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
480 if (old_pe_calc)
then
482 dt_k_t2 = (t0(k-1)-t0(k))
483 ds_k_t2 = (s0(k-1)-s0(k))
484 dte_t2 = 0.0 ; dse_t2 = 0.0
486 dte_t2 = kddt_h(k+1) * ((t0(k+1) - t0(k)) + dte(k+1))
487 dse_t2 = kddt_h(k+1) * ((s0(k+1) - s0(k)) + dse(k+1))
488 dt_k_t2 = (t0(k-1)-t0(k)) - &
489 (kddt_h(k+1)/ hp_b(k)) * ((t0(k+1) - t0(k)) + dte(k+1))
490 ds_k_t2 = (s0(k-1)-s0(k)) - &
491 (kddt_h(k+1)/ hp_b(k)) * ((s0(k+1) - s0(k)) + dse(k+1))
493 dte_term = dte_t2 + hp_b(k) * (t0(k)-t0(k-1))
494 dse_term = dse_t2 + hp_b(k) * (s0(k)-s0(k-1))
496 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
498 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
500 th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te(k+1)
501 sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se(k+1)
506 kddt_h_guess = kddt_h(k)
508 if (old_pe_calc)
then
509 call find_pe_chg_orig(kddt_h_guess, h_tr(k-1), hp_b(k), &
510 dte_term, dse_term, dt_k_t2, ds_k_t2, &
511 dt_to_dpe(k-1), ds_to_dpe(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
512 pres_z(k), dt_to_dcolht(k-1), ds_to_dcolht(k-1), &
513 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
514 pe_chg=pe_chg_k(k,2), dpec_dkd=dpeb_dkd(k))
516 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
517 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
518 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
519 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
520 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
521 pe_chg=pe_chg_k(k,2), dpec_dkd=dpeb_dkd(k), &
522 colht_cor=colht_cor_k(k,2))
528 kddt_h_guess = (1.0+0.01*(itt-3))*kddt_h(k)
530 if (old_pe_calc)
then
531 call find_pe_chg_orig(kddt_h_guess, h_tr(k-1), hp_b(k), &
532 dte_term, dse_term, dt_k_t2, ds_k_t2, &
533 dt_to_dpe(k-1), ds_to_dpe(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
534 pres_z(k), dt_to_dcolht(k-1), ds_to_dcolht(k-1), &
535 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
538 call find_pe_chg(0.0, kddt_h_guess, hp_a(k-1), hp_b(k), &
539 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
540 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
541 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
542 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
547 dpeb_dkd_est(k) = (4.0*(pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
548 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))) / 3.0
549 dpeb_dkd_trunc(k) = (pe_chg(4)-pe_chg(2))/(0.02*kddt_h(k)) - &
550 (pe_chg(5)-pe_chg(1))/(0.04*kddt_h(k))
551 dpeb_dkd_err(k) = (dpeb_dkd_est(k) - dpeb_dkd(k))
552 dpeb_dkd_err_norm(k) = (dpeb_dkd_est(k) - dpeb_dkd(k)) / &
553 (abs(dpeb_dkd_est(k)) + abs(dpeb_dkd(k)) + 1e-100)
559 b1 = 1.0 / (hp_b(k) + kddt_h(k))
560 c1_b(k) = kddt_h(k) * b1
562 te(nz) = b1* (h_tr(nz)*t0(nz))
563 se(nz) = b1* (h_tr(nz)*s0(nz))
565 te(k) = b1 * (h_tr(k) * t0(k) + kddt_h(k+1) * te(k+1))
566 se(k) = b1 * (h_tr(k) * s0(k) + kddt_h(k+1) * se(k+1))
568 if (old_pe_calc)
then
569 dte(k) = b1 * ( kddt_h(k)*(t0(k-1)-t0(k)) + dte_t2 )
570 dse(k) = b1 * ( kddt_h(k)*(s0(k-1)-s0(k)) + dse_t2 )
573 hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kddt_h(k)
574 dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
575 ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
576 dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
577 ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
582 tf(1) = b1 * (h_tr(1) * t0(1) + kddt_h(2) * te(2))
583 sf(1) = b1 * (h_tr(1) * s0(1) + kddt_h(2) * se(2))
584 if (old_pe_calc)
then
585 dte(1) = b1 * kddt_h(2) * ((t0(2)-t0(1)) + dte(2))
586 dse(1) = b1 * kddt_h(2) * ((s0(2)-s0(1)) + dse(2))
590 tf(k) = te(k) + c1_b(k)*tf(k-1)
591 sf(k) = se(k) + c1_b(k)*sf(k-1)
595 do k=1,nz ; tb(k) = tf(k) ; sb(k) = sf(k) ;
enddo
596 pe_chg_tot1b = 0.0 ; pe_chg_tot2b = 0.0 ; t_chg_totb = 0.0
598 pe_chg_tot1b = pe_chg_tot1b + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
599 ds_to_dpe(k) * (sf(k) - s0(k)))
600 t_chg_totb = t_chg_totb + h_tr(k) * (tf(k) - t0(k))
603 pe_chg_tot2b = pe_chg_tot2b + (pe_chg_k(k,2) - colht_cor_k(k,2))
613 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
614 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
615 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
616 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
617 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
626 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
628 th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te_a(k-2)
629 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se_a(k-2)
631 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
633 kddt_h_a(k) = 0.0 ;
if (k < k_cent) kddt_h_a(k) = kddt_h(k)
636 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
637 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
638 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
639 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
640 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
641 pe_chg=pe_change, colht_cor=colht_cor)
642 pe_chg_k(k,3) = pe_change
643 colht_cor_k(k,3) = colht_cor
645 b1 = 1.0 / (hp_a(k-1) + kddt_h_a(k))
646 c1_a(k) = kddt_h_a(k) * b1
648 te_a(1) = b1*(h_tr(1)*t0(1))
649 se_a(1) = b1*(h_tr(1)*s0(1))
651 te_a(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kddt_h_a(k-1) * te_a(k-2))
652 se_a(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kddt_h_a(k-1) * se_a(k-2))
655 hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kddt_h_a(k)
656 dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
657 ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
658 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
659 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
668 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
674 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
676 th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te_b(k+1)
677 sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se_b(k+1)
680 kddt_h_b(k) = 0.0 ;
if (k > k_cent) kddt_h_b(k) = kddt_h(k)
683 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
684 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
685 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
686 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
687 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
688 pe_chg=pe_change, colht_cor=colht_cor)
689 pe_chg_k(k,3) = pe_chg_k(k,3) + pe_change
690 colht_cor_k(k,3) = colht_cor_k(k,3) + colht_cor
692 b1 = 1.0 / (hp_b(k) + kddt_h_b(k))
693 c1_b(k) = kddt_h_b(k) * b1
695 te_b(k) = b1 * (h_tr(k)*t0(k))
696 se_b(k) = b1 * (h_tr(k)*s0(k))
698 te_b(k) = b1 * (h_tr(k) * t0(k) + kddt_h_b(k+1) * te_b(k+1))
699 se_b(k) = b1 * (h_tr(k) * s0(k) + kddt_h_b(k+1) * se_b(k+1))
702 hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kddt_h_b(k)
703 dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
704 ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
705 dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
706 ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
716 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
718 th_a(k-1) = h_tr(k-1) * t0(k-1) + kddt_h(k-1) * te_a(k-2)
719 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kddt_h(k-1) * se_a(k-2)
722 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
724 th_b(k) = h_tr(k) * t0(k) + kddt_h(k+1) * te_b(k+1)
725 sh_b(k) = h_tr(k) * s0(k) + kddt_h(k+1) * se_b(k+1)
730 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
731 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
732 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
733 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
734 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
735 pe_chg=pe_change, colht_cor=colht_cor)
736 pe_chg_k(k,3) = pe_chg_k(k,3) + pe_change
737 colht_cor_k(k,3) = colht_cor_k(k,3) + colht_cor
752 b1 = 1.0 / (hp_a(k-1)*hp_b(k) + kddt_h(k)*(hp_a(k-1) + hp_b(k)))
753 tf(k-1) = ((hp_b(k) + kddt_h(k)) * th_a(k-1) + kddt_h(k) * th_b(k) ) * b1
754 sf(k-1) = ((hp_b(k) + kddt_h(k)) * sh_a(k-1) + kddt_h(k) * sh_b(k) ) * b1
755 tf(k) = (kddt_h(k) * th_a(k-1) + (hp_a(k-1) + kddt_h(k)) * th_b(k) ) * b1
756 sf(k) = (kddt_h(k) * sh_a(k-1) + (hp_a(k-1) + kddt_h(k)) * sh_b(k) ) * b1
758 c1_a(k) = kddt_h(k) / (hp_a(k-1) + kddt_h(k))
759 c1_b(k) = kddt_h(k) / (hp_b(k) + kddt_h(k))
764 tf(k) = te_a(k) + c1_a(k+1)*tf(k+1)
765 sf(k) = se_a(k) + c1_a(k+1)*sf(k+1)
768 tf(k) = te_b(k) + c1_b(k)*tf(k-1)
769 sf(k) = se_b(k) + c1_b(k)*sf(k-1)
773 pe_chg_tot1c = 0.0 ; pe_chg_tot2c = 0.0 ; t_chg_totc = 0.0
775 pe_chg_tot1c = pe_chg_tot1c + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
776 ds_to_dpe(k) * (sf(k) - s0(k)))
777 t_chg_totc = t_chg_totc + h_tr(k) * (tf(k) - t0(k))
780 pe_chg_tot2c = pe_chg_tot2c + (pe_chg_k(k,3) - colht_cor_k(k,3))
790 hp_a(k) = h_tr(k) ; hp_b(k) = h_tr(k)
791 dt_to_dpe_a(k) = dt_to_dpe(k) ; ds_to_dpe_a(k) = ds_to_dpe(k)
792 dt_to_dpe_b(k) = dt_to_dpe(k) ; ds_to_dpe_b(k) = ds_to_dpe(k)
793 dt_to_dcolht_a(k) = dt_to_dcolht(k) ; ds_to_dcolht_a(k) = ds_to_dcolht(k)
794 dt_to_dcolht_b(k) = dt_to_dcolht(k) ; ds_to_dcolht_b(k) = ds_to_dcolht(k)
806 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
808 th_a(k-1) = h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2)
809 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2)
811 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
813 dkd = 0.5 * kddt_h(k) - kd_so_far(k)
815 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
816 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
817 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
818 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
819 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
820 pe_chg=pe_change, colht_cor=colht_cor)
822 pe_chg_k(k,4) = pe_change
823 colht_cor_k(k,4) = colht_cor
825 kd_so_far(k) = kd_so_far(k) + dkd
827 b1 = 1.0 / (hp_a(k-1) + kd_so_far(k))
828 c1_a(k) = kd_so_far(k) * b1
830 te(1) = b1*(h_tr(1)*t0(1))
831 se(1) = b1*(h_tr(1)*s0(1))
833 te(k-1) = b1 * (h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2))
834 se(k-1) = b1 * (h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2))
837 hp_a(k) = h_tr(k) + (hp_a(k-1) * b1) * kd_so_far(k)
838 dt_to_dpe_a(k) = dt_to_dpe(k) + c1_a(k)*dt_to_dpe_a(k-1)
839 ds_to_dpe_a(k) = ds_to_dpe(k) + c1_a(k)*ds_to_dpe_a(k-1)
840 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1_a(k)*dt_to_dcolht_a(k-1)
841 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1_a(k)*ds_to_dcolht_a(k-1)
849 th_a(k-1) = h_tr(k-1) * t0(k-1) ; sh_a(k-1) = h_tr(k-1) * s0(k-1)
851 th_a(k-1) = h_tr(k-1) * t0(k-1) + kd_so_far(k-1) * te(k-2)
852 sh_a(k-1) = h_tr(k-1) * s0(k-1) + kd_so_far(k-1) * se(k-2)
855 th_b(k) = h_tr(k) * t0(k) ; sh_b(k) = h_tr(k) * s0(k)
857 th_b(k) = h_tr(k) * t0(k) + kd_so_far(k+1) * te(k+1)
858 sh_b(k) = h_tr(k) * s0(k) + kd_so_far(k+1) * se(k+1)
861 dkd = kddt_h(k) - kd_so_far(k)
863 call find_pe_chg(kd0, dkd, hp_a(k-1), hp_b(k), &
864 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
865 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe_b(k), ds_to_dpe_b(k), &
866 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
867 dt_to_dcolht_b(k), ds_to_dcolht_b(k), &
868 pe_chg=pe_change, colht_cor=colht_cor)
870 pe_chg_k(k,4) = pe_chg_k(k,4) + pe_change
871 colht_cor_k(k,4) = colht_cor_k(k,4) + colht_cor
874 kd_so_far(k) = kd_so_far(k) + dkd
876 b1 = 1.0 / (hp_b(k) + kd_so_far(k))
877 c1_b(k) = kd_so_far(k) * b1
879 te(k) = b1 * (h_tr(k)*t0(k))
880 se(k) = b1 * (h_tr(k)*s0(k))
882 te(k) = b1 * (h_tr(k) * t0(k) + kd_so_far(k+1) * te(k+1))
883 se(k) = b1 * (h_tr(k) * s0(k) + kd_so_far(k+1) * se(k+1))
886 hp_b(k-1) = h_tr(k-1) + (hp_b(k) * b1) * kd_so_far(k)
887 dt_to_dpe_b(k-1) = dt_to_dpe(k-1) + c1_b(k)*dt_to_dpe_b(k)
888 ds_to_dpe_b(k-1) = ds_to_dpe(k-1) + c1_b(k)*ds_to_dpe_b(k)
889 dt_to_dcolht_b(k-1) = dt_to_dcolht(k-1) + c1_b(k)*dt_to_dcolht_b(k)
890 ds_to_dcolht_b(k-1) = ds_to_dcolht(k-1) + c1_b(k)*ds_to_dcolht_b(k)
897 tf(1) = b1 * (h_tr(1) * t0(1) + kddt_h(2) * te(2))
898 sf(1) = b1 * (h_tr(1) * s0(1) + kddt_h(2) * se(2))
900 tf(k) = te(k) + c1_b(k)*tf(k-1)
901 sf(k) = se(k) + c1_b(k)*sf(k-1)
905 pe_chg_tot1d = 0.0 ; pe_chg_tot2d = 0.0 ; t_chg_totd = 0.0
907 pe_chg_tot1d = pe_chg_tot1d + (dt_to_dpe(k) * (tf(k) - t0(k)) + &
908 ds_to_dpe(k) * (sf(k) - s0(k)))
909 t_chg_totd = t_chg_totd + h_tr(k) * (tf(k) - t0(k))
912 pe_chg_tot2d = pe_chg_tot2d + (pe_chg_k(k,4) - colht_cor_k(k,4))
918 energy_kd = 0.0 ;
do k=2,nz ; energy_kd = energy_kd + pe_chg_k(k,1) ;
enddo
919 energy_kd = energy_kd / dt
923 if (cs%id_ERt>0)
call post_data(cs%id_ERt, pe_chg_k(:,1), cs%diag)
924 if (cs%id_ERb>0)
call post_data(cs%id_ERb, pe_chg_k(:,2), cs%diag)
925 if (cs%id_ERc>0)
call post_data(cs%id_ERc, pe_chg_k(:,3), cs%diag)
926 if (cs%id_ERh>0)
call post_data(cs%id_ERh, pe_chg_k(:,4), cs%diag)
927 if (cs%id_Kddt>0)
call post_data(cs%id_Kddt, kddt_h, cs%diag)
928 if (cs%id_Kd>0)
call post_data(cs%id_Kd, kd, cs%diag)
929 if (cs%id_h>0)
call post_data(cs%id_h, h_tr, cs%diag)
930 if (cs%id_zInt>0)
call post_data(cs%id_zInt, z_int, cs%diag)
931 if (cs%id_CHCt>0)
call post_data(cs%id_CHCt, colht_cor_k(:,1), cs%diag)
932 if (cs%id_CHCb>0)
call post_data(cs%id_CHCb, colht_cor_k(:,2), cs%diag)
933 if (cs%id_CHCc>0)
call post_data(cs%id_CHCc, colht_cor_k(:,3), cs%diag)
934 if (cs%id_CHCh>0)
call post_data(cs%id_CHCh, colht_cor_k(:,4), cs%diag)
935 if (cs%id_T0>0)
call post_data(cs%id_T0, t0, cs%diag)
936 if (cs%id_Tf>0)
call post_data(cs%id_Tf, tf, cs%diag)
937 if (cs%id_S0>0)
call post_data(cs%id_S0, s0, cs%diag)
938 if (cs%id_Sf>0)
call post_data(cs%id_Sf, sf, cs%diag)
939 if (cs%id_N2_0>0)
then
940 n2(1) = 0.0 ; n2(nz+1) = 0.0
942 call calculate_density(0.5*(t0(k-1) + t0(k)), 0.5*(s0(k-1) + s0(k)), &
943 pres(k), rho_here, tv%eqn_of_state)
944 n2(k) = ((us%L_to_Z**2*gv%g_Earth) * rho_here / (0.5*gv%H_to_Z*(h_tr(k-1) + h_tr(k)))) * &
945 ( 0.5*(dsv_dt(k-1) + dsv_dt(k)) * (t0(k-1) - t0(k)) + &
946 0.5*(dsv_ds(k-1) + dsv_ds(k)) * (s0(k-1) - s0(k)) )
948 call post_data(cs%id_N2_0, n2, cs%diag)
950 if (cs%id_N2_f>0)
then
951 n2(1) = 0.0 ; n2(nz+1) = 0.0
953 call calculate_density(0.5*(tf(k-1) + tf(k)), 0.5*(sf(k-1) + sf(k)), &
954 pres(k), rho_here, tv%eqn_of_state)
955 n2(k) = ((us%L_to_Z**2*gv%g_Earth) * rho_here / (0.5*gv%H_to_Z*(h_tr(k-1) + h_tr(k)))) * &
956 ( 0.5*(dsv_dt(k-1) + dsv_dt(k)) * (tf(k-1) - tf(k)) + &
957 0.5*(dsv_ds(k-1) + dsv_ds(k)) * (sf(k-1) - sf(k)) )
959 call post_data(cs%id_N2_f, n2, cs%diag)