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!
Specifically, we have seen:
I ran a 1 year standard job on Ruby (xccab), with Ocean Temperature included in the Ocean STASH:
Some observations from the model output:
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.
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:
The model output showed the following:
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.
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.
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.