Open ars599 opened 2 years ago
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
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
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