[MITgcm-support] Heat Budget in tutorial (Global Ocean Simulation at 4 deg)

[Martin Losch]

Hi Krishnakumar,

so you use tutorial_global_oce_latlon for this test. I did not realize that (o:

I cannot see any obvious problems in your configurations, except for maybe this one:
I do not understand, what diagnostics you are actually saving. in data.diagnostics, there is only the RHOanom(aly) and in data you have dumpFreq=1month, which gives you a lot of snap shot output in 800 years (you may have noticed).

In order to check the heat budget I you need to compare the cumulative (or mean*avaraging interval) heatflux into the ocean to the change in heat content (potential temperature) over the averaging interval. In order to get these diagnostics, you need to save (in data.diagnostics)
something like oceQnet (surface heat flux, but without relaxing contribution, as far as I remember), or TFLUX (which contains everything) with a "frequency" of, say, 1year to get yearly averages, then THETA with the same negative frequency (-1year) and a "timePhase" of 0 to get snapshots at the beginning and the end of the this interval, and then compute your spatial averages/integrals from that. 
Alternatively you can use the spatial averages of the diagnostics package like this:
> &DIAG_STATIS_PARMS
> #--- to check conservation of Volume, Heat & Salt (with z*+realFW+AB_2):
> stat_fields(1:,1)= 'TRELAX  ','SRELAX  ',
>                   'surForcT','surForcS','TFLUX   ','SFLUX   ',
>                   'oceQnet ','oceSflux','oceFWflx',
>    stat_fname(1)= 'surfStDiag',
>     stat_freq(1)= 31104000.,
> stat_fields(1:,2)= 'ETAN    ','THETA   ','SALT    ',
>    stat_fname(2)= 'instStDiag',
>     stat_freq(2)= -31104000.,
>    stat_phase(2)= 0.,
> &
and use the corresponding variables to check your balance.

Good luck,

Martin

PS. Completely unrelated to your question are  the following general comments on the tutorial, which is meant to illustrate some of the capabilities of the model, and thus has turned on things that are not necessarily useful for your purposes:
- SIZE.h, on a sinlge processor I'd use sNx =  90 and nSx = 1, that reduces the number of output files.
- I'd turn off the tave-package (packages.conf). You do not use it (according to your "data") and all of its capabilities have been superseeded by the diagnostics package.
- although very effective, the cd_code is not very popular, even with it's developer, because it introduces viscosity of unknow size. "We" now prefer bi-harmonic viscosity to control noise (viscA4).
- turn of the ptracers (usePtracers=.false.), unless you really want to use it (speeds up the code).


On May 20, 2011, at 12:49 AM, Krishnakumar Rajagopalan wrote:

> Hi Martin,
> 
> Thanks a lot for your suggestions.
> 
> I did make some errors in the calculation of the mass_cell. Earlier I considered the thickness of the cell as rf(i,j,k), but now I multiply  rf(i,j,k) with hFacC(i,j,k)  to find the thicknes of the cells. The volume of the  cell is now calculated as rAc(i,k)*rf(i,j,k)*hFac(i,j,k), where rAc(i,j) is the area of tracer cell.  As the temperature and salinity fields in the cells change with time, there are some changes in the  potential density of the cells. I have averaged the densities at time=0 and time = 800 yrs to find the mean density. Now mass of the cell is calculated as volume_cell*mean_density_cell.  I have also changed the value of  Cp from 4000 to 3994 J/kg/k , which is the value the solver is using. These changes have reduced del H by a small amount. Please take a look at attached figure 1. I have also revised the calculation of the surface heat flux. It is now calculated as surfFacT.*rAc and integrated over the ocean surface as sum(sum(surfFacT.*rAc )), at any time, which is then integrated over time, to find Q. ( Q is practically unchanged). Figure 2 shows the ratio of Q/del_H and the ratio increases linearly.
> 
> The attached config.tar file contains CPP_OPTIONS.h (this tutorial uses the default CPP_OPTIONS.h  and is located only in /build) ,all files in /code directory and the data* files.
> 
> It would be great if you could give us some comments.
> 
> Best regards
> 
> Krishnakumar
> 
> ( Sorry about the delay in replying, I was carefully checking my post processing)
> 
> From: Martin Losch <Martin.Losch at awi.de>
> To: mitgcm-support at mitgcm.org
> Sent: Mon, May 16, 2011 10:12:00 PM
> Subject: Re: [MITgcm-support] Heat Budget in tutorial (Global Ocean Simulation at 4 deg)
> 
> Krishnakumar,
> 
> what would be interesting in this context are details of your configuration, e.g. the content of your "code" directory (main the cpp-flags that you use), and your namelists (data* files).
> 
> Wild guess:
> Did you ever plot the ratio of the two curves? From the figure it looks like this ratio might be constant, which would then point to an inconsistency in your analysis (e.g. Cp is no the same as in the model, or there is a problem/systematic inaccuracy with the computation of the "mass_cell", etc.)
> 
> Martin
> 
> On May 16, 2011, at 9:25 PM, Krishnakumar Rajagopalan wrote:
> 
>> Hi All,
>> 
>> 
>> It seems that the heat globally added at the ocean surface does not match the heat absorbed by the ocean.
>> 
>> 
>> Let Ho and Ht denote the initial  enthalpy and enthalpy at time, t,  of the ocean.  Then (Ht – Ho) represents the change in enthalpy , del H, of the ocean which is calculated from  summing  over all the cells, mass_cell*Cp*( T(t) - T_initial  ),  where T_initial  and T(t) represent the initial temperature and temperature at time, t, of the cell.  The quantity, del H,  at any time , must  equal the heat supplied to the ocean till that time,  which in the case of this tutorial comes from the restoring term { -lamba*(theta-SST) }  and Qnet.  I have integrated the diagnostic,  surfForcT (model surface forcing for Temperature), over the area of the ocean and with time to calculate the heat, Q entering the ocean. (The surfForcT is output every month). In (attached) figure, the red line shows del  H as a function of time and the blue line shows Q supplied to the ocean  upto that time, t .  Ideally the two lines must coincide.  What could be the reason for this apparent discrepency which actually increases with time ?
>> 
>> Thanks in advance
>> 
>> Krishnakumar