The mom output of the Zonal/meridional/vertical advections

[ars599]

Dear all,
I have tried to follow up with Fabio's heat budget analysis and know more details about the tropical advection process. So far, in Fabio's script, he is using the model output temp_sweby_advec, and that is the total advection. However, I would like to get zonal/meridional/vertical components.

Would anyone please let me know how I can calculate or output the zonal/meridional/vertical advections?

Regards,

Arnold

[StephenGriffies]

Assuming the model uses advect_sweby_all=.true., then you can enable the diagnostic table entries temp_x_adv, temp_y_adv, and temp_z_adv to get the three separate components. However, no other advection scheme has that coded.

The reason we reluctantly coded the separated advection direction diagnostics is that most modern advection schemes do not allow for a clear separation between the three advection directions when considering regional budgets. We coded the separate contributions for advect_sweby_all=.true. only after some debate among the GFDL developers, and only for that scheme since it was used in the GFDL CM2.1 model.

Besides the numerical niceties mentioned above, there are fundamental questions about how to interpret separate components of an advection operator for a non-divergent velocity field. Splitting directional contributions in regional budgets takes some care to render a physically meaningful diagnostic. Examples of the issues I am referring to can be found here:

Appendix B in Gregory (2000): "Vertical heat transports in the ocean and their effect on time-dependent climate change"
Section 3 in Ray et al (2018): DOI: 10.1175/JCLI-D-18-0152.1

[ars599]

Dear Steve,

Many thanks for your help. I will go through the papers you referred. Here I would like to share a script that I got from Jaison Kurian more than ten years ago. It's an offline calculation. However, I haven't got chance to clearify all the calculations, and would like to have somebody's help. It was the method described in Vialard and Delecluse (1998a [1]) to find the terms of mixed layer heat budget. I guess after the correction of the script, it will be very helpful for the COSIMA community.

[1] https://journals.ametsoc.org/view/journals/phoc/28/6/1520-0485_1998_028_1071_aosftt_2.0.co_2.xml

\ cancel mode verify
! Original copy from
! http://web.atmos.ucla.edu/~jaison/ctrl_ml_heatb_const_vdiff_corr.jnl
!
! No vertical diffusion values from model outputs are used here.
! Then the entraintment term E is calculated with K/(rhocp)
! = 510^-5 m2/s as in prange_et_al_2003.pdf
!
! Please note that the value for constant vertical diffusivity (m^2/sec)
! "kappa_h" from /mom4p0c_dec102004/mom4/src/mom4/ocean_param/mixing/vert/const/
! ocean_vert_mix_coeff.f90 is 0.1e-4 --> 110^-5
!
!
! Use the method described in Vialard and Delecluse (1998a) to find the
! terms of mixed layer heat budget.
!
! . --> subscript
! ' --> integral ov. ML
! dT'
! --- = Ax' + Ay' + Az' + Dh' + F + E
! dt
!
! -1 | 0 dT -1 | 0 dT
! Ax = --- | u---dz ; Ay = --- | v---dz
! h |-h dx h |-h dy
!
! 1 |0 dw
! Az = --- ( W.-h (T.-h - T') + | (T - T')* -- * dz
! h |-h dz
!
! 1 |0
! Dh = ---* | D
! h |-h
!
! Q* + Qs[1-f(-h)]
! F = ---------------- ; f(-h) = R * exp(-h/L1) + (1-R) * exp(-h/L2)
! p.0 * Cp * h
!
! Qnet = sfc_hflux_coupler + sfc_hflux_from_runoff
! + sfc_hflux_pme + frazil_2d + sw_heat[summed in ML]
!
! 1 [ dh [ dT ] ]
! E = ---* [----(T.-h - T') + [K---- ] / p.0Cp ]
! h [ dt [ dz }.-h ] ]
!
!
! NOTE : sfc_hflux is contained in the from k=1 level of vertical diffusion
! but this should be added to k=1 level of sw_heat according the
! term balance method for ACCESS-OM. Hence the Eqn.A9 in Jerome's
! paper is modified as follows:
!
! Q + Qs[1-f(-h)] --> net heat gained by the mixed layer
!
! From the term balance jnl file
!
! let source1= sw_heat[d=1] !short wave heating, without correction.
! let source =
! IF K[g=source1] EQ 1 THEN source1+sfc_hflux[d=5] ELSE source1 ! corrected
!
! ==> now add source from k=1 to ML base. --> net heat gained by the mixed layer
!
! We can use swflx from ocean_diag.nc and sfc_hflux to use the default
! method. But let us believe on the sw pen of ACCESS-OM.
!
! MLD : from ACCESS-OM or ACCESS-OMd outputs
!
! Note : This is a copy of original script at
! ~/analysis/ACCESS-OM/clim_omip/ml_heatbudg/mom_mld_extn/
! ctrl_ml_heatb_const_vdiff.jnl
! Use Ferret version 5.8 to get proper variable attributes.
!
! Created On : 11/Dec/2006
! Created By : JK
!
! Modifications : 1. On 03/Jul/2007 : Minimal changes to suit the writing part.
! However, the code remains exactly the same as the one
! used for JGR paper (results confirmed manually).
!
!-----------------------------------------------------------------------
!
!----USER INPUTS
!
define symbol do_write = YES
def symb expName = $1 !! hPIC2C !!
define symbol ncfl = ($expName)_ml_heatb_const_vdiff_fsg599.nc
! set reg/t=16-JAN-1900:16-DEC-2005
! set reg/x=0:360/y=-90:90/t=01-JAN-0001:31-DEC-0050
!
! cache enough memory size
!
set memory/size=6000
!
! open datasets (keep the order as it is)
!
define symbol mer = /data/sul086/($expName)/month
! define symbol mer = /short/p66/ars599/($expName)/month/
set data "($mer)/($expName).ocean.pot_temp.nc" ! d=1
set data "($mer)/($expName).ocean.mld.nc" ! d=2
set data "($mer)/($expName).ocean.u.nc" ! d=3
set data "($mer)/($expName).ocean.v.nc" ! d=4
set data "($mer)/($expName).ocean.wt.nc" ! d=5
set data "($mer)/($expName).ocean.sw_heat.nc" ! d=6
set data "($mer)/($expName).ocean.sfc_hflux_coupler.nc" ! d=7
set data "($mer)/($expName).ocean.sfc_hflux_pme.nc" ! d=8
set data "($mer)/($expName).ocean.sfc_hflux_from_runoff.nc" ! d=9
! set data "($mer)/($expName).ocean.rho.nc" ! d=10
set data "($mer)/($expName).ocean.frazil_2d.nc" ! d=10
sh/bri da
let l_start=pot_temp[d=1],r=lstart
let l_end =pot_temp[d=1],r=lend

 let t_ctrl   = pot_temp[d=1]         ! unit: 30  C
 let h_ctrl   = mld[d=2]              ! unit: 100 m
 let u_ctrl   = u[d=3,g=t_ctrl@ASN]   ! unit: 0.5 m/s
 let v_ctrl   = v[d=4,g=t_ctrl@ASN]   ! unit: 0.5 m/s
 let w_ctrl   = wt[d=5]               ! unit: 1e5 m/s

! let r_ctrl = rho[d=10] ! unit: 1000

!----define the depth of cell tops (cell_top) such that--------
!----if MLD is 5, it should choose only k=1 but if MLD is 5.01, then
! it should choose k=1,2. and so on

! let zgap = zbox[gz=pot_temp[d=1]] ! suggested by sjm599
! let ztop = z[gz=pot_temp[d=1]] ! suggested by sjm599
let ztop = ZSEQUENCE({5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125, 135, 145,
155, 165, 175, 185, 195, 205, 216.846755981445, 241.349014282227,
280.780731201172, 343.250457763672, 427.315551757812, 536.715637207031,
665.414123535156, 812.781616210938, 969.065124511719, 1130.93493652344,
1289.60461425781, 1455.77014160156, 1622.92565917969, 1801.55810546875,
1984.85461425781, 2182.90478515625, 2388.41748046875, 2610.93505859375,
2842.564453125, 3092.20483398438, 3351.29467773438, 3628.0576171875,
3913.26440429688, 4214.4951171875, 4521.91796875, 4842.56591796875,
5166.1298828125, 5499.2451171875, 5831.29443359375})
let cell_top = ztop[GZ=t_ctrl@ASN]
let cell_bot = cell_top[k=@shf:1] !the bottom depth corresponding to cell_top

! define constants used

 let rho     = 1035.00000       ! Kg/m^3
 let Cp      = 3989.24495292815 ! J/kg/deg
 let convert = 365.25/12*24*60*60      ! sec to month

!----extend MLD to base of that cell--------------------------------------
let mld_mask = IF cell_top LT h_ctrl THEN 1 ! to mask vars on T-cell
let h_extn_1 = cell_bot * mld_mask ! because "K" will be the same
let h_kvals = K[g=h_extn_1] + h_extn_1 * 0
let h_kzero = h_kvals - h_kvals[k=1:h_kvals,return=kend@max]
let h_extn_z = h_extn_1 * h_kzero[k=@weq:0]
let h_extn = h_extn_z[k=@sum]
let r_mask = r_ctrl * mld_mask

! Find Temp over mixed layer --> T' --> use the bottom values of that level
! define symbol do_write = NO

 let t_mld    = t_ctrl * mld_mask
 let T_mix    = t_mld[Z=@AVE]
 let dT_dt    = T_mix[t=@DDC]*convert
 set var/title="SST Tendency : Av. Ov. MLD (ACCESS-OM Ctrl)"/\
                 units="^oC/month"  dT_dt

 IF ($do_write) THEN
   sp rm -f ($ncfl)
   REPEAT/L=`l_start`:`l_end`:1 (;\
      SAVE/file=($ncfl)/append/quiet dT_dt;\
   )
 ENDIF

! Find net heat gained by the mixed layer : our own method based on ACCESS-OM outputs
! define symbol do_write = YES

! sw_heat[k=1] = sw_heat[k=2:50@sum] * (-1)
! sw_into_layer1 = swflx - sw_heat[k=2:50@sum] = swflx + sw_heat[k=1]
! sw_into_mixed_layer = swflx - sw_heat[k=@mld:50@sum]
!
let sw_ctrl = sw_heat[d=6,g=t_ctrl@ASN] ! unit: -100 W/m^2
let sfc_c_ctrl = sfc_hflux_coupler[d=7] ! unit: -300 W/m^2
let sfc_rf_ctrl = sfc_hflux_from_runoff[d=9] ! unit: 1 W/m^2
let sfc_pme_ctrl = sfc_hflux_pme[d=8] ! unit: 20 W/m^2
let fr_ctrl = frazil_2d[d=10] ! unit: 20 W/m^2

     let source1   = sw_ctrl * mld_mask !! sw pen in mld
     let sfc_hflux = sfc_c_ctrl + sfc_rf_ctrl + sfc_pme_ctrl + fr_ctrl
 let Qnet      = IF K[g=source1] EQ 1 THEN source1+sfc_hflux ELSE source1

! let SEF = (Qnet[z=@sum]/(r_mask[k=@ave]Cph_extn) ) * convert
let SEF = (Qnet[z=@sum]/(rhoCph_extn) ) * convert

     set var/title="Atm. Forcing : Av. Ov. MLD(ACCESS-OM Ctrl)"/\
                     units="^oC/month"  SEF

     IF ($do_write) THEN
         REPEAT/L=`l_start`:`l_end`:1 (;\
            save/file=($ncfl)/append/quiet SEF;\
         )
     ENDIF

! Find advection fields --> no explicit masking for u/v is needed since we use t_mld
! define symbol do_write = NO

 let dT_dx = t_mld[x=@DDC] 
 let u_adv = u_ctrl * dT_dx                    
 let Ax    = ( (-1/h_extn)* u_adv[z=@DIN] ) * convert
 set var/title="U Adv. : Av. Ov. MLD (ACCESS-OM Ctrl)"/\
                 units="^oC/month"  Ax

 IF ($do_write) THEN
    REPEAT/L=`l_start`:`l_end`:1 (;\
        save/file=($ncfl)/append/quiet Ax;\
    )
 ENDIF

 let dT_dy = t_mld[y=@DDC] 
 let v_adv = v_ctrl * dT_dy
 let Ay    = ( (-1/h_extn)* v_adv[z=@DIN] ) * convert
 set var/title="V Adv. : Av. Ov. MLD (ACCESS-OM Ctrl)"/\
                 units="^oC/month"  Ay

 IF ($do_write) THEN
    REPEAT/L=`l_start`:`l_end`:1 (;\
        save/file=($ncfl)/append/quiet Ay;\
    )
 ENDIF

 ! find T at base of ML ie -h. This part has been modified on 12/Dec/2006
 ! do it for exact depth of -h instead of going to next level..That is at
 ! the depth of extended MLD. First we need to get the temperature values 
 ! interpolated at "ztop" values so that k values will be matching exactly 
 ! with original temp variable. In essence, only the temperature variable
 ! will change all other lines remains the same (except the line which
 ! assigns the orginal depth grid before @SUM).

     define axis/z/units=meters/DEPTH/from_data zcelltop=ztop
     let t_ztop       = t_ctrl[gz=zcelltop]
     let bot_mask     = IF cell_top GE h_ctrl THEN 1
     let t_bot        = t_ztop * bot_mask
     let kvals        = k[g=t_bot] + t_bot * 0
     let kzero        = kvals - kvals[k=1:`kvals,return=kend`@MIN]
     let T_mix_base_z = t_bot * kzero[k=@WEQ:0]
     let T_mix_base_z1= T_mix_base_z[gz=t_ctrl@ASN]
     let T_mix_base   = T_mix_base_z1[k=@SUM]

 ! find w at base of ML

     let mld_mask_w1  = IF cell_top LE h_ctrl THEN 1 ! it is "cell_top" itself
     let w1_mld       = w_ctrl * mld_mask_w1         !  then only we will get the wvel
     let wkvals       = K[g=w1_mld] + w1_mld *0      !  exactly at the base of ML
     let wkzero       = wkvals - wkvals[k=1:`wkvals,return=kend`@MAX]
     let w_mix_base_z = w1_mld * wkzero[k=@WEQ:0]
     let w_mix_base   = w_mix_base_z[k=@SUM]

 ! mask w for integrating within ML (we should exclude w at ML base)

     let mld_mask_w   = IF Z[GZ=w_ctrl] LE h_ctrl THEN 1
     let w_mld        = w_ctrl * mld_mask_w 
     
 ! do the integration part
 ! NOTE: this term will give missing values if mld is too shallow 
 !       because of w_mld[z=@DDC]. According to Jerome (1999a) this
 !       Az_t2_in is a negligiable term. Due to these reasons this
 !       term is not include in the current calculation
      
   ! let t_wgrd   = t_ctrl[g=w_ctrl@ASN] + w_mld*0    ! Temp. on w grid 
   ! let Az_t2    = (t_wgrd-T_mix)*w_mld[z=@DDC]      ! 2nd term in Az
   ! let Az_t2_in = Az_t2[z=@DIN]                     ! integral of 2nd term

 ! do first term in Az

     let Az_t1    = w_mix_base * (T_mix_base - T_mix)  

                 ! let Az = (1/h_extn) * ( Az_t1 + az_t2_in ) 
     let Az = (1/h_extn) * Az_t1 * convert
     set var/title="W Adv. : Av. Ov. MLD (ACCESS-OM Ctrl)"/\
                 units="^oC/month"  Az

     IF ($do_write) THEN
         REPEAT/L=`l_start`:`l_end`:1 (;\
            save/file=($ncfl)/append/quiet Az;\
         )
     ENDIF

! Find the horizontal diffusion integrated over ML (not used now)

! let hdiff_mld = temp_hdiff_neu[d=7]* mld_mask ! within ML
! let Dh = (1/h_extn)*hdiff_mld[z=@din]

! Find Entraintment using new constant value of K/(rhoCp) 510^-5 m2/s

 ! find term 1 of E
     let dh_dt           = h_extn[t=@DDC]
     let gamma           = IF dh_dt GT 0 THEN 1 ELSE 0 ! step function
     let E_term1         = dh_dt*gamma*(T_mix_base - T_mix)
     let E_t1            = (1/h_extn) * E_term1 * convert
     set var/title="Entraintment term 1: At Base of ML (ACCESS-OM Ctrl)"/\
                 units="^oC/month"  E_t1

    IF ($do_write) THEN
         REPEAT/L=`l_start`:`l_end`:1 (;\
            save/file=($ncfl)/append/quiet E_t1;\
         )
     ENDIF

 ! find term 2 of E

! define symbol do_write =YES

 let dT_dz           = t_ctrl[k=@DDC]
     let dT_dz_mld       = dT_dz*bot_mask
     let vd_kvals        = k[g=dT_dz_mld] + dT_dz_mld * 0
     let vd_kzero        = vd_kvals - vd_kvals[k=1:`vd_kvals,return=kend`@MIN]
     let dT_dz_base_z    = dT_dz_mld * vd_kzero[k=@WEQ:0]
     let dT_dz_base      = dT_dz_base_z[k=@SUM]

     let E_term2 = 1e-5 * dT_dz_base
     let E_t2    = (1/h_extn) * E_term2 * convert
     set var/title="Entraintment term 2: At Base of ML (ACCESS-OM Ctrl)"/\
                 units="^oC/month"  E_t2

     IF ($do_write) THEN
         REPEAT/L=`l_start`:`l_end`:1 (;\
            save/file=($ncfl)/append/quiet E_t2;\
         )
     ENDIF