Jump to latest

Introducing the coupled Hines scheme at version 4.5 of the portable Unified Model

This is a continuation of the work of Bryan Lawrence, documented in Bryan's gravity wave notes.

19/09/03
Run xacvw has completed successfully. This uses UM 3A gravity wave drag code up to 50 km and Rayleigh friction between 50 km and the model top at 80 km. Experiment was run for one month and will be a useful control run for comparison with the Hines scheme. This run may be extended at a later date.

Run xacvx submitted to run over the weekend. This calls the Hines code but does not allow updating of the model u and v fields. This allows the code to be tested for bugs, but because the rest of the model does not see the Hines accelerations it is equivalent to running without any gravity wave drag. Rayleigh friction is turned on above 50 km. The results of this run will be compared with xacvw.

26/09/03
Run xacvx completed one month successfully. However, realised that Hines orographic waves were turned off in xacvx (i.e., surface stresses were set to zero before calling Hines) so this run only tests the spectral part for bugs.

Set up experiment xacvy which is identical to xacvx except that orographic waves are turned back on, so this should test the full coupled Hines scheme for bugs in the code. Submitted to run over the weekend.

01/10/03
Run xacvy completed one month successfully. The success of runs xacvx and xacvy is good evidence that there are no more bugs in the Hines code itself which would cause the model to crash. However, currently the model crashes when Hines is allowed to update the winds (experiments xacvb-e). According to the model output the symptoms of all the crashes are:

Model completed with the following :
Error Code : 1
Message : QT_POS : MASS-WEIGHTED QT SUMMED OVER A LEVEL WAS NEGATIVE.

and the "fort6" output files for some processors contain the following:

*** SETFIL : WARNING *** Filtering entire model domain.
Timestep should be reduced.
Filtering has been switched off.
FILTERING SWITCHED OFF.
See earlier in output for error message.

Looking at the diagnostics from runs xacvb-d the most obvious explanation for the error messages is the enormous accelerations at some grid points in the top model level. For example, see the plots xacvd_u_map and xacvd_u_profile. There are a number of ways in which the top level accelerations could be reduced:

However, all these will affect the circulation lower down which may be undesirable. Also, setting tflag=.false. was tried in run xacve with the result that the model still crashed with the same error messages and very large accelerations occurred in at least the top two levels.

The next steps are to:
(a) check whether the convection scheme is switched on in the mesosphere;
(b) set up a run to try Claudio Piani's suggestion of having a very small source term of 0.1 or even 0.01 to test for numerical problems.

03/10/03
Jamie Kettleborough informs me that the convection scheme loops around the number of wet levels specified in the model. In the standard L64 runs there are 60 wet levels, thus there is no convection occurring in the top 4 levels. (This was done because Alan and Jamie found that with 64 wet levels the model crashes every 10 years or so for no obvious reason, which is a serious problem for long climate runs. However, for the purposes of assessing the Hines scheme my runs will be approx. 1 year so hopefully this issue won't arise).

Experiment xacvf was submitted with 64 wet levels. Otherwise parallel to xacvd, hence no orographic component. Critical relative humidity and moisture diffusion parameters for levels 61-64 were copied from level 60, since they were constant anyway from level 22 upwards. Crashed on timestep 616 (12 days and 20 hours) with usual QTPOS and filtering error messages.

Spoke to Alan re. correct values for diffusion coefficients - apparently top layer should be more spongy than layers below (i.e., 1st order diffusion and smaller coefficient). In all runs to date this has not been the case. Not altering the diffusion parameters on the top level can lead to instabilities and model crashes, although usually after a few months (according to Alan). This may not fix the problem with the Hines experiments, but it's still better to have the right values.

Resubmitted xacvf to run over weekend with diffusion parameters copied from one of Alan's runs with 64 wet levels.

08/10/03
Some progress! The run of xacvf with corrected diffusion parameters ran for 1188 timesteps (24 days and 18 hours) so almost doubling the run length of the previous attempt. However, this is not the final cure as the model still crashed with QTPOS and filtering errors.

10/10/03
Due to the comparative success of xacvf this was taken as the basis for two runs to test Claudio's idea of a very small source term. Experiment xacvg had SOURCE =0.1 and xacvh had SOURCE=0.01 but were otherwise identical to xacvf. Both these experiments ran for 698 time steps (14 days and 13 hours) before falling over with the usual errors.

Discussed latest results with Bryan who points out that the same errors keep arising because the CFL criterion is being broken in all the experiments, albeit for different reasons. If the source term is large and the gravity waves deposit a lot of momentum at the top of the model this results in large accelerations and large values of zonal wind, hence the CFL criterion is broken. If the source term is small and the gravity waves provide little or no drag at upper levels then the polar vortex may become too intense and the CFL criterion will be broken.

The most dramatic recent improvement in run length resulted from changing the moisture diffusion parameters in the top wet level, therefore the next logical step is to try tuning the diffusion rather than the gravity wave source term. Bryan has suggested two ways of doing this:

Two experiments were submitted to run over the weekend. Xacvi is identical to xacvf (i.e., only level 64 diffusion has been tuned) but the orographic gravity wave component has been included so that the full coupled Hines scheme can be tested. Xacvj is a copy of xacvi but the level 63 diffusion has been tuned to look similar to the theta diffusion profile.

15/10/03
Experiment xacvi crashed after 544 timesteps (11 days and 8 hours) and xacvj crashed after 435 timesteps (9 days, 1.5 hours) both with QTPOS and filtering error messages. It should be noted that the previous longest run obtained with the full coupled Hines scheme was xacva which crashed after only 2 days, so these latest runs represent progress!

17/10/03
Experiment xacvk was submitted with del squared moisture diffusion in level 63 and a coefficient that was midway between the values used in levels 62 and 64. This ran for 554 timesteps (11 days and 13 hours) so is marginally the longest integration so far.

Experimented with changing the level 63 u, v and theta diffusion coefficient in job xacvl. Tried a variety of values and found that they all resulted in earlier crashes than xacvk, so this remains the best run to date. Learned that it is very easy to make the model performance worse and very difficult to make it better when changing these parameters. From this I conclude that the u, v, theta diffusion parameters are already optimised for this model configuration and there is little to be gained from adjusting them. However, it may be possible to obtain some further improvement by continued tuning of the moisture diffusion.

24/10/03
Used experiment xacvm to explore further the tuning of the moisture diffusion parameters in levels 63 and 64. Started with the same values as xacvk and tried a number of small adjustments to the level 63 diffusion coefficient only, since the value used in xacvk was chosen somewhat arbitrarily, but found that none of the values produced runs as long as xacvk. Then tried tuning the coefficient in levels 63 and 64 by incrementing them both gradually so as to maintain a similar profile to xacvk. Found the values of kappa=2.8e7 in level 63 and 4.1e6 in level 64 allowed the model to run for 657 time steps (13 days, 16.5 hours) before crashing.

Experimented again briefly with tuning the u, v, theta diffusion parameters in levels 63 and 64 while retaining the q diffusion values used in xacvm. As before, found that even small increments make the model crash much sooner. Job xacvm therefore remains the longest run to date.

14/11/03
Produced modset /home/tolkien/alisonp/mods/japmods/rayleigh_l64_3a to add Rayleigh friction to the model between 60 and 70 km in conjunction with the coupled Hines gravity wave drag. The gravity wave drag is still allowed to penetrate to the top of the model. The Rayleigh profile is based on Bryan's code in v18a_sr96.upd (in /home/tornado/lawrence/smm/src/v18). The actual function used is
ARAYL=1E-7+3.0E-6*(Z-70.)**2/100.
which gives the same drag as Bryan's code at the top of the model but a faster increase in drag over the layer in which Rayleigh friction is applied.

Experiment xacvo (parallel to xacvm but including the Rayleigh friction mod) was left to run over the weekend.

13/02/04

  1. Trying to get to the bottom of "segmentation faults" in some of my recent runs.

    N.B. It seems that this problem only affects runs that use Mike Palmer's Hines code. Only need to solve this problem if it becomes necessary to do long runs with Mike's code.
    The time profile for diagnostics tagged for climate meaning is set up to produce 5 day means at 5 day intervals. This is consistent with restart dumps every 5 days and should produce a correct 30 day mean at the end of the month. However, it produces a segmentation fault. Need to check again my profiles, but there could be a problem with STASH processing. Mike's code reuses some of the UM 3A GWD scheme diagnostic storage, so the problem may well occur in STASH.

  2. Modset hines-call10 fully turns off the orographic waves in the coupled Hines scheme by setting the amplitude, argument and the individual x and y components of the surface stress to zero before calling the spectral code. (Hines-call8 sets only the amplitude and argument of the surface stress to zero so doesn't properly turn off the orographic waves).

    Set up experiment xadwu. This is a copy of xadwq (in turn a copy of xacvo). It uses the coupled Hines scheme with the orographic waves COMPLETELY off (i.e., Hines-call10) and tuned Rayleigh friction. This run should now be regarded as our "best ever" (xacvo) minus the orographic waves. It also has corrected climate meaning. Left to run over the weekend.

18/2/04
Checked again the circumstances in which Mike Palmer's code crashes when the climate meaning is corrected. The problem seems to have disappeared! Reran xadws (Mike's code, stress diagnostics, corrected climate meaning). Last time, it crashed with a segmentation fault in the first timestep but now seems to be working perfectly. Suspect possible tolkien problems may have been responsible for earlier crashes. Left to run overnight.

Experiment xadwu completed 30 days satisfactorily. Still need to check the diagnostics to see if it looks sensible. Bryan suggests trying to run this one on.

19/2/04
Run xadws completed 30 days without problem - correct climate meaning now works with both versions of the Hines code.

Checked the ancillary files used in xadwu before submitting longer run. Seems to be using AMIPII zonal mean ozone climatology and AMIPII SSTs and sea ice. Most other ancillaries seem to have sensible settings (the majority are "not used"). For soil moisture and snow depth an ancillary file is named but the fields are marked "not used", so they must be coming from the start dump and then evolving freely. Not sure if this is correct for an AMIPII run or whether there should be some updating. Almost certainly not important for looking at gravity waves in the stratosphere but might be as well to check the correct setup of AMIP style runs for future reference. An additional moisture ancillary file is named as a multi-level user ancillary. Looks as though it was prepared by Jamie and probably goes with his stratospheric moisture mod in tolkien:/soft/model/um/local/vn4.5/mods/l64/qrad_old_4.5.

Continuation run of xadwu submitted - target run length is 3 years and 3 months which will give us four DJF seasons and 3 JJAs. Ran until time step 1837 (38 days, 6.5 hours) and before crashing with normal QT_POS and filtering error messages.

Decided to try running on Palmer code (xadws) to 3 years and 3 months - rather than just take it on trust that this code is stable!

20/2/04
Experiment xadws ran until time step 1714 (35 days, 17 hours) then crashed with QT_POS and filtering error messages, so the Palmer code seems to be suffering from similar stability problems to our version of Hines.

Checked orographic stresses and accelerations from xadwu at end of days 1 and 10 by plotting zonal means and a few individual levels. Both quantities seem to be zero everywhere which confirms that the orographic waves are switched off correctly by modset hines-call10.

Went back to my original control run, xacvw (3A GWD plus standard L64 UM Rayleigh friction above 50 hPa, no Hines scheme). Corrected diagnostic list to get all quantities from 3A scheme and corrected climate meaning.  Left to run over the weekend.

25/02/04

  1. Control experiment xacvw completed one month successfully.  Submitted CRUN with target end date of 3 years and 3 months.
  2. When either version of the Hines scheme is included the model crashes with filtering and QT_POS error messages.  However, we know that Mike Palmer completed runs of 24 years at vn4.4 using his code.  The implication is that the precise details of the Hines scheme may not be responsible for the failures and that they are caused by another part of the UM.  Modsets directly affecting the filtering are an obvious candidate.  Sought advice from Alan Iwi on appropriate modsets.
  3. Alan compared his  L64 run, xabke, which ran for 10 years, with my xadwu.  There were very few significant differences between the jobs, and none that affected the filtering, but Alan suggested three specific measures that may improve the stability of my runs.
    1. Change back to 60 wet levels OR include Alan's mod conv60 (in ~iwi/um/vn4.5/mods on tolkien) which allows 64 wet levels but limits convection to the lowest 60.  Alan has never been able to produce a stable run with convection on 64 levels.
    2. Switch on the modset lux_32bit.mod to fix some very small numbers that are hard coded in the model.
    3. Add the hand edit file mkghganc (in window SubModelIndep --> PostProc --> LocalPostProc) to alter the way the greenhouse gases are used from the ancillary files.
  4. Decided to try Alan's suggested changes in a job with Mike Palmer's code, since we know that Mike's code should be capable of producing a stable run.  Decided to use 60 wet levels, rather than conv60, as then the run will be parallel to my xacvw control.  Job id is xadwt.  Left to run overnight.
26/2/04

  1. Control run xacvw has completed 7 months so far and is still running.  Therefore, my control run  is stable.
  2. xadwt crashed on timestep 1272 (26½ days).  Look at differences between my control and xadwt:

  3. Job xacvw Title 3A orographic GWD + Rayleigh, NO Hines call
    Job xadwt Title As xadws, changes to q, ghgs near top

    Difference in window subindep_Compile_Mods
     -> Model Selection
       -> Sub-Model Independent
         -> Compilation and Modifications
           -> Modifications for the model
    Differences in Table Fortran mods
     4c4
    <  $MOD64/gsdev20_l64_3a Y
    ---
    >  $JAP_MODS/gsdev20_l64_hines_mpalmer Y
    19c19
    <  $LOCAL_MODS/lux_32bit.mod N
    ---
    >  $LOCAL_MODS/lux_32bit.mod Y
    23,24c23,24
    <  $JAP_MODS/hines-call8 N
    <  $JAP_MODS/hines-cbnl5 N
    ---
    >  $JAP_MODS/gsdev44_l64_4_norf Y
    >  $JAP_MODS/hines-um44-v8-xacsh-stress Y

    Difference in window subindep_PostProc_ScrLocal
     -> Model Selection
       -> Sub-Model Independent
         -> Post Processing
           -> Local post-processing scripts
    Differences in Table Local post-processing scripts
     1a2
    >  /usr/local/um/local/vn4.5/hand_edits/mkghganc Y

    Apart from Hines and Rayleigh mods, the only differences are in the additional hand edit and the lux_32bit mod suggested by Alan.
  4. Decided to remove the last two changes from xadwt so that the experiment is exactly parallel to my control with the only remaining differences being the removal of Rayleigh friction and the addition of the Hines spectral code:

    Job xacvw Title 3A orographic GWD + Rayleigh, NO Hines call
    Job xadwt Title As xadws, changes to q, ghgs near top

    Difference in window subindep_Compile_Mods
     -> Model Selection
       -> Sub-Model Independent
         -> Compilation and Modifications
           -> Modifications for the model
    Differences in Table Fortran mods
     4c4
    <  $MOD64/gsdev20_l64_3a Y
    ---
    >  $JAP_MODS/gsdev20_l64_hines_mpalmer Y
    23,24c23,24
    <  $JAP_MODS/hines-call8 N
    <  $JAP_MODS/hines-cbnl5 N
    ---
    >  $JAP_MODS/gsdev44_l64_4_norf Y
    >  $JAP_MODS/hines-um44-v8-xacsh-stress Y

    Resubmitted xadwt overnight.

27/2/04
  1. Control run xacvw has completed 1 year so far and is continuing.
  2. Run xadwt crashed again after 1271 timesteps (26½) days exactly as before, so the lux_32bit.mod and mkghganc changes have made no difference to the stability.
  3. So, we have an experimental setup that can run stably, as demonstrated by the xacvw control, and we have Mike Palmer's Hines scheme mods that Mike was able to use in runs of 24 years in length.  However, when they are combined the model becomes unstable.
  4. What else is different about my runs and those done by Mike?
  5. Following discussion with Bryan and Alan it was agreed that the following actions, listed in order of priority, should be pursued:

21/5/04
Notes from L64 meeting held at Reading University 20/5/04

Present: Lesley Gray (meeting organiser) (CGAM), Lois Steenman-Clark (CGAM), Peter Braesicke (DAMTP, Cambridge), Joanna Haigh (Imperial College), Andrew ? (Imperial College), Wenshou ? (Leeds University), Alison Pamment.

  1. Lois explained the process of migrating L64 from Turing to Newton. On Turing there were at least three distinct flavours of L64:
  2. Newton (SGI Altix) is a little-endian machine so it requires UM4.5+fixes in order to run. Current issues are:
  3. One of the main aims of the meeting was to agree a "base version" of the model that Lois will concentrate on building. The idea is to have one, rather than three, basic flavours on Newton that everyone can use as the starting point for building their own experiments. The following points were agreed:
  4. It was agreed that, once this model has been built, a 130 year control integration will be run which will hopefully save everyone having to do separate long control runs. Not all the details of this run have yet been agreed, but the points that have been agreed are as follows:
  5. At Lesley's suggestion, the same group will meet monthly to discuss progress with building the model and (in future) some science! The next meeting will be held at Reading on Wednesday 23rd June, 1330-1530.
28/5/04
Next steps agreed today with Bryan: 01/07/04
Plots for experiments xafba and xafbb.

17/09/04

  1. Plots of xafba show very large u accelerations at 0.03 hPa (third model level from top), in excess of 10,000 ms-1day-1, in one or two grid points on the west coast of Norway at 0Z 04/12/1978, i.e., after the model has run for 72 hours.
    Hines coupled U Acc. at 0.03 hPa 04/12/1978 00:00 Z
    This is caused by the orographic component as shown below.
    Hines Orographic U Acc. at 0.03 hPa 04/12/1978 00:00 Z
    There is very little contribution from the spectral acceleration in these grid points (see main results page for plot). The model does not crash at this stage in the integration - it falls over part way through the eighth day. However, large accelerations of the kind seen here may be responsible for some of the crashes we have seen and in any case this behaviour is unphysical so it needs to be understood.

  2. UM values for the grid point where the U acceleration is highest can be fed into the testbed software to try to see exactly what is going on. The testbed requires values for standard deviation of the orography, xx, xy and yy gradients of the orography, pstar and profiles of U and T for the point in question. The maximum U acceleration occurs at the grid point located at 65 degrees N, 7.5 degrees East. A zoomed plot is shown below.
    Zoomed Hines coupled U Acc. at 0.03 hPa 04/12/1978 00:00 Z
    However, the grid point of interest appears to be entirely over the sea. The land-sea mask (not shown) agrees that this grid point is sea, as does the orography. The following plot shows the standard deviation of the orography for the zoomed region.
    Zoomed sd of orography
    .

  3. I need to check where in the code the winds, stresses and accelerations are regridded between the U-grid and the P-grid to see if I can explain this discrepancy between the location of the orography and the location of the maximum acceleration. It may also be worth reviewing the way the grid is handled for the parallel processing - maybe I am a grid point out somewhere. This needs to be resolved so that the correct values can be fed into the testbed to diagnose the reason for the large acceleration.
26/10/04
Plots for experiment xafbe and xafba-xafbe differences.

The main conclusion to draw from the plots is that the U winds at 0.03 hPa are far more zonal in appearance in the control run (xafbe) than in the Hines run (xafba). In the control run the winds are westerly throughout the northern hemisphere. In the Hines run "bullseyes" of strong easterly winds are clearly visible over northern hemisphere land; they are embedded within the generally westerly flow leading to regions of strong horizontal wind shear. The easterlies tend to occur at 60-70 degrees N and may well be the cause of the filtering error messages which accompany the model crashes.

02/11/04
Profiles for experiment xafba for a single grid point over northern hemisphere orography show the development of what appears to be 2 grid length noise at upper levels in the coupled Hines zonal acceleration and its subsequent effect on the zonal wind. The profiles are plotted every 2 timesteps (i.e., hourly) for the first six hours of the run. The problem occurs initially at around 65 km but quickly affects both higher and lower levels. By timestep 12 this gives rise to noise in the zonal wind from approximately 55 km to the top of the model. Clearly there is strong vertical wind shear in this grid point in addition to the strong horizontal shear noted on 26/10/04. It is no wonder that the model cannot survive this!

A more detailed analysis of timesteps 2 , 4 and 6 (the first three hours of the run) shows that it is the orographic component that determines the nature of the total acceleration and the contribution of the spectral component is negligible by comparison. One point to note is that the total zonal momentum flux is roughly three orders of magnitude less than the orographic flux. This should not be the case, so need to check where in the code this is happening (the orographic flux seems very big).

01/12/04
As noted on 17/9/04, Bryan's gravity wave testbed can be used to investigate the reason for the unrealistic orographic accelerations over the northern hemisphere. Also noted was the apparent mismatch in location between the maximum orographic acceleration and the orography itself which makes it difficult to know precisely which model profiles of U and T correspond to a particular set of orographic values. Another manifestation of this problem is the band of zero zonal accelerations around the equator in the first of the 17/9/04 plots. The problem appears to be that in the part of the model domain being processed in the top row of the logical processor array the atmospheric quantities are being shifted two grid rows to the north of their correct position. The problem is visible in Hines diagnostics output on both the P/T and U/V grids.


Diagnostics on P/T grid

Hines zonal orographic acceleration
Hines meridional orographic acceleration
Hines zonal spectral acceleration
Hines meridional spectral acceleration
Hines zonal momentum flux
Hines meridional momentum flux
Hines zonal orographic momentum flux
Hines meridional orographic momentum flux

Diagnostics on U/V grid

Hines zonal acceleration
Hines meridional acceleration
U component of wind after GWD
V component of wind after GWD





To diagnose the cause of the location error I performed a number of single timestep runs on 1, 4 and 8 processors under experiment id XAFBJ.
Total no. of processors

1
4
8

No. of processors in NS direction

1
2
4

No. of processors in EW direction

1
2
2

The stripe of zero values in the coupled Hines zonal acceleration can clearly be seen around the equator in the 2x2 processor version of XAFBJ and around the South Pole in the 1x1 processor version (see results page). This supports the view that the location error in the accelerations is occurring in the top row of processors, rather than in the remaining rows. The most compelling evidence, however, comes from comparing the accelerations from the 2x2 and 4x2 processor runs with the orography. The first plot below shows the Hines zonal acceleration (U/V grid) over Europe for the 2x2 run and the second shows the same for the 4x2 run. The third plot shows the standard deviation of the orography (P/T grid) for the same area.

XAFBJ 2x2: Hines coupled U Acc. at 0.03 hPa
XAFBJ 4x2: Hines coupled U Acc. at 0.03 hPa
UM Standard Deviation of Orography
The most important area to note is the Black Sea (42 - 45 degrees N, 30 - 38 degrees E) where there are clearly no orographic data in a block of 6 adjacent grid boxes. The 2x2 run shows a block of 6 grid points having appreciable gravity wave acceleration apparently centred over the Black Sea. In the 4x2 run these same 6 acceleration values are shifted two grid rows to the south and correspond much more closely with the location of the orography. The slight northeastwards offset of the accelerations compared to the orography is due to the staggering of the U/V and P/T grids. Also visible in the 4x2 run are two grid rows of zero acceleration in the vicinity of 45 degrees N. Everything to the south of this latitude is correctly located with respect to the orography while everything to the north is shifted northwards by two grid rows. The northward shift can most clearly be seen near Denmark where accelerations occur over the sea despite there being no orographic values. In the 4x2 configuration 45 degrees N is the southern extreme of the domain handled by the top row of processors. Therefore the location error is definitely associated with the top processor row, rather than with all the other rows.

There are two possible causes of the location error - regridding or domain splitting between the processors.

21/12/2004

The location error is cured!

A detailed analysis of the way in which the gravity wave code handles the model grid was carried out. The location error was caused by logic pertaining to the top row of processors (i.e., the domain splitting) and the way this relates to the arrays used to output diagnostics from the Hines scheme. For parts of the model domain being handled at the top of the logical processor group the first two rows of u, v, theta, exner pressure and q (i.e., the prognostic variables) are "removed" in the top level control routine GWAV_CTL. The logic is as follows:

  SUBROUTINE GWAV_CTL...
  .
  .
  .

  FIRST_POINT=START_POINT_INC_HALO ! This is the first non-polar point of the
                                   ! field.  It includes the EW halo.
                                   ! If (at_top_of_lpg) this excludes the north
                                   ! halo row (garbage) and the North Pole,
                                   ! otherwise it includes all rows of the array

  LAST_POINT=END_P_POINT_INC_HALO  ! This is the last non-polar point of the
                                   ! field.  It includes the EW halo.
                                   ! If (at_base_of_lpg) this excludes the South
                                   ! Pole and the south halo row (garbage),
                                   ! otherwise it includes all rows of the array

  POINTS = LAST_POINT - FIRST_POINT
  ROWS = POINTS / ROW_LENGTH
  JS = FIRST_POINT - 1
  .
  .
  .

  CALL GWAV_INTCTL(
 &     D1(JPSTAR+JS),D1(JP_EXNER(1)+JS),D1(JTHETA(1)+JS),D1(JQ(1)+JS),
 &     D1(JU(1)+JS-ROW_LENGTH),D1(JV(1)+JS-ROW_LENGTH),
 &     P_FIELD,U_FIELD,                    ! P_FIELD and U_FIELD are always
                                           ! dimensions of full arrays handled
                                           ! by an individual processor

 &     ROWS,                               ! If (at_top_of_lpg)
                                           ! or (at_base_of_lpg) ROWS will be
                                           ! full array minus two rows,
                                           ! otherwise it will be total rows
 &     ROW_LENGTH, ... , ARGFLDPT data,
 &     D1(JOROG_SD+JSL), more orographic data, ...)

In contrast to the prognostic variable arrays, the arrays used for diagnostics are always passed between code levels as the whole of the domain being handled by a particular processor. This means that, when copying to the diagnostic arrays on processors in the top row, an offset of two rows must be applied before any diagnostic data are written. Failure to do so will cause the last two rows of the diagnostic array to be left unfilled, resulting in the two rows of zeros seen in many of the plots.

The offset problem is avoided at the South Pole because for processors on the bottom row the whole prognostic array is passed to the lower level routines. The code doesn't process the South Pole or the unused row beyond because of the way the ROWS variable is calculated (see above). Although the whole array is passed only ROWS number of rows are actually processed. For processors in between the top and bottom rows the whole array is both passed and processed.

The plot below shows an eight processor (4 NS by 2 EW) run, XAFBK, which is identical to XAFBJ except that a 2 row offset has been applied to the diagnostics in the top row of processors. Comparing this to the 1/12/04 plots shows that the locations of the accelerations and orography are now in agreement at all latitudes. Look particularly at the Black Sea and the area between Denmark and Norway.
Hines coupled U acceleration 0.03 hPa

17/01/2005

The offset of two rows needs to be applied to all Hines scheme diagnostics in both the gwave and gwvert subroutines. (For the updated u and v winds the offset is already taken into account in the call to p_to_uv so it doesn't have to be applied when filling the output array). Results for a 4 processor run, XAFBL, show that the stripe of zero values around the equator has been eliminated from all the diagnostic fields.

On 2/11/04 the question was raised as to why the total momentum flux is apparently several orders of magnitude smaller than the orographic flux. In fact the diagnostic labelled as "momentum flux" in the list of STASH diagnostics is the spectral momentum flux and the total value would be obtained by adding the two components. The STASHmaster file has been amended to label the diagnostics more clearly.

26/01/2005

A version of the testbed software that reads orography and pstar data from an input file (rather than being hard coded) is working and tested. It will work with Bryan's original test cases aswell as UM data.

The single timestep run XAFBL shows an area of large easterly U accelerations close to the St. Lawrence River in Canada (see below).

Hines coupled U acceleration 0.03 hPa

The maximum acceleration has a magnitude of 10,949.32 m/s/day and several grid points have accelerations in excess of 5000 m/s/day. The orography at the points of greatest acceleration has a height of 200 - 600 m and standard deviation of 100 - 200 m, so this is not an area of unusually high or rough orography.

Orographic height

Standard deviation of orography

The testbed was run for the point of maximum acceleration (294.375 degrees longitude, 56.25 degrees N latitude on the U/V grid and 292.5 degrees longitude, 55 degrees N latitude on the P/T grid). See the testbed results. The main point to note is that the orographic momentum flux is between 10 and 100 times larger than would be expected. The testbed code needs to be checked to see exactly how these values are being calculated. I need to check the literature for observed values.

The results were compared to UM profiles for the same grid point. The orographic fluxes look similar in the testbed and the model but the spectral flux profiles seem very different. Why is this? I have imported hines-cbnl5 from my UM mods directory back into the testbed, so everything below the level of HNSWRAPC should be running identical code. The top control level of the testbed is very simple and does not do any processing to the fluxes or accelerations. Therefore the differences must occur in the GWAVE3A or GWVERT3A subroutines. GWAVE3A does no processing of the fluxes, so GWVERT3A seems the most likely place to start to look for differences between the codes.

14/02/2005

The model and testbed would not be expected to produce identical results if only a single u profile is fed into the testbed. In the model, the wind is first regridded from the U/V to the P/T grid which requires averaging the four wind profiles surrounding the orography grid point. I have not been doing this with the wind profiles input to the testbed. Also, in the model the total (i.e. orographic + spectral) U and V accelerations are regridded back to the U/V grid before being applied to the winds. This means the testbed would have to be run four times to produce accelerations at four orography grid points so that these could be interpolated onto the U/V grid to calculate a final acceleration that would be consistent with the UM. I have not been doing this with the testbed. See below for a schematic representation of what is done in the UM to calculate the acceleration at a single U/V grid point such as the one in the centre of the diagram. Note that the consequence of the interpolations backwards and forwards between the two horizontal grids is that the acceleration at the central U/V point depends not only on the wind profile at that point but also on the wind at the 8 nearest neighbours. The gravity wave code is therefore not strictly working on a single column.

            X                        X                        X


                        O                        O

                                      (u/v)
            X                        X                        X

                         (p/t)
                        O                        O


            X                        X                        X

01/03/2005

Analysis of XAFBL (1 timestep run with latest Hines version) compared with XAFBG (1 timestep run of control (3A + Hines)) shows that the orographic flux profiles for a grid point near the St. Lawrence River are very different.

XAFBL:
Orographic u mom flux
XAFBG:
Orographic u mom flux

Apart from the differences above 60 km, the fluxes near the surface are different. At the lowest level the surface stresses are very similar (about 0.13 N m-2). They are not identical because the diagnostic from the control run has been regridded to the U/V grid and has therefore undergone some horizontal interpolation, whereas the Hines diagnostic is output on the P/T grid.

Over the next three levels the surface stress is reduced by approximately one third in the control run, whereas it actually increases in the Hines run. The control run includes the parametrization of trapped lee waves, but the Hines run does not. This is one possible source of additional drag near the surface in the control run and perhaps it should be reintroduced in the Hines run. The fact that the orographic stress actually grows with height in the Hines run seems odd - this may well be a bug.

Further analysis of the control run (see below) shows that while the lee wave code does indeed supply some drag in the lowest two or three layers, the greatest drag in fact comes from the hydraulic jump component. The linear saturation component makes no contribution in this grid point.

XAFBG u lee acceleration XAFBG u hydraulic jump acceleration

This raises the following question: if the hydraulic jump regime is producing significant drag in the control run, why is it not doing so in the Hines run? This, together with the question of why the orographic momentum flux grows in Hines at low levels, makes me think it is worth looking at exactly what is happening in that part of the BNL code.

There seem to be two possible next steps:

  1. Put the lee wave code back into the Hines run:
  2. Look carefully at what the hydraulic jump code is doing in the Hines run. How does this differ from the control run and why?
On balance it seems unnecessary to reintroduce the lee code. The hydraulic jump code is a far more likely candidate for reducing the low level stresses.

A further point has come to light from an analysis of the code. Both the control and Hines runs call the subroutine GW_SCOR as one of the first tasks in GW_VERT. GW_SCOR calculates the vertical profile of the Scorer parameter and uses it to (a) decide whether lee waves are present, indicated by L_LEE=.TRUE. and (b) determine the value of TRANS, the "transmission coefficient" which is the scheme's method of accounting for the reflection of hydrostatic waves at the tropopause (the reflection of such waves reduces the surface stress). The value of TRANS varies between 0 and 1 where 1 indicates no reflection of waves from the tropopause. On returning from GW_SCOR the code adjusts the surface stresses as follows:

      DO I=1,POINTS
        S_X_STRESS=S_X_STRESS*TRANS(I)
        S_Y_STRESS=S_Y_STRESS*TRANS(I)
      ENDDO
where S_X_STRESS and S_Y_STRESS are the zonal and meridional components of the surface stress. The 3A code uses these as the initial values in all subsequent calculations of momentum flux in the hydraulic jump and/or linear saturation regimes. (The lee wave code determines its own additional surface stress and calculates how that would be distributed in the vertical). The Hines code also uses the TRANS adjusted S_X_STRESS and S_Y_STRESS in its hydraulic jump calculation. However, prior to calling GW_SCOR and adjusting the surface stresses the Hines code does the following:
      S_STRESS_SQ = S_X_STRESS(I)**2 + S_Y_STRESS(I)**2
      .
      .
      .
      S_STRESS=SQRT(S_STRESS_SQ)
      .
      .
      .
      S_STRESS_AMP(I) = S_STRESS
      !NB UM stress direction is opposite to that in Hines
      S_STRESS_ARG(I) = ATAN2(-S_Y_STRESS(I),-S_X_STRESS(I))
and later calls the spectral code, via subroutine HINESWRAPC, using S_STRESS_AMP and S_STRESS_ARG. Thus the Hines code is using inconsistent surface stresses for the hydraulic jump and the coupled orographic/spectral code. This is a bug. The stresses fed to the Hines code should be adjusted by TRANS in the same way as those used elsewhere.

16/03/2005

The adjustment of S_STRESS_AMP by TRANS has now been added to the Hines code:

      S_STRESS_AMP(I)=S_STRESS_AMP(I)*TRANS(I)
and TRANS_D has been reintroduced as an available diagnostic (see below). The changes are incorporated in the latest Hines run, XAFBM.
XAFBM vertical transmission coefficient

Somewhat surprisingly, there are differences in the transfer coefficient between the Hines and the control runs. I would have expected them to be identical. Understanding the reason for the differences is not a high priority currently, but is noted here because some of the values are quite large (see below).
XAFBM minus XAFBG vertical transmission coefficient

The adjustment of S_STRESS_AMP by TRANS cures the problem of the inreasing momentum flux between levels 1 and 2 in the Hines run. Previously the code was using the adjusted value of surface stress in the lowest level, but passing the unadjusted stress into the bottom of the Hines calculations, hence the increase. The corrected profile is shown below (cf the same plot from XAFBL above):
Orographic u mom flux

The question of why the hydraulic jump code seems to have no effect in the Hines run has also been answered: it does have an effect, just not in the grid point that I was studying! I was looking at a grid point with co-ordinates 292.50 degrees E, 55.00 degrees N in the Hines run. A neighbouring point (296.25 degrees E, 57.5 degrees N) does show hydraulic jump activity:
Orographic u mom flux

However, the above profile is alarming as it actually allows the momentum flux to change sign between levels 5 and 6. I suspect that level 5 is the top of the hydraulic jump and level 6 is calculated by the coupled Hines code and perhaps an inadvertent sign change is introduced here, but more investigation is needed to confirm this. I would also like to investigate the diagnosis of hydraulic jumps - if one had occurred in the grid point at 292.50 degrees E, 55.00 degrees N it would have reduced the surface stress to either one third of its initial value, or even zero if a critical level was also diagnosed. If hydraulic jumps are not being diagnosed when they should be, this could explain why the orographic accelerations are sometimes excessive at upper levels in the model.

29/03/2005
The BNL code introduces a sign change in the subroutine HNSWRAPC by multiplying the momentum flux profiles (both orographic and spectral) by minus one. This is done to maintain compatibility between the inner Hines routine and the UM while calculating the fluxes in the correct physical sense (see Bryan's notes on the UM sign convention). However, the hydraulic jump fluxes are calculated in GWVERT, outside the main Hines routine and are already in the correct sense for the UM. These fluxes are placed in the lowest few levels of the profile in grid points where a hydraulic jump has been diagnosed, then the Hines routine is called to calculate fluxes for the remaining levels. The sign inversion should only be applied to those levels in which the fluxes are calculated by the Hines subroutine, HNS_UM8C. Adding an IF statement so that only levels above any hydraulic jump have their sign changed seems to cure the problem with the flux profile:

c
c following because of warped sign convention in UM!
c
      DO K=1,LEVELS
         do j=1,rows_p
            do jj=1,row_length
               i=(j-1)*row_length+jj
               EWF(I,K)=-EWF(I,K)
               NSF(I,K)=-NSF(I,K)
               if (k .ge. qlaunch(i)) then
                 EWFQ(I,K)=-EWFQ(I,K)
                 NSFQ(I,K)=-NSFQ(I,K)
               endif
            enddo
         ENDDO
      ENDDO
The plot below (for the same grid point as the last plot) is taken from experiment XAFBN in which the correction has been applied:
Orographic u mom flux

However, when the corresponding accelerations are plotted, there is no sign of the hydraulic jump! Why? Once again we are back to asking why the jumps are not being applied in the Hines run.
Orographic u acceleration

Below is another pair of plots of momentum flux and acceleration from XAFBN, this time for grid point 270 degrees E, 50 degrees N.
Orographic u mom flux Orographic u acceleration

In this case there seems to be a critical level at about 18km above which the momentum flux is zero. A corresponding peak in the accelerations can be seen. However, there is also a phantom acceleration in the top level which does not correspond to any feature in the momentum fluxes. Where is this coming from? These two examples serve to illustrate that there are some issues surrounding the way the accelerations are calculated from the fluxes that need to be resolved.

29/07/2005
Many extra diagnostics have been added to the gravity wave scheme in both the Hines and control versions. It appears that all these diagnostics are identical except those calculated inside the Scorer routine, GW_SCOR (gwscor3a.f). The diagnostics that differ are TRANS (as noted above) and L_LEE which is a logical flag to indicate where lee waves are diagnosed in the model. This is a puzzle given that the Hines scheme mods do not alter the Scorer routine in any way. U, V, exner pressure, theta and P* are all identical on entering the Scorer routine and the code should be identical between the control and Hines runs. However, the outputs from the routine differ! How can this be? Below is the difference in the TRANS diagnostic between runs with the latest modsets for the Hines (XAFBZ) and control (XAFBY) versions. This is a new version of the diagnostic on the P-grid, i.e., the native grid of the GWD scheme.

TRANS on P-grid

The earlier plot of the differences in TRANS (see 16/3/2005 entry) used the diagnostic coded in the 3A scheme which is regridded onto the U/V grid before being output. The fact that differences remain when TRANS is plotted on the P-grid shows two things:

It turns out that, although the Hines scheme does not modify the Scorer routine, the two model versions call the subroutine differently. One of the arguments to the Scorer routine is LEVELS which, in the Hines experiments, is just the number of levels in the model, namely 64. However, in the control version of the 64 level model Rayleigh friction is applied above 20 hPa (50 km) and orographic gravity waves are turned off just below this height. This is achieved by calculating a quantity referred to as TOP_GWD_LEVEL which is the model level immediately below 20 hPa. This quantity is passed to GW_SCOR as its LEVELS argument and effectively limits the 3A orographic waves to the lowest 35 levels of the model. The LEVELS argument is used within the Scorer routine to calculate values within the 1D array K_LEE_LIM. This in turn is used to define the levels between which the average Brunt-Vaisala frequency, N, is calculated. These average values of N are used ultimately in the calculation of TRANS. Thus, simply passing a different number of levels to GW_SCOR in the Hines and control runs results in the differences noted in the TRANS diagnostic and therefore affects the surface stresses at the grid points where TRANS differs. A second version of XAFBZ (Hines experiment) was run in which the values in K_LEE_LIM were artificially set equal to those in the control run. The differences in TRANS were eliminated as shown below.

TRANS on P-grid

It is probably unphysical to limit the levels in Scorer in the Hines experiments in this way, but now at least the differences between the control and Hines runs at the beginning of the first time step are understood.

In the process of investigating the values in K_LEE_LIM it was found that this array is never initialized in the code. In the control run the array contained large random values from the 13th element onwards. In the Hines run, fortuitously, it contained zeroes. Although the uninitialized elements of the array do not appear to be used, and therefore probably do not affect the experiments, it seems advisable to correct this problem. A mod to initialize K_LEE_LIM was written and added to modset init_gwd which should now be used in both control and Hines experiments.

The current state of play regarding which modsets to select is as follows:

  1. CONTROL - without debug diagnostics
    $MODS64/gsdev20_l64_3a
    $JAP_MODS/init_gwd
  2. CONTROL - with debug diagnostics
    As item (1) plus:
    $JAP_MODS/gwdiags_3a
    $JAP_MODS/gwscor_3a_call_debug
  3. HINES - without debug diagnostics
    $BNL_MODS/gsdev20_l64_hines
    $JAP_MODS/hines-call20-diag $JAP_MODS/hines-cbnl6 $JAP_MODS/init_gwd
  4. HINES - with debug diagnostics
    As item (3) plus:
    $JAP_MODS/gwdiags_hines
N.B. $JAP_MODS=/home/tolkien/alisonp/mods/japmods; $MODS64=$UM_HOME/local/vn4.5/ mods/l64; $BNL_MODS=/home/tolkien/alisonp/mods/bnlmods.

01/08/2005
At a meeting with Bryan and Scott Osprey, at which the problems with the Hines experiments were discussed, Neal Butchart (Met Office, Stratospheric Group) suggested checking whether second order finite differencing is being used in my runs. Within the Met Office this is standard practice for stratospheric versions of the UM in order to achieve stable running (operationally, 4th order differencing is used). Second order finite differencing is turned on by setting NU_BASIC=0.000000 in window atmos_Science_Section_Advec (i.e., the Primary field advection window underneath Scientific Parameters and Sections -> Section by section choices) in the UMUI. In fact, all my experiments are already using second order differencing. However, Neal emailed a job basis library for a recent UM 4.7 run on the Met Office NECs and suggested comparing all the parameter settings with my own. 4.7 is similar to 4.5 and is currently being used in the Met Office to run stratospheric chemistry experiments. The job supplied by Neal has not been used in long runs and so is not guaranteed to be stable enough for climate experiments.

Another issue discussed with Neal was time stepping. The 4.7 runs in the Met Office use 4 dynamics sweeps per physics time step and the test for winds exceeding a threshold is switched on (Atmosphere -> Scientific Parameters and Sections -> Time-Stepping window in the UMUI). The maximum wind check further halves the dynamics time step if the specified threshold, usually 225.0 m/s, is exceeded. The result is a dynamics timestep of 3.75 minutes when winds faster than 225 m/s are encountered.

03/08/2005
The results from the latest runs which experiment with the parameter settings for time stepping, q diffusion and K_LEE_LIM in the Scorer routine are summarized in the table below.

Experiment name

XAGZB




XAGZC



XAGZD



XAGZE



Experiment description

Coupled Hines. K_LEE_LIM (11) and (12) set to 32 in Scorer (as per control). 4 dynamics sweeps per physics time step + max wind test for halving dynamics time step.

As XAGZB but K_LEE_LIM set from LEVELS argument to Scorer (i.e., elements 11 and 12 set to 61).

As XAGZC but 5 dynamics sweeps per physics time step + max wind test for halving dynamics timestep.

As XAGZC but q diffusion coefficients and RHcrit copied from Neal Butchart's experiment (xahba).

No. of time steps completed

865 (18 days 30 mins)




1279 (26 days 15 hours 30 mins)



1440 (30 days = 1 month)
RUN COMPLETED
WITHOUT CRASHING

1393 (29 days 30 mins)



From this I conclude that:

  1. Neal's q diffusion parameters improve model stability from 26.5 days to 29 days and therefore should be used in future runs;
  2. it is necessary to use 5 dynamics sweeps per physics time step to achieve a full month's run;
  3. changing elements 11 and 12 of the K_LEE_LIM array in the Scorer routine from 32 to 61 improves model stability from 18 days to 26.5 days. Clearly the Hines runs are very sensitive to the levels used in the Scorer calculations and the next step is to decide the correct values for K_LEE_LIM.

23/08/2005
A detailed analysis of the Scorer code has been carried out to determine the correct values to use in elements 11 and 12 of K_LEE_LIM. In the 64 level model these are the elements that determine the full vertical extent of the non-hydrostatic trapped lee waves and any hydrostatic waves reflected at the tropopause. The documentation of the 3A UM gravity wave parametrization is very clear that both the reflected waves and the trapped lee waves are confined to the troposphere. Therefore, an upper limit of level 61 (near the top of the stratosphere in this model version) in the Hines runs is clearly wrong and means that the Scorer routine is not working as intended.

For both lee waves and reflected hydrostatic waves the algorithm requires an atmospheric column which is divided into two parts of equal height. The mean Scorer parameter is calculated over those model levels occurring within the lower and upper sections, respectively. The mean parameters are used within the Scorer routine to calculate the transmission coefficient (TRANS) and, in the control experiment, are also passed to subroutine GW_LEE to allow the lee wave surface fluxes and vertical profiles to be calculated. GW_LEE is not called in the Hines experiment.

The code's implementation of the algorithm leaves a lot to be desired as regards splitting the tropospheric column in half with respect to height. The relevant code section is as follows (comments in italics are mine and do not appear in the UM code):

! Subroutine arguments:
      INTEGER
     * LEVELS              ! IN   NUMBER OF MODEL LEVELS
      REAL
     * AKH(LEVELS+1)       ! IN   VALUE AT LAYER BOUNDARY
     *,BKH(LEVELS+1)       ! IN   VALUE AT LAYER BOUNDARY

        .
        .
        .

! Local scalars
      REAL    PU
      LOGICAL FLAG
      INTEGER K,M

        .
        .
        .

! Local dynamic arrays
! The next three variables are all MAXIMA used to limit lee waves and/or TRANS
! in the vertical, so they don't necessarily influence every lee grid point
! directly
      INTEGER
     * K_LEE_LIM(LEVELS-3) ! TOP MODEL LEVEL HEIGHT FOR CALCULATION
     *                     ! OF L_SQ(I,2), FOR EACH POSSIBLE LEE LEVEL
     *,K_LEE_MAX           ! MAXIMUM LEVEL FOR FINDING LEE WAVE HEIGHT
                           ! I.e., max level for lee INTERFACE
     *,K_TRANS             ! LEVEL OF L_SQ SPLIT FOR CALCULATION OF
                           ! TRANSMISSION COEFFICIENT
      REAL
     * P(LEVELS-3)         ! PRESSURE AT LEVEL TOP BOUNDARIES
     *                     ! FOR STANDARD POINT
     *,P_LIM(LEVELS-3)     ! PRESSURE AT TOP BOUNDARIES FOR
     *                     ! K_LEE_LIM
        .
        .
        .

!---------------------------------------------------------------------
! Find approximate level of 0.6 tropopause height - maximum lee limit
! Find approx half tropopause height - transmittion coefficient top
!---------------------------------------------------------------------
      FLAG = .TRUE.
! N.B. LEVELS = 35 in 3A + Rayleigh, LEVELS=64 in Hines, therefore
! K_LEE_MAX = 31 or 60, respectively
      K_LEE_MAX = LEVELS-4
      K_TRANS = K_LEE_MAX
      DO K= 3,LEVELS-4
         IF (FLAG) THEN
! The use of 100000 Pascals in the calculation of PU means that a surface
! pressure of 1000 mb is assumed.  This is almost certainly a poor assumption
! near orography!
! This point is noted in the UM documentation, but the code makes the assumption
! anyway!!!
           PU=100000.*BKH(K+1) + AKH(K+1)
           IF ( PU .GT. 50000. ) THEN
             K_TRANS = K
! K_TRANS = 11 in my 64 level runs
           ELSE IF ( PU .LT. 45000. ) THEN
             K_LEE_MAX = K
! K_LEE_MAX = 12 in my 64 level runs
             FLAG = .FALSE.
           ENDIF
         ENDIF
      ENDDO
!---------------------------------------------------------------------
! Find array of top level scorer calculation limits K_LEE_LIM, for each
! possible lee level by using the same standard point as above.
!---------------------------------------------------------------------
      DO K= 2,LEVELS-3
        P(K)=100000.*BKH(K+1) + AKH(K+1)
! This next bit is horrible!!!
! Note that for P(K)=500 mb or less, P_LIM (a pressure!) <= 0
! Also, contrary to the UM documentation, the tropospheric column is being split
! by pressure and not by height!!!
! N.B. The Scorer routine was originally coded in 1994 when the climate version
! of the UM (and maybe the operational version) had 19 levels in the vertical.
! The levels were positioned differently to those eventually used when higher
! vertical resolution was introduced.  This horrible code may have worked
! (fortuitously) at 19 levels but it certainly doesn't work at 64!
!
! Initially, all elements of K_LEE_LIM = LEVELS-3 for K <= lee interface level
! (L12)
        IF ( K .LE. K_LEE_MAX ) THEN
          P_LIM(K) = 2*P(K) - 100000.
          K_LEE_LIM(K) = LEVELS-3
        ENDIF
      ENDDO

      DO M= 2,K_LEE_MAX
        DO K= M+1,LEVELS-3
! Now for the bit where it all comes unstuck!!!
! P_LIM (which is negative in places!) is used to define the maximum vertical
! extent of lee waves and TRANS by setting the elements of K_LEE_LIM.  Because
! a model level pressure < 0 is never found the loop simply continues to iterate
! to LEVELS-3 with the result that K_LEE_LIM elements (11) and (12) are set to
! 32 in the control run and 61 in the Hines version.  (The addition of Rayleigh
! friction to the control run and the accompanying lid on orographic gravity
! waves is all that stops K_LEE_LIM being set to 61 in that version also).
          IF ( P(K) .LT. P_LIM(M) .AND.
     *       K_LEE_LIM(M).EQ.LEVELS-3  ) THEN
             K_LEE_LIM(M) = K
          ENDIF
        ENDDO
      ENDDO
      K_LEE_LIM(1)=2
        .
        .
        .
06/09/2005
Experiment name

XAGZF





XAGZI



XAGZJ



XAGZK

Experiment description

q diffusion and RHcrit as Neal Butchart's experiment, 5 dynamics sweeps per physics time step + max wind test for halving dynamics time step, K_LEE_LIM(11) and(12) set to 32 in Scorer routine.

As XAGZF but K_LEE_LIM(11) and (12) set to 26 in Scorer.


As XAGZI but 4 dynamics sweeps per physics time step.


As XAGZI but 2 dynamics sweeps per physics time step.
No. of time steps completed

1153 (24 days 30 mins)





1440 (30 days = 1 month)
RUN COMPLETED
WITHOUT CRASHING

1440 (30 days = 1 month)
RUN COMPLETED
WITHOUT CRASHING

77 (1 day 14 hours 30 mins)

The above results from XAGZF show that, even using moisture diffusion and RH parameters as per Neal Butchart's experiment and 5 dynamics sweeps per time step, setting K_LEE_LIM(11) and (12) equal to 32 (as in the control run) is not sufficient to obtain stable running for a full month. The value of 32 is somewhat arbitrary as it is based on the vertical limit imposed on orographic gravity waves when the Rayleigh friction parametrization is included. There is no reason to believe that this is physically consistent with the actual position of the model levels. Examination of the global mean pressure on model levels suggests that level 24 is at around 150 hPa which might be considered a reasonable mean height for the tropopause. However, K_LEE_LIM(10) is set to level 26 by the code, so if elements (11) and (12) were artificially set to 24 this would have the effect of limiting lee and reflected waves with the highest layer interfaces to have a smaller total vertical extent than waves with a lower interface. This seems inconsistent with the way in which the Scorer routine treats all the other low level waves. Therefore, a level of 26 for K_LEE_LIM(11) and (12) was chosen. Apart from maintaining consistency with lower level waves this also limits the TRANS and lee wave calculations at around the 100 hPa level and thus should include all of the tropopause, even in the tropics. This more physically based limit was used in Hines runs with both 5 and 4 dynamics sweeps per time step (runs XAGZI and XAGZJ, respectively). As the table shows, both were stable for a month. However, an attempt to lengthen the dynamics time step further (XAGZK) was unsuccessful.

14/09/2005
Based on the latest results (see 06/09/2005) two climate runs, each on 16 processors, were submitted using Neal's moisture parameters, a level limit of 26 in the Scorer routine and 4 dynamics sweeps per physics time step. The runs were XAHEA (coupled Hines) and XAHEB (uncoupled Hines). Apart from the addition of some climate meaning diagnostics XAHEA was identical to XAGZJ. XAHEB was identical to XAHEA except that the coupling flag within the Hines code was set to .false. Neither run completed the first month. This was particularly puzzling in the case of XAHEA (which ran for 29 days) as it was scientifically identical to the 8 processor job XAGZJ which did complete a month.

XAHEA and XAHEB were resubmitted on 16 processors, this time using 5 dynamics sweeps per time step. They ran until crashing with QTPOS and filtering error messages as follows:

Examination of the job log files showed that they were teetering on the edge of breaking the CFL criterion almost from the outset with filtering being intermittently switched off in both jobs in almost every month of integration. The frequency and duration of periods without filtering seemed to increase with time in XAHEB before the job finally crashed, but in XAHEA the crash came almost without warning.

Comparison of the log files of XAHEA and XAGZI, which are scientifically identical except for the addition of diagnostics and which both used 5 dynamics sweeps per time step, revealed that filtering messages were appearing in different time steps in the two runs. This suggests a difference in the model's behaviour on different numbers of processors (XAHEA: 16, XAGZI: 8).

10/10/2005
The next priority is to examine the apparent differences in model behaviour between different numbers of processors as this is very worrying. The first plot below shows the difference in zonal mean monthly mean U wind for the first month of integration of XAHEB (uncoupled Hines) and XAHEA (coupled Hines). The second plot shows the same thing but for the difference between XAHEA (coupled Hines, 16 processors) and XAGZI (coupled Hines, 8 processors).
Zonal Mean Monthly Mean U
Zonal Mean Monthly Mean U

The main point to note in both plots is the dipole at upper levels in the polar jet, indicating a shift in position. Uncoupling the Hines scheme moves the jet polewards in the mesosphere and equatorwards in the stratosphere. However, changing from 8 to 16 processors produces an effect of roughly equal magnitude in the opposite direction, particularly in the mesosphere. Similar differences can be seen between the various runs in plots of instantaneous U at the end of the first day of integration. The fact that changing the number of processors has an effect equal in magnitude to turning off the Hines coupling indicates that the results of these runs are unreliable and that the cause of the differences between the 8 and 16 processor runs needs further investigation.

12/12/2005
Zonal Mean U DJF 1978-79

30/01/2006
Zonal mean U ensemble mean

Back to top


Email: J.A.Pamment@rl.ac.uk
Tel: +44 1235 778065
Fax: +44 1235 445848

Last modified 10/10/2005