Description

A slight 'decentering' in the calculation of the divergence damping coefficient was found which caused a model crash

over steep orography in a 7km simulation over the Irakian mountains.

A more symmetric calculation cures this problem (and even enables one to reduce the coefficient over steep terrain).

This action has a certain, however small, influence to the results.

(this work has finished in 2015 and came into COSMO version 5.1)

NEW:

Recently (Dec. 2015), a real bug in the subroutine has been identified (missing metric coefficient)

which should be fixed after some testing.

In accordance with this bug fix, there is a certain hope that the namelist-parameter divdamp_slope is not longer necessary.

Documentation

This is a very technical change in the calculation of the divergence damping coefficient in tilted terrain.

The description given in Baldauf (2013) COSMO Techn. Report is still correct (i.e. no additional documentation seems necessary.)

User's guide: no adaptation necessary.

Process documentation (history/changes log file): is updated

NEW:

no change in documentation necessary

At DWD some unrealistic spikes have been observed lately in the meteograph output of the latent heat flux at grid point Frankfurt. An investigation showed, that this was due to a dry soil moisture layer 5, where soil moisture nearly was at the plant wilting point (PWP). The plants are fully evaporating, until no soil moisture is available any more, then they stop evaporating. And if some moisture is added, the evaporation starts again fully.

A more detailed description of the fix and the visualization of the latent heat fluxes can be found here.

grib2 input leads to wrong internal level index "kflat" for those PE domains in a parallel run which are really flat (i.e., PE domains over sea).

kflat is the uppermost level where hhl has horizontal variations, right below height vcflat.

Because in grib2, kflat is no longer stored in the grib2 grid meta data, it has to be recomputed during model initialization.

Here, the bug is a missing parallel communication (global maximum) of kfalt, after it has been computed from the local criterion

"max(hhl(:,:,k)) /= min(hhl(:,:,k))" in a vertical loop on each single PE. Now, for all "flat" PE domains, this criterion is never fulfilled and kflat remains at its initial value 1.

Then, the initialization loop for the upper Rayleigh damping coeficients, which runs from 1 to kflat, does not compute the damping coefficients,

which means they are 0.0 and the Rayleigh damping is deactivated.

Description:

In the case, that density is transported additionally in the finite volume advection schemes (as in the Bott-scheme), it can occur, that a grid box is completely emptied due to the direction splitting. This would result in an explosion of the specific mass of cloud variables. Therefore a limiter is introduced, which sets the density to at least 10% of the density of the reference state.

This changes the results.

Description:

When fixing the bug with using the variable tgrlat, which was calculated on u-grid points, also on v-grid points, another problem was introduced:  Now the values are calculated on the v-grid points, but at one point it was forgotten to remove the averaging factor 0.5. This has been eliminated now.

The bug fix has been tested in DWD's parallel suite with neutral verification results.

Description:

Some time ago, the initialization of several coefficients in the microphysics has been changed. But there are several other constants, derived from these coefficients, which are not computed in the code, but also initialized with a certain value. And these values have not been updated accordingly.

This has been corrected now in the cloude ice scheme hydci_pp and the graupel scheme hydci_pp_gr.

The correction changes the results. The effect has been tested at DWD for several weeks with neutral impact on verification.

In spring 2012, occasional crashes of the COSMO model occurred at DWD. The errors were traced back to the routines of the lake parameterization scheme FLake. To clear the trouble, changes have been made in the COSMO-model files data_flake.f90 and src_flake.f90.

In SUBROUTINE flake_radflux (file src_flake.f90), the solar radiation flux at various levels (at the ice-water or air-water, mixed layer-thermocline, and water column-bottom sediment interfaces) and the mean fluxes over the mixed layer and over the thermocline are computed. An exponential decay law for the solar radiation flux is used. Provision for a two-band approximation of the decay law is made. The simplest one-band approximation is used at the moment, where the fraction of the total solar radiation flux for the “unused band†is set to zero and the extinction coefficient is set to a very large value. A very large extinction coefficient is also used for the ice that effectively prevents the penetration of solar radiation into the ice interior. The use of a large extinction coefficient results in a large argument of the EXP function (large magnitude as the argument is negative). This usually causes no problem. However, some compilers (with certain combinations of the compiler options) report an “invalid operation†when the argument of the EXP function is large. To eliminate the problem, the magnitude of the EXP function argument is now limited by 1.E+02. A new constant (PARAMETER) c_maxearg_flk is defined in MODULE data_flake (file data_flake.f90) that is used in SUBROUTINE flake_radflux (file src_flake.f90).

Since the lake parameterization scheme FLake is based on the self-similar representation (assumed shape) of the evolving temperature profile, certain relations between FLake prognostic variables (temperature and depth of the mixed layer, bottom temperature, mean temperature of the water column, and shape factor for the temperature profile in the thermocline) must hold at any time. This is ensured during the COSMO model runs. However, a loss of accuracy due to GRIB encoding and decoding may result in a violation of the above relations. Consider the following example. Suppose the mixed-layer temperature, the bottom temperature, and the mean temperature of the water column proved to be equal to each other at the end of a COSMO run, i.e. the water column is fully mixed and neutrally stratified. If these temperatures are GRIB encoded and then decoded to be used as initial conditions for the next COSMO run, a loss of accuracy during GRIB IO may result in slightly unstable stratification (e.g. bottom temperature exceeds mixed-layer temperature where both temperature are larger than the fresh-water temperature of maximum density) which is forbidden within FLake. In the overwhelming majority of situations, the violations of the FLake rules are harmless due to security constraints in the FLake code. However, in some very rare “exotic†cases fatal errors may occur, e.g. an attempt to take a square root of a negative number, leading to the COSMO model crash during the first time step (just after the model initialization). To remedy the trouble, changes are made in SUBROUTINE flake_init (file src_flake.f90). Consistency between the FLake prognostic variables is now enforced during the initialization, i.e. prior to the first COSMO-model time step.

Notes

There are some fields which are used only for special settings of Namelist switches (for example depth_lk is only allocated if llake=.TRUE.) Some of these fields are used in interfaces and have to be allocated anyhow, otherwise there could be problems in the corresponding subroutines.

Affected fields are: depth_lk, tket_conv, tket_hshr, tket_sso, ut_sso, vt_sso

Notes

If no microphysics are computed (lphys=.FALSE. or lgsp=.FALSE.), the program will crash, if DBZ should be written to disk.

The variable DBZ has to be eliminated from the output list, if lphys=.FALSE or lgsp=.FALSE.

Notes

A verification of cloud cover with the Lindenberg cloud radar showed, that the PDF of high ice clouds is too low in the model and should be changed.

Notes

The mpe_io-module (as is) crashes on the new Cray XE6. There is a problem in the initialization phase.

The problem was a MPI_Send from one processor to itself. The MPI Library on the Cray seems not to like that. A workaround at MeteoSwiss was to use Buffered Send, with which this problem did not occur.

The solution now is not to send data within a processor, but just to copy it.

Notes

A strange behaviour of the Lahey Fujitsu Compiler has been reported:
For the routines advection_ef_xyz or advection_ef_zyx, the field for the tracer that should be advected, is given twice to the subroutine: Once with INTENT(IN) and once with INTENT(OUT). Within the subroutine, these are two different fields then, called "s" and "su".

The Lahey compiler did not modify the field, most probably it only accepted the INTENT(IN) field.

Is this a compiler-bug or a program-bug?

According to the Fortran Standard, this kind of argument aliasing is explicitely forebidden. So it really is a program bug.

Notes

Some more diagnostic debug output in double precision is desired
Mail by Xavier Lapillonne from 03. November 2010 to Uli Schaettler.

Notes

Passing non-allocated arrays in parameter lists is not accepted by all compilers and can lead to model crashes.

Notes

Just for completeness. These fluxes can only be written as averaged values up to now.

Notes

The extrapolation of pressure to the uppermost half level is wrong in src_gridpoints (only for diagnostic print outs) and for interpolation to pressure level (subroutine p_int in src_output.f90). This could even lead to a model crash

Notes

The variable nsat_steps is only allocated, if output of synthetic satellite images is switched on. But in the init- and compute-phase of organize_satellites it is accessed in all cases. Therefore allocate this variable also in all cases.

Notes

Correction of size distribution for rain droplets (in hydci_pp, hydci_pp_gr). Generalization of implementation for arbitrary mu_rain.

Bugfix in density correction of fall speeds (by Thorsten Reinhardt): Added missing factor zrho1o2 for consideration of density in the sedimentation of rain and snow (in hydci_pp)

Notes

In the computation of the fall speeds a factor is missing for consideration of density in the sedimentation of rain and snow (in hydci_pp)

Notes

The variable govcp is initialized only for itheta_adv>1, but in case of irefatm=3, it is also used in other cases.

Notes

tgrlat (tangens of transformed latitude) is computed only on u-latitudes (same latitudes as mass grid points), but has also to be considered on v-latitudes, which are different from the u-latitudes because of the staggering.

Notes

Mail by Schneider (Uni Bonn) to Michael Baldauf.

Notes

There are stripes in the output, when running CCLM on COSMO-DE resolution.

Mail by Erwan Brisson to Uli Schaettler.

Notes

The 2m temperature T_2M is used in the soil model and therefore influences the prognostic results. This is a problem, because there are different schemes to compute the 2m temperature. The whole COSMO results therefore dependent on the choice of the special scheme.

Notes

Problem with advection for confluent horizontal flows

Adapted horizontal advection of T and pp to handle flow situations with confluent horizontal flows in neighbouring grid boxes (opposite signs of U and/or V) better. This avoids a spurious and potentially detrimental heat source in such grid boxes, which lead to extreme Theta-values, to "grid point storms" and other unwanted effects.

New subroutines "horiz_adv_driver_stagmix", "adv_upwind?_???_stagmix", new namelist switch "itype_tppadv";

Notes

The SR potential_vorticity_rho is called with variable zvarlev as argument, but with wrong limits. The vertical dimension must be between (1:ke)

Notes

The interface of SR calc_Theta_Tppp expects reference temperature and temperature deviation from the calling routine. But in the output, only full temperature is available.

In src_advection_rk.f90 there is a out-of-bounds array access in case ldyn_bbc=.true. and nboundlines<3 for the computation of zw_ke(i,j,1) where the following is done...

    DO   j = jstart-3, jend+3
      DO i = istart-3, iend+3
        zw_ke(i,j,1) = 0.5_wp * ( w (i,j,ke,nadv) + w (i,j,ke+1,nadv) )
      ENDDO
    ENDDO

If really required, this should probably be checked at model setup time and the model should abort in case of these settings.

Notes

The results of the COSMO simulations are not reproducible on MeteoSwiss's Cray, if luseobs=.TRUE., i.e. the results differ, if a different processor configuration is used.

Notes

Depending on tuning constants q_crit and clc_diag, negative values of ql could occur in SR cloud_diag in meteo_utilities.f90 and have to be avoided.

Notes

There are potential inconsistencies in computing the density RHO in src_advection_rk.f90, regarding the timelevel.

Notes

The subroutines of module pp_utilities are written for the full threedimensional arrays. For the parallelized cosmo-model this is no problem. But for sequential fieldextra, this is a memory problem.

The routines should be re-written for two-dimensional slices, as is done with the physical routines for COSMO-ICON Physics.

Notes

In connection with the revision of the Runge-Kutta core, NEC proposed some more optimizations for the SX-9.

Notes

Mail by Andreas Will from 30. November 2010 to Michael Baldauf, Uli Schaettler

Notes

When using nboundlines=4, the Latent Heat Nudging crashes

Mail by Xavier Lapillonne from 12. November 2010 to Uli Schaettler et al.

Notes

The computation of radiation on a coarser grid is not working with nboundlines>3

Mail by Beate Saendig from 11. August 2010 to Uli Schaettler.

Notes

There are differences between grib- and meteograph-output.

Mail by Gerd Vogel to Uli Schaettler.

Corrected wrong vertical index range in final computation of pp_lm for non-hydrostatic input models (COSMO, ICON, UM), possible segfault

Corrected wrong index bound when computing domain decomposition, possible crash

Corrected problems with T_SO at coastal points for very coarse input model. A previous fix for this problem had erroneously lead to T_SO = 0 K in areas with 0 soil moisture (deserts)

no change of results

Previously, the PBL height was determined by the minimum layer index k over all PE domains for which the PE domain max P is just above 850 hPa. This is problematic for elevated domains without sea points and surely fails if the entire domain is a high plain. Now this is has been changed: Now we search for the highest level below the pressure ps-150 hPa*(1-min(hsurf)/vcflat_icon) at the point(s) with the domain minimum orography. If more than one such points exists, the one(s) with the maximum surface pressure over these points are chosen. If again more than one such points exists, the point with the lowest level index wins.

For this: new namelist parameter vcflat_icon (default 16000 m). This changes the results because of the max(ps)-150 hPa instead of fixed 850 hPa value.

Fixed subroutine make_fn(), so that now also ICON analyses can be interpolated and used as initial and boundary conditions for COSMO.

Note: Currently does not work for ICON analyses from DWD, because some fields are missing there (W, QV_S, T_G).
Here, the forecast fields from vv=0 should be used instead: rename manually to igafYYYYMMDDhh and set yinput_type=‘analysis’ and
yinput_prefix='igaf'.

Does not change results!

For nonhydrostatic input data (ICON and COSMO): updated vertical hydrostatic pressure computation (lbalance_pp=.true.) to be consistent with the RK-dynamics with its two different fast-waves solvers. Have overtaken routine calc_p_hydrostat_psts() from COSMO 5.3 src_artifdata.f90.

New namlist parameters itype_balance_pp (1=old, 2=new) and itype_fast_waves_lm (same as itype_fastwaves in COSMO).

Changes the results slightly for p respectively pp, qv and qc.

Fixed copying of yvarbd_d(:) to youtlist_bd(:), where an erroneous nzmxini_d instead of nzmxbd_d was used to specify the copy-range.

Does not change results!

Fixed problems with nlevskip if some of the skipped levels do appear in the coarse input files (added some „return“ in src_read_coarse_grid.f90).

Does not change results!

Introduced some tolerance in refatm check for COSMO2COSMO, because there could be a loss of accuracy in grib coding.

Does not change results!

For ICON2COSMO, ncdf dimension name of „cell“ in grid file was „ncells“ for older grid files. Now both names are possible.

Does not change results!

Surface temperature T_SO (initial data) should be adapted to the new orography by way of conservation of vertical differences to near-surface atmospheric temperature T_ke.

But for ICON input, this was done wrong so that effectively T_SO was not adapted but only horizontally interpolated.

Serious, if COSMO model is started from interpolated analysis.

See: Slides by U. Blahak (starting page 13 ff)

src_coarse_interpol.f90: rh_s_gl, the rel. moisture at the surface, is derived from the input field qv_s_in by using t_so_in(0) instead of t_g_in. This is wrong over land in case of snow cover, where

t_so_in(0) instead of t_g_in are different. t_g_in itself would have to be computed by using the subroutine tgcom(). What about lake- or sea ice points?

src_2d_fields.f90: Here, qv_s_lm is computed from rh_s_gl, and here the correct zt_g_lm this used.

Notes

The variable czmls_in, which holds the depths of the soil layers in cm, i only known in processor 0, not in the other ones, because in src_namelists.f90, erroneously the variable czml_soil_in is distributed.

Notes

Hi Uli,
Andrea is trying to make some COSMO runs on the Sochi area and he ran into a bug with reading grib2 data: startlon does not match between grib and namelist.
I have looked at io_utilities.f90, without doing any test, but I think the problem is very simple:

      IF (igrbed == 1) THEN
        zstartlon_g = REAL(rlo1, ireals)
        zendlon_g   = REAL(rlo2, ireals)
      ELSEIF (igrbed == 2) THEN
        IF (rlo1 >= 180_ireals) THEN
          zstartlon_g = REAL(rlo1, ireals) - 360_ireal
        ENDIF
        IF (rlo2 >  180_ireals) THEN
          zendlon_g   = REAL(rlo2, ireals) - 360_ireal
        ENDIF
      ENDIF
      

i.e. these 2 lines:

       zstartlon_g = REAL(rlo1, ireals)
       zendlon_g   = REAL(rlo2, ireals)
      

should be repeated also for grib2 otherwise zstartlon_g and zendlon_g are undefined if they start east of Greenwich.
best regards, Davide

Notes

Hi Uli, I tested int2lm 1.18 with IFS data on frames and there are some problems since that situation is apparently not handled by the new "isolated points" procedure, so wrong data appear, e.g. in T interpolated at the lowest hybrid level.

Up to now I applied some quick and dirty survival patches, reverting to the previous algorithm in case of frames, but I will try to implement something more clean in the next weeks, I think that, when running with frames, there will have to be 2 different isolated point structures, one for full domain (interpolation of analysis) and one for framed data (interpolation of forecasts), or maybe we should just disable Match interpolation when using frames? What do you think about it?
bye, Davide

Notes

It is not checked (not in NetCDF and not in Grib) that the grid points of the COSMO-Model and the grid points of an external parameter data set really do match. So it could be possible that the User specifies a grid, where the external parameters can be shifted up to half of the mesh width.

Notes

The following exception is not handled in the INT2LM:

If there is an isolated point in the grid of the COSMO-Model but no matching grid point (with the same land-sea-mask) in the driving model, the program stops. This could really happen e.g. if the COSMO-grid is only over land, and lake points of that grid do not appear on the coarser grid.

An exception should be implemented, that initializes the variables for the isolated points also in that case.

Notes

Values for isolated points have been initialized by the nearest available point with the same land-sea-mask. Depending on the domain decomposition, the nearest point could be different.

The new method determines the nearest available point with the same land-sea-mask in the total domain, regardless of the decomposition. Additional communication is necessary then to distribute the information to the corresponding processors.

Notes

There were problems with the calculation of the interpolation weights in module src_grids.f90. If the coarse grid point and the fine grid point are identical but the precision of calculation gives a minimal difference, huge interpolation weights occured, that produced wrong results. Now "identical" is checked against an epsilon.

Notes

The above variables have been allocated with a wrong number of vertical levels. Now they are allocated with ke_in and ke1in.

Notes

The MARS archive of ECMWF has an option to retrieve IFS data on a rotated lat-lon grid. But the INT2LM cannot handle such data. From the IFS it expects only a regular (not rotated) lat-lon grid. When reading a rotated grid, INT2LM crashes.

Short term action: Check, whether the IFS data are rotated, then do a proper model abort of INT2LM with an appropriate error message

Long term action: Evaluate, whether reading of rotated IFS data would be beneficial