Ocean Temperature Bug
Details of investigation strange results observed for the Ocean Temperature diagnostic with stash code (30,301). This turned out simply to be a problem with packing, and was solved by switching it off!
1. Problem
Specifically, we have seen:
A. Exceptionally large values in the Red Sea (approximately 0.3×1011 and 0.3×1010) and in the Arabian Gulf (around 0.7×108).
B. Vertical bands of one value extending across the whole grid, even over land points, on rows corresponding to latitudes 15 - 22.5, 27.5 and 30.
2. Observations
I ran a 1 year standard job on Ruby (xccab), with Ocean Temperature included in the Ocean STASH:
* TDAYMNEV: Time mean, meaning period - 1 dump period, sampling frequency - every 2 timesteps with offset 1, and output every dump period
* DO20: Variables on full levels, range of model levels - 1 to 20
* UPMEAN:Dump store with climate mean tag, period 1 and period 2
Some observations from the model output:
1. The large sea values (A) occur in the dump files only.
2. The vertical bands (B) occur in the annual mean files only.
3. There are also exceptionally large values of salinity in the Red Sea and Arabian Gulf.
4. The large sea values extend through the top 5 levels, i.e., the levels where the Red Sea and Arabian Gulf exist.
5. The horizontal bands are in the top 5 levels only and cover only the exact rows on which the Red Sea and Arabian Gulf stand.
3. Salinity
Salinity is of interest since it is used in the UM's STATE_T routine to calculate Ocean Temperature. The point raised in remark 3. is a known issue with FAMOUS and the way it deals with Salinity in the inland seas. However, it is possible to remove the error in the top layer.
In Ocean GCM -> scientific params and sections -> salinity controls, switch on "apply limits to range of surface salinity permitted". This should restrain salinity to lie within the values specified of 0 and 0.045.
I applied this to job xccad and found that the unphysical Salinity and Ocean Temperature values vanished in the top layer, confirming the connection between the two quantities.
4. Further Observations
To investigate Remarks 1 and 2 in more detail, I set up xccad to run for 10 days, producing daily Ocean dumps, and from STASH, daily instantaneous values and daily means of Ocean Temperature. Note that this was with the top layer Salinity correction in place.
To set up the ocean STASH in xccad:
* Daily mean:
o TDAYMNA: time mean, meaning period - 1 day, sampling frequency - every timestep, output - every day
o DO20: Levels 1-20
o UPO: PP file, output unit 65
* Instantaneous daily values:
o TDAILY: no time processing, output - every day
o DO20
o UPO
The model output showed the following:
*
Daily dumps:
After the first day, there are large values in the Red Sea and Arabian Gulf (A) across all top 5 layers. The following 9 days show normal values on the top layer and large values on the 4 layers below.
*
Daily instantaneous values:
All 10 days show a normal top-layer and vertical bands (B) on the following 4 layers.
*
Daily means:
The first day shows vertical bands on all top 5 layers and a normal top layer for the next 9 days.
We see that the salinity correction fixes the vertical bands in the top layer, confirming the connection between large values of Salinity and Ocean Temperature in the Red Sea and Arabian Gulf, and vertical bands of Ocean Temperature extending across continents.
The only exceptions are the first dump and first daily mean which could be explained by large salinity values in the start file. This would cause the results of the first timestep to be flawed, and as there are two Ocean timesteps per day, this would affect the first daily mean, but the instantaneous value after 1 day would be correct. The first daily dump is puzzling as I expected this to be the same as the instantaneous values but then I do not know exactly what the dump represents.
5. Debugging
In order to see where the unexpected land values first appear, I put printed out parts of the Ocean Temperature array at various points in the code.
The calling tree of the UM code with storage of the Ocean Temperature array in that routine:
ROW_CTL stashwork(si(tempitem,30,im_index))
-> BLOKCALC Temperature(*)
-> BLOKCNTL Temperature(*)
-> ROWCALC Temperature(imt_stash,jmt,km) and Temp_row(imt,km)
-> STATE_T Temperature(imt,km)
The actual calculation of temperatures is done in STATE_T, one row at a time. I put PRINT statements in ROWCALC just after the main J loop where the diagnostic array rows are calculated.
The diagram below illustrates the part of the Ocean Temperature array that I was interested in, with Fortran and IDL indexes (in fact, I printed these entire rows). The green boxes represent the Arabian Gulf and the purple boxes represent the Red Sea.
30.0 48 50
27.5 47 49
25.0 46 48
22.5 45 47
20.0 44 46
17.5 43 45
15.0 42 44
latitude idl fortran 11 12 13 14 15
idl 10 11 12 13 14
longitude 37.50 41.25 45.00 48.75 52.50
After the 1st timestep there are large values in the Red Sea and Arabian Gulf over all top 5 layers, whereas as before, after the 2nd timestep the top layer looked normal and the lower 4 layers contained .
There were, in fact, no missing data values at all, as these are added in ROW_CTL by MASKODIAGN. Therefore next, I printed the sub-arrays after the Ocean mask was applied in ROW_CTL and found that the expected land points were marked as missing data.
Nothing further is done to the data in the stashwork array before it is processed by the STASH routine. Consequently it seems the vertical bands are created by some error in STASH, which would explain why they show up in the instantaneous output and mean files (defined in the STASH panel), but not in the restart dumps.
Appendix
A. Ocean Diagnostics
From an email written by Jonathan Gregory explaining how to add an ocean diagnostic to the UM code. It usefully explains where the diagnostics are calculated and how diagnostic arrays are passed between routines and into the stash workspace.
The stash workspace for section 30 is dimensioned as local array in ROW_CTL, at line ROW_CTL.133. (Section 30 is all the ocean model apart from sea ice, which is section 32. Stash codes are section*1000+item.) It is passed into lower-level routines in separate pieces for each diagnostic. For instance, I added the diagnostic for total 3D ocean velocity (barotropic+baroclinic), with the following lines in ROW_CTL
integer utotitem ! Stash item number for total velocity OJG5F401.37
parameter(utotitem=320) OJG5F401.38
This is so that the stash code (30320 for u-velocity, 30321 for v) isn't hard-coded in more than one place.
&,stashwork(si(utotitem,30,im_index)) OJG5F401.40
&,stashwork(si(utotitem+1,30,im_index)) OJG5F401.41
&,sf(utotitem,30),sf(utotitem+1,30) OJG5F401.42
in the CALL BLOKCALC, which is the next routine down. This passes BLOKCALC the stash workspace and stash flags for the diagnostics
if (sf(utotitem,30)) OJG5F401.43
& call maskodiagn(imt,uv_j_dim,km,.true.,0.,fkmq ORH0F404.36
&,stashwork(si(utotitem,30,im_index))) OJG5F401.45
if (sf(utotitem+1,30)) OJG5F401.46
& call maskodiagn(imt,uv_j_dim,km,.true.,0.,fkmq ORH0F404.37
&,stashwork(si(utotitem+1,30,im_index))) OJG5F401.48
to put missing data at the inactive points. You could avoid this if you put missing data at those points in the diagnostic as you calculated it. If you want to use this procedure, you would use jmt and fkmp for the TS-grid, whereas my diags are on the UV-grid.
In the SUBROUTINE BLOKCALC statement, we have
&,utot,vtot,sfutot,sfvtot,Temperature,SFTemp OMB1F404.160
and they are dimensioned
real OJG5F401.31
& utot(*),vtot(*) ! Total 3D velocity diagnostics OJG5F401.32
&,Temperature(*) ! actual temperature diagnostic OMB1F404.161
logical OJG5F401.33
& sfutot,sfvtot ! Stash flags for 3D velocity diagnostics OJG5F401.34
They are passed to BLOKCNTL
&,utot,vtot,sfutot,sfvtot,Temperature,SFTemp OMB1F404.163
in which they are dimensioned and passed in just the same way into ROWCALC. In ROWCALC they are dimensioned as 3D arrays
& utot(imt_stash,jmtm1,km) ! Total u-velocity diagnostic OJG5F401.4
&,vtot(imt_stash,jmtm1,km) ! Total v-velocity diagnostic OJG5F401.5
Yours would be jmt rather than jmtm1. The diags are actually calculated in ROWCALC and that would be appropriate too for the speed of sound. The tracer evolution is done by TRACER, which is called from ROWCALC, but a purely diagnostic calculation doesn't have to be in there. My diags don't really have to be worked out, since the numbers are already available, but you would need a similar loop:
IF (SFUTOT) THEN OJG5F401.8
DO K=1,KM OJG5F401.9
DO I=1,IMT_STASH OJG5F401.10
UTOT(I,J,K)=U(I,K) OJG5F401.11
ENDDO OJG5F401.12
ENDDO OJG5F401.13
ENDIF OJG5F401.14
IF (SFVTOT) THEN OJG5F401.15
DO K=1,KM OJG5F401.16
DO I=1,IMT_STASH OJG5F401.17
VTOT(I,J,K)=V(I,K) OJG5F401.18
ENDDO OJG5F401.19
ENDDO OJG5F401.20
ENDIF OJG5F401.21
ROWCALC deals with just one row (row J), so these loops copy the values from the 2D lon-depth "slab" for row J into the 3D diagnostic array.
Finally, you would need to add the diags to the call to OSWAPDIAGS near the end of ROWCALC. I think you can see what to do there; it's the same for all diags. This routine propagates information in the "halo" of each processor to the other processors.
The temperature "slab" at the present time-level is in T(1:IMT,1:KM,1) and the salinity in T(1:IMT,1:KM,2). "T" is for "tracer". The model uses a leap-frog timestep. After TRACER has been executed, the tracers at the next time level are in TA and those in the previous one in TB, but you probably don't need those. I expect you need the depths too, don't you? Perhaps you can locate those arrays, whose names I can't remember, by looking in TRACER and its subroutines, but if not I can have a look. There are arrays for the mid-layer, layer bottoms and layer thicknesses, I think.