Project

General

Profile

Bug #11341

emilnox produces a floating point exception

Added by Tommi Bergman about 3 years ago. Updated about 3 years ago.

Status:
Closed
Priority:
High
Assignee:
-
Category:
-
Target version:
Start date:
09/07/2018
Due date:
% Done:

100%


Description

During the CRESCENDO runs with EC-Earth I found that emission_nox.f90 causes a floating point exception due to division by zero or nan or even nan/nan (which one it is probably varies between occurrences).

itop = min(itop,itropo-1)
nlay = itop - itmin15  
if (nlay.le.0) then
   itmin15=itop-1
   nlay = 1
   if (lightning_output_2)      write(6,*) 'WARNING noxlight_cvp: itmin15>=itropo: ',i,j,itropo,itmin15,itop
endif !nlay le 0
! distributing according to airmass

airmass=0.
do l=itmin15+1,itop
   airmass = airmass + m(i,j,l)
enddo

do l=itmin15+1,itop
   emilig(i,j,l) = dsum*m(i,j,l)/airmass
enddo

Here the airmass is not well defined if calculation of troposphere height (index of the level) produces default value which is 1 :

Because it is assumed that lightning does not occur above tropopause, there is min(itop,itropo-1) which is to make sure that cloud top (itop) is 1 level below tropopause. Now when itropo=1 itop=0. Furthermore since part of the lightning nox is divided between level with the temperature less than 258K (itmin15) and itop we calculate the separation between these layers nlay.

In the special case when itropo=1 the difference between itop (=0) and itmin15 becomes zero or negative and it is used to decide that all lighting happens on one level.

Later in the loop level indices from itmin15+1 to itop, so from 0 to 0 is used with airmass(i,j,l) l being the iterator going from 0 to 0 here, which is of course not defined for normal fortran arrays.

It is rare, but still a special case for itropo=1 must be added, or decide if the itropo=1 is a good default value.

History

#1 Updated by Philippe Le Sager about 3 years ago

  • % Done changed from 0 to 10

The assumption (as coded) is that the tropopause is at itropo GT 2 to accommodate the two domains:

    do l=itmin15+1,itop   ! 70% of GC + 100% of IC
and
    do l=1,itmin15   ! 10% of GC
with itop=itropo-1 and itmin15=itropo-2 if they are too high. That means using:
itropo=max(2,itropo)
is not enough although it would not crash, because there is 10% of GC NOX not distributed.

Several ways to tackle that:
  1. set a minimum tropopause level of 3: itropo=max(3,itropo), which distributes 70% of GC in layer 2 and 10% in layer 1
  2. if tropopause level is below 0.8*ground pressure (defining the PBL where the last 20% of GC is distributed), skip the two domains above and distribute all GC and IC in the PBL
  3. if tropopause level is below 0.8*ground pressure (defining the PBL where the last 20% of GC is distributed), set itrop=level_pblh+1. With extra check: itrop=max(level_pblh+1,3)

The 3rd solution, compared to the 2nd one, adds the distinction between IC/GC.

#2 Updated by Twan van Noije about 3 years ago

I propose to set the default to 100 hPa as in IFS (see line 218 in CMIP6_STRATO_AERO_PROCESS).

#3 Updated by Philippe Le Sager about 3 years ago

  • Status changed from New to In Progress

We already set the default tropopause to the highest level where the pressure is above 100hPa (l.185 of base/toolbox.F90):

    DO l=1,lmx

       ! register highest troposheric level
       if (papf(l) .gt. 10000.) ltropo_ifs_default = max(ltropo_ifs_default,l)

That means, to get a tropopause in level 1, a lapse rate of 0.002 K/m between level 1 and level 1 + 2km was found, and not above.

#4 Updated by Twan van Noije about 3 years ago

You're right. However, I don't see how this routine can ever return a tropopause level of 1 given that it should be located between 550 and 75 hPa.

#5 Updated by Tommi Bergman about 3 years ago

At my current runs all instances with errors are located in grid boxes around [87,62] lon=82.5 lat=35 ( In Himalaya) and in the hourly output for 2012 the surface pressure at this location is usually less than 550 hPa during first two months of the year. Even after that it goes regularly below that.

#6 Updated by Twan van Noije about 3 years ago

OK. So, this brings us back to finding a solution to deal with these cases in the LNOx routine.

#7 Updated by Philippe Le Sager about 3 years ago

  • Status changed from In Progress to Resolved
  • % Done changed from 10 to 100

Resolved by setting itropo to 3 (if lower than 3) in r891.

#8 Updated by Philippe Le Sager about 3 years ago

  • Status changed from Resolved to Closed

#9 Updated by Philippe Le Sager about 3 years ago

  • Target version set to TM5-MP 3.0

Also available in: Atom PDF