Documentation of the Changes in the COSMO-Model
Version 5.3

Version 5.3 is the next version, which should become a COSMO release.

Contents:

1. Changed Calling Sequence of Assimilation and Relaxation
2. Technical Changes in the Assimilation
3. Redesign of the 3D diffusion to improve stability
4. Implemented interface to C++ dynamical core and serialization
5. Further Changes in the Dynamics
6. COSMO-ICON Physics
7. Technical Changes
8. Further Bug Fixes
9. Changes to the Namelists
10. Changes of Results

1. Changed Calling Sequence of Assimilation and Relaxation

(by O.Fuhrer, POMPA)

The new C++ dynamical core treats the dynamics and the relaxation. In the Fortran code the assimilation is done in between the dynamics and the relaxation. The reason for this is historical, since the Asselin time filtering - which is only used for the Leapfrog dynamical core - is done in the same subroutine sardass as the relaxation.

To call the relaxation before the assimilation, the two tasks of subroutine sardass in module src_relaxation.f90 have been splitted to the new subroutine relaxation (only doing the relaxation) and the new subroutine timefilter (doing the Asselin time filter), which has to be called in any case after the assimilation.

It is not expected that switching the assimilation and relaxation leads to huge differences in the results, since both are simply relaxing the fields against a reference on each timestep and switching the order of operations is only a second order effect. The version with switched assimilation and relaxation is being used for the regular COSMO-1 runs at MeteoSwiss.

There is an important change regarding the latent heat nudging: the resetting of tt_lheat to 0.0 must not be done in the relaxation any more, but has been moved to the end of the latent heat nudging module src_lheat_nudge.

Back to Contents

2. Technical Changes in the Assimilation

(by C. Schraff)

The following technical changes have been implemented:

• Adaptations to read and process AMDAR data according to WIGOS template and scatterometer (METOP) data according to Eumetsat template.
• Additional observed quantities written to the feedback file for verification (model equivalent calculator: MEC) purposes: wind speed + direction, dewpoint temperature, radiation, 3-hourly precip and wind gusts, max. and min. T-2m over arbitrary period instead of fixed 12 hours. These quantities are either read and processed or computed from other observed variables.
• minor changes related to feedback files:
  • Observations for which the obs operator is not implemented in COSMO but exists in MEC (i.e. temporally non-local obs), are set to st_passive (instead of st_obs_only, even though they have missing value in veri-data!) and flagged as fl_obstype.
  • Feedback file body entries dlat, dlon added (according to 3DVAR V1_42).
  • Entries introduced for GPS Slant Path Delay.
  • Operator flag added for verification runs.
  • index_x, z_station declared as integer instead of short.
• Enlarged horizontal, vertical and temporal check limits in redundancy check of SYNOP to account for differences in BUFR and TAC (alphanumeric) reports. (this can change the results slightly).
• Use of obs_fof_read conditioned on ifdef __COSMO__ to make it compatible with LETKF / MEC code.
• For DWD: prefer SYNOP with small KZ (TAC) over SYNOP with large KZ (BUFR) in redundancy check.

Back to Contents

3. Redesign of the 3D diffusion to improve stability

(by M. Baldauf)

Now all the terms of the 3D-diffusion, which can be treated with a tridiagonal solver, are indead handled so. This strongly improves stability in steep terrain. In fact, the 3D diffusion alone now shouldn't become unstable at all.

To this purpose, the subroutine explicit_horizontal_diffusion has been completely rewritten and a suitable subroutine complete_tendencies_uvwt_3D has been added (complete_tendencies_trcr and complete_tendencies_tke have been adapted, too).

Note:

• The changes are only active, if l3dturb=.TRUE. or/and l3dturb_metr=.TRUE. (in this case, results are changed).
• If l3dturb=.FALSE., also l3dturb_metr is set to .FALSE.
• The default value for l3dturb_metr is now .FALSE.

Back to Contents

4. Implemented interface to C++ dynamical core and serialization

(by P. Spoerri, MCH colleagues)

The calls to the C++ dynamical core have been implemented in this version of the COSMO-Model using the pragma CPP_DYCORE. In the main program, the initialization of this dynamical core is called, while the core itself is called from organize_dynamics. A new module src_cpp_dycore.f90 has been implemented, which is the interface between the Fortran code and the C++ dynamical core. Also, the module time_utilities.f90 has been adapted, to measure the wrapper routines for the C++ core.

The C++ dynamical core can be activated with the new namelist switch lcpp_dycore (default: .FALSE.)

Together with the C++ dynamical core, also the serialization has been implemented. Serialization is a method to dump certain fields from the Fortran dycore to disk, to compare the results of the C++ core with these fields for unit testing. It is implemented using the construct "!$ser", which appears in the Fortran code as a comment.

Back to Contents

5. Further Changes in the Dynamics

(by M. Baldauf)

Implemented possibility to switch on/off the Euler dynamics

A new Namelist-Parameter l_euler_dynamics allows to switch off/on the Euler dynamical core, i.e. advection, pressure gradient, divergence and buoyancy-terms for the u, v, w, T and p-equation.

This switch does not affect tracers or Coriolis force, i.e. even if l_euler_dynamics=.FALSE., tracer advection and diffusion can take place and (if lcori=.TRUE.) Coriolis force can act. (However, there is an exception for the leapfrog or semi-implicit dynamical cores: if l_euler_dynamics=.FALSE. then vertical advection of the dynamical variables still is active).

Implemented possibility to switch off the tracer advection

The namelist variable y_scalar_advect now accepts the value NONE, to switch off all tracer advection.

Bug Fixes in the Dynamics

• src_runge_kutta.f90: (by A. Arteaga)
  Corrections in the computation of q_cond.
• hori_diffusion.f90: (by A. Arteaga, O. Fuhrer)
  Bug fix in the calculation of ztha (for targeted diffusion): due to the bug, the targeted diffusion approximately acts twice as strong as intended. This has been corrected now. (This changes the results)
  Additionally, not only excessive cold pools but also excessive warm pools are removed by the targeted diffusion for cold pools (this potentially changes the results slightly)
• fast_waves_sc.f90: (by M. Baldauf, G. deMorsier)
  In the fast waves solver, again a modification in the calculation of the slope-dependent divergence damping coefficient took place to enhance stability in steep terrain.

Back to Contents

6. COSMO-ICON Physics

Microphysics

The possibility to run the microphysics before all other parameterizations has been introduced. A new switch lgsp_first has been introduced for that. The default setting is .FALSE., and the microphysics are still computed after the dynamical time stepping.

Radiation

A blocked version of the Ritter-Geleyn (RG) radiation (in module radiation_rg.f90) and the corresponding interface (radiation_interface.f90) has been implemented. The RG radiation scheme needs special input variables (similar to other radiation schemes as RRTM), which are computed in the interface. The different computations are grouped in subroutines, which have been put to a new module (radiation_utilities.f90).

This version also supports the possibility to work on a coarser grid, and because of that we need a special data structure for the variables in the interface, which is for the example of a variable xvar:

The interface module computes all input variables for the radiation scheme on this new data structure. If a coarse grid is used, an additional subroutine (average_to_coarse_grid) is called, which computes the averaging of nradcoarse*nradcoarse to one coarse grid point. The results of these averaging are given to the driving routine fesft of the RG radiation.

The decision that we want to support the coarse grid for the radiation and also the fact that the input variables for the radiation scheme are different ones than for all other parameterizations lead to the decision, that we do not use the "copy-in/copy-out" that we use for the other blocked versions of the physics, but the method described above for the computations.

A more detailed documentation of the "Radiation in blocked data structure" is in preparation.

The code has been modified in a way that the radiation can now run on GPUs, if it is compiled with -D_OPENACC. The changes are:

• radiation_interface.f90 and radiation_utilities.f90:
  inserted necessary !$acc directives which were missing up to now.
• src_block_fields_org.f90:
  Only the fields that are used for the radiation are defined and updated on host / device. This list has to be updated when more parameterizations of the COSMO-ICON physics are added.
• src_tracer.f90:
  Included new routines trcr_acc_update_device, trcr_acc_update_host for GPUs.
• organize_physics.f90:
  Implemented updates for host and device, resp., when using OpenACC Implementing possibility to run all parameterizations in one block by using -D ONE_BLOCK_PHY.
• lmorg.f90:
  Some of the already present acc-statements and calls to GPU-routines are deactivated (m_device_management, acc_utilities), because only the radiation can now run on GPUs.
  These statements have to be activated again in a later stage.
• acc_global_data.f90 and acc_utilities.f90:
  These modules have been introduced.

Turbulence

A first version of the COSMO-ICON physics version of the turbulence code has been implemented, but is not activated yet. There are 2 new modules:

• turbulence_data.f90: defines data necessary for the turbulence scheme
• turbulence_turbdiff.f90: contains all necessary subroutines (among them: turbtran and turbdiff)

All routines are using already the data structure (nproma,ke), but the input data is just taken out of the COSMO ijk-structure and not yet in blocked format. Also, there still is a consolidation process between this current COSMO version and the ICON version. For this consolidation, much more tests have to be performed with the COSMO-Model, before the version can be taken as official version.

The interface routine (turbulence_organize) for the COSMO-ICON physics version has been implemented in the existing module turbulence_interface.f90. The calls to this routine (and some initialization routines) are already implemented in organize_physics, but commented. To activate them, the comments have to be deleted, and the calls to the turbulence scheme in ijk-data format have to be commented out.

2-Moment Scheme

3D Mixing has been switched off for NCCLOUD and NCICE, which caused a crash when testing with the new 3D turbulence (l3dturb=.TRUE.)

Back to Contents

7. Technical Changes

Modification of grib_api implementation for centers /= DWD

With the earlier implementation of grib_api (use of local shortname concept, local use sections) it was not possible for other centers than DWD (WMO Code 78) to use GRIB2. Modifications have been implemented, which now allow all centers to run with grib_api (in io_metadata.f90):

• remove local section before changing the center and set it again afterwards
• set localInformationNumber in local section also for other centers
• use key backgroundProcess instead of backgroundGeneratingProcessIdentifier

Also there have been changes to the grib2 definition files that have to be used together with grib_api. These changes are first implemented for the DWD definition files (definitions.edzw) of grib_api-1.13.1.

Some of these definition files are especially written for center=78 (DWD), e.g. local.78.def or grib2LocalSectionNumber.78.table. For other centers we provide a script, which creates links with the centers WMO number to the DWD files, e.g. for the center COSMO=250:

In this way all centers can use INT2LM (and COSMO-Model) together with grib_api and the DWD definition files.

Computation of pure diabatic temperature tendencies

A new output variable TTENS_DIAB has been added.

The variable tinc_lh does not contain the correct temperature increment due to latent heating. First, it misses the main contribution due to the saturation adjustment after the dynamics (RK and leapfrog), and second, for the leapfrog dynamics it is not correctly reset and initialized in each timestep.

Therefore, a new output variable TTENS_DIAB has been introduced, which is coupled to the new internal model field ttens_diab and contains all pure diabatic tendencies of the model:

• heating/cooling by radiative flux divergences
• heating/cooling by latent heat release due to phase changes in:
  • microphycis parameterization and satad after dynamics
  • parameterized convection (deep and shallow)
    NOTE:
    the adding of another argument in Subroutines cu_tied / cu_shallow changes the results using the Cray Compiler. This is not the case with gfortran. It must be an optimization issue
• if ltmpcor = .TRUE.: associated heating/cooling by TKE sources in the turbulence (itype_turb=3 only). Not sure however, if this exclusively contains phase changes or also some SGS adiabatic mixing, which is strictly speaking not a diabatic process!

Computation of Lightning Potential Index (LPI) after Lynn et al. (2010)

The LPI is based on the coexistence of updrafts, supercooled water, graupel/hail and other ice species (cloud ice, snow) in the same grid box in a temperature range of 0 to -20 deg C. It is basically an integral over the square of the grid scale vertical speed over the height interval between 0 to -20 deg C, weighted by a function depending on the coexistence of graupel/hail with other ice species.

We have added two additional filter steps to

1. eliminate spurious signals in the vicinity of updraft cores (inactive regions of convective clouds) and
2. in graupel regions of intense orographic wave related clouds:
   • Majority of grid columns in a neighbourhood (± 10 × 10 km²) must have wmax ≥ 1.1 m/s.
   • The total buoyancy (between ps-50hPa and ps-550 hPa) must be ≥ -1500 J/Kg. Total buoyancy is defined as buo = r_d ⋅ ∫_{{ps-50 hPa}}^{{ps-550 hPa}} (Tv_parcel - Tv_surroundings) d ln(p) The above threshold has been determined by U. Blahak based on summer season data from 2014 of COSMO-DE and lightning observations. The input buoyancy field is smoothed over a certain neighbourhood region (± 20 x 20 km²) to obtain spatially more stable results.

NOTE:

• Spatial filtering does not work correctly in case of periodic BCs!
• If the interior domain is smaller than the filtering/smoothing neighbourhood(s), the neighbourhoods are reduced accordingly (could be the case in small-domain idealized cases).

There are new subroutines in pp_utilities.f90 that are called from src_output.f90, if LPI is added to yvarml:

• SUBROUTINE lightning_potential_index(...)
  For computing the LPI after Lynn et al. (2010) with the possibility of an additional filtering of stable conditions to inhibit spurious flash signals from deep orographic wave clouds. The LPI is based on the coexistence of updrafts, supercooled water, graupel/hail and other ice species (cloud ice, snow) in the same grid box in a temperature range of 0 to -20 deg C.
  This subroutine also computes the unsmoothed total buoyancy as a 2D field. But the actual buoyancy filtering is not done in this routine and the internal switch elim_stable_conditions is set to .FALSE. It is done later by a call to lpi_spatial_filter(...) after a spatial smoothing of buo.
• SUBROUTINE lpi_spatial_filter(...)
  Spatial filtering of the LPI, where the filtering works on total domain fields, not on processor subdomains:
  • set LPI to 0 if the majority of vertical columns in a 10x10 km² neighbourhood have MAX(W) < 1.1 m/s.
  • set LPI to 0 if the smoothed (20x20 km²) (or optionally unsmoothed) total buoyancy between -50 and -550 hPa above the surface is < -1500 J/kg.
    For this, we hard-code the internal switches elim_stable_conditions=.TRUE. and filter_buoyancy=.TRUE.

Back to Contents

8. Further Bug Fixes

• data_parameters.f90:
  Editorial changes
• .f90:
  Corrected settings of kzdims (for a boundary exchange in case of lperi) Allocate array wrdfac only, if it is not allocated (can be the case when running DFI) (LT)
• gscp_interface.f90:
  Also put tt_lheat_new_b and qrsflux_b to ifdef NUDGING
• mo_fdbk.f90:
  Adapted call to add_history in SR create_fdbk
• mo_fdbk.f90, src_obs_rad.f90, src_sat_rttov.f90:
  Corrected interpolation of input for RTTOV (linterp=.TRUE.). Changes results for SYNSAT fields
• numeric_utilities_rk.f90, parallel_utilities.f90:
  Replace KIND=8 for INTEGER variables by generic definition i8 from kind_parameters
• organize_assimilation.f90, lmorg.f90, src_allocation.f90:
  
  • Implemented allocation of extra fields in init-phase of organize_assimilation (from src_allocation)
  • Implemented an extra action for cleanup in organize_assimilation: deallocation of memory at the end of the program Call that action during cleanup of lmorg
• organize_data.f90:
  
  • Write HHL to restart files: this is necessary when starting from GRIB2 files and does not hurt otherwise
  • Moved computation of ndiff_bd_ini out of src_input.f90, so it can be computed once here for all cases (including restarts).
  • Bug fix for setting next output step in restart runs: If an output group finishes output before a restart begins, but constant data should be written, there was a wrong memory access in src_output.f90, because the last output step was set to nstop-1. But the corresponding array outblock%ngrib is only allocated up to outblock%outsteps < nstop-1. The last output step is now set to outblock%outsteps.
• organize_data.f90 and organize_dynamics.f90:
  Corrections for ASCII output in YUSPECIF.
• src_artifdata.f90:
  Bugfix calc_p_hydrostat_xxx: corrected initialization of qv and qc for RH>1.0.
• src_input.f90:
  
  • From now on, HHL is written to restart files and taken from there, if available (is necessary for GRIB2 based simulations, because HHL cannot be reconstructed then)
  • If HHL is not available, it is re-computed
  • Moved computation of ndiff_ini_bd to organize_data, to do it for all cases once (including restarts)
  • Fix computation of vcflat when reading GRIB2 data (had to initialize izexvc)
  • Fix in SR create_file_name in case of lbdana=.TRUE. and a dt, which does not always fit to a full hour
• src_output.f90:
  Take care that TKETENS output variable has the correct unit
• src_radiation.f90:
  
  • Do not use data types for exchg_boundaries, because local memory is involved
  • Optimizations in coe_th and coe_so for better vectorization
• src_slow_tendencies_rk.f90:
  For computation of lhfl_s use the correct lh_v or lh_s depending on interfacial temperature t_g. This changes results for variables LHFL_S, ALHFL_S
• src_soil_multlay.f90:
  
  • Modification of the maximum infiltration rate of surface water by hydraulic conductivity at saturation.
  • Implementation of the soil water transport by macro pores in the root zone with a vertical profile of the soil hydraulic conductivity at saturation.
  • Limit evapotranspiration if soil moisture is below 1.05*zpwp.
These modifications will change the results.

Back to Contents

9. Changes to the Namelists

Back to Contents

10. Changes of Results

The following changes do have an influence on the results:

• The change of the calling sequence of Assimilation and Relaxation changes the results.
• The enlarged horizontal, vertical and temporal check limits in redundancy check of SYNOP can change the results slightly.
• The bug fixes in the dynamics do influence the results.
• The fix in src_soil_multlay (limit evapotranspiration if soil moisture is below 1.05*zpwp) changes the results.
• The fix in src_slow_tendencies_rk.f90 changes the values of the fields LHFL_S, ALHFL_S
• The fix in src_radiation.f90 might influence the results, but depending on compiler
• The fix for the computation of the synthetic satellite images influences the values for the SYNSAT-products.
• The changes in src_conv_tiedtke, src_conv_shallow (adding an additional argument in subroutines cu_tied, cu_shallow) changes the results using the Cray compiler with option -hflex_mp=conservative. With -hflex=intolerant, the results are not changed. There are also no changes with gfortran. For the Cray compiler it must be an optimization issue

Back to Contents