[MITgcm-support] Heat budget in MITGCM
Abhisek Chakraborty
abhisek.sac at gmail.com
Fri Nov 21 08:43:31 EST 2014
Dear Jean-Michel,
To close the surface heat budget, I followed your suggestion:
TsurfCor = SUM( WTHMASS(i,j,1)*RAC(i,j) ) / globalArea
and finally:
tend_T = (...) - ( WTHMASS(i,j,1) - TsurfCor )/(drF(1)*hFacC(i,j,1))
Also, the shortwave tendency terms is computed as:
tend_Qsw = oceQsw/(rhoConst*Cp)/(drF(k)*hFacC(I,j,k))
* ( swfrac(k) - swfrac(k+1) )
where for surface, k=1
But in this way, the surface tendency term is becoming one order higher in
magnitude than other terms (advection, diffusion, surface correction. Pls.
see the figure in my earlier mail).
So, please guide me how to close the budget at surface.
regards,
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 Thu, Nov 20, 2014 at 6:59 PM, Abhisek Chakraborty <abhisek.sac at gmail.com>
wrote:
> Hi Jean-Michel,
>
> Sorry for delayed reply.
>
> Following your suggestions, I think I am now able to close the heat budget
> for the subsurface levels.
>
> However, in spite of my different attempts I am not able to close the
> budget at surface where I am finding the surface correction term (as
> specified in your mail) along with other advection-diffusion terms are too
> small to balance the shortwave heating term.
>
> I have attached two figures (one for subsurface at 20 m depth and another
> for surface).
>
> Pls. let me know if I am missing something.
>
> 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 Wed, Nov 19, 2014 at 1:19 AM, Jean-Michel Campin <jmc at ocean.mit.edu>
> wrote:
>
>> Hi Abhisek,
>>
>> This budget looks better. Just few remarks (see below).
>>
>> On Tue, Nov 18, 2014 at 11:20:46PM +0530, Abhisek Chakraborty wrote:
>> > Dear Jean-Michel,
>> >
>> > Thanks for your prompt suggestion.
>> >
>> >
>> > (1) the advection term for k>1 will be
>> >
>> > - [ (ADVx_TH(i+1,j,k) – ADVx_TH(I,j,k))/CV +
>> >
>> > (ADVy_TH(I,j+1,k) – ADVy_TH(I,j,k))/CV +
>> >
>> > (ADVr_TH(I,j,k+1) – ADVr_TH(I,j,k))/CV ]
>>
>> The vertical advective & diffusive transport/flux are (I think)
>> positive in the direction of the vertical coordinate, i.e.,
>> when using z-coords, positive upward. So you will have to reverse
>> (or add a minus sign) in front of:
>> (ADVr_TH(I,j,k+1) – ADVr_TH(I,j,k))/CV
>>
>> >
>> > Where CV=Ac * del RF * Hc = RAC(I,j) * (RF(K+1)-RF(K)) * hFacC(I,j,k)
>> >
>> > (
>> http://mitgcm.org/public/r2_manual/latest/online_documents/node73.html)
>> >
>> >
>> >
>> > And similarly the diffusion and KPP terms. Correct?
>> >
>> >
>> > (2) The source term (shortwave) for k>1 will be the following
>> >
>> > depth=RF(k)
>> >
>> > swfrac=0.62*exp(depth/0.6) + (1.0-0.62) * exp(depth/20)
>> >
>> > if(depth < -200.0) swfrac=0.0
>> >
>> > source term = (oceQsw/(rhoConst*Cp)) * swfrac * hFacC(I,j,k)
>> >
>> > Is this correct for k>1 ?
>>
>> This does not look right:
>> from depth=RF(k) , we compute swfrac(k) as above
>> from depth=RF(k+1), we compute swfrac(k+1) as above except that
>> at the bottom (deepest level of this column) swfrac(k+1)=0.
>> and the tendency contribution from shortwave heating should be something
>> like:
>> tend_Qsw = oceQsw/(rhoConst*Cp)/(drF(k)*hFacC(I,j,k))
>> * ( swfrac(k) - swfrac(k+1) )
>> >
>> > (3) I am using the tendency as TOTTTEND/86400. Is it correct?
>>
>> yes, providing you are using linear free-surface.
>>
>> >
>> > (4) I am using linFSConserveTr=TRUE. Yet I have to use
>> > WTHMASS(i,j,1)*RAC(i,j) to correct the heat budget at the surface?
>>
>> Yes, but since you are using linFSConserveTr=TRUE., the global mean
>> value of WTHMASS has to be removed in the budget (like it is in the
>> model):
>> TsurfCor = SUM( WTHMASS(i,j,1)*RAC(i,j) ) / globalArea
>> and finally:
>> tend_T = (...) - ( WTHMASS(i,j,1) - TsurfCor )/(drF(1)*hFacC(i,j,1))
>>
>> Cheers,
>> Jean-Michel
>>
>> >
>> > 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 10:42 PM, Jean-Michel Campin <jmc at ocean.mit.edu
>> >
>> > wrote:
>> >
>> > > Hi Abhisek,
>> > >
>> > > I would recommand to concentrate first on any level except the surface
>> > > (k=1),
>> > > and try to close the budget ; then, once you get a clean closed budget
>> > > for levels at depth, you can revisit the surface level problem.
>> > >
>> > > at depth (k>1):
>> > > the way you estimate the tendency is not what is done in the model,
>> > > see e.g.:
>> > >
>> http://mitgcm.org/public/r2_manual/latest/online_documents/node73.html
>> > >
>> > > The diagnostics for the KPP non-local term is like a vertical flux
>> > > (added to fVerT), so it can be treated (for budget purpose) the same
>> way
>> > > as DIFrI_TH.
>> > >
>> > > at the surface, you can search the mitgcm-support archive and find
>> some
>> > > information there, e.g.:
>> > > http://mitgcm.org/pipermail/mitgcm-support/2014-April/009093.html
>> > > but it's likely that you would need additional diagnostics output,i
>> > > e.g. WTHMASS & WSLTMASS at k=1.
>> > >
>> > > Cheers,
>> > > Jean-Michel
>> > >
>> > > On Tue, Nov 18, 2014 at 09:49:34PM +0530, Abhisek Chakraborty wrote:
>> > > > 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
>> > > > >
>> > > > >
>> > >
>> > > > _______________________________________________
>> > > > 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
>> > >
>>
>> > _______________________________________________
>> > 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/20141121/1265200f/attachment-0001.htm>
More information about the MITgcm-support
mailing list