[MITgcm-support] Heat budget in MITGCM

Tue Nov 18 11:19:34 EST 2014

Many thanks Ryan, Dimitris, John and Gael for your suggestions.

I am not much familiar with MATLAB, so for me it was too much difficult to
understand the code shared by Gael, though I tried to look into it.

As a matter of fact is that I am using is the ECCO1 version (1 deg global,
50 levels, excluding poles; i=1:360; j=1:160; k=1:50) of MITGCM with linear
free surface and linFSConserveTR=True. The diagnostic parameters are
written in every day with diagnostic timephase=0. Model parameterizations
are KPP and GM-Redi.

I’d like to analyze the heat budget at a particular location, say, (i,j,k).
So, from the suggestions by Ryan, Dimitris and looking into the Gael’s
code, I am trying the following:

!!Cell volumes (center difference)

Cvx = (RAC(i+1,j) - RAC(i-1,j)) * RF(k) * hFacC(i,j,k)

Cvy = (RAC(i,j+1) - RAC(i,j-1)) * RF(k) * hFacC(i,j,k)

Cvz = RAC(i,j) * (RF(k+1) – RF(k)) * hFacC(i,j,k)

!!Divergence of advection & diffusion by center difference

if (k>1) then

         adc_dif_vert_k = (ADVr_TH(i,j,k+1) – ADVr_TH(I,j,k-1))/cvz +

                              (DIFrE_TH(I,j,k+1) – DIFrE_TH(I,j,k-1))/cvz +

                              (DIFrI_TH(I,j,k+1) – DIFrI_TH(I,j,k-1))/cvz

else

         adc_dif_vert_k = 0.0

endif

adc_dif_hori_k = (ADVx_TH(i+1,j,k) – ADVx_TH(i-1,j,k))/cvx +

                         (ADVy_TH(i,j+1,k) – ADVy_TH(i,j-1,k))/cvy +

                         (DIFxE_TH(i+1,j,k) – DIFxE_TH(i-1,j,k))/cvx +

                         (DIFyE_TH(i,j+1,k) – DIFyE_TH(i,j-1,k))/cvy

!!To account for the  shortwave heating

depth=RF(k)

swfrac=0.62*exp(depth/0.6) + (1.0-0.62) * exp(depth/20)

if(depth < -200.0) then

swfrac=0.0

endif

if (k=1)then

               swfrac=1.0

endif

!!Finally balance

TOTTTEND/86400.0 = - adc_dif_vert_k -  adc_dif_hori_k +
(oceQsw/(rhoConst*Cp)) * swfrac * hFacC(I,j,k) + KPPg_TH(i,j,k)

Please let me know whether I am doing properly or not.

Thanking in advance,

Abhisek

*--------------------------------------------------------------------------------ABHISEK
CHAKRABORTY*

*Scientist - "SD"Oceanic Sciences Division (*
*AOSG/EPSA)*

*Space Applications Centre*

*Indian Space Research Organization*

*Ahmedabad - 380 015*

*Gujarat, INDIA*
*Contact: +91-79-2691-6054 (O), +91-79-2686-1929(R)*

On Tue, Nov 18, 2014 at 5:06 AM, gael forget <gforget at mit.edu> wrote:

> Hi Abhisek, et al,
>
> the following matlab code may provide some guidance :
>
> http://mitgcm.org/viewvc/MITgcm/MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_D.m
>
> It does not cover every possible combination of model options
> (e.g., no kpp or thsice) of course, but treats a few. In that code,
> it is assumed that
> 1) the setting for the diagnostics package was according to
>
> http://mitgcm.org/viewvc/MITgcm/MITgcm_contrib/gael/verification/global_oce_llc90/input/data.diagnostics
>      which includes vertical integration. The relevant diagnostics
>      for budget are those with filename like 'diags/budg2d_snap_set1',
> 2) tendencies have been pre-computed by differencing between
>     consecutive (monthly) snap shots — so that ETAN
>     really means dETAN/dt in this context, and similarly for
>     THETA, SALT, etc. The other variable names used in
>     diags_set_D.m
> <http://mitgcm.org/viewvc/MITgcm/MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_D.m> match
> those used in pkg/diagnostics
>
> I should also mention that diags_set_D.m
> <http://mitgcm.org/viewvc/MITgcm/MITgcm_contrib/gael/matlab_class/gcmfaces_diags/diags_set_D.m> is
> part of a
> broader framework (gcmfaces) which is documented in
>
> http://mitgcm.org/viewvc/*checkout*/MITgcm/MITgcm_contrib/gael/matlab_class/gcmfaces.pd
> <http://mitgcm.org/viewvc/MITgcm/MITgcm_contrib/gael/matlab_class/gcmfaces.pd>
> f
> … although I am afraid this documentation is also
> lacking wrt the budget diagnostic codes.
>
> Cheers,
> Gael
>
> On Nov 17, 2014, at 5:53 PM, Dimitris Menemenlis <dmenemenlis at gmail.com>
> wrote:
>
> … and you probably already know this, but not clear from your email:
> shortwave is not dumped in surface level only, but rather is distributed
> with exponentially
> decaying profile in top 200-m, with subtle modifications if this heating
> is occurring within
> the KPP mixing layer depth.
>
> On Nov 17, 2014, at 2:44 PM, Ryan Abernathey <ryan.abernathey at gmail.com>
> wrote:
>
> Abhisek,
>
> This is a very common question that unfortunately is not documented very
> well. A couple of suggestions...
>
> For the transport terms, (ADV* and DF*), keep in mind that these are FLUX
> terms, defined on the cell boundaries. In order to calculate a heat budget,
> you need to calculate the DIVERGENCE of those fluxes.
> http://mitgcm.org/public/r2_manual/latest/online_documents/node71.html
> Maybe this is obvious, but your email did not make it clear whether you
> were doing this or not.
>
> Likewise, TFLUX is the downward flux at the surface, in W/m^2. In order to
> convert to temperature units, you need to divide by
> HeatCapacity_Cp*rUnit2mass. For z coordinates, rUnit2mass is equal
> to rhoConst.
>
> There are some subtleties regarding tracer budgets under different free
> surface treatments. Maybe someone else on the list understands this issue
> better and is willing to explain. Anyway, those should be second-order
> effects.
>
> Good luck with your analysis.
>
> -Ryan
>
>
> On Sat, Nov 15, 2014 at 2:29 AM, Abhisek Chakraborty <
> abhisek.sac at gmail.com> wrote:
>
>> Dear Users,
>>
>> I am trying to analyze heat budget in MITGCM. From the diagnostic
>> outputs I have the following terms
>>
>> TOTTTEND , ADVr_TH, ADVx_TH, ADVy_TH, DFrE_TH, DFxE_TH, DFyE_TH,
>> DFrI_TH, KPPg_TH, TFLUX (i.e. all terms corresponding to tendency,
>> advection, diffusion, KPP and total heat flux).
>>
>> Apart from these, I have other usual outputs like THETA, SALT etc.
>>
>> The tendency term is in degC/day, the other terms are in degC*m^3/s
>> and TFLUX is in W/m^2. Thus I have converted the tendency term into
>> degC/s by dividing by 86400.0 and the advection-diffusion terms are
>> divided by cell volume (=RAC*RF*HFACC).
>>
>> I am using "linFSconserveTr=True".
>>
>> For the surface layer I have to consider TFLUX (watt/m^2), but how to
>> convert it into degC/s unit?
>>
>> For subsurface levels, I am trying to equate tendency term with
>> (advection + diffusion + KPP term). But there is a mismatch.
>>
>> Can somebody please guide me how to achieve the exact heat budget for
>> both surface and subsurface levels ?
>>
>> Thanking in advance,
>> Abhisek
>>
>> --
>>
>>
>> *ABHISEK CHAKRABORTY*
>>
>> *Scientist - "SD"Oceanic Sciences Division (*AOSG/EPSA)*
>>
>> *Space Applications Centre*
>>
>> *Indian Space Research Organization*
>>
>> *Ahmedabad - 380 015*
>>
>> *Gujarat, INDIA*
>> *Contact: +91-79-2691-6054 (O), +91-79-2686-1929(R)*
>>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
>
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20141118/5ee7331d/attachment-0001.htm>