emilnox produces a floating point exception
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.
#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 ICand
do l=1,itmin15 ! 10% of GCwith
itmin15=itropo-2if 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:
- set a minimum tropopause level of 3:
itropo=max(3,itropo), which distributes 70% of GC in layer 2 and 10% in layer 1
- 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
- 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:
The 3rd solution, compared to the 2nd one, adds the distinction between IC/GC.
#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
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.
#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.