diff --git a/A-INSTALL b/A-INSTALL
index 1a1a6d61efc83d7bb5d825e445cd638ca14b5d49..706b184e120f08545f02dd0d2703dd195fd41e01 100644
--- a/A-INSTALL
+++ b/A-INSTALL
@@ -1084,7 +1084,19 @@ etc ...
# c) MNH_ECRAD for optional compilation of new ECRAD radiative library from ECMWF
# --------------------------------------
#
-# The full ECRAD package was not included into the open source version of Meso-NH
+# The default version of ECRAD is 1.4.0 (open-source)
+#
+# Configure & Compilation
+export MNH_ECRAD=1
+./configure
+
+etc ...
+# The version of ECRAD is set by (by default):
+# export VER_ECRAD=140
+#
+# To use the previous version 1.0.1:
+#
+# The full ECRAD package 1.0.1 was not included into the open source version of Meso-NH
# because it needs a licence agrement.
#
# See here to get the licence & full sources : https://software.ecmwf.int/wiki/display/ECRAD/ECMWF+Radiation+Scheme+Home
@@ -1098,6 +1110,7 @@ tar xvfz ecrad-1.0.1.tar.gz
# Configure & Compilation
export MNH_ECRAD=1
+export VER_ECRAD=101
./configure
etc ...
@@ -1107,8 +1120,8 @@ etc ...
#
# Usage :
# 1) In namelist replace RAD='ECMW' by RAD='ECRA'
-# 2) Add link to all 'ecrad-1.0.1/data' files in your mesonh run directory
-ln -sf ${SRC_MESONH}/src/LIB/RAD/ecrad-1.0.1/data/* .
+# 2) Add link to all 'ecrad-1.X.X/data' files in your mesonh run directory
+ln -sf ${SRC_MESONH}/src/LIB/RAD/ecrad-1.X.X/data/* .
#
# REM : you can replace CDATADIR = "." by CDATADIR = "data" of ini_radiations_ecrad.f90 to link only the data folder instead of all the files one by one
#
diff --git a/README_MNH_CONDA b/README_MNH_CONDA
index d8b0bc0d7f56205e27f10717bcdd302162e0b205..5aa86aa5914c7678202a888b28399675e4c6d7fc 100644
--- a/README_MNH_CONDA
+++ b/README_MNH_CONDA
@@ -67,7 +67,7 @@ source ${MNH_MINICONDA}/miniconda3/etc/profile.d/conda.sh
conda create -n mnh_conda_cartopy_offlinedata
conda activate mnh_conda_cartopy_offlinedata
-conda install -c conda-forge netcdf4 cartopy matlibplot
+conda install -c conda-forge netcdf4 cartopy matplotlib
#
# this is the minimum packages needed for the MesoNH script of the test cases
diff --git a/src/LIB/MPIvide/mpi.h b/src/LIB/MPIvide/mpi.h
index 72a1b9a677de36d2a4923e7d4a9c2ae22108ab70..0676c6749831d4b43a2c5bd299881de5fe02efc7 100644
--- a/src/LIB/MPIvide/mpi.h
+++ b/src/LIB/MPIvide/mpi.h
@@ -4,6 +4,10 @@
* (C) 1993 by Argonne National Laboratory and Mississipi State University.
* All rights reserved. See COPYRIGHT in top-level directory.
*/
+/* Modifications:
+ P.Wautelet 19/11/2021: add MPI_LOGICAL4 and MPI_LOGICAL8 optional types
+! + modify MPI_REAL4/8 and MPI_INTEGER4/8
+*/
/* user include file for MPI programs */
@@ -94,6 +98,13 @@ typedef int MPI_Datatype;
#define MPI_2DOUBLE_PRECISION ((MPI_Datatype)33)
#define MPI_CHARACTER ((MPI_Datatype)1)
+#define MPI_INTEGER4 ((MPI_Datatype)34)
+#define MPI_INTEGER8 ((MPI_Datatype)35)
+#define MPI_LOGICAL4 ((MPI_Datatype)36)
+#define MPI_LOGICAL8 ((MPI_Datatype)37)
+#define MPI_REAL4 ((MPI_Datatype)38)
+#define MPI_REAL8 ((MPI_Datatype)39)
+
/* Communicators */
typedef int MPI_Comm;
#define MPI_COMM_WORLD 91
diff --git a/src/LIB/MPIvide/mpif.h b/src/LIB/MPIvide/mpif.h
index 97f91109c93e21ebb5e4b2f2131e79459448cb28..3557e3f48a6dc9f63d4bae9675f239d208534d7c 100644
--- a/src/LIB/MPIvide/mpif.h
+++ b/src/LIB/MPIvide/mpif.h
@@ -1,3 +1,6 @@
+!Modifications:
+! P.Wautelet 19/11/2021: add MPI_LOGICAL4 and MPI_LOGICAL8 optional types
+! + modify MPI_REAL4/8 and MPI_INTEGER4/8
!
!
! (C) 1993 by Argonne National Laboratory and Mississipi State University.
@@ -125,7 +128,7 @@
PARAMETER (MPI_2REAL=32,MPI_2DOUBLE_PRECISION=33,MPI_CHARACTER=1)
PARAMETER (MPI_BYTE=3,MPI_UB=16,MPI_LB=15,MPI_PACKED=14)
- INTEGER MPI_ORDER_C, MPI_ORDER_FORTRAN
+ INTEGER MPI_ORDER_C, MPI_ORDER_FORTRAN
PARAMETER (MPI_ORDER_C=56, MPI_ORDER_FORTRAN=57)
INTEGER MPI_DISTRIBUTE_BLOCK, MPI_DISTRIBUTE_CYCLIC
INTEGER MPI_DISTRIBUTE_NONE, MPI_DISTRIBUTE_DFLT_DARG
@@ -144,16 +147,19 @@
INTEGER MPI_INTEGER16
INTEGER MPI_REAL4, MPI_REAL8, MPI_REAL16
INTEGER MPI_COMPLEX8, MPI_COMPLEX16, MPI_COMPLEX32
+ INTEGER MPI_LOGICAL4, MPI_LOGICAL8
PARAMETER (MPI_INTEGER1=1,MPI_INTEGER2=4)
- PARAMETER (MPI_INTEGER4=6)
- PARAMETER (MPI_INTEGER8=13)
+ PARAMETER (MPI_INTEGER4=34)
+ PARAMETER (MPI_INTEGER8=35)
PARAMETER (MPI_INTEGER16=0)
- PARAMETER (MPI_REAL4=10)
- PARAMETER (MPI_REAL8=11)
+ PARAMETER (MPI_REAL4=38)
+ PARAMETER (MPI_REAL8=39)
PARAMETER (MPI_REAL16=0)
PARAMETER (MPI_COMPLEX8=23)
PARAMETER (MPI_COMPLEX16=24)
PARAMETER (MPI_COMPLEX32=0)
+ PARAMETER (MPI_LOGICAL4=36)
+ PARAMETER (MPI_LOGICAL8=37)
COMMON /MPIPRIV/ MPI_BOTTOM,MPI_STATUS_IGNORE,MPI_STATUSES_IGNORE
!
diff --git a/src/LIB/MPIvide/mpivide.c b/src/LIB/MPIvide/mpivide.c
index 313f162a8e4d4717add7bf80eeb03ddd0a144f20..eebe91b71a8e3178634fb98688f9118c8325037b 100644
--- a/src/LIB/MPIvide/mpivide.c
+++ b/src/LIB/MPIvide/mpivide.c
@@ -1,51 +1,39 @@
-/*
-MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+/*
+MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
MNH_LIC for details. version 1.
-*/
+*/
+/* Modifications :
+ P. Wautelet 19/11/2021: add function findtypesize + improve/add support for 32 and 64 bits variables + add mpi_reduce
+*/
+
+#include <stdio.h>
#include <string.h>
#include "mpi.h"
-/* Variables defined in meso-nh code */
-#ifdef FUJI
-#if MNH_REAL == 4
- #define MPI_PRECISION MPI_REAL
- #define MPI_2PRECISION MPI_2REAL
-#else
- #define MPI_PRECISION MPI_DOUBLE_PRECISION
- #define MPI_2PRECISION MPI_2DOUBLE_PRECISION
-#endif
-#else
- #define MPI_PRECISION MPI_REAL
- #define MPI_2PRECISION MPI_2REAL
-#endif
+/* MPI_INTEGER is defined in mpi.h */
+#define SIZEINTEGER4 4
+#define SIZEINTEGER8 8
-/* MPI_INTEGER is defined in mpi.h */
+#define SIZELOGICAL4 4
+#define SIZELOGICAL8 8
-#define MPI_INTEGER8 MPI_LONG_LONG_INT
+#define SIZEREAL4 4
+#define SIZEREAL8 8
-#ifdef FUJI
-#if MNH_INT == 8
+#if MNH_INT == 8
#define SIZEINTEGER 8
-#define SIZEINTEGER8 8
#define SIZELOGICAL 8
#else
#define SIZEINTEGER 4
-#define SIZEINTEGER8 8
#define SIZELOGICAL 4
#endif
#if MNH_REAL == 4
#define SIZEPRECISION 4
#define SIZE2PRECISION 8
#else
-#define SIZEPRECISION 8
-#define SIZE2PRECISION 16
-#endif
-#else
-#define SIZEINTEGER 8
-#define SIZEINTEGER8 8
#define SIZEPRECISION 8
#define SIZE2PRECISION 16
#endif
@@ -65,6 +53,57 @@ char *fct;
/* printf("MPIVIDE::Passage dans %s \n", fct); */
}
+int findtypesize(int type)
+{
+ int size;
+
+ switch(type)
+ {
+ case MPI_INTEGER4:
+ size = SIZEINTEGER4 ;
+ break;
+ case MPI_INTEGER8:
+ size = SIZEINTEGER8 ;
+ break;
+ case MPI_LOGICAL4:
+ size = SIZELOGICAL4 ;
+ break;
+ case MPI_LOGICAL8:
+ size = SIZELOGICAL8 ;
+ break;
+ case MPI_REAL4:
+ size = SIZEREAL4 ;
+ break;
+ case MPI_REAL8:
+ size = SIZEREAL8 ;
+ break;
+ case MPI_2REAL:
+ size = 2*SIZEPRECISION ;
+ break;
+ case MPI_2DOUBLE_PRECISION:
+ size = 2*SIZE2PRECISION ;
+ break;
+ case MPI_INTEGER:
+ size = SIZEINTEGER;
+ break;
+ case MPI_REAL:
+ size = SIZEPRECISION;
+ break;
+ case MPI_DOUBLEDOUBLE:
+ size = SIZE_DOUBLEDOUBLE;
+ break;
+ case MPI_LOGICAL:
+ size = SIZELOGICAL ;
+ break;
+ default:
+ printf("ERROR : unknown precision in findtypesize (MPIVIDE library)\n");
+ size = SIZEPRECISION;
+ break;
+ }
+
+ return size;
+}
+
#pragma weak mpi_cart_sub__ = mpi_cart_sub
#pragma weak mpi_cart_sub_ = mpi_cart_sub
void mpi_cart_sub
@@ -127,26 +166,10 @@ void mpi_alltoallv(void *sendbuf, int *sendcounts,
void *recvbuf, int *recvcounts,
int *rdispls, int *recvtype, int *comm, int *__ierr)
{
- int size = SIZE2PRECISION;
+ int size;
+
disppass("alltoallv");
- switch(*sendtype)
- {
- case MPI_INTEGER:
- size = SIZEINTEGER;
- break;
- case MPI_PRECISION:
- size = SIZEPRECISION;
- break;
- case MPI_2PRECISION:
- size = SIZE2PRECISION;
- break;
- case MPI_DOUBLEDOUBLE:
- size = SIZE_DOUBLEDOUBLE;
- break;
- case MPI_LOGICAL:
- size = SIZELOGICAL ;
- break;
- }
+ size = findtypesize(*sendtype);
memcpy(recvbuf, sendbuf, (*recvcounts)*size);
*__ierr = 0;
@@ -176,30 +199,10 @@ int *recvtype;
int *comm;
int *__ierr;
{
- int size = SIZE2PRECISION;
- disppass("allgatherv");
-
- switch(*sendtype)
- {
- case MPI_INTEGER:
- size = SIZEINTEGER;
- break;
- case MPI_INTEGER8:
- size = SIZEINTEGER8;
- break;
- case MPI_PRECISION:
- size = SIZEPRECISION;
- break;
- case MPI_2PRECISION:
- size = SIZE2PRECISION;
- break;
- case MPI_DOUBLEDOUBLE:
- size = SIZE_DOUBLEDOUBLE;
- break;
- case MPI_LOGICAL:
- size = SIZELOGICAL ;
- break;
- }
+ int size;
+
+ disppass("allgatherv");
+ size = findtypesize(*sendtype);
memcpy(recvbuf, sendbuf, (*recvcounts)*size);
*__ierr = 0;
}
@@ -219,29 +222,10 @@ int *root;
int *comm;
int *__ierr;
{
- int size = SIZE2PRECISION;
+ int size;
+
disppass("gather");
- switch(*sendtype)
- {
- case MPI_INTEGER:
- size = SIZEINTEGER;
- break;
- case MPI_PRECISION:
- size = SIZEPRECISION;
- break;
- case MPI_2PRECISION:
- size = SIZE2PRECISION;
- break;
- case MPI_DOUBLEDOUBLE:
- size = SIZE_DOUBLEDOUBLE;
- break;
- case MPI_DOUBLE:
- size = 8 ;
- break;
- case MPI_LOGICAL:
- size = SIZELOGICAL ;
- break;
- }
+ size = findtypesize(*sendtype);
memcpy(recvbuf, sendbuf, (*recvcount)*size);
*__ierr = 0;
@@ -263,26 +247,10 @@ int *root;
int *comm;
int *__ierr;
{
- int size = SIZE2PRECISION;
+ int size;
+
disppass("gatherv");
- switch(*sendtype)
- {
- case MPI_INTEGER:
- size = SIZEINTEGER;
- break;
- case MPI_PRECISION:
- size = SIZEPRECISION;
- break;
- case MPI_2PRECISION:
- size = SIZE2PRECISION;
- break;
- case MPI_DOUBLEDOUBLE:
- size = SIZE_DOUBLEDOUBLE;
- break;
- case MPI_LOGICAL:
- size = SIZELOGICAL ;
- break;
- }
+ size = findtypesize(*sendtype);
memcpy(recvbuf, sendbuf, (*recvcounts)*size);
*__ierr = 0;
@@ -352,26 +320,10 @@ int *recvtype;
int *comm;
int *__ierr;
{
- int size = SIZE2PRECISION;
+ int size;
+
disppass("allgather");
- switch(*sendtype)
- {
- case MPI_INTEGER:
- size = SIZEINTEGER;
- break;
- case MPI_PRECISION:
- size = SIZEPRECISION;
- break;
- case MPI_2PRECISION:
- size = SIZE2PRECISION;
- break;
- case MPI_DOUBLEDOUBLE:
- size = SIZE_DOUBLEDOUBLE;
- break;
- case MPI_LOGICAL:
- size = SIZELOGICAL ;
- break;
- }
+ size = findtypesize(*sendtype);
memcpy(recvbuf, sendbuf, (*recvcount)*size);
*__ierr = 0;
@@ -425,28 +377,32 @@ int *comm;
int *op;
int *__ierr;
{
- int size = SIZE2PRECISION;
+ int size;
disppass("allreduce");
- switch(*datatype)
- {
- case MPI_INTEGER:
- size = SIZEINTEGER;
- break;
- case MPI_PRECISION:
- size = SIZEPRECISION;
- break;
- case MPI_2PRECISION:
- size = SIZE2PRECISION;
- break;
- case MPI_DOUBLE:
- size = 8 ;
- break;
- case MPI_LOGICAL:
- size = SIZELOGICAL ;
- break;
- }
- memcpy(recvbuf, sendbuf, (*count)*size);
+ size = findtypesize(*datatype);
+ memcpy(recvbuf, sendbuf, (*count)*size);
+ *__ierr = 0;
+}
+
+#pragma weak mpi_reduce__ = mpi_reduce
+#pragma weak mpi_reduce_ = mpi_reduce
+void mpi_reduce
+( sendbuf, recvbuf, count, datatype, op, root, comm, __ierr )
+void *sendbuf;
+void *recvbuf;
+int *count;
+int *datatype;
+int *root;
+int *comm;
+int *op;
+int *__ierr;
+{
+ int size;
+
+ disppass("reduce");
+ size = findtypesize(*datatype);
+ memcpy(recvbuf, sendbuf, (*count)*size);
*__ierr = 0;
}
diff --git a/src/LIB/Python/Panel_Plot.py b/src/LIB/Python/Panel_Plot.py
index 8e095faf07fd46b6dcd4f01df5dddca6b0b5759a..328d26264ad1f209de65a0ecdc36f563e07dde70 100644
--- a/src/LIB/Python/Panel_Plot.py
+++ b/src/LIB/Python/Panel_Plot.py
@@ -20,33 +20,54 @@ import cartopy.feature as cfeature
class PanelPlot():
def __init__(self, nb_l, nb_c, Lfigsize, bigtitle, titlepad=40, minmaxpad=1.03, timepad=-0.06, lateralminmaxpad=0.86,
- labelcolorbarpad=6.0, colorbaraspect=20, colorbarpad=0.04 ):
- self.bigtitle = bigtitle
- self.Lfigsize = Lfigsize
- self.nb_l = nb_l
- self.nb_c = nb_c
+ labelcolorbarpad=6.0, colorbaraspect=20, colorbarpad=0.04, tickspad=0.8,
+ minmaxTextSize=10, bigtitleSize=13, titleSize=12,
+ xlabelSize=11, ylabelSize=11, timeSize=11, cbTicksLabelSize=11, cbTitleSize=11, xyTicksLabelSize=10, figBoxLinewidth=1,
+ xyTicksWidth=1, xyTicksLength=6):
+
+ self.bigtitle = bigtitle # Panel title
+ self.Lfigsize = Lfigsize # Panel size
+ self.nb_l = nb_l # Panel number of lines
+ self.nb_c = nb_c # Panel number of rows
self.nb_graph = 0 # New independent graph within the subplot
+
self.titlepad = titlepad # Title pad (vertical shift) from graph
+ self.tickspad = tickspad # Ticks pad (between ticks and axis label)
self.minmaxpad = minmaxpad # Min/Max print pad (vertical shift)
self.timepad = timepad # Time print pad (vertical shift)
self.colorbarpad = colorbarpad # Colorbar pad (horizontal shift from graph)
self.lateralminmaxpad = lateralminmaxpad
- self.labelcolorbarpad = labelcolorbarpad # Vertical colorbal label pad
- self.colorbaraspect = colorbaraspect # Ratio of long to short dimensions of colorbar w.r.t. the figure
-
+ self.labelcolorbarpad = labelcolorbarpad # Vertical colorbal label pad
+ self.colorbaraspect = colorbaraspect # Ratio of long to short dimensions of colorbar w.r.t. the figure
+
+ self.minmaxTextSize = minmaxTextSize # min/max text fontsize
+ self.bigtitleSize = bigtitleSize # Panel title fontsize
+ self.titleSize = titleSize # Graph title fontsize
+ self.xlabelSize = xlabelSize # X-label fontsize
+ self.ylabelSize = ylabelSize # Y-label fontsize
+ self.timeSize = timeSize # Time attribute of the graphs fontsize
+ self.cbTicksLabelSize = cbTicksLabelSize # Colorbar ticks label fontsize
+ self.cbTitleSize = cbTitleSize # Colorbar title fontsize
+ self.xyTicksLabelSize = xyTicksLabelSize # X/Y ticks label fontsize
+ self.figBoxLinewidth = figBoxLinewidth # Figure Box contour line width
+
+ self.xyTicksWidth = xyTicksWidth # Ticks width
+ self.xyTicksLength = xyTicksLength # Ticks length
+
# Initialization of the panel plots
self.fig = plt.figure(figsize=(self.Lfigsize[0],self.Lfigsize[1]))
self.fig.set_dpi(125)
self.fig.suptitle(self.bigtitle,fontsize=16)
- def save_graph(self, iplt, fig):
+ def save_graph(self, iplt, fig, fig_name='tempgraph'):
"""
- Create a temporary png file of (sub)-plot(s) which can be converted to PDF
+ Create a temporary png file of the panel plot which can be converted to PDF
"""
self.iplt = iplt
self.fig = fig
+ self.fig_name=fig_name # .png figure prefix name
- self.fig.savefig('tempgraph'+str(self.iplt)) #TODO possibility to change the default value of .png file name
+ self.fig.savefig(self.fig_name+str(self.iplt))
self.iplt+=1
return self.iplt
@@ -57,7 +78,7 @@ class PanelPlot():
self.drawCoastLines = drawCoastLines
self.projo = projo
- # Grid lines and labels
+ # Grid lines and labels
if 'PlateCarree' in str(projo):
gl = ax.gridlines(crs=self.projo, draw_labels=True, linewidth=1, color='gray')
if float(cartopy.__version__[:4]) >= 0.18:
@@ -67,7 +88,7 @@ class PanelPlot():
gl.xlabels_top = False
gl.ylabels_right = False
- # Coastlines
+ # Coastlines
if self.drawCoastLines and 'GeoAxes' in str(type(ax)):
ax.coastlines(resolution='10m')
@@ -75,14 +96,19 @@ class PanelPlot():
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.LAKES, alpha=0.7)
- def addWhitecm(self, Laddwhite, colormap_in, nb_level):
- if Laddwhite:
+ def addWhitecm(self, colormap_in, nb_level,whiteTop=False):
+ """
+ Add a white color at the top (whiteTop=True) or bottom of the colormap w.r.t. the number of independent colors used
+ """
color_map = cm.get_cmap(colormap_in, 256)
newcolor_map = color_map(np.linspace(0, 1, 256))
whites = np.array([1, 1, 1, 1]) #RBG code + opacity
- for i in range(int(256/nb_level)): newcolor_map[:i, :] = whites
+ if whiteTop:
+ for i in range(int(256/nb_level)): newcolor_map[-i, :] = whites
+ else:
+ for i in range(int(256/nb_level)): newcolor_map[:i, :] = whites
newcmp = ListedColormap(newcolor_map)
- return newcmp
+ return newcmp
def set_Title(self, ax, i, title, Lid_overlap, xlab, ylab):
@@ -97,10 +123,10 @@ class PanelPlot():
self.i = i
#self.ax[self.i].set_xlabel("test", fontweight='bold')
if not Lid_overlap:
- self.ax[self.i].set_title(title, pad=self.titlepad)
+ self.ax[self.i].set_title(title, pad=self.titlepad, fontsize=self.titleSize)
else: # If graph overlap, title is concatenated
new_title = self.ax[self.i].get_title() + ' and ' + title
- self.ax[self.i].set_title(new_title, pad=self.titlepad)
+ self.ax[self.i].set_title(new_title, pad=self.titlepad, fontsize=self.titleSize)
def set_xlim(self, ax, i, xlim):
"""
@@ -110,7 +136,7 @@ class PanelPlot():
self.xlim = xlim
self.i = i
- self.ax[self.i].set_xlim(xlim[0],xlim[1])
+ self.ax[self.i].set_xlim(xlim[0],xlim[1])#, fontsize=self.xlabelSize)
def set_ylim(self, ax, i, ylim):
"""
@@ -120,7 +146,7 @@ class PanelPlot():
self.ylim = ylim
self.i = i
- self.ax[self.i].set_ylim(ylim[0],ylim[1])
+ self.ax[self.i].set_ylim(ylim[0],ylim[1])#, fontsize=self.ylabelSize)
def set_XYaxislab(self, ax, i, xlab, ylab):
"""
@@ -136,14 +162,14 @@ class PanelPlot():
# https://github.com/SciTools/cartopy/issues/1332
if 'GeoAxes' in str(type(self.ax[self.i])):
self.ax[self.i].text(-0.11, 0.45, ylab, verticalalignment='top', horizontalalignment='left',
- rotation='vertical', rotation_mode='anchor', transform=self.ax[self.i].transAxes, color='black', fontsize=11)
+ rotation='vertical', rotation_mode='anchor', transform=self.ax[self.i].transAxes, color='black', fontsize=self.ylabelSize)
self.ax[self.i].text(0.45, -0.06, xlab, verticalalignment='top', horizontalalignment='left',
- rotation='horizontal', rotation_mode='anchor', transform=self.ax[self.i].transAxes, color='black', fontsize=11)
+ rotation='horizontal', rotation_mode='anchor', transform=self.ax[self.i].transAxes, color='black', fontsize=self.xlabelSize)
else:
- self.ax[self.i].set_xlabel(xlab)
- self.ax[self.i].set_ylabel(ylab)
+ self.ax[self.i].set_xlabel(xlab, fontsize=self.xlabelSize, labelpad=0.1)
+ self.ax[self.i].set_ylabel(ylab, fontsize=self.ylabelSize, labelpad=0.1)
- def addLine(self, ax, beg_coord, end_coord, color='black', linewidth=1):
+ def addLine(self, ax, beg_coord, end_coord, color='black', linewidth=0.2):
self.ax = ax
self.beg_coord = beg_coord
self.end_coord = end_coord
@@ -171,10 +197,10 @@ class PanelPlot():
strtext = " min = " + "{:.3e}".format(np.nanmin(var*facconv)) + " max = " + "{:.3e}".format(np.nanmax(var*facconv))
if not Lid_overlap:
self.ax[self.i].text(0.01, self.minmaxpad, strtext, verticalalignment='top', horizontalalignment='left',
- transform=self.ax[self.i].transAxes, color='black', fontsize=10)
+ transform=self.ax[self.i].transAxes, color='black', fontsize=self.minmaxTextSize)
else:
self.ax[self.i].text(self.lateralminmaxpad, self.minmaxpad, strtext, verticalalignment='top', horizontalalignment='right',
- transform=self.ax[self.i].transAxes, color='black', fontsize=10)
+ transform=self.ax[self.i].transAxes, color='black', fontsize=self.minmaxTextSize)
# Print to help choose min/max value for ticks
self.print_minmax(var*facconv, title)
@@ -188,12 +214,12 @@ class PanelPlot():
strtext = "Time = " + timetxt
self.ax[self.i].text(0.0, self.timepad, strtext, verticalalignment='top', horizontalalignment='left',
- transform=self.ax[self.i].transAxes, color='black', fontsize=10)
+ transform=self.ax[self.i].transAxes, color='black', fontsize=self.timeSize)
def psectionV(self, Lxx=[], Lzz=[], Lvar=[], Lxlab=[], Lylab=[], Ltitle=[], Lminval=[], Lmaxval=[],
Lstep=[], Lstepticks=[], Lcolormap=[], Lcbarlabel=[], LcolorLine=[],
- Lfacconv=[], ax=[], Lid_overlap=[], colorbar=True, orog=[], Lxlim=[], Lylim=[], Ltime=[], Lpltype=[], LaddWhite_cm=[]):
+ Lfacconv=[], ax=[], Lid_overlap=[], colorbar=True, orog=[], Lxlim=[], Lylim=[], Ltime=[], Lpltype=[], LaddWhite_cm=[], LwhiteTop=[]):
"""
Horizontal cross section plot
Parameters :
@@ -216,9 +242,10 @@ class PanelPlot():
- Lfacconv : List of factors for unit conversion of each variables
- ax : List of fig.axes for ploting multiple different types of plots in a subplot panel
- Lid_overlap: List of number index of plot to overlap current variables
- - Lpltype : List of types of plot 'cf' or 'c'. cf=contourf, c=contour (lines only)
+ - Lpltype : List of types of plot 'cf' or 'c'. cf=contourf, c=contour (lines only)
- colorbar : show colorbar or not
- - LaddWhite_cm : List of boolean to add white color to a colormap at the first (low value) tick colorbar
+ - LaddWhite_cm : List of boolean to add white color to a colormap at the last bottom (low value) tick colorbar
+ - LwhiteTop : List of boolean to add the white color at the first top (high value). If false, the white is added at the bottom if Laddwhite_cm=T
"""
self.ax = ax
firstCall = (len(self.ax) == 0)
@@ -235,6 +262,7 @@ class PanelPlot():
if not Lcolormap: LcolorLine=['black']*len(Lvar)
if not Lpltype: Lpltype=['cf']*len(Lvar)
if not LaddWhite_cm: LaddWhite_cm=[False]*len(Lvar)
+ if not LwhiteTop: LwhiteTop=[False]*len(Lvar)
if not Lylab: Lylab = ['']*len(Lvar)
if not Lxlab: Lxlab = ['']*len(Lvar)
# Add an extra percentage of the top max value for forcing the colorbar show the true user maximum value (correct a bug)
@@ -261,7 +289,7 @@ class PanelPlot():
# Print time validity
if Ltime: self.showTimeText(self.ax, iax, str(Ltime[i]))
-
+
# Number of contours level
if not Lstep[i]: # Default value of number of steps is 20
Lstep[i] = (Lmaxval[i] - Lminval[i])/20
@@ -270,16 +298,16 @@ class PanelPlot():
levels_contour = np.arange(Lminval[i],Lmaxval[i],step=Lstep[i])
# Add White to colormap
- if LaddWhite_cm[i] and Lcolormap: Lcolormap[i]=self.addWhitecm(LaddWhite_cm[i], Lcolormap[i], len(levels_contour))
+ if LaddWhite_cm[i] and Lcolormap: Lcolormap[i]=self.addWhitecm(Lcolormap[i], len(levels_contour), LwhiteTop[i])
# Plot
if Lpltype[i]=='c': # Contour
if LcolorLine:
cf = self.ax[iax].contour(Lxx[i], Lzz[i], var*Lfacconv[i], levels=levels_contour,
- norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], colors=LcolorLine[i])
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], colors=LcolorLine[i], linewidths=0.1)
else:
cf = self.ax[iax].contour(Lxx[i], Lzz[i], var*Lfacconv[i], levels=levels_contour,
- norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i])
+ norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i], linewidths=0.1)
else: # Contourf
cf = self.ax[iax].contourf(Lxx[i], Lzz[i], var*Lfacconv[i], levels=levels_contour,
norm=norm, vmin=Lminval[i], vmax=Lmaxval[i], cmap=Lcolormap[i])
@@ -289,6 +317,13 @@ class PanelPlot():
# X/Y Axis label
self.set_XYaxislab(self.ax, iax, Lxlab[i], Lylab[i])
+
+ # Ticks label
+ self.ax[iax].tick_params(axis='both', labelsize=self.xyTicksLabelSize, width=self.xyTicksWidth, length=self.xyTicksLength, pad=self.tickspad)
+
+ # Bounding box of the plot line width
+ for axis in ['top','bottom','left','right']:
+ self.ax[iax].spines[axis].set_linewidth(self.figBoxLinewidth)
# X/Y Axis limits value
if Lxlim:
@@ -304,17 +339,17 @@ class PanelPlot():
# Color label on contour-line
if Lpltype[i]=='c': # Contour
- self.ax[iax].clabel(cf)
- #self.ax[iax].clabel(cf, levels=np.arange(Lminval[i],Lmaxval[i],step=Lstep[i])) #TODO bug, levels not recognized
+ self.ax[iax].clabel(cf, fontsize=self.cbTicksLabelSize)
+ #self.ax[iax].clabel(cf, levels=np.arange(Lminval[i],Lmaxval[i],step=Lstep[i]), fontsize=self.cbTicksLabelSize) #TODO bug, levels not recognized
#Filling area under topography
if not orog==[]:
- self.ax[iax].fill_between(Lxx[i][0,:], orog, color='black')
+ self.ax[iax].fill_between(Lxx[i][0,:], orog, color='black', linewidth=0.2)
# Colorbar
if colorbar:
cb=plt.colorbar(cf, ax=self.ax[iax], fraction=0.031, pad=self.colorbarpad, ticks=np.arange(Lminval[i],Lmaxval[i],Lstepticks[i]), aspect=self.colorbaraspect)
- cb.ax.set_title(Lcbarlabel[i], pad = self.labelcolorbarpad, loc='left') #This creates a new AxesSubplot only for the colorbar y=0 ==> location at the bottom
+ cb.ax.set_title(Lcbarlabel[i], pad=self.labelcolorbarpad, loc='left', fontsize=self.cbTitleSize) # This creates a new AxesSubplot only for the colorbar y=0 ==> location at the bottom
return self.fig
@@ -324,7 +359,7 @@ class PanelPlot():
"""
XY (multiple)-lines plot
Parameters :
- - Lxx : List of variables to plot or coordinates along the X axis #TODO : ajouter Lfacconv pour les deux axes. Impact tous les cas test avec lignes X/Y
+ - Lxx : List of variables to plot or coordinates along the X axis
- Lyy : List of variables to plot or coordinates along the Y axis
- Lxlab : List of x-axis label
- Lylab : List of y-axis label
@@ -410,7 +445,7 @@ class PanelPlot():
def psectionH(self, lon=[],lat=[], Lvar=[], Lcarte=[], Llevel=[], Lxlab=[], Lylab=[], Ltitle=[], Lminval=[], Lmaxval=[],
Lstep=[], Lstepticks=[], Lcolormap=[], LcolorLine=[], Lcbarlabel=[], Lproj=[], Lfacconv=[], coastLines=True, ax=[],
- Lid_overlap=[], colorbar=True, Ltime=[], LaddWhite_cm=[], Lpltype=[], Lcbformatlabel=[]):
+ Lid_overlap=[], colorbar=True, Ltime=[], LaddWhite_cm=[], LwhiteTop=[], Lpltype=[], Lcbformatlabel=[]):
"""
Horizontal cross section plot
Parameters :
@@ -438,6 +473,7 @@ class PanelPlot():
- colorbar : show colorbar or not
- Lpltype : List of types of plot 'cf' or 'c'. cf=contourf, c=contour (lines only)
- LaddWhite_cm : List of boolean to add white color to a colormap at the first (low value) tick colorbar
+ - LwhiteTop : List of boolean to add the white color at the first top (high value). If false, the white is added at the bottom if Laddwhite_cm=T
- Lcbformatlabel: List of boolean to reduce the format to exponential 1.1E+02 format colorbar label
"""
self.ax = ax
@@ -457,6 +493,7 @@ class PanelPlot():
if not Lcolormap: LcolorLine=['black']*len(Lvar)
if not Lpltype: Lpltype=['cf']*len(Lvar)
if not LaddWhite_cm: LaddWhite_cm=[False]*len(Lvar)
+ if not LwhiteTop: LwhiteTop=[False]*len(Lvar)
if not Lcbformatlabel: Lcbformatlabel=[False]*len(Lvar)
# Add an extra percentage of the top max value for forcing the colorbar show the true user maximum value (correct a bug)
if Lstep: Lmaxval = list(map(lambda x, y: x + 1E-6*y, Lmaxval, Lstep) ) #The extra value is 1E-6 times the step ticks of the colorbar
@@ -511,7 +548,7 @@ class PanelPlot():
levels_contour = np.arange(Lminval[i],Lmaxval[i],step=Lstep[i])
# Add White to colormap
- if LaddWhite_cm[i] and Lcolormap: Lcolormap[i]=self.addWhitecm(LaddWhite_cm[i], Lcolormap[i], len(levels_contour))
+ if LaddWhite_cm[i] and Lcolormap: Lcolormap[i]=self.addWhitecm(Lcolormap[i], len(levels_contour), LwhiteTop[i])
# Plot
if Lproj:
@@ -544,14 +581,14 @@ class PanelPlot():
# Color label on contour-line
if Lpltype[i]=='c': # Contour
if 'GeoAxes' in str(type(self.ax[self.i])): # cartopy does not like the levels arguments in clabel, known issue
- self.ax[iax].clabel(cf)
+ self.ax[iax].clabel(cf, fontsize=self.cbTicksLabelSize)
else:
- self.ax[iax].clabel(cf, levels=np.arange(Lminval[i],Lmaxval[i],step=Lstep[i]))
+ self.ax[iax].clabel(cf, levels=np.arange(Lminval[i],Lmaxval[i],step=Lstep[i]), fontsize=self.cbTicksLabelSize)
# Colorbar
if colorbar:
cb=plt.colorbar(cf, ax=self.ax[iax], fraction=0.031, pad=self.colorbarpad, ticks=np.arange(Lminval[i],Lmaxval[i],Lstepticks[i]), aspect=self.colorbaraspect)
- cb.ax.set_title(Lcbarlabel[i], pad = self.labelcolorbarpad, loc='left') #This creates a new AxesSubplot only for the colorbar y=0 ==> location at the bottom
+ cb.ax.set_title(Lcbarlabel[i], pad = self.labelcolorbarpad, loc='left', fontsize=self.cbTitleSize) #This creates a new AxesSubplot only for the colorbar y=0 ==> location at the bottom
if Lcbformatlabel[i]: cb.ax.set_yticklabels(["{:.1E}".format(i) for i in cb.get_ticks()])
return self.fig
diff --git a/src/LIB/Python/misc_functions.py b/src/LIB/Python/misc_functions.py
index aa91aa7a311dc721c17f3f28cc7a8fd4f78df5b1..efc0916e04a60a7ee3cd2a276d1bb1eba3de0182 100644
--- a/src/LIB/Python/misc_functions.py
+++ b/src/LIB/Python/misc_functions.py
@@ -67,12 +67,12 @@ def windvec_verti_proj(u, v, level, angle):
return projected_wind
def oblique_proj(var, ni, nj, lvl, i_beg, j_beg, i_end, j_end):
- """Compute an oblique projection of a 3D variable w.r.t. its axes
+ """Compute an oblique projection of a variable w.r.t. its axes
Parameters
----------
- var : array 3D
- the 3D variable to project (e.g. THT)
+ var : array 3D or 2D
+ the variable to project (e.g. THT, ZS)
ni : array 1D
1D x-axis of the 3D dimension
@@ -94,27 +94,37 @@ def oblique_proj(var, ni, nj, lvl, i_beg, j_beg, i_end, j_end):
angle_proj : float
the angle (radian) of the new axe w.r.t the x/ni axes (West-East)
- out_var : array 2D
- a 2D (z,m) variable projected on the oblique axe
+ out_var : array 2D or 1D
+ a 2D (z,m) or 1D (m) variable projected on the oblique axe
axe_m : array 1D
a 1D m new axe (distance from the beggining point)
"""
- dist_seg=np.sqrt((i_end-i_beg)**2.0 + (j_end-j_beg)**2.0) # Distance de la section oblique m
- out_var = np.zeros((len(lvl),int(dist_seg)+1)) # Initialisation du nouveau champs projeté dans la coupe (z,m)
- axe_m = np.zeros(int(dist_seg)+1) #Axe des abscisses qui sera tracé selon la coupe
- axe_m_coord = [] #Coordonnées x,y des points qui composent l'axe
- axe_m_coord.append( (ni[i_beg],nj[j_beg]) ) #Le premier point est celui donné par l'utilisateur
- for m in range(int(dist_seg)): #Discrétisation selon distance de la coupe / int(distance_de_la_coupe)
+ dist_seg=np.sqrt((i_end-i_beg)**2.0 + (j_end-j_beg)**2.0) # Distance de la section oblique m
+ if var.ndim ==3:
+ out_var = np.zeros((len(lvl),int(dist_seg)+1)) # Initialisation du nouveau champs projeté dans la coupe (z,m)
+ else: # 2D
+ out_var = np.zeros(int(dist_seg)+1) # Initialisation du nouveau champs projeté dans la coupe (m)
+
+ axe_m = np.zeros(int(dist_seg)+1) # Axe des abscisses qui sera tracé selon la coupe
+ axe_m_coord = [] # Coordonnées x,y des points qui composent l'axe
+ axe_m_coord.append( (ni[i_beg],nj[j_beg]) ) # Le premier point est celui donné par l'utilisateur
+
+ for m in range(int(dist_seg)): # Discrétisation selon distance de la coupe / int(distance_de_la_coupe)
axe_m_coord.append( (axe_m_coord[0][0] + (ni[i_end]-ni[i_beg])/(int(dist_seg))*(m+1),
axe_m_coord[0][1] + (nj[j_end]-nj[j_beg])/(int(dist_seg))*(m+1) ))
axe_m[m+1] = np.sqrt((ni[i_beg]-axe_m_coord[m+1][0])**2 + (nj[j_beg]-axe_m_coord[m+1][1])**2)
-
- for k in range(len(lvl)):
- a=RectBivariateSpline(ni, nj,var[k,:,:],kx=1,ky=1) #Interpolation par niveau à l'ordre 1 pour éviter des valeurs négatives de champs strictement > 0
+
+ if var.ndim ==3: # 3D variables to project
+ for k in range(len(lvl)):
+ a=RectBivariateSpline(ni, nj,var[k,:,:],kx=1,ky=1) # Interpolation par niveau à l'ordre 1 pour éviter des valeurs négatives de champs strictement > 0
+ for m in range(int(dist_seg)+1):
+ out_var[k,m] = a.ev(axe_m_coord[m][0],axe_m_coord[m][1]) # La fonction ev de RectBivariate retourne la valeur la plus proche du point considéré
+ else: # 2D variables to project
+ a=RectBivariateSpline(ni, nj,var[:,:],kx=1,ky=1)
for m in range(int(dist_seg)+1):
- out_var[k,m] = a.ev(axe_m_coord[m][0],axe_m_coord[m][1]) # La fonction ev de RectBivariate retourne la valeur la plus proche du point considéré
+ out_var[m] = a.ev(axe_m_coord[m][0],axe_m_coord[m][1])
angle_proj = math.acos((ni[i_end]-ni[i_beg])/axe_m[-1])
return angle_proj, out_var, axe_m
@@ -206,4 +216,4 @@ def comp_altitude2DVar(oneVar3D, orography, ztop, level, n_y, n_x):
else:
for k in range(len(level)):
altitude[k,j,i] = orography[j,i] + level[k]*((ztop-orography[j,i])/ztop)
- return altitude, n_x3D, n_y3D
\ No newline at end of file
+ return altitude, n_x3D, n_y3D
diff --git a/src/LIB/Python/read_MNHfile.py b/src/LIB/Python/read_MNHfile.py
index 6c1311c1730eddb20446c428ddedc1de2d683cf2..fba716062817184b910fdab3d1bcf2d7bd414082 100644
--- a/src/LIB/Python/read_MNHfile.py
+++ b/src/LIB/Python/read_MNHfile.py
@@ -11,6 +11,47 @@ MNH_LIC for details. version 1.
import netCDF4 as nc
import numpy as np
+def read_withEPY(LnameFiles,Dvar_input, Dvar_output={}, path='.'):
+ import epygram
+ epygram.init_env()
+ for i,keyFiles in enumerate(Dvar_input.keys()):
+ print('Reading file ' + keyFiles)
+ theFile = epygram.formats.resource(LnameFiles[i],'r')
+ Dvar_output[keyFiles] = {} #initialize dic for each files
+ for var in Dvar_input[keyFiles]: #For each files
+ # Read variables
+ if(theFile.format == 'FA'):
+ Dvar_output[keyFiles][var] = theFile.readfield(var)
+ elif(theFile.format == 'LFI'):
+ if(var[1]==None or var[1]==0): # 2D Field
+ Dvar_output[keyFiles][var[0]] = theFile.readfield(var)
+ else: # 3D Field
+ Dvar_output[keyFiles][var[0]+str(var[1])] = theFile.readfield(var).getlevel(k=var[1])
+ elif(theFile.format == 'netCDFMNH'):
+ if(var[1]==None or var[1]==0): # 2D Field
+ Dvar_output[keyFiles][var[0]] = theFile.readfield(var[0])
+ else:
+ Dvar_output[keyFiles][var[0]+str(var[1])] = theFile.readfield(var[0]).getlevel(k=var[1])
+ elif(theFile.format == 'GRIB'):
+ if len(var)==6: # GRIB2
+ Dvar_output[keyFiles][var[5]] = theFile.readfield({'discipline': var[0], 'parameterCategory': var[1], 'typeOfFirstFixedSurface': var[2],'parameterNumber': var[3], 'level': var[4]})
+ elif len(var)==5: # GRIB1
+ Dvar_output[keyFiles][var[4]] = theFile.readfield({'indicatorOfParameter': var[0], 'paramId': var[1], 'indicatorOfTypeOfLevel': var[2], 'level': var[3]})
+ else: epygramError("GRIB format error. GRIB1 expects 4 values : [indicatorOfParameter, paramId, indicatorOfTypeOfLevel, level, 'casual name'], GRIB2 expects 5 values [discipline, parameterCategory, typeOfFirstFixedSurface, parameterNumber, level, casual name]")
+ else:
+ raise epygramError("Unknown format file, please use FA, LFI, GRIB or MNH NetCDF")
+ theFile.close()
+
+ # Transform spectral data to physics space (for AROME and ARPEGE)
+ for f in Dvar_output:
+ for var in Dvar_output[f]:
+ try:
+ if(Dvar_output[f][var].spectral):
+ Dvar_output[f][var].sp2gp()
+ except:
+ break
+ return Dvar_output
+
def read_netcdf(LnameFiles, Dvar_input, path='.', get_data_only=True, del_empty_dim=True, removeHALO=True):
"""Read a netCDF4 Meso-NH file
For each file, call functions to read diachronic or synchronous file
@@ -20,9 +61,9 @@ def read_netcdf(LnameFiles, Dvar_input, path='.', get_data_only=True, del_empty_
LnameFiles : list of str
list of Meso-NH netCDF4 file (diachronic or synchronous)
- Dvar_input : Dict{'fileNumber' : 'var_name',('group_name','var_name')}
+ Dvar_input : Dict{'keyFile' : 'var_name',('group_name','var_name')}
where
- 'fileNumber' is a str corresponding to 'f' + the file number in LnameFiles (by order)
+ 'keyFile' is a str corresponding to a key for the file number in LnameFiles (by order)
'var_name' is the exact str of the netCDF4 variable name
('group_name','var_name') is the exact tuple of the (sub-)groups name and the netCDF4 variable name
e.g. : {'f1':['ZS', 'WT','ni', 'level'],
@@ -51,19 +92,18 @@ def read_netcdf(LnameFiles, Dvar_input, path='.', get_data_only=True, del_empty_
Dvar[ifile][('group_name','var_name')] if the group contains more than one variable
"""
Dvar = {}
- for i,nameFiles in enumerate(LnameFiles):
- f_nb = 'f' + str(i+1)
- print('Reading file ' + f_nb)
- print(path + nameFiles)
- theFile = nc.Dataset(path + nameFiles,'r')
- Dvar[f_nb] = {}
- if '000' in nameFiles[-6:-3]:
+ for i,keyFiles in enumerate(Dvar_input.keys()):
+ print('Reading file ' + keyFiles)
+ print(path + LnameFiles[i])
+ theFile = nc.Dataset(path + LnameFiles[i],'r')
+ Dvar[keyFiles] = {}
+ if '000' in LnameFiles[i][-6:-3]:
if theFile['MASDEV'][0] <= 54:
raise TypeError('The python lib is available for MNH >= 5.5')
else:
- Dvar[f_nb] = read_TIMESfiles_55(theFile, Dvar_input[f_nb], Dvar[f_nb], get_data_only, del_empty_dim, removeHALO)
+ Dvar[keyFiles] = read_TIMESfiles_55(theFile, Dvar_input[keyFiles], Dvar[keyFiles], get_data_only, del_empty_dim, removeHALO)
else:
- Dvar[f_nb]= read_BACKUPfile(theFile, Dvar_input[f_nb], Dvar[f_nb], get_data_only, del_empty_dim, removeHALO)
+ Dvar[keyFiles]= read_BACKUPfile(theFile, Dvar_input[keyFiles], Dvar[keyFiles], get_data_only, del_empty_dim, removeHALO)
theFile.close()
return Dvar
diff --git a/src/LIB/RAD/data_mnh/spectral_albedo.nc b/src/LIB/RAD/data_mnh/spectral_albedo.nc
new file mode 100644
index 0000000000000000000000000000000000000000..d6e6335e8a63baeeafd7fdf57fac27eff43bfbb7
Binary files /dev/null and b/src/LIB/RAD/data_mnh/spectral_albedo.nc differ
diff --git a/src/LIB/RAD/data_mnh/spectral_emissivity.nc b/src/LIB/RAD/data_mnh/spectral_emissivity.nc
new file mode 100644
index 0000000000000000000000000000000000000000..14ac4ae328034fb4ef41b4850e05b80faa51dd2d
Binary files /dev/null and b/src/LIB/RAD/data_mnh/spectral_emissivity.nc differ
diff --git a/src/LIB/SURCOUCHE/src/mode_field.f90 b/src/LIB/SURCOUCHE/src/mode_field.f90
index 97466a414e1dbb88e5fb8746cb10d3e65706df1c..ce06ed50bf7d80102600d3423136c97d754915b4 100644
--- a/src/LIB/SURCOUCHE/src/mode_field.f90
+++ b/src/LIB/SURCOUCHE/src/mode_field.f90
@@ -2002,6 +2002,20 @@ ALLOCATE(TFIELDLIST(IDX)%TFIELD_X0D(IMODEL))
IDX = IDX+1
!
IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
+TFIELDLIST(IDX)%CMNHNAME = 'DRYMASSS'
+TFIELDLIST(IDX)%CSTDNAME = ''
+TFIELDLIST(IDX)%CLONGNAME = 'DRYMASSS'
+TFIELDLIST(IDX)%CUNITS = 'kg'
+TFIELDLIST(IDX)%CDIR = '--'
+TFIELDLIST(IDX)%CCOMMENT = 'Total Dry Mass Source'
+TFIELDLIST(IDX)%NGRID = 0
+TFIELDLIST(IDX)%NTYPE = TYPEREAL
+TFIELDLIST(IDX)%NDIMS = 0
+TFIELDLIST(IDX)%LTIMEDEP = .TRUE.
+ALLOCATE(TFIELDLIST(IDX)%TFIELD_X0D(IMODEL))
+IDX = IDX+1
+!
+IF(IDX>MAXFIELDS) CALL ERR_INI_FIELD_LIST()
TFIELDLIST(IDX)%CMNHNAME = 'BL_DEPTH'
TFIELDLIST(IDX)%CSTDNAME = ''
TFIELDLIST(IDX)%CLONGNAME = 'BL_DEPTH'
@@ -3740,6 +3754,13 @@ IF (.NOT.ASSOCIATED(XDRYMASST)) THEN
TFIELDLIST(IID)%TFIELD_X0D(1)%DATA=>XDRYMASST
END IF
!
+IF (.NOT.ASSOCIATED(XDRYMASSS)) THEN
+ CALL PRINT_MSG(NVERB_DEBUG,'GEN','INI_FIELD_SCALARS',' XDRYMASSS was not associated')
+ ALLOCATE(XDRYMASSS)
+ CALL FIND_FIELD_ID_FROM_MNHNAME('DRYMASSS',IID,IRESP)
+ TFIELDLIST(IID)%TFIELD_X0D(1)%DATA=>XDRYMASSS
+END IF
+!
IF (.NOT.ASSOCIATED(NRIMX)) THEN
CALL PRINT_MSG(NVERB_DEBUG,'GEN','INI_FIELD_SCALARS',' NRIMX was not associated')
ALLOCATE(NRIMX)
@@ -4341,6 +4362,14 @@ IF (.NOT.ASSOCIATED(TFIELDLIST(IID)%TFIELD_X0D(KTO)%DATA)) THEN
END IF
XDRYMASST => TFIELDLIST(IID)%TFIELD_X0D(KTO)%DATA
!
+CALL FIND_FIELD_ID_FROM_MNHNAME('DRYMASSS',IID,IRESP)
+IF (.NOT.ASSOCIATED(TFIELDLIST(IID)%TFIELD_X0D(KTO)%DATA)) THEN
+ CALL PRINT_MSG(NVERB_DEBUG,'GEN','FIELDLIST_GOTO_MODEL',&
+ 'TFIELDLIST(IID)%TFIELD_X0D(KTO)%DATA was not associated for '//TFIELDLIST(IID)%CMNHNAME)
+ ALLOCATE(TFIELDLIST(IID)%TFIELD_X0D(KTO)%DATA)
+END IF
+XDRYMASSS => TFIELDLIST(IID)%TFIELD_X0D(KTO)%DATA
+!
! MODD_DYN_n variables
!
CALL FIND_FIELD_ID_FROM_MNHNAME('RIMX',IID,IRESP)
diff --git a/src/LIB/SURCOUCHE/src/mode_mppdb.f90 b/src/LIB/SURCOUCHE/src/mode_mppdb.f90
index 28b3d7cbbb229b182a0f30d753e5def5ab16b47e..3fb6e15422f287e14ebde48979ed3b19b3d1e276 100644
--- a/src/LIB/SURCOUCHE/src/mode_mppdb.f90
+++ b/src/LIB/SURCOUCHE/src/mode_mppdb.f90
@@ -70,6 +70,7 @@ MODULE MODE_MPPDB
INTERFACE MPPDB_CHECK1D
MODULE PROCEDURE MPPDB_CHECK1D_REAL
END INTERFACE
+ LOGICAL :: MPPDB_ACTIVED = .FALSE.
INTERFACE MPPDB_CHECK2D
MODULE PROCEDURE MPPDB_CHECK2D_REAL
@@ -125,6 +126,7 @@ MODULE MODE_MPPDB
IF (MPPDB_INITIALIZED) RETURN
!
MPPDB_INITIALIZED = .TRUE.
+ MPPDB_ACTIVED = .TRUE.
!
! Init MPI
!
@@ -284,8 +286,28 @@ MODULE MODE_MPPDB
END SUBROUTINE MPPDB_INIT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ SUBROUTINE MPPDB_START_DEBUG()
+ IMPLICIT NONE
+ MPPDB_ACTIVED = .TRUE.
+ END SUBROUTINE MPPDB_START_DEBUG
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+ SUBROUTINE MPPDB_STOP_DEBUG()
+ IMPLICIT NONE
+ MPPDB_ACTIVED = .FALSE.
+ END SUBROUTINE MPPDB_STOP_DEBUG
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ SUBROUTINE MPPDB_GET_ACTIVED(OACTIVE)
+ IMPLICIT NONE
+ LOGICAL , INTENT(OUT) :: OACTIVE
+ OACTIVE = MPPDB_ACTIVED
+ END SUBROUTINE MPPDB_GET_ACTIVED
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ SUBROUTINE MPPDB_SET_ACTIVED(OACTIVE)
+ IMPLICIT NONE
+ LOGICAL , INTENT(IN) :: OACTIVE
+ MPPDB_ACTIVED = OACTIVE
+ END SUBROUTINE MPPDB_SET_ACTIVED
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE MPPDB_BARRIER()
#ifdef MNH_SP4
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
@@ -296,7 +318,7 @@ MODULE MODE_MPPDB
!
! synchronize all father & sons
!
- IF ( .NOT. MPPDB_INITIALIZED ) RETURN
+ IF (( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED )) RETURN
!
CALL MPI_BARRIER(MPPDB_INTRA_COMM,IERR)
!
@@ -334,7 +356,7 @@ MODULE MODE_MPPDB
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
RETURN
#else
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
!get the global size of KTAB
CALL MPI_ALLREDUCE(SIZE(KTAB), IGLBSIZEPTAB, 1, MNHINT_MPI, MPI_SUM, MPPDB_INTER_COMM, IINFO_ll)
IF ( IGLBSIZEPTAB == 0 ) RETURN
@@ -511,7 +533,7 @@ MODULE MODE_MPPDB
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
RETURN
#else
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
!get the global size of OTAB
CALL MPI_ALLREDUCE(SIZE(OTAB), IGLBSIZEPTAB, 1, MNHINT_MPI, MPI_SUM, MPPDB_INTER_COMM, IINFO_ll)
IF ( IGLBSIZEPTAB == 0 ) RETURN
@@ -691,7 +713,7 @@ MODULE MODE_MPPDB
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
RETURN
#else
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
!
IF ( PRESENT(PPRECISION) ) THEN
ZPRECISION = PPRECISION
@@ -897,7 +919,7 @@ MODULE MODE_MPPDB
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
RETURN
#else
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
!get the global size of OTAB
CALL MPI_ALLREDUCE(SIZE(OTAB), IGLBSIZEOTAB, 1, MNHINT_MPI, MPI_SUM, MPPDB_INTER_COMM, IINFO_ll)
IF ( IGLBSIZEOTAB == 0 ) RETURN
@@ -1171,7 +1193,7 @@ MODULE MODE_MPPDB
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
RETURN
#else
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
!
IF ( PRESENT(PPRECISION) ) THEN
ZPRECISION = PPRECISION
@@ -1511,7 +1533,7 @@ MODULE MODE_MPPDB
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
RETURN
#else
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
!
IF ( PRESENT(PPRECISION) ) THEN
ZPRECISION = PPRECISION
@@ -1786,14 +1808,16 @@ MODULE MODE_MPPDB
SUBROUTINE MPPDB_CHECKLB(PLB,MESSAGE,PPRECISION,HLBTYPE,KRIM)
- USE MODD_IO, ONLY: GSMONOPROC, ISP, ISNPROC, L2D, LPACK
- USE MODD_MPIF, ONLY: MPI_STATUS_IGNORE
- USE MODD_PARAMETERS_ll, ONLY: JPHEXT
- USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD
- use modd_precision, only: MNHINT_MPI, MNHREAL_MPI
+ USE MODD_PARAMETERS_ll, ONLY : JPHEXT
+ USE MODD_VAR_ll , ONLY : NMNH_COMM_WORLD
+ USE MODD_IO , ONLY : ISP,ISNPROC,GSMONOPROC,LPACK,L2D
+ USE MODD_MPIF , ONLY : MPI_STATUS_IGNORE, MPI_SUM
USE MODE_DISTRIB_LB
-
+ USE MODE_TOOLS_ll, ONLY : GET_GLOBALDIMS_ll
+ USE MODE_MODELN_HANDLER, ONLY : GET_CURRENT_MODEL_INDEX
+ use modd_precision, only: MNHINT_MPI, MNHREAL_MPI
+
IMPLICIT NONE
REAL, DIMENSION(:,:,:) , TARGET :: PLB
@@ -1822,13 +1846,18 @@ MODULE MODE_MPPDB
INTEGER :: IIU_SON_ll,IJU_SON_ll,IKU_SON_ll
INTEGER :: IIB_SON_ll,IIE_SON_ll,IJB_SON_ll,IJE_SON_ll
INTEGER :: IHEXT_SON_ll , IDIFF_HEXT , IRIM_ll , IRIM_SON_ll
+ INTEGER :: IMI , IGLBSIZEPTAB
+
#ifdef MNH_SP4
!pas de mpi_spawn sur IBM-SP ni MPI_ARGV_NULL etc ...
RETURN
#else
- IF ( .NOT. MPPDB_INITIALIZED ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
+ !get the global size of PLB
+ CALL MPI_ALLREDUCE(SIZE(PLB), IGLBSIZEPTAB, 1,MNHINT_MPI, MPI_SUM, MPPDB_INTER_COMM, IINFO_ll)
+ IF ( IGLBSIZEPTAB == 0 ) RETURN
!
- CALL MPPDB_BARRIER()
+ CALL MPPDB_BARRIER()
!
IF(MPPDB_FATHER_WORLD) THEN
!
@@ -1838,26 +1867,30 @@ MODULE MODE_MPPDB
IIU_ll = IIMAX_ll+2*JPHEXT
IJU_ll = IJMAX_ll+2*JPHEXT
IKU_ll = SIZE(PLB,3)
- IRIM_ll = (KRIM+JPHEXT)*2
-
- IF (HLBTYPE == 'LBX' .OR. HLBTYPE == 'LBXU') THEN
- IIU_ll = IRIM_ll
- ELSE
- IJU_ll = IRIM_ll
+ IRIM_ll = MAX(1,KRIM)
+
+ IF (HLBTYPE == 'LBX' ) THEN
+ IIU_ll = JPHEXT*2
+ ELSE IF ( HLBTYPE == 'LBXU') THEN
+ IIU_ll = (IRIM_ll+JPHEXT)*2
+ ELSE IF ( HLBTYPE == 'LBY') THEN
+ IJU_ll = JPHEXT*2
+ ELSE IF ( HLBTYPE == 'LBYV') THEN
+ IJU_ll = (IRIM_ll+JPHEXT)*2
END IF
IF (MPPDB_IRANK_WORLD.EQ.0) THEN
! I/O proc case
ALLOCATE(Z3D(IIU_ll,IJU_ll,SIZE(PLB,3)))
DO JI = 1,ISNPROC
- CALL GET_DISTRIB_LB(HLBTYPE,JI,'FM','WRITE',KRIM,IIB,IIE,IJB,IJE)
+ CALL GET_DISTRIB_LB(HLBTYPE,JI,'FM','WRITE',IRIM_ll,IIB,IIE,IJB,IJE)
IF (IIB /= 0) THEN
TX3DP=>Z3D(IIB:IIE,IJB:IJE,:)
IF (ISP /= JI) THEN
CALL MPI_RECV(TX3DP,SIZE(TX3DP),MNHREAL_MPI,JI-1 &
,99,NMNH_COMM_WORLD,MPI_STATUS_IGNORE,IINFO_ll)
ELSE
- CALL GET_DISTRIB_LB(HLBTYPE,JI,'LOC','WRITE',KRIM,IIB,IIE,IJB,IJE)
+ CALL GET_DISTRIB_LB(HLBTYPE,JI,'LOC','WRITE',IRIM_ll,IIB,IIE,IJB,IJE)
TX3DP = PLB(IIB:IIE,IJB:IJE,:)
END IF
END IF
@@ -1867,7 +1900,7 @@ MODULE MODE_MPPDB
ELSE
! Other processors
- CALL GET_DISTRIB_LB(HLBTYPE,ISP,'LOC','WRITE',KRIM,IIB,IIE,IJB,IJE)
+ CALL GET_DISTRIB_LB(HLBTYPE,ISP,'LOC','WRITE',IRIM_ll,IIB,IIE,IJB,IJE)
IF (IIB /= 0) THEN
TX3DP=>PLB(IIB:IIE,IJB:IJE,:)
CALL MPI_BSEND(TX3DP,SIZE(TX3DP),MNHREAL_MPI,0,99,NMNH_COMM_WORLD,IINFO_ll)
@@ -1891,13 +1924,17 @@ MODULE MODE_MPPDB
IIU_SON_ll = IIMAX_ll+2*IHEXT_SON_ll
IJU_SON_ll = IJMAX_ll+2*IHEXT_SON_ll
IKU_SON_ll = SIZE(PLB,3)
- IRIM_SON_ll = (KRIM+IHEXT_SON_ll)*2
+ IRIM_SON_ll = MAX(1,KRIM)
!
- IF (HLBTYPE == 'LBX' .OR. HLBTYPE == 'LBXU') THEN
- IIU_SON_ll = IRIM_SON_ll
- ELSE
- IJU_SON_ll = IRIM_SON_ll
- END IF
+ IF (HLBTYPE == 'LBX' ) THEN
+ IIU_SON_ll = IHEXT_SON_ll*2
+ ELSE IF ( HLBTYPE == 'LBXU') THEN
+ IIU_SON_ll = (IRIM_SON_ll+IHEXT_SON_ll)*2
+ ELSE IF ( HLBTYPE == 'LBY') THEN
+ IJU_SON_ll = IHEXT_SON_ll*2
+ ELSE IF ( HLBTYPE == 'LBYV') THEN
+ IJU_SON_ll = (IRIM_SON_ll+IHEXT_SON_ll)*2
+ END IF
!
ALLOCATE(TAB_SON_ll(IIU_SON_ll,IJU_SON_ll,IKU_SON_ll))
!
@@ -1929,12 +1966,13 @@ MODULE MODE_MPPDB
IJB_SON_ll-IDIFF_HEXT:IJE_SON_ll+IDIFF_HEXT,1:IKU_SON_ll) ) )
IF ( MAX_VAL .EQ. 0.0 ) MAX_VAL = 1.0
!
+ IMI = GET_CURRENT_MODEL_INDEX()
MAX_DIFF=MAXVAL(Z3D(IIB_ll-IDIFF_HEXT:IIE_ll+IDIFF_HEXT,IJB_ll-IDIFF_HEXT:IJE_ll+IDIFF_HEXT,1:IKU_ll)/MAX_VAL)
!
IF (MAX_DIFF > PPRECISION ) THEN
- print*," MPPDB_CHECKLB :: KO MPPDB_CHECKLB =", MESSAGE ," Error=",MAX_DIFF , MAX_VAL
+ print*," MPPDB_CHECKLB :: KO MPPDB_CHECKLB =", MESSAGE ," Error=",MAX_DIFF , MAX_VAL, IMI
ELSE
- print*," MPPDB_CHECKLB :: OK MPPDB_CHECKLB =", MESSAGE ," Error=",MAX_DIFF , MAX_VAL
+ print*," MPPDB_CHECKLB :: OK MPPDB_CHECKLB =", MESSAGE ," Error=",MAX_DIFF , MAX_VAL, IMI
END IF
flush(unit=OUTPUT_UNIT)
!
@@ -1994,7 +2032,7 @@ MODULE MODE_MPPDB
INTEGER :: KSIZE
INTEGER :: KSIZE_FULL
!
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
!
! IF ( SIZE(PTAB) == 0 ) THEN
! ALLOCATE(ZFIELD2D(0,0))
@@ -2098,7 +2136,7 @@ MODULE MODE_MPPDB
INTEGER :: IINFO_ll
INTEGER :: IKSIZE_ll
!
- IF ( ( .NOT. MPPDB_INITIALIZED ) ) RETURN
+ IF ( ( .NOT. MPPDB_INITIALIZED ) .OR. ( .NOT. MPPDB_ACTIVED ) ) RETURN
CALL MPI_ALLREDUCE(SIZE(PTAB), IGLBSIZEPTAB, 1,MNHINT_MPI, MPI_SUM, MPPDB_INTRA_COMM, IINFO_ll)
IF ( IGLBSIZEPTAB == 0 ) RETURN
CALL MPI_ALLREDUCE(SIZE(PTAB,2),IKSIZE_ll, 1, MNHINT_MPI, MPI_MAX, MPPDB_INTRA_COMM, IINFO_ll)
diff --git a/src/MNH/ares.f b/src/MNH/ares.fx90
similarity index 98%
rename from src/MNH/ares.f
rename to src/MNH/ares.fx90
index 2a300286f4320976d6edfc3de07d6ab24a574337..88f5e5d29428bce1433119715c84ce0bc7ff8e58 100644
--- a/src/MNH/ares.f
+++ b/src/MNH/ares.fx90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1987-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1987-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -833,11 +833,8 @@ C....................................................................
C...........PARAMETERS and their descriptions:
- INTEGER NCAT ! number of cations
- PARAMETER ( NCAT = 2 )
-
- INTEGER NAN ! number of anions
- PARAMETER ( NAN = 3 )
+ INTEGER, PARAMETER :: NCAT = 2 ! number of cations
+ INTEGER, PARAMETER :: NAN = 3 ! number of anions
C...........ARGUMENTS and their descriptions
@@ -1374,16 +1371,32 @@ c and the excess ammonium forms ammonum nitrate
end if
c
return
- end
+ end subroutine awater
c23456789012345678901234567890123456789012345678901234567890123456789012
+ MODULE MODI_poly4
+ INTERFACE
+ function poly4(A,X)
+ real A(4), X
+ end function poly4
+ END INTERFACE
+ END MODULE MODI_poly4
+
function poly4(A,X)
real poly4
real A(4), X
poly4 = A(1) + X * ( A(2) + X * ( A(3) + X * ( A(4) )))
- return
- end
+ return
+ end function poly4
+
+ MODULE MODI_poly6
+ INTERFACE
+ function poly6(A,X)
+ real A(6), X
+ end function poly6
+ END INTERFACE
+ END MODULE MODI_poly6
function poly6(A,X)
real poly6
@@ -1391,5 +1404,5 @@ c23456789012345678901234567890123456789012345678901234567890123456789012
poly6 = A(1) + X * ( A(2) + X * ( A(3) + X * ( A(4) +
& X * ( A(5) + X * (A(6) )))))
return
- end ! awater
+ end function poly6
c //////////////////////////////////////////////////////////////////
diff --git a/src/MNH/boundaries.f90 b/src/MNH/boundaries.f90
index 046a28af303cc6401540b799ad74289b7794fb3d..1b73df6c2d508cc3eb67c3af1c1ea3f1314660b8 100644
--- a/src/MNH/boundaries.f90
+++ b/src/MNH/boundaries.f90
@@ -199,6 +199,7 @@ USE MODD_PARAMETERS
USE MODD_PARAM_LIMA, ONLY : NMOD_CCN, NMOD_IFN, LBOUND, LWARM, LCOLD
USE MODD_PARAM_n, ONLY : CELEC,CCLOUD
USE MODD_PASPOL, ONLY : LPASPOL
+USE MODD_PRECISION, ONLY: MNHREAL32
USE MODD_REF_n
USE MODD_SALT, ONLY : LSALT
@@ -420,6 +421,24 @@ ELSE
!
END IF
!
+! ============================================================
+!
+! Reproductibility for RSTART -> truncate ZLB to real(knd=4) to have reproductible result
+!
+ZLBXVT(:,:,:) = real(ZLBXVT(:,:,:),kind=MNHREAL32)
+ZLBXWT(:,:,:) = real(ZLBXWT(:,:,:),kind=MNHREAL32)
+ZLBXTHT(:,:,:) = real(ZLBXTHT(:,:,:),kind=MNHREAL32)
+IF ( SIZE(PTKET,1) /= 0 ) THEN
+ ZLBXTKET(:,:,:) = real(ZLBXTKET(:,:,:),kind=MNHREAL32)
+END IF
+IF ( KRR > 0) THEN
+ ZLBXRT(:,:,:,:) = real(ZLBXRT(:,:,:,:),kind=MNHREAL32)
+END IF
+IF ( KSV > 0) THEN
+ ZLBXSVT(:,:,:,:) = real(ZLBXSVT(:,:,:,:),kind=MNHREAL32)
+END IF
+! ============================================================
+!
IF ( SIZE(PLBYTHS,1) /= 0 .AND. &
( HLBCY(1)=='OPEN' .OR. HLBCY(2)=='OPEN' )) THEN
ZLBYUT(:,:,:) = PLBYUM(:,:,:) + ZTSTEP * PLBYUS(:,:,:)
@@ -453,6 +472,24 @@ ELSE
END IF
!
!
+! ============================================================
+!
+! Reproductibility for RSTART -> truncate ZLB to real(knd=4) to have reproductible result
+!
+ZLBYUT(:,:,:) = real(ZLBYUT(:,:,:),kind=MNHREAL32)
+ZLBYWT(:,:,:) = real(ZLBYWT(:,:,:),kind=MNHREAL32)
+ZLBYTHT(:,:,:) = real(ZLBYTHT(:,:,:),kind=MNHREAL32)
+IF ( SIZE(PTKET,1) /= 0 ) THEN
+ ZLBYTKET(:,:,:) = real(ZLBYTKET(:,:,:),kind=MNHREAL32)
+END IF
+IF ( KRR > 0) THEN
+ ZLBYRT(:,:,:,:) = real(ZLBYRT(:,:,:,:),kind=MNHREAL32)
+END IF
+IF ( KSV > 0) THEN
+ ZLBYSVT(:,:,:,:) = real(ZLBYSVT(:,:,:,:),kind=MNHREAL32)
+END IF
+! ============================================================
+!
!-------------------------------------------------------------------------------
! PONDERATION COEFF for Non-Normal velocities and pot temperature
!
diff --git a/src/MNH/budget.f90 b/src/MNH/budget.f90
index 13c746c4344cd996220e6905df3e2020bbf2fab9..5e7dd0809c8f6296aca79452016cafe437950254 100644
--- a/src/MNH/budget.f90
+++ b/src/MNH/budget.f90
@@ -7,6 +7,7 @@
! P. Wautelet 28/01/2020: new subroutines: Budget_store_init, Budget_store_end and Budget_source_id_find in new module mode_budget
! P. Wautelet 17/08/2020: treat LES budgets correctly
! P. Wautelet 05/03/2021: measure cpu_time for budgets
+! J.Escobar : 06/10/2021 :for bit reproductiblity use MPPDB_CHECK if LCHECK=T
!-----------------------------------------------------------------
!#################
@@ -36,6 +37,8 @@ contains
subroutine Budget_store_init( tpbudget, hsource, pvars )
use modd_les, only: lles_call
+ USE MODE_MPPDB
+ USE MODD_CONF, ONLY : LCHECK
type(tbudgetdata), intent(inout) :: tpbudget ! Budget datastructure
character(len=*), intent(in) :: hsource ! Name of the source term
@@ -43,7 +46,15 @@ subroutine Budget_store_init( tpbudget, hsource, pvars )
integer :: iid ! Reference number of the current source term
- call Print_msg( NVERB_DEBUG, 'BUD', 'Budget_store_init', trim( tpbudget%cname )//':'//trim( hsource ) )
+ character(len=:),allocatable :: hbudget
+
+ hbudget = trim( tpbudget%cname )//':'//trim( hsource )
+
+ IF (LCHECK) THEN
+ CALL MPPDB_CHECK3D(PVARS,'BUD_INI::'//hbudget,PRECISION)
+ END IF
+
+ call Print_msg( NVERB_DEBUG, 'BUD', 'Budget_store_init', hbudget )
if ( lles_call ) then
call Second_mnh( ztime1 )
@@ -112,6 +123,8 @@ subroutine Budget_store_init( tpbudget, hsource, pvars )
subroutine Budget_store_end( tpbudget, hsource, pvars )
use modd_les, only: lles_call
+ USE MODE_MPPDB
+ USE MODD_CONF, ONLY : LCHECK
use modi_les_budget, only: Les_budget
@@ -123,7 +136,14 @@ subroutine Budget_store_end( tpbudget, hsource, pvars )
integer :: igroup ! Number of the group where to store the source term
real, dimension(:,:,:), allocatable :: zvars_add
- call Print_msg( NVERB_DEBUG, 'BUD', 'Budget_store_end', trim( tpbudget%cname )//':'//trim( hsource ) )
+ character(len=:),allocatable :: hbudget
+
+ hbudget = trim( tpbudget%cname )//':'//trim( hsource )
+
+ IF (LCHECK) THEN
+ CALL MPPDB_CHECK3D(PVARS,'BUD_END::'//hbudget,PRECISION)
+ END IF
+ call Print_msg( NVERB_DEBUG, 'BUD', 'Budget_store_end', hbudget )
if ( lles_call ) then
if ( hsource /= tpbudget%clessource ) &
diff --git a/src/MNH/compute_bl89_ml.f90 b/src/MNH/compute_bl89_ml.f90
index 20c9a078db5e550dbf3bb03cd2fcc6b13ddc89a5..b2df24bd9b4d00d0e37550e347af3dcef9fc3abf 100644
--- a/src/MNH/compute_bl89_ml.f90
+++ b/src/MNH/compute_bl89_ml.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2006-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2006-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -57,6 +57,7 @@ END MODULE MODI_COMPUTE_BL89_ML
!! R.Honnert Oct 2016 : Update with AROME
!! Q.Rodier 01/2019 : support RM17 mixing length as in bl89.f90
! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+! P. Wautelet 17/12/2021: bugfix: KK instead of JKK
!!
!-------------------------------------------------------------------------------
!
@@ -162,9 +163,9 @@ IF (OUPORDN.EQV..TRUE.) THEN
! Lenght travelled by parcel to nullify energy
ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
( ZHLVPT(J1D,KK) - ZVPT_DEP(J1D) ) &
- - XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
+ - XRM17*PSHEAR(J1D,KK)*sqrt(abs(PTKEM_DEP(J1D))) &
+ SQRT (ABS( &
- (XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) + &
+ (XRM17*PSHEAR(J1D,KK)*sqrt(abs(PTKEM_DEP(J1D))) + &
PG_O_THVREF(J1D) * (ZHLVPT(J1D,KK) - ZVPT_DEP(J1D)) )**2 &
+ 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
* ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )) ) / &
diff --git a/src/MNH/ini_lb.f90 b/src/MNH/ini_lb.f90
index b4d44b50adff30ff2ee57769d3b6be3bbbceba63..79eac4a6e755e814bc908f200675fb50b5707275 100644
--- a/src/MNH/ini_lb.f90
+++ b/src/MNH/ini_lb.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1998-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1998-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -689,7 +689,14 @@ IF (CCLOUD=='LIMA' ) THEN
CASE ('READ')
WRITE(INDICE,'(I2.2)')(JSV - NSV_LIMA_CCN_FREE + 1)
IF ( KSIZELBXSV_ll /= 0 ) THEN
- TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3))//INDICE)
+ IF ( TPINIFILE%NMNHVERSION(1) < 5 &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) < 5 ) &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) == 5 &
+ .AND. TPINIFILE%NMNHVERSION(3) < 1 ) ) THEN
+ TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3))//INDICE)
+ ELSE
+ TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3)))//INDICE
+ END IF
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBX'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'2_Y_Z_','LBXSVM',JSV
@@ -708,7 +715,14 @@ IF (CCLOUD=='LIMA' ) THEN
END IF
!
IF (KSIZELBYSV_ll /= 0 ) THEN
- TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3))//INDICE)
+ IF ( TPINIFILE%NMNHVERSION(1) < 5 &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) < 5 ) &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) == 5 &
+ .AND. TPINIFILE%NMNHVERSION(3) < 1 ) ) THEN
+ TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3))//INDICE)
+ ELSE
+ TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3)))//INDICE
+ END IF
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBY'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'X_2_Z_','LBYSVM',JSV
@@ -746,7 +760,14 @@ IF (CCLOUD=='LIMA' ) THEN
CASE ('READ')
WRITE(INDICE,'(I2.2)')(JSV - NSV_LIMA_IFN_FREE + 1)
IF ( KSIZELBXSV_ll /= 0 ) THEN
- TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2))//INDICE)
+ IF ( TPINIFILE%NMNHVERSION(1) < 5 &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) < 5 ) &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) == 5 &
+ .AND. TPINIFILE%NMNHVERSION(3) < 1 ) ) THEN
+ TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2))//INDICE)
+ ELSE
+ TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2)))//INDICE
+ END IF
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBX'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'2_Y_Z_','LBXSVM',JSV
@@ -765,7 +786,14 @@ IF (CCLOUD=='LIMA' ) THEN
END IF
!
IF (KSIZELBYSV_ll /= 0 ) THEN
- TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2))//INDICE)
+ IF ( TPINIFILE%NMNHVERSION(1) < 5 &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) < 5 ) &
+ .OR. ( TPINIFILE%NMNHVERSION(1) == 5 .AND. TPINIFILE%NMNHVERSION(2) == 5 &
+ .AND. TPINIFILE%NMNHVERSION(3) < 1 ) ) THEN
+ TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2))//INDICE)
+ ELSE
+ TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2)))//INDICE
+ END IF
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBY'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'X_2_Z_','LBYSVM',JSV
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 3cf3c9525806ce866d1510a45a83de703a3ee133..d46bb33959f348f1f85efc006769d4cdaf8259b3 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -1802,7 +1802,7 @@ gles = lles_mean .or. lles_resolved .or. lles_subgrid .or. lles_updraft &
.or. lles_downdraft .or. lles_spectra
!Called if budgets are enabled via NAM_BUDGET
!or if LES budgets are enabled via NAM_LES (condition on kmi==1 to call it max once)
-if ( ( cbutype /= "NONE" .and. nbumod == kmi ) .or. ( gles .and. kmi == 1 ) ) THEN
+if ( ( cbutype /= "NONE" .and. nbumod == kmi ) .or. ( gles .and. kmi == 1 ) .or. LCHECK ) THEN
call Budget_preallocate()
end if
@@ -1917,7 +1917,7 @@ CALL READ_FIELD(KMI,TPINIFILE,IIU,IJU,IKU, &
NSIZELBXR_ll,NSIZELBYR_ll,NSIZELBXSV_ll,NSIZELBYSV_ll, &
XUM,XVM,XWM,XDUM,XDVM,XDWM, &
XUT,XVT,XWT,XTHT,XPABST,XTKET,XRTKEMS, &
- XRT,XSVT,XZWS,XCIT,XDRYMASST, &
+ XRT,XSVT,XZWS,XCIT,XDRYMASST,XDRYMASSS, &
XSIGS,XSRCT,XCLDFR,XBL_DEPTH,XSBL_DEPTH,XWTHVMF,XPHC,XPHR, &
XLSUM,XLSVM,XLSWM,XLSTHM,XLSRVM,XLSZWSM, &
XLBXUM,XLBXVM,XLBXWM,XLBXTHM,XLBXTKEM, &
@@ -1972,8 +1972,8 @@ CALL INI_LES_n
!* 11. INITIALIZE THE SOURCE OF TOTAL DRY MASS Md
! ------------------------------------------
!
-IF((KMI==1).AND.LSTEADYLS) THEN
- XDRYMASSS = 0.
+IF((KMI==1).AND.LSTEADYLS .AND. (CCONF=='START') ) THEN
+ XDRYMASSS = 0.
END IF
!
!-------------------------------------------------------------------------------
@@ -2167,20 +2167,22 @@ IF ( KMI > 1) THEN
DPTR_XLBYRM=>XLBYRM
DPTR_XLBXSVM=>XLBXSVM
DPTR_XLBYSVM=>XLBYSVM
- CALL INI_ONE_WAY_n(NDAD(KMI),KMI, &
- DPTR_XBMX1,DPTR_XBMX2,DPTR_XBMX3,DPTR_XBMX4,DPTR_XBMY1,DPTR_XBMY2,DPTR_XBMY3,DPTR_XBMY4, &
- DPTR_XBFX1,DPTR_XBFX2,DPTR_XBFX3,DPTR_XBFX4,DPTR_XBFY1,DPTR_XBFY2,DPTR_XBFY3,DPTR_XBFY4, &
- NDXRATIO_ALL(KMI),NDYRATIO_ALL(KMI), &
- DPTR_CLBCX,DPTR_CLBCY,NRIMX,NRIMY, &
- DPTR_NKLIN_LBXU,DPTR_XCOEFLIN_LBXU,DPTR_NKLIN_LBYU,DPTR_XCOEFLIN_LBYU, &
- DPTR_NKLIN_LBXV,DPTR_XCOEFLIN_LBXV,DPTR_NKLIN_LBYV,DPTR_XCOEFLIN_LBYV, &
- DPTR_NKLIN_LBXW,DPTR_XCOEFLIN_LBXW,DPTR_NKLIN_LBYW,DPTR_XCOEFLIN_LBYW, &
- DPTR_NKLIN_LBXM,DPTR_XCOEFLIN_LBXM,DPTR_NKLIN_LBYM,DPTR_XCOEFLIN_LBYM, &
- CCLOUD, LUSECHAQ, LUSECHIC, &
- DPTR_XLBXUM,DPTR_XLBYUM,DPTR_XLBXVM,DPTR_XLBYVM,DPTR_XLBXWM,DPTR_XLBYWM, &
- DPTR_XLBXTHM,DPTR_XLBYTHM, &
- DPTR_XLBXTKEM,DPTR_XLBYTKEM, &
- DPTR_XLBXRM,DPTR_XLBYRM,DPTR_XLBXSVM,DPTR_XLBYSVM )
+ IF (CCONF=='START') THEN
+ CALL INI_ONE_WAY_n(NDAD(KMI),KMI, &
+ DPTR_XBMX1,DPTR_XBMX2,DPTR_XBMX3,DPTR_XBMX4,DPTR_XBMY1,DPTR_XBMY2,DPTR_XBMY3,DPTR_XBMY4, &
+ DPTR_XBFX1,DPTR_XBFX2,DPTR_XBFX3,DPTR_XBFX4,DPTR_XBFY1,DPTR_XBFY2,DPTR_XBFY3,DPTR_XBFY4, &
+ NDXRATIO_ALL(KMI),NDYRATIO_ALL(KMI), &
+ DPTR_CLBCX,DPTR_CLBCY,NRIMX,NRIMY, &
+ DPTR_NKLIN_LBXU,DPTR_XCOEFLIN_LBXU,DPTR_NKLIN_LBYU,DPTR_XCOEFLIN_LBYU, &
+ DPTR_NKLIN_LBXV,DPTR_XCOEFLIN_LBXV,DPTR_NKLIN_LBYV,DPTR_XCOEFLIN_LBYV, &
+ DPTR_NKLIN_LBXW,DPTR_XCOEFLIN_LBXW,DPTR_NKLIN_LBYW,DPTR_XCOEFLIN_LBYW, &
+ DPTR_NKLIN_LBXM,DPTR_XCOEFLIN_LBXM,DPTR_NKLIN_LBYM,DPTR_XCOEFLIN_LBYM, &
+ CCLOUD, LUSECHAQ, LUSECHIC, &
+ DPTR_XLBXUM,DPTR_XLBYUM,DPTR_XLBXVM,DPTR_XLBYVM,DPTR_XLBXWM,DPTR_XLBYWM, &
+ DPTR_XLBXTHM,DPTR_XLBYTHM, &
+ DPTR_XLBXTKEM,DPTR_XLBYTKEM, &
+ DPTR_XLBXRM,DPTR_XLBYRM,DPTR_XLBXSVM,DPTR_XLBYSVM )
+ ENDIF
END IF
!
!
diff --git a/src/MNH/isocom.f b/src/MNH/isocom.f
index 451c1b485c342af59e48c962566e12ed682dc631..42481f83a3f489686c70fa8a696f9ffe8966ff04 100644
--- a/src/MNH/isocom.f
+++ b/src/MNH/isocom.f
@@ -3,6 +3,13 @@ CMNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
CMNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
CMNH_LIC for details. version 1.
C=======================================================================
+C Modifications:
+C P. Wautelet 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
+C P. Wautelet 22/01/2019: replace obsolete SNGL intrinsics by REAL intrinsics
+C P. Wautelet 19/04/2019: use kind(0.0d0) instead of kind=8
+C P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+C P. Wautelet 17/12/2021: add missing definitions of parameter ONE (in POLY3 and POLY3B)
+C=======================================================================
C
C *** ISORROPIA CODE
C *** SUBROUTINE ISOROPIA
@@ -123,11 +130,6 @@ C
C *** COPYRIGHT 1996-2000, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY
C *** WRITTEN BY ATHANASIOS NENES
C
-C Modifications:
-C P. Wautelet 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
-C P. Wautelet 22/01/2019: replace obsolete SNGL intrinsics by REAL intrinsics
-C P. Wautelet 19/04/2019: use kind(0.0d0) instead of kind=8
-C P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
C=======================================================================
C
SUBROUTINE ISOROPIA (WI, RHI, TEMPI, CNTRL,
@@ -3882,7 +3884,7 @@ C
SUBROUTINE POLY3B (A1, A2, A3, RTLW, RTHI, ROOT, ISLV)
C
IMPLICIT REAL(kind(0.0d0)) (A-H, O-Z)
- PARAMETER (ZERO=0.D0, EPS=1D-15, MAXIT=100, NDIV=5)
+ PARAMETER (ZERO=0.D0, ONE=1.D0, EPS=1D-15, MAXIT=100, NDIV=5)
C
FUNC(X) = X**3.d0 + A1*X**2.0 + A2*X + A3
C
diff --git a/src/MNH/mesonh.f90 b/src/MNH/mesonh.f90
index 84017a6b994ac7be3f3681429bfd3a61771b67ed..ac5f3cdef04db6c7ffa837bbbd1f849b8974e9a8 100644
--- a/src/MNH/mesonh.f90
+++ b/src/MNH/mesonh.f90
@@ -134,6 +134,7 @@ INTEGER :: IINFO_ll ! return code of // routines
! Switch to model 1 variables
#ifndef CPLOASIS
CALL MPPDB_INIT()
+CALL MPPDB_STOP_DEBUG()
#endif
!
CALL GOTO_MODEL(1,ONOFIELDLIST=.TRUE.)
diff --git a/src/MNH/mode_RBK90_Integrator.f90 b/src/MNH/mode_RBK90_Integrator.f90
index ae11720e3cd8fdc1ce7ab60d602dd20a5e4ce0d5..536726def5d7ba4cca6db9fab1f67c5e867db171 100644
--- a/src/MNH/mode_RBK90_Integrator.f90
+++ b/src/MNH/mode_RBK90_Integrator.f90
@@ -1,10 +1,11 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! Modifications:
-! Philippe 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
+! P. Wautelet 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
+! P. Wautelet 17/12/2021: remove unused time variables (and remove use of non initialised time values)
!-----------------------------------------------------------------
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
@@ -1316,9 +1317,6 @@ USE MODI_CH_FCN
REAL :: T, Y(N)
!~~~> Output variables
REAL :: Ydot(N)
-!~~~> Local variables
- REAL :: Told
- REAL :: TIME
!
INTEGER, INTENT(IN) :: KMI ! model number
INTEGER, INTENT(IN) :: KVECNPT
@@ -1330,8 +1328,6 @@ INTEGER, INTENT(IN) :: KVECNPT
INTEGER :: ISPECIES ! Ancillary variable
- Told = TIME
- TIME = T
!JPP CALL Update_SUN()
!JPP CALL Update_RCONST()
!JPP CALL Fun( Y, FIX, RCONST, Ydot )
@@ -1352,7 +1348,6 @@ INTEGER, INTENT(IN) :: KVECNPT
Ydot(JLSHIFT+1:JLSHIFT+ISPECIES) = ZFCN(JL,1:ISPECIES)
JLSHIFT = JLSHIFT + ISPECIES
END DO
- TIME = Told
END SUBROUTINE FunTemplate
@@ -1389,8 +1384,6 @@ INTEGER, INTENT(IN) :: KMI ! model number
INTEGER, INTENT(IN) :: KVECNPT
!
!~~~> Local variables
- REAL :: Told
- REAL :: TIME
!#ifdef FULL_ALGEBRA
!## INTEGER :: i, j
!#endif
@@ -1399,8 +1392,6 @@ INTEGER, INTENT(IN) :: KVECNPT
INTEGER :: JL, JLL, JLSHIFT ! Loop indexes
INTEGER :: ISPECIES ! Ancillary variable
- Told = TIME
- TIME = T
!JPP CALL Update_SUN()
!JPP CALL Update_RCONST()
!#ifdef FULL_ALGEBRA
@@ -1441,7 +1432,6 @@ INTEGER, INTENT(IN) :: KVECNPT
JLSHIFT = JLSHIFT + NSPARSEDIM_NAQ
END IF
END DO
- TIME = Told
END SUBROUTINE JacTemplate
diff --git a/src/MNH/mode_RBK90_linearalgebra.f90 b/src/MNH/mode_RBK90_linearalgebra.f90
index 8fd23aa139e89b8c906915ccda98c36e80e102e0..de787d50d3e04f7b43d59230b4f92f7d525bf730 100644
--- a/src/MNH/mode_RBK90_linearalgebra.f90
+++ b/src/MNH/mode_RBK90_linearalgebra.f90
@@ -1,10 +1,11 @@
-!MNH_LIC Copyright 1994-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! Modifications:
! P. Wautelet 22/02/2019: DOUBLE COMPLEX -> COMPLEX(kind(0.0d0)) to respect Fortran standard
+! P. Wautelet 17/12/2021: convert arithmetic if to classic if (deleted from Fortran 2018 standard)
!-----------------------------------------------------------------
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
@@ -1261,46 +1262,49 @@ END SUBROUTINE KppSolveTR
WDOT = 0.0D0
IF (N .LE. 0) RETURN
- IF (incX .EQ. incY) IF (incX-1) 5,20,60
+ IF (incX .EQ. incY) THEN
!
! Code for unequal or nonpositive increments.
!
- 5 IX = 1
- IY = 1
- IF (incX .LT. 0) IX = (-N+1)*incX + 1
- IF (incY .LT. 0) IY = (-N+1)*incY + 1
- DO i = 1,N
- WDOT = WDOT + DX(IX)*DY(IY)
- IX = IX + incX
- IY = IY + incY
- END DO
- RETURN
+ IF ( incX < 1 ) THEN
+ IX = 1
+ IY = 1
+ IF (incX .LT. 0) IX = (-N+1)*incX + 1
+ IF (incY .LT. 0) IY = (-N+1)*incY + 1
+ DO i = 1,N
+ WDOT = WDOT + DX(IX)*DY(IY)
+ IX = IX + incX
+ IY = IY + incY
+ END DO
+ ELSE IF ( incX == 1 ) THEN
!
! Code for both increments equal to 1.
!
! Clean-up loop so remaining vector length is a multiple of 5.
!
- 20 M = MOD(N,5)
- IF (M .EQ. 0) GO TO 40
- DO i = 1,M
- WDOT = WDOT + DX(i)*DY(i)
- END DO
- IF (N .LT. 5) RETURN
- 40 MP1 = M + 1
- DO i = MP1,N,5
- WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) + &
- DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)
- END DO
- RETURN
+ M = MOD(N,5)
+ IF (M .EQ. 0) GO TO 40
+ DO i = 1,M
+ WDOT = WDOT + DX(i)*DY(i)
+ END DO
+ IF (N .LT. 5) RETURN
+ 40 MP1 = M + 1
+ DO i = MP1,N,5
+ WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) + &
+ DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)
+ END DO
+ ELSE
!
! Code for equal, positive, non-unit increments.
!
- 60 NS = N*incX
- DO i = 1,NS,incX
- WDOT = WDOT + DX(i)*DY(i)
- END DO
+ NS = N*incX
+ DO i = 1,NS,incX
+ WDOT = WDOT + DX(i)*DY(i)
+ END DO
+ END IF
+ END IF
- END FUNCTION WDOT
+ END FUNCTION WDOT
!--------------------------------------------------------------
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index 5c14265f0c3324c64bf283ff5032eb8b7e010f7d..146c49c400b15a8e2287a111f56ad299aca1c502 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -951,6 +951,14 @@ CALL SECOND_MNH2(ZTIME2)
!
XT_BOUND = XT_BOUND + ZTIME2 - ZTIME1
!
+!
+! For START/RESTART MPPDB_CHECK use
+!IF ( (IMI==1) .AND. (CCONF == "START") .AND. (KTCOUNT == 2) ) THEN
+! CALL MPPDB_START_DEBUG()
+!ENDIF
+!IF ( (IMI==1) .AND. (CCONF == "RESTA") .AND. (KTCOUNT == 1) ) THEN
+! CALL MPPDB_START_DEBUG()
+!ENDIF
!-------------------------------------------------------------------------------
!* initializes surface number
IF (CSURF=='EXTE') CALL GOTO_SURFEX(IMI)
@@ -979,6 +987,10 @@ IF ( nfile_backup_current < NBAK_NUMB ) THEN
TFILE_SURFEX => TZBAKFILE
CALL GOTO_SURFEX(IMI)
CALL WRITE_SURF_ATM_n(YSURF_CUR,'MESONH','ALL',.FALSE.)
+ IF ( KTCOUNT > 1) THEN
+ CALL DIAG_SURF_ATM_n(YSURF_CUR,'MESONH')
+ CALL WRITE_DIAG_SURF_ATM_n(YSURF_CUR,'MESONH','ALL')
+ END IF
NULLIFY(TFILE_SURFEX)
END IF
!
@@ -1433,17 +1445,6 @@ IF (CDCONV/='NONE') THEN
END IF
END IF
!
-IF ( nfile_backup_current > 0 .AND. nfile_backup_current <= NBAK_NUMB ) THEN
- IF ( KTCOUNT == TBACKUPN(nfile_backup_current)%NSTEP ) THEN
- IF (CSURF=='EXTE') THEN
- CALL GOTO_SURFEX(IMI)
- CALL DIAG_SURF_ATM_n(YSURF_CUR,'MESONH')
- TFILE_SURFEX => TZBAKFILE
- CALL WRITE_DIAG_SURF_ATM_n(YSURF_CUR,'MESONH','ALL')
- NULLIFY(TFILE_SURFEX)
- END IF
- END IF
-END IF
!
CALL SECOND_MNH2(ZTIME2)
!
@@ -1744,14 +1745,20 @@ ZTIME1 = ZTIME2
ZRUS=XRUS
ZRVS=XRVS
ZRWS=XRWS
-
+!
if ( .not. l1d ) then
if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'PRES', xrus(:, :, :) )
if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'PRES', xrvs(:, :, :) )
if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'PRES', xrws(:, :, :) )
end if
-
-CALL RAD_BOUND (CLBCX,CLBCY,CTURB,XCARPKMAX, &
+!
+CALL MPPDB_CHECK3DM("before RAD_BOUND : other var",PRECISION,XUT,XVT,XRHODJ,XTKET)
+CALL MPPDB_CHECKLB(XLBXUM,"modeln XLBXUM",PRECISION,'LBXU',NRIMX)
+CALL MPPDB_CHECKLB(XLBYVM,"modeln XLBYVM",PRECISION,'LBYV',NRIMY)
+CALL MPPDB_CHECKLB(XLBXUS,"modeln XLBXUS",PRECISION,'LBXU',NRIMX)
+CALL MPPDB_CHECKLB(XLBYVS,"modeln XLBYVS",PRECISION,'LBYV',NRIMY)
+!
+ CALL RAD_BOUND (CLBCX,CLBCY,CTURB,XCARPKMAX, &
XTSTEP, &
XDXHAT, XDYHAT, XZHAT, &
XUT, XVT, &
@@ -1797,6 +1804,7 @@ IF(.NOT. L1D) THEN
XRUS_PRES = XRUS - XRUS_PRES + ZRUS
XRVS_PRES = XRVS - XRVS_PRES + ZRVS
XRWS_PRES = XRWS - XRWS_PRES + ZRWS
+ CALL MPPDB_CHECK3DM("after pressurez:XRU/V/WS",PRECISION,XRUS,XRVS,XRWS)
!
END IF
!
diff --git a/src/MNH/nband_model.fx90 b/src/MNH/nband_model.fx90
index e0aface0e240a114c19048fbf5ddf99cb6c969f5..6005b838759cf188856ece7a479f72945cf033a2 100644
--- a/src/MNH/nband_model.fx90
+++ b/src/MNH/nband_model.fx90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1996-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1996-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -13,6 +13,7 @@
* named nband_model.f90 and compiled with -Fixed
* J.Escobar (1/12/2017) bug => intialized all ZV=0.0 in spectr
* P. Wautelet 21/11/2019: replace several CONTINUE (workaround of problems with gfortran OpenACC)
+* P. Wautelet 17/12/2021: comment ZBSUI variable (not used and was not initialized)
*
SUBROUTINE NBMVEC
I ( KIDIA ,KFDIA ,KLON,KLEV,KGL,KCABS,KNG1,KUABS
@@ -768,7 +769,7 @@ C
DO 255 JL=KIDIA,KFDIA
ZBLEV(JL,KLEV+1)=ZRES(JL)
ZBINT(JL,KLEV+1)=ZBINT(JL,KLEV+1)+ZRES(JL)
- ZBSUI(JL)=ZBSUI(JL)+ZBSUR(JL)
+C ZBSUI(JL)=ZBSUI(JL)+ZBSUR(JL)
255 CONTINUE
C IF (NIMP.EQ.0) THEN
C JL=KIDIA
diff --git a/src/MNH/nhoa_coupln.f90 b/src/MNH/nhoa_coupln.f90
index f8ea30639b021065de23652e151303524f0644f4..072dfbeb3ca7aa64c133a8d0b2855e4a3d1fdea8 100644
--- a/src/MNH/nhoa_coupln.f90
+++ b/src/MNH/nhoa_coupln.f90
@@ -49,6 +49,7 @@ SUBROUTINE NHOA_COUPL_n(KDAD,PTSTEP,KMI,KTCOUNT,KKU)
!! JL Redelsperger 03/2021 Version 0
!! MODIFICATIONS
!! -------------
+! P. Wautelet 17/12/2021: bugfix: set IIU and IJU values
!!-----------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -98,17 +99,20 @@ CHARACTER(LEN=4) :: ZINIT_TYPE
!---Coupled OA MesoNH----------------------------------------------------------------------------
!* 0. INITIALISATION
! --------------
+IIU = SIZE(XRHODJ,1)
+IJU = SIZE(XRHODJ,2)
+
! allocate flux local array
-ALLOCATE(ZCOUPTFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPUFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPVFL(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPTFL(IIU,IJU))
+ALLOCATE(ZCOUPUFL(IIU,IJU))
+ALLOCATE(ZCOUPVFL(IIU,IJU))
! allocate sfc variable local array
-ALLOCATE(ZCOUPUA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPVA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPTA(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPUO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPVO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
-ALLOCATE(ZCOUPTO(SIZE(XRHODJ,1),SIZE(XRHODJ,2)))
+ALLOCATE(ZCOUPUA(IIU,IJU))
+ALLOCATE(ZCOUPVA(IIU,IJU))
+ALLOCATE(ZCOUPTA(IIU,IJU))
+ALLOCATE(ZCOUPUO(IIU,IJU))
+ALLOCATE(ZCOUPVO(IIU,IJU))
+ALLOCATE(ZCOUPTO(IIU,IJU))
! values in ocean sfc
IKE=KKU-JPVEXT
ZCOUPUO(:,:)= XUT(:,:,IKE)
diff --git a/src/MNH/num_diff.f90 b/src/MNH/num_diff.f90
index 65c221c739d179b1116501ac1ce8fe9700be733a..688304a1a2bce2db56e2bc8aadbe9774692ae689 100644
--- a/src/MNH/num_diff.f90
+++ b/src/MNH/num_diff.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -227,6 +227,7 @@ use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, lbudget_th, lbudg
tbudgets
USE MODD_CONF
USE MODD_PARAMETERS
+use modd_precision, only: MNHREAL32
use mode_budget, only: Budget_store_init, Budget_store_end
USE MODE_ll
@@ -1120,18 +1121,18 @@ CASE ('CYCL') ! In that case one must have HLBCX(1) == HLBCX(2)
TPHALO2%WEST(:,:) + PFIELDM(IW+1,:,:) &
-4.*( PFIELDM(IW-2,:,:) + PFIELDM(IW,:,:) ) &
+6.* PFIELDM(IW-1,:,:) &
- - TPHALO2LS%WEST(:,:) - PLSFIELD(IW+1,:,:) &
- +4.*( PLSFIELD(IW-2,:,:) + PLSFIELD(IW,:,:) ) &
- -6.* PLSFIELD(IW-1,:,:) )
+ - real(TPHALO2LS%WEST(:,:),kind=MNHREAL32) - real(PLSFIELD(IW+1,:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(IW-2,:,:),kind=MNHREAL32) + real(PLSFIELD(IW,:,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(IW-1,:,:),kind=MNHREAL32) )
!
PRFIELDS(IE+1,:,:) = PRFIELDS(IE+1,:,:) - PRHODJ(IE+1,:,:) * &
PDK4*( &
PFIELDM(IE-1,:,:) + TPHALO2%EAST(:,:) &
-4.*( PFIELDM(IE,:,:) + PFIELDM(IE+2,:,:) ) &
+6.* PFIELDM(IE+1,:,:) &
- - PLSFIELD(IE-1,:,:) - TPHALO2LS%EAST(:,:) &
- +4.*( PLSFIELD(IE,:,:) + PLSFIELD(IE+2,:,:) ) &
- -6.* PLSFIELD(IE+1,:,:) )
+ - real(PLSFIELD(IE-1,:,:),kind=MNHREAL32) - real(TPHALO2LS%EAST(:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(IE,:,:),kind=MNHREAL32) + real(PLSFIELD(IE+2,:,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(IE+1,:,:),kind=MNHREAL32) )
!
!!$ ENDIF
!
@@ -1142,10 +1143,10 @@ CASE ('CYCL') ! In that case one must have HLBCX(1) == HLBCX(2)
PDK4*( &
PFIELDM(IW-2:IE-2,:,:) + PFIELDM(IW+2:IE+2,:,:) &
-4.*( PFIELDM(IW-1:IE-1,:,:) + PFIELDM(IW+1:IE+1,:,:) ) &
- +6.* PFIELDM(IW:IE,:,:) &
- - PLSFIELD(IW-2:IE-2,:,:) - PLSFIELD(IW+2:IE+2,:,:) &
- +4.*( PLSFIELD(IW-1:IE-1,:,:) + PLSFIELD(IW+1:IE+1,:,:) ) &
- -6.* PLSFIELD(IW:IE,:,:))
+ +6.* PFIELDM(IW:IE,:,:) &
+ - real(PLSFIELD(IW-2:IE-2,:,:),kind=MNHREAL32) - real(PLSFIELD(IW+2:IE+2,:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(IW-1:IE-1,:,:),kind=MNHREAL32) + real(PLSFIELD(IW+1:IE+1,:,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(IW:IE,:,:),kind=MNHREAL32) )
!
ELSE
!
@@ -1215,7 +1216,8 @@ CASE ('OPEN','WALL','NEST')
PRFIELDS(IW-1,:,:) = PRFIELDS(IW-1,:,:) + PRHODJ(IW-1,:,:) * &
PDK2*( &
PFIELDM(IW-2,:,:) -2.*PFIELDM(IW-1,:,:) + PFIELDM(IW,:,:) &
- -PLSFIELD(IW-2,:,:) +2.*PLSFIELD(IW-1,:,:) - PLSFIELD(IW,:,:) )
+ - real(PLSFIELD(IW-2,:,:),kind=MNHREAL32) +2.*real(PLSFIELD(IW-1,:,:),kind=MNHREAL32) &
+ - real(PLSFIELD(IW,:,:),kind=MNHREAL32) )
!
!!$ ELSEIF (NHALO == 1) THEN
ELSE
@@ -1225,9 +1227,9 @@ CASE ('OPEN','WALL','NEST')
TPHALO2%WEST(:,:) + PFIELDM(IW+1,:,:) &
-4.*( PFIELDM(IW-2,:,:) + PFIELDM(IW,:,:) ) &
+6.* PFIELDM(IW-1,:,:) &
- - TPHALO2LS%WEST(:,:) - PLSFIELD(IW+1,:,:) &
- +4.*( PLSFIELD(IW-2,:,:) + PLSFIELD(IW,:,:) ) &
- -6.* PLSFIELD(IW-1,:,:) )
+ - real(TPHALO2LS%WEST(:,:),kind=MNHREAL32) - real(PLSFIELD(IW+1,:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(IW-2,:,:),kind=MNHREAL32) + real(PLSFIELD(IW,:,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(IW-1,:,:),kind=MNHREAL32) )
!
ENDIF
!
@@ -1236,7 +1238,8 @@ CASE ('OPEN','WALL','NEST')
PRFIELDS(IE+1,:,:) = PRFIELDS(IE+1,:,:) + PRHODJ(IE+1,:,:) * &
PDK2*( &
PFIELDM(IE,:,:) -2.*PFIELDM(IE+1,:,:) + PFIELDM(IE+2,:,:) &
- - PLSFIELD(IE,:,:) +2.*PLSFIELD(IE+1,:,:) - PLSFIELD(IE+2,:,:) )
+ - real(PLSFIELD(IE,:,:),kind=MNHREAL32) +2.*real(PLSFIELD(IE+1,:,:),kind=MNHREAL32) &
+ - real(PLSFIELD(IE+2,:,:),kind=MNHREAL32) )
!
!!$ ELSEIF (NHALO == 1) THEN
ELSE
@@ -1246,9 +1249,9 @@ CASE ('OPEN','WALL','NEST')
PFIELDM(IE-1,:,:) + TPHALO2%EAST(:,:) &
-4.*( PFIELDM(IE ,:,:) + PFIELDM(IE+2,:,:) ) &
+6.* PFIELDM(IE+1,:,:) &
- - PLSFIELD(IE-1,:,:) - TPHALO2LS%EAST(:,:) &
- +4.*( PLSFIELD(IE ,:,:) + PLSFIELD(IE+2,:,:)) &
- -6.* PLSFIELD(IE+1,:,:))
+ - real(PLSFIELD(IE-1,:,:),kind=MNHREAL32) - real(TPHALO2LS%EAST(:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(IE ,:,:),kind=MNHREAL32) + real(PLSFIELD(IE+2,:,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(IE+1,:,:),kind=MNHREAL32) )
!
ENDIF
@@ -1262,9 +1265,9 @@ CASE ('OPEN','WALL','NEST')
PFIELDM(IW-2:IE-2,:,:) + PFIELDM(IW+2:IE+2,:,:) &
-4.*( PFIELDM(IW-1:IE-1,:,:) + PFIELDM(IW+1:IE+1,:,:) ) &
+6.* PFIELDM(IW:IE,:,:) &
- - PLSFIELD(IW-2:IE-2,:,:) - PLSFIELD(IW+2:IE+2,:,:) &
- +4.*( PLSFIELD(IW-1:IE-1,:,:) + PLSFIELD(IW+1:IE+1,:,:) ) &
- -6.* PLSFIELD(IW:IE,:,:))
+ - real(PLSFIELD(IW-2:IE-2,:,:),kind=MNHREAL32) - real(PLSFIELD(IW+2:IE+2,:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(IW-1:IE-1,:,:),kind=MNHREAL32) + real(PLSFIELD(IW+1:IE+1,:,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(IW:IE,:,:),kind=MNHREAL32) )
!
ELSE
!
@@ -1353,18 +1356,18 @@ IF ( .NOT. L2D ) THEN
TPHALO2%SOUTH(:,:) + PFIELDM(:,IS+1,:) &
-4.*( PFIELDM(:,IS-2,:) + PFIELDM(:,IS,:) ) &
+6.* PFIELDM(:,IS-1,:) &
- - TPHALO2LS%SOUTH(:,:) - PLSFIELD(:,IS+1,:) &
- +4.*( PLSFIELD(:,IS-2,:) + PLSFIELD(:,IS,:) ) &
- -6.* PLSFIELD(:,IS-1,:) )
+ - real(TPHALO2LS%SOUTH(:,:),kind=MNHREAL32) - real(PLSFIELD(:,IS+1,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(:,IS-2,:),kind=MNHREAL32) + real(PLSFIELD(:,IS,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(:,IS-1,:),kind=MNHREAL32) )
!
PRFIELDS(:,IN+1,:) = PRFIELDS(:,IN+1,:) - PRHODJ(:,IN+1,:) * &
PDK4*( &
PFIELDM(:,IN-1,:) + TPHALO2%NORTH(:,:) &
-4.*( PFIELDM(:,IN,:) + PFIELDM(:,IN+2,:) ) &
+6.* PFIELDM(:,IN+1,:) &
- - PLSFIELD(:,IN-1,:) - TPHALO2LS%NORTH(:,:) &
- +4.*( PLSFIELD(:,IN,:) + PLSFIELD(:,IN+2,:) ) &
- -6.* PLSFIELD(:,IN+1,:) )
+ - real(PLSFIELD(:,IN-1,:),kind=MNHREAL32) - real(TPHALO2LS%NORTH(:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(:,IN,:),kind=MNHREAL32) + real(PLSFIELD(:,IN+2,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(:,IN+1,:),kind=MNHREAL32) )
!
!!$ ENDIF
!
@@ -1376,9 +1379,9 @@ IF ( .NOT. L2D ) THEN
PFIELDM(:,IS-2:IN-2,:) + PFIELDM(:,IS+2:IN+2,:) &
-4.*( PFIELDM(:,IS-1:IN-1,:) + PFIELDM(:,IS+1:IN+1,:) ) &
+6.* PFIELDM(:,IS:IN,:) &
- - PLSFIELD(:,IS-2:IN-2,:) - PLSFIELD(:,IS+2:IN+2,:) &
- +4.*( PLSFIELD(:,IS-1:IN-1,:) + PLSFIELD(:,IS+1:IN+1,:) ) &
- -6.* PLSFIELD(:,IS:IN,:) )
+ - real(PLSFIELD(:,IS-2:IN-2,:),kind=MNHREAL32) - real(PLSFIELD(:,IS+2:IN+2,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(:,IS-1:IN-1,:),kind=MNHREAL32) + real(PLSFIELD(:,IS+1:IN+1,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(:,IS:IN,:),kind=MNHREAL32) )
!
ELSE
!
@@ -1448,8 +1451,9 @@ IF ( .NOT. L2D ) THEN
!
PRFIELDS(:,IS-1,:) = PRFIELDS(:,IS-1,:) + PRHODJ(:,IS-1,:) * &
PDK2*( &
- PFIELDM(:,IS-2,:) -2.*PFIELDM(:,IS-1,:) + PFIELDM(:,IS,:) &
- -PLSFIELD(:,IS-2,:) +2.*PLSFIELD(:,IS-1,:) - PLSFIELD(:,IS,:) )
+ PFIELDM(:,IS-2,:) -2.*PFIELDM(:,IS-1,:) + PFIELDM(:,IS,:) &
+ - real(PLSFIELD(:,IS-2,:),kind=MNHREAL32) +2.*real(PLSFIELD(:,IS-1,:),kind=MNHREAL32) &
+ - real(PLSFIELD(:,IS,:),kind=MNHREAL32) )
!
!!$ ELSEIF (NHALO == 1) THEN
ELSE
@@ -1459,9 +1463,9 @@ IF ( .NOT. L2D ) THEN
TPHALO2%SOUTH(:,:) + PFIELDM(:,IS+1,:) &
-4.*( PFIELDM(:,IS-2,:) + PFIELDM(:,IS,:) ) &
+6.* PFIELDM(:,IS-1,:) &
- - TPHALO2LS%SOUTH(:,:) - PLSFIELD(:,IS+1,:) &
- +4.*( PLSFIELD(:,IS-2,:) + PLSFIELD(:,IS,:) ) &
- -6.* PLSFIELD(:,IS-1,:) )
+ - real(TPHALO2LS%SOUTH(:,:),kind=MNHREAL32) - real(PLSFIELD(:,IS+1,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(:,IS-2,:),kind=MNHREAL32) + real(PLSFIELD(:,IS,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(:,IS-1,:),kind=MNHREAL32) )
!
ENDIF
!
@@ -1470,7 +1474,8 @@ IF ( .NOT. L2D ) THEN
PRFIELDS(:,IN+1,:) = PRFIELDS(:,IN+1,:) + PRHODJ(:,IN+1,:) * &
PDK2*( &
PFIELDM(:,IN,:) -2.*PFIELDM(:,IN+1,:) + PFIELDM(:,IN+2,:) &
- -PLSFIELD(:,IN,:) +2.*PLSFIELD(:,IN+1,:) - PLSFIELD(:,IN+2,:) )
+ - real(PLSFIELD(:,IN,:),kind=MNHREAL32) +2.*real(PLSFIELD(:,IN+1,:),kind=MNHREAL32) &
+ - real(PLSFIELD(:,IN+2,:),kind=MNHREAL32) )
!
!!$ ELSEIF (NHALO == 1) THEN
ELSE
@@ -1480,9 +1485,9 @@ IF ( .NOT. L2D ) THEN
PFIELDM(:,IN-1,:) + TPHALO2%NORTH(:,:) &
-4.*( PFIELDM(:,IN,:) + PFIELDM(:,IN+2,:) ) &
+6.* PFIELDM(:,IN+1,:) &
- - PLSFIELD(:,IN-1,:) - TPHALO2LS%NORTH(:,:) &
- +4.*( PLSFIELD(:,IN,:) + PLSFIELD(:,IN+2,:) ) &
- -6.* PLSFIELD(:,IN+1,:) )
+ - real(PLSFIELD(:,IN-1,:),kind=MNHREAL32) - real(TPHALO2LS%NORTH(:,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(:,IN,:),kind=MNHREAL32) + real(PLSFIELD(:,IN+2,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(:,IN+1,:),kind=MNHREAL32) )
!
ENDIF
!
@@ -1496,9 +1501,9 @@ IF ( .NOT. L2D ) THEN
PFIELDM(:,IS-2:IN-2,:) + PFIELDM(:,IS+2:IN+2,:) &
-4.*( PFIELDM(:,IS-1:IN-1,:) + PFIELDM(:,IS+1:IN+1,:) ) &
+6.* PFIELDM(:,IS:IN,:) &
- - PLSFIELD(:,IS-2:IN-2,:) - PLSFIELD(:,IS+2:IN+2,:) &
- +4.*( PLSFIELD(:,IS-1:IN-1,:) + PLSFIELD(:,IS+1:IN+1,:) ) &
- -6.* PLSFIELD(:,IS:IN,:) )
+ - real(PLSFIELD(:,IS-2:IN-2,:),kind=MNHREAL32) - real(PLSFIELD(:,IS+2:IN+2,:),kind=MNHREAL32) &
+ +4.*( real(PLSFIELD(:,IS-1:IN-1,:),kind=MNHREAL32) + real(PLSFIELD(:,IS+1:IN+1,:),kind=MNHREAL32) ) &
+ -6.* real(PLSFIELD(:,IS:IN,:),kind=MNHREAL32) )
!
ELSE
diff --git a/src/MNH/one_wayn.f90 b/src/MNH/one_wayn.f90
index 035f9498615835453bd56cbdf1462d6dabeeeb5c..eac7238d63ed4e9ee1b641155964a222bdc1273e 100644
--- a/src/MNH/one_wayn.f90
+++ b/src/MNH/one_wayn.f90
@@ -126,7 +126,7 @@ SUBROUTINE ONE_WAY_n(KDAD,PTSTEP,KMI,KTCOUNT, &
!* 0. DECLARATIONS
! ------------
USE MODD_CH_MNHC_n, only: LUSECHAQ, LUSECHIC
-USE MODD_CONF, only: CEQNSYS
+USE MODD_CONF, only: CEQNSYS,CCONF
USE MODD_CST, only: XCPD, XP00, XRD, XRV, XTH00
USE MODD_DYN_n, ONLY: LOCEAN
USE MODD_FIELD_n, only: XPABST, XRT, XSVT, XUT, XVT, XWT, XTHT, XTKET
@@ -683,22 +683,8 @@ IF(.NOT. OSTEADY_DMASS) THEN
!
!* 4.5 segment beginning (we have first to recover the dry mass at T-DT)
!
- IF(SIZE(XRT,4) == 0) THEN
- ! dry air case
-! ------------
- ZRHOD(:,:,:) = XPABST(:,:,:)/(XPABST(:,:,:)/XP00)**ZRD_O_CPD/(XRD*XTHT(:,:,:))
- ELSE ! moist air case
-! --------------
- ZRHOD(:,:,:) = XPABST(:,:,:)/(XPABST(:,:,:)/XP00)**ZRD_O_CPD/(XRD*XTHT(:,:,:) &
- *(1.+ZRV_O_RD*XRT(:,:,:,1)))
- ENDIF
-!
-!
- ZDRYMASSM = SUM3D_ll (ZJ(:,:,:)*ZRHOD(:,:,:),IINFO_ll,NXOR_ALL(KMI)+JPHEXT,NYOR_ALL(KMI)+JPHEXT, &
- 1+JPVEXT,NXEND_ALL(KMI)-JPHEXT,NYEND_ALL(KMI)-JPHEXT,SIZE(XRHODJ,3)-JPVEXT)
-!
- PDRYMASST = ZDRYMASST
- PDRYMASSS = (PDRYMASST - ZDRYMASSM) / (PTSTEP*KDTRATIO)
+ PDRYMASST = ZDRYMASST
+ IF ( CCONF /= 'RESTA' ) PDRYMASSS = 0.
ENDIF
!
END IF
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index 9c5525716dbf86a3953e636a14f29abc0cf781d4..4bfda374f9c9bbfe451b4564d89a245aa157afa1 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -679,6 +679,7 @@ IF (.NOT. OCLOUD_ONLY .AND. KTCOUNT /= 1) THEN
NDLON,NFLEV,CAER,NAER,NSTATM, &
XSINDEL,XCOSDEL,XTSIDER,XCORSOL, &
XSTATM,XOZON, XAER)
+ XAER_CLIM = XAER
END IF
END IF
!
diff --git a/src/MNH/rad_bound.f90 b/src/MNH/rad_bound.f90
index 296d476b44d32be4d2814e4614cfe0c7cfc1e1b2..22c423583eabe9d1a038f11cd51ef3fefd047ec0 100644
--- a/src/MNH/rad_bound.f90
+++ b/src/MNH/rad_bound.f90
@@ -159,9 +159,10 @@ END MODULE MODI_RAD_BOUND
!* 0. DECLARATIONS
! ------------
!
-USE MODD_CONF
+USE MODD_CONF
USE MODD_CTURB
USE MODD_PARAMETERS
+USE MODD_PRECISION, ONLY: MNHREAL32
USE MODD_RECYCL_PARAM_n, ONLY: LRECYCL, XRCOEFF
!
USE MODE_ll
@@ -317,6 +318,14 @@ SELECT CASE ( HLBCX(1) )
! + ZKTSTEP*( ZLBXU(:,:) ) ) &
! ) * ZINVTSTEP / (1.+ ZCPHASX(:,:) +ZKTSTEP)
!
+! ============================================================
+!
+! Reproductibility for RSTART -> truncate ZLB to real(knd=4) to have reproductible result
+!
+ ZLBEU = real(ZLBEU,kind=MNHREAL32)
+ ZLBGU = real(ZLBGU,kind=MNHREAL32)
+ ZLBXU = real(ZLBXU,kind=MNHREAL32)
+! ============================================================
PRUS (IIB,:,:) =(PRHODJ(IIB-1,:,:) + PRHODJ(IIB,:,:)) * 0.5 * &
ZINVTSTEP / (1.+ ZKTSTEP * ZALPHA2 ) * &
( (1. - ZCPHASX(:,:) - ZKTSTEP * (1. - ZALPHA2)) * PUT(IIB,:,:) &
@@ -374,7 +383,7 @@ SELECT CASE ( HLBCX(2) )
ZLBXU(:,:) = PLBXUM(ILBX-JPHEXT+1,:,:)
END IF
ELSE
- ZLBEU (:,:) = PLBXUS(ILBX-JPHEXT+1,:,:)
+ ZLBEU (:,:) = PLBXUS(ILBX-JPHEXT+1,:,:)
ZLBGU (:,:) = PLBXUM(ILBX-JPHEXT+1,:,:) - PLBXUM(ILBX-JPHEXT,:,:) + &
PTSTEP * (PLBXUS(ILBX-JPHEXT+1,:,:) - PLBXUS(ILBX-JPHEXT,:,:))
IF ( LRECYCL ) THEN
@@ -393,7 +402,15 @@ SELECT CASE ( HLBCX(2) )
! + ZLBGU (:,:) * ZCPHASX(:,:) &
! + ZKTSTEP*ZLBXU(:,:) ) &
! ) * ZINVTSTEP / (1.+ZCPHASX(:,:) +ZKTSTEP)
-!
+!
+! ============================================================
+!
+! Reproductibility for RSTART -> truncate ZLB to real(knd=4)
+!
+ ZLBEU = real(ZLBEU,kind=MNHREAL32)
+ ZLBGU = real(ZLBGU,kind=MNHREAL32)
+ ZLBXU = real(ZLBXU,kind=MNHREAL32)
+! ============================================================
PRUS (IIE+1,:,:) =(PRHODJ(IIE+1,:,:) + PRHODJ(IIE,:,:)) * 0.5 * &
ZINVTSTEP / (1.+ ZKTSTEP * ZALPHA2 ) * &
( (1. - ZCPHASX(:,:) - ZKTSTEP * (1. - ZALPHA2) ) * PUT(IIE+1,:,:) &
@@ -469,6 +486,15 @@ SELECT CASE ( HLBCY(1) )
! - ZLBGV (:,:) * ZCPHASY(:,:) &
! + ZKTSTEP*ZLBYV(:,:) ) &
! ) * ZINVTSTEP / (1.+ ZCPHASY(:,:) +ZKTSTEP)
+!
+! ============================================================
+!
+! Reproductibility for RSTART -> truncate ZLB to real(knd=4) to have reproductible result
+!
+ ZLBEV = real(ZLBEV,kind=MNHREAL32)
+ ZLBGV = real(ZLBGV,kind=MNHREAL32)
+ ZLBYV = real(ZLBYV,kind=MNHREAL32)
+! ============================================================
PRVS (:,IJB,:) =(PRHODJ(:,IJB-1,:) + PRHODJ(:,IJB,:)) * 0.5 * &
ZINVTSTEP / (1.+ ZKTSTEP * ZALPHA2 ) * &
( (1. - ZCPHASY(:,:) - ZKTSTEP * (1. - ZALPHA2) ) * PVT(:,IJB,:)&
@@ -546,6 +572,14 @@ SELECT CASE ( HLBCY(2) )
! + ZKTSTEP* ZLBYV(:,:) ) &
! ) * ZINVTSTEP / (1.+ ZCPHASY(:,:) +ZKTSTEP)
!
+! ============================================================
+!
+! Reproductibility for RSTART -> truncate ZLB to real(knd=4) to have reproductible result
+!
+ ZLBEV = real(ZLBEV,kind=MNHREAL32)
+ ZLBGV = real(ZLBGV,kind=MNHREAL32)
+ ZLBYV = real(ZLBYV,kind=MNHREAL32)
+! ============================================================
PRVS (:,IJE+1,:) =(PRHODJ(:,IJE+1,:) + PRHODJ(:,IJE,:)) * 0.5 * &
ZINVTSTEP / (1.+ ZKTSTEP * ZALPHA2 ) * &
( (1. - ZCPHASY(:,:) - ZKTSTEP * (1. - ZALPHA2) ) * PVT(:,IJE+1,:)&
diff --git a/src/MNH/read_all_data_grib_case.f90 b/src/MNH/read_all_data_grib_case.f90
index 84c13799185ffe6e85d837934b7348cf51c022f6..841ae62b0c834902a5fec42fcbbac2c947e6d8c3 100644
--- a/src/MNH/read_all_data_grib_case.f90
+++ b/src/MNH/read_all_data_grib_case.f90
@@ -136,6 +136,7 @@ END MODULE MODI_READ_ALL_DATA_GRIB_CASE
! Q. Rodier 21/04/2020: correction GFS u and v wind component written in the right vertical order
! Q. Rodier 02/09/2020: Read and interpol geopotential height for interpolation on isobaric surface Grid of NCEP
! P. Wautelet 09/03/2021: move some chemistry initializations to ini_nsv
+!JP Chaboureau 02/08/2021: add ERA5 reanlysis in pressure levels
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -202,10 +203,11 @@ REAL, INTENT(INOUT) :: PTIME_HORI ! time spent in hor. interpolati
! ------------------------------
! General purpose variables
INTEGER :: ILUOUT0 ! Unit used for output msg.
-INTEGER :: IRESP ! Return code of FM-routines
+INTEGER :: IRESP ! Return code of FM-routines
INTEGER :: IRET ! Return code from subroutines
-INTEGER(KIND=kindOfInt) :: IRET_GRIB ! Return code from subroutines
+INTEGER(KIND=kindOfInt) :: IRET_GRIB ! Return code from subroutines
INTEGER, PARAMETER :: JP_GFS=31 ! number of pressure levels for GFS model
+INTEGER, PARAMETER :: JP_ERA=37 ! number of pressure levels for ERA5 reanalysis
REAL :: ZA,ZB,ZC ! Dummy variables
REAL :: ZD,ZE,ZF ! |
REAL :: ZTEMP ! |
@@ -230,7 +232,7 @@ REAL, DIMENSION(:,:), ALLOCATABLE :: ZYM ! Y of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZLATM ! Lat of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZLONM ! Lon of PGD mass points
! Variable involved in the task of reading the grib file
-INTEGER(KIND=kindOfInt) :: IUNIT ! unit of the grib file
+INTEGER(KIND=kindOfInt) :: IUNIT ! unit of the grib file
CHARACTER(LEN=50) :: HGRID ! type of grid
INTEGER :: IPAR ! Parameter identifier
INTEGER :: ITYP ! type of level (Grib code table 3)
@@ -246,13 +248,13 @@ INTEGER :: IMODEL ! Type of Grib file :
! 10 -> NCEP - GFS
INTEGER :: ICENTER ! number of center
INTEGER :: ISIZE ! size of grib message
-INTEGER(KIND=kindOfInt) :: ICOUNT ! number of messages in the file
+INTEGER(KIND=kindOfInt) :: ICOUNT ! number of messages in the file
INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: IGRIB ! number of the grib in memory
INTEGER :: INUM ,INUM_ZS ! number of a grib message
-REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM ! parameter of girb grid
+REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM ! parameter of grib grid
INTEGER,DIMENSION(:),ALLOCATABLE :: IINLO ! longitude of grib grid
INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: IINLO_GRIB ! longitude of grib grid
-REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM_ZS ! parameter of girb grid for ZS
+REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM_ZS ! parameter of grib grid for ZS
INTEGER,DIMENSION(:),ALLOCATABLE :: IINLO_ZS ! longitude of grib grid for ZS
REAL,DIMENSION(:),ALLOCATABLE :: ZVALUE ! Intermediate array
REAL,DIMENSION(:),ALLOCATABLE :: ZOUT ! Intermediate arrays
@@ -264,11 +266,11 @@ TYPE(DATE_TIME) :: TPTCUR ! Date & time of the grib fi
INTEGER :: ITWOZS
! surface pressure
REAL, DIMENSION(:), ALLOCATABLE :: ZPS_G ! Grib data : Ps
-REAL, DIMENSION(:), ALLOCATABLE :: ZLNPS_G ! Grib data : ln(Ps)
-REAL, DIMENSION(:), ALLOCATABLE :: ZWORK_LNPS ! Grib data on zs grid: ln(Ps)
+REAL, DIMENSION(:), ALLOCATABLE :: ZLNPS_G ! Grib data : ln(Ps)
+REAL, DIMENSION(:), ALLOCATABLE :: ZWORK_LNPS ! Grib data on zs grid: ln(Ps)
INTEGER :: INJ,INJ_ZS
! orography
-CHARACTER(LEN=50) :: HGRID_ZS ! type of grid
+CHARACTER(LEN=50) :: HGRID_ZS ! type of grid
!
! Reading and projection of the wind vectors u, v
REAL :: ZALPHA ! Angle of rotation
@@ -330,14 +332,17 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTEV_LS ! T V
REAL, DIMENSION(:), ALLOCATABLE :: ZPV ! vertical level in grib file
INTEGER :: IPVPRESENT ,IPV
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZR_DUM
-INTEGER :: IMI
+INTEGER :: IMI
TYPE(TFILEDATA),POINTER :: TZFILE
INTEGER, DIMENSION(JP_GFS) :: IP_GFS ! list of pressure levels for GFS model
+INTEGER, DIMENSION(JP_ERA) :: IP_ERA ! list of pressure levels for ERA5 reanalysis
INTEGER :: IVERSION,ILEVTYPE
-LOGICAL :: GFIND ! to test if sea wave height is found
+LOGICAL :: GFIND ! to test if sea wave height is found
!---------------------------------------------------------------------------------------
IP_GFS=(/1000,975,950,925,900,850,800,750,700,650,600,550,500,450,400,350,300,&
- 250,200,150,100,70,50,30,20,10,7,5,3,2,1/)!
+ 250,200,150,100,70,50,30,20,10,7,5,3,2,1/)
+IP_ERA=(/1000,975,950,925,900,875,850,825,800,775,750,700,650,600,550,500,450,&
+ 400,350,300,250,225,200,175,150,125,100,70,50,30,20,10,7,5,3,2,1/)
!
TZFILE => NULL()
!
@@ -567,7 +572,6 @@ ELSE IF (HFILE=='CHEM') THEN
END IF
DEALLOCATE (ZOUT)
!
-! *** BEGIN MODIF SB ADD HS ***
!---------------------------------------------------------------------------------------
!* 2.3 bis Read and interpol Sea Wave significant height
!---------------------------------------------------------------------------------------
@@ -593,7 +597,7 @@ SELECT CASE (IMODEL)
GFIND=.TRUE.
END IF
!
- IF(GFIND) THEN
+ IF (GFIND) THEN
!!!!!!!!!!! Faire en sorte de le faire que pour le CASE(0)
! Sea wave significant height disponible uniquement pour ECMWF
! recuperation du tableau de valeurs
@@ -615,7 +619,6 @@ SELECT CASE (IMODEL)
DEALLOCATE (ZOUT)
END IF
END SELECT
- ! *** END MODIF SB ADD HS ***
!
!---------------------------------------------------------------------------------------
!* 2.4 Interpolation surface pressure
@@ -628,6 +631,12 @@ WRITE (ILUOUT0,'(A)') ' | Searching pressure'
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=152)
+ IF( INUM < 0 ) THEN
+ WRITE (ILUOUT0,'(A)') ' | Logarithm of surface pressure is missing. It is then supposed that'
+ WRITE (ILUOUT0,'(A)') ' | this ECMWF file has atmospheric fields on pressure levels (e.g. ERA5)'
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=134)
+ IMODEL = 11
+ END IF
CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=1)
CASE(6,7) ! NEW AROME,ARPEGE
@@ -647,7 +656,7 @@ SELECT CASE (IMODEL)
ALLOCATE (ZLNPS_G(ISIZE))
ZLNPS_G(:) = ZVALUE(1:ISIZE)
ZPS_G (:) = EXP(ZVALUE(1:ISIZE))
- CASE(1,2,3,4,5,10) ! arpege mocage aladin aladin-reunion NCEP
+ CASE(1,2,3,4,5,10,11) ! arpege mocage aladin aladin-reunion NCEP ERA5
ALLOCATE (ZPS_G (ISIZE))
ALLOCATE (ZLNPS_G(ISIZE))
ZPS_G (:) = ZVALUE(1:ISIZE)
@@ -708,7 +717,7 @@ DEALLOCATE (ZLNPS_G)
!
WRITE (ILUOUT0,'(A)') ' | Reading T and Q fields'
!
-IF (IMODEL/=10) THEN
+IF (IMODEL/=10.AND.IMODEL/=11) THEN
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
ISTARTLEVEL=1
@@ -732,7 +741,7 @@ IF (IMODEL/=10) THEN
IF(INUM < 0) call Print_msg( NVERB_FATAL, 'IO', 'READ_ALL_DATA_GRIB_CASE', 'air temperature is missing' )
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ISTARTLEVEL)
IF(INUM < 0) call Print_msg( NVERB_FATAL, 'IO', 'READ_ALL_DATA_GRIB_CASE', 'atmospheric specific humidity is missing' )
-ELSE ! NCEP
+ELSEIF (IMODEL==10) THEN ! NCEP
ISTARTLEVEL=1000
IT=130
IQ=157
@@ -740,20 +749,32 @@ ELSE ! NCEP
IF(INUM < 0) call Print_msg( NVERB_FATAL, 'IO', 'READ_ALL_DATA_GRIB_CASE', 'air temperature is missing' )
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ISTARTLEVEL)
IF(INUM < 0) call Print_msg( NVERB_FATAL, 'IO', 'READ_ALL_DATA_GRIB_CASE', 'atmospheric relative humidity is missing' )
+ELSE ! ERA5
+ ISTARTLEVEL=1000
+ IT=130
+ IQ=133
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ISTARTLEVEL)
+ IF(INUM < 0) call Print_msg( NVERB_FATAL, 'IO', 'READ_ALL_DATA_GRIB_CASE', 'air temperature is missing' )
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ISTARTLEVEL)
+ IF(INUM < 0) call Print_msg( NVERB_FATAL, 'IO', 'READ_ALL_DATA_GRIB_CASE', 'atmospheric specific humidity is missing' )
ENDIF
!
-IF (IMODEL/=10) THEN ! others than NCEP
+IF (IMODEL/=10.AND.IMODEL/=11) THEN ! others than NCEP AND ERA5
CALL GRIB_GET(IGRIB(INUM),'NV',INLEVEL)
INLEVEL = NINT(INLEVEL / 2.) - 1
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
ELSE
- INLEVEL=JP_GFS
+ IF (IMODEL==10) THEN
+ INLEVEL=JP_GFS
+ ELSE
+ INLEVEL=JP_ERA
+ END IF
END IF
!
ALLOCATE (ZT_G(ISIZE,INLEVEL))
ALLOCATE (ZQ_G(ISIZE,INLEVEL))
!
-IF (IMODEL/=10) THEN ! others than NCEP
+IF (IMODEL/=10.AND.IMODEL/=11) THEN ! others than NCEP AND ERA5
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ILEV1)
@@ -770,7 +791,7 @@ IF (IMODEL/=10) THEN ! others than NCEP
CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,INLEVEL-JLOOP1+1))
CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
END DO
-ELSE ! NCEP
+ELSEIF (IMODEL==10) THEN ! NCEP
DO JLOOP1=1, INLEVEL
ILEV1 = IP_GFS(JLOOP1)
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ILEV1)
@@ -779,14 +800,29 @@ ELSE ! NCEP
CALL PRINT_MSG(NVERB_FATAL,'IO','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
CALL GRIB_GET(IGRIB(INUM),'values',ZQ_G(:,JLOOP1),IRET_GRIB)
- WRITE (ILUOUT0,*) 'Q ',ILEV1,IRET_GRIB
CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=0,KNUMBER=0,KLEV1=ILEV1,KTFFS=100)
IF (INUM< 0) THEN
WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
CALL PRINT_MSG(NVERB_FATAL,'IO','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,JLOOP1),IRET_GRIB)
- WRITE (ILUOUT0,*) 'T ',ILEV1,IRET_GRIB
+ CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
+ END DO
+ELSE ! ERA5
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = IP_ERA(JLOOP1)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ILEV1)
+ IF (INUM< 0) THEN
+ WRITE(YMSG,*) 'atmospheric humidity level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'IO','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET(IGRIB(INUM),'values',ZQ_G(:,JLOOP1),IRET_GRIB)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ILEV1)
+ IF (INUM< 0) THEN
+ WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'IO','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,JLOOP1),IRET_GRIB)
CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
END DO
END IF
@@ -797,7 +833,7 @@ CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
!
!* 2.5.2 Load level definition parameters A and B
!
-IF (IMODEL/=10) THEN ! others than NCEP
+IF (IMODEL/=10.AND.IMODEL/=11) THEN ! others than NCEP AND ERA5
IF (HFILE(1:3)=='ATM') THEN
XP00_LS = 101325.
@@ -853,14 +889,19 @@ IF (IMODEL/=10) THEN ! others than NCEP
ELSE
ALLOCATE (XA_LS(INLEVEL))
ALLOCATE (XB_LS(0))
- XA_LS = 100.*IP_GFS
+ IF (IMODEL==10) THEN
+ XA_LS = 100.*IP_GFS
+ ELSE
+ XA_LS = 100.*IP_ERA
+ END IF
END IF
!
!* 2.5.3 Compute atmospheric pressure on grib grid
!
WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on Grib grid is being computed'
+
ALLOCATE (ZPF_G(INI,INLEVEL))
-IF (IMODEL/=10) THEN ! others than NCEP
+IF (IMODEL/=10.AND.IMODEL/=11) THEN ! others than NCEP and ERA5
IF (HFILE(1:3)=='ATM') THEN
ZPF_G(:,:) = SPREAD(XA_LS,1,INI)*XP00_LS + &
SPREAD(XB_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
@@ -869,7 +910,11 @@ IF (IMODEL/=10) THEN ! others than NCEP
SPREAD(XB_SV_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
END IF
ELSE
- ZPF_G(:,:) = 100.*SPREAD(IP_GFS,1,INI)
+ IF (IMODEL==10) THEN
+ ZPF_G(:,:) = 100.*SPREAD(IP_GFS,1,INI)
+ ELSE
+ ZPF_G(:,:) = 100.*SPREAD(IP_ERA,1,INI)
+ END IF
END IF
DEALLOCATE (ZPS_G)
!
@@ -965,6 +1010,8 @@ ALLOCATE (ZRV_G(INI))
ALLOCATE (ZOUT(INO))
IF (IMODEL/=10) THEN ! others than NCEP
DO JLOOP1=1, INLEVEL
+ !WRITE (ILUOUT0,*) 'JLOOP1=',JLOOP1,MINVAL(ZPM_G(:,JLOOP1)),MINVAL(ZT_G(:,JLOOP1)),MINVAL(ZQ_G(:,JLOOP1))
+ !WRITE (ILUOUT0,*) ' ',MAXVAL(ZPM_G(:,JLOOP1)),MAXVAL(ZT_G(:,JLOOP1)),MAXVAL(ZQ_G(:,JLOOP1))
!
! Compute Theta V and relative humidity on grib grid
!
@@ -994,15 +1041,15 @@ IF (IMODEL/=10) THEN ! others than NCEP
END DO
ELSE !NCEP
DO JLOOP1=1, INLEVEL
- WRITE (ILUOUT0,*) 'JLOOP1=',JLOOP1,MINVAL(ZPM_G(:,JLOOP1)),MINVAL(ZT_G(:,JLOOP1)),MINVAL(ZQ_G(:,JLOOP1))
- WRITE (ILUOUT0,*) ' ',MAXVAL(ZPM_G(:,JLOOP1)),MAXVAL(ZT_G(:,JLOOP1)),MAXVAL(ZQ_G(:,JLOOP1))
+ !WRITE (ILUOUT0,*) 'JLOOP1=',JLOOP1,MINVAL(ZPM_G(:,JLOOP1)),MINVAL(ZT_G(:,JLOOP1)),MINVAL(ZQ_G(:,JLOOP1))
+ !WRITE (ILUOUT0,*) ' ',MAXVAL(ZPM_G(:,JLOOP1)),MAXVAL(ZT_G(:,JLOOP1)),MAXVAL(ZQ_G(:,JLOOP1))
ZH_G(:) =ZQ_G(:,JLOOP1)
ZRV_G(:) = (XRD/XRV)*SM_FOES(ZT_G(:,JLOOP1))*0.01*ZH_G(:) &
/(ZPM_G(:,JLOOP1) -SM_FOES(ZT_G(:,JLOOP1))*0.01*ZH_G(:))
- WRITE (ILUOUT0,*) ' ',MINVAL(ZRV_G(:)),MAXVAL(ZRV_G(:))
+ !WRITE (ILUOUT0,*) ' ',MINVAL(ZRV_G(:)),MAXVAL(ZRV_G(:))
ZTHV_G(:)=ZT_G(:,JLOOP1) * ((XP00/ZPM_G(:,JLOOP1))**(XRD/XCPD)) * &
((1. + XRV*ZRV_G(:)/XRD) / (1. + ZRV_G(:)))
- WRITE (ILUOUT0,*) ' ',MINVAL(ZTHV_G(:)),MAXVAL(ZTHV_G(:))
+ !WRITE (ILUOUT0,*) ' ',MINVAL(ZTHV_G(:)),MAXVAL(ZTHV_G(:))
!
! Interpolation : H
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
@@ -1027,20 +1074,29 @@ DEALLOCATE (ZOUT)
!
ALLOCATE (ZGH_G(ISIZE,INLEVEL))
!
-IF(IMODEL==10) THEN !NCEP with pressure grid only
+IF (IMODEL==10.OR.IMODEL==11) THEN !NCEP or ERA5 with pressure grid only
DO JLOOP1=1, INLEVEL
- ILEV1 = IP_GFS(JLOOP1)
- WRITE (ILUOUT0,'(A)') ' | Searching geopotential height'
- CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=3,KNUMBER=5,KLEV1=ILEV1)
+ IF (IMODEL==10) THEN
+ ILEV1 = IP_GFS(JLOOP1)
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=3,KNUMBER=5,KLEV1=ILEV1)
+ ELSE
+ ILEV1 = IP_ERA(JLOOP1)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=129,KLEV1=ILEV1)
+ END IF
IF (INUM< 0) THEN
!callabortstop
- WRITE(YMSG,*) 'Geopoential height level ',JLOOP1,' is missing'
+ WRITE(YMSG,*) 'Geopotential height level ',JLOOP1,' is missing'
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
!
CALL GRIB_GET(IGRIB(INUM),'values',ZGH_G(:,JLOOP1),IRET_GRIB)
CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
!
+ IF (IMODEL/=10) THEN
+ ! Data given in archives are multiplied by the gravity acceleration
+ ZGH_G(:,JLOOP1) = ZGH_G(:,JLOOP1) / XG
+ END IF
+ !
END DO
!
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM_ZS),IIU,IJU,ZLONOUT,ZLATOUT,&
@@ -1059,9 +1115,9 @@ END IF
!
!* 2.5.5 Compute atmospheric pressure on MESO-NH grid
!
-WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on MesoNH grid is being computed'
+WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on the Meso-NH grid is being computed'
ALLOCATE (ZPF_LS(IIU,IJU,INLEVEL))
-IF (IMODEL/=10) THEN ! others than NCEP
+IF (IMODEL/=10.AND.IMODEL/=11) THEN ! others than NCEP and ERA5
IF (HFILE(1:3)=='ATM') THEN
ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_LS,1,IIU),2,IJU)*XP00_LS + &
SPREAD(SPREAD(XB_LS,1,IIU),2,IJU)*SPREAD(XPS_LS,3,INLEVEL)
@@ -1070,7 +1126,11 @@ IF (IMODEL/=10) THEN ! others than NCEP
SPREAD(SPREAD(XB_SV_LS,1,IIU),2,IJU)*SPREAD(XPS_SV_LS,3,INLEVEL)
END IF
ELSE
- ZPF_LS(:,:,:) = 100.*SPREAD(SPREAD(IP_GFS,1,IIU),2,IJU)
+ IF(IMODEL==10) THEN
+ ZPF_LS(:,:,:) = 100.*SPREAD(SPREAD(IP_GFS,1,IIU),2,IJU)
+ ELSE
+ ZPF_LS(:,:,:) = 100.*SPREAD(SPREAD(IP_ERA,1,IIU),2,IJU)
+ END IF
END IF
!
ALLOCATE (ZEXNF_LS(IIU,IJU,INLEVEL))
@@ -1420,11 +1480,6 @@ END IF
! After this projection, the file is simil
!
IF (HFILE(1:3)=='ATM') THEN
-IF (IMODEL/=10) THEN ! others than NCEP
- ISTARTLEVEL = 1
-ELSE
- ISTARTLEVEL = 10
-END IF
ALLOCATE (XU_LS(IIU,IJU,INLEVEL))
ALLOCATE (XV_LS(IIU,IJU,INLEVEL))
ALLOCATE (ZTU_LS(INO))
@@ -1445,16 +1500,20 @@ SELECT CASE (IMODEL)
ISTARTLEVEL = 0
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
END IF
- CASE (10)
+ CASE (10,11)
IPAR = 131
ISTARTLEVEL = 1
END SELECT
DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
- IF (IMODEL/=10) THEN ! others than NCEP
+ IF (IMODEL/=10.AND.IMODEL/=11) THEN ! others than NCEP and ERA5
ILEV1 = JLOOP1
ELSE
- ILEV1 = IP_GFS(INLEVEL+ISTARTLEVEL-JLOOP1)
+ IF(IMODEL==10) THEN
+ ILEV1 = IP_GFS(INLEVEL+ISTARTLEVEL-JLOOP1)
+ ELSE
+ ILEV1 = IP_ERA(INLEVEL+ISTARTLEVEL-JLOOP1)
+ END IF
END IF
! read component u
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ILEV1)
@@ -1483,10 +1542,14 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
END IF
DEALLOCATE (ZVALUE)
! read component v and perform interpolation if not Arpege grid
- IF (IMODEL/=10) THEN ! others than NCEP
+ IF (IMODEL/=10.AND.IMODEL/=11) THEN ! others than NCEP and ERA5
ILEV1 = JLOOP1
ELSE
- ILEV1 = IP_GFS(INLEVEL+ISTARTLEVEL-JLOOP1)
+ IF(IMODEL==10) THEN
+ ILEV1 = IP_GFS(INLEVEL+ISTARTLEVEL-JLOOP1)
+ ELSE
+ ILEV1 = IP_ERA(INLEVEL+ISTARTLEVEL-JLOOP1)
+ END IF
END IF
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR+1,KLEV1=ILEV1)
IF (INUM < 0) THEN
@@ -1661,7 +1724,7 @@ TPTCUR%nmonth=INT((IDATE-TPTCUR%nyear*10000)/100)
TPTCUR%nday=IDATE-TPTCUR%nyear*10000-TPTCUR%nmonth*100
CALL GRIB_GET(IGRIB(INUM),'startStep',ITIMESTEP,IRET_GRIB)
CALL GRIB_GET(IGRIB(INUM),'stepUnits',CSTEPUNIT,IRET_GRIB)
-IF (IMODEL==0) THEN
+IF (IMODEL==0.OR.IMODEL==11) THEN
ITWOZS=0
IF ((TPTCUR%nyear ==2000).AND.(TPTCUR%nmonth >11)) ITWOZS=1
IF ((TPTCUR%nyear ==2000).AND.(TPTCUR%nmonth ==11)) THEN
@@ -2135,7 +2198,7 @@ ALLOCATE(INLO_GRIB(SIZE(KINLO)))
INO= KNOLON*KNOLARG
SELECT CASE (KMODEL)
!
-CASE(0,5) ! CEP/MOCAGE
+CASE(0,5,11) ! CEP/MOCAGE/ERA5
! en theorie il faut ces 4 lignes
! CALL GRIB_GET(KGRIB,'latitudeOfFirstGridPointInDegrees',ZILA1)
! CALL GRIB_GET(KGRIB,'longitudeOfFirstGridPointInDegrees',ZILO1)
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
index a8c847c3f280ae1e4f7e1a27561252b676a3c574..9db15fc163ad08c775913431fff5c97226d3839e 100644
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -2149,6 +2149,13 @@ IF ( CRAD /= 'NONE' .AND. CPROGRAM=='MESONH' ) THEN
IF(CCONF=='START') THEN
CGETRAD='INIT'
END IF
+ IF(CCONF=='RESTA' .AND. (.NOT. LAERO_FT) .AND. (.NOT. LORILAM) &
+ .AND. (.NOT. LSALT) .AND. (.NOT. LDUST)) THEN
+ WRITE(UNIT=ILUOUT,FMT=9001) KMI
+ WRITE(UNIT=ILUOUT,FMT=*) '!!! WARNING !!! FOR REPRODUCTIBILITY BETWEEN START and START+RESTART,'
+ WRITE(UNIT=ILUOUT,FMT=*) 'YOU MUST USE LAERO_FT=T WITH CAER=TEGE IF CCONF=RESTA IN ALL SEGMENTS'
+ WRITE(UNIT=ILUOUT,FMT=*) 'TO UPDATE THE OZONE AND AEROSOLS CLIMATOLOGY USED BY THE RADIATION CODE;'
+ END IF
END IF
!
! 3.6 check the initialization of the deep convection scheme
diff --git a/src/MNH/read_field.f90 b/src/MNH/read_field.f90
index 1f8d4b3cab8e98abd0698ec2c34249c956ec75ba..f7ccb114e6c605c3aad6f3585688b88e4f2b6b8b 100644
--- a/src/MNH/read_field.f90
+++ b/src/MNH/read_field.f90
@@ -18,7 +18,7 @@ INTERFACE
KSIZELBXR_ll,KSIZELBYR_ll,KSIZELBXSV_ll,KSIZELBYSV_ll, &
PUM,PVM,PWM,PDUM,PDVM,PDWM, &
PUT,PVT,PWT,PTHT,PPABST,PTKET,PRTKEMS, &
- PRT,PSVT,PZWS,PCIT,PDRYMASST, &
+ PRT,PSVT,PZWS,PCIT,PDRYMASST,PDRYMASSS, &
PSIGS,PSRCT,PCLDFR,PBL_DEPTH,PSBL_DEPTH,PWTHVMF,PPHC,PPHR, &
PLSUM,PLSVM,PLSWM,PLSTHM,PLSRVM, PLSZWSM, &
PLBXUM,PLBXVM,PLBXWM,PLBXTHM,PLBXTKEM,PLBXRM,PLBXSVM, &
@@ -81,6 +81,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRCT ! turbulent flux
! <s'Rc'> at t
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCIT ! ice conc. at t
REAL, INTENT(OUT) :: PDRYMASST ! Md(t)
+REAL, INTENT(OUT) :: PDRYMASSS ! d Md(t) / dt
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSIGS ! =sqrt(<s's'>) for the
! Subgrid Condensation
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCLDFR ! cloud fraction
@@ -141,7 +142,7 @@ END MODULE MODI_READ_FIELD
KSIZELBXR_ll,KSIZELBYR_ll,KSIZELBXSV_ll,KSIZELBYSV_ll, &
PUM,PVM,PWM,PDUM,PDVM,PDWM, &
PUT,PVT,PWT,PTHT,PPABST,PTKET,PRTKEMS, &
- PRT,PSVT,PZWS,PCIT,PDRYMASST, &
+ PRT,PSVT,PZWS,PCIT,PDRYMASST,PDRYMASSS, &
PSIGS,PSRCT,PCLDFR,PBL_DEPTH,PSBL_DEPTH,PWTHVMF,PPHC,PPHR, &
PLSUM,PLSVM,PLSWM,PLSTHM,PLSRVM,PLSZWSM, &
PLBXUM,PLBXVM,PLBXWM,PLBXTHM,PLBXTKEM,PLBXRM,PLBXSVM, &
@@ -361,6 +362,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRCT ! turbulent flux
! <s'Rc'> at t
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCIT ! ice conc. at t
REAL, INTENT(OUT) :: PDRYMASST ! Md(t)
+REAL, INTENT(OUT) :: PDRYMASSS ! d Md(t) / dt
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSIGS ! =sqrt(<s's'>) for the
! Subgrid Condensation
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCLDFR ! cloud fraction
@@ -1464,6 +1466,11 @@ CALL INI_LB(TPINIFILE,GLSOURCE,ISV, &
!* 2.3 Some special variables:
!
CALL IO_Field_read(TPINIFILE,'DRYMASST',PDRYMASST) ! dry mass
+IF (CCONF=='RESTA') THEN
+ CALL IO_Field_read(TPINIFILE,'DRYMASSS',PDRYMASSS) ! dry mass tendency
+ELSE
+ PDRYMASSS=XUNDEF ! should not be used
+END IF
!
SELECT CASE(HGETSRCT) ! turbulent flux SRC at time t
CASE('READ')
diff --git a/src/MNH/recycling.f90 b/src/MNH/recycling.f90
index 9734eebc00e2f5bc16e8e923bfdc04cba2f120e7..18e0956cd05cb710e26718962a054e801d79636b 100644
--- a/src/MNH/recycling.f90
+++ b/src/MNH/recycling.f90
@@ -65,35 +65,17 @@ END MODULE MODI_RECYCLING
!**** 0. DECLARATIONS
! ---------------
!
-! module
-USE MODE_POS
-USE MODE_ll
-USE MODE_IO
-!USE MODI_SHUMAN
-!
-USE MODD_PARAMETERS
-USE MODD_CONF
-!
-USE MODD_CST
-!
-USE MODD_DIM_n
-USE MODD_CONF
-USE MODD_CONF_n
-USE MODD_GRID
-USE MODD_GRID_n
-USE MODD_METRICS_n
-USE MODD_TIME
-USE MODD_TIME_n
-USE MODD_DYN_n
-USE MODD_FIELD_n
-USE MODD_CURVCOR_n
-USE MODD_REF
-!
-USE MODD_VAR_ll, ONLY: IP, NPROC
+USE MODD_CONF, ONLY: CCONF
+USE MODD_FIELD_n, ONLY: XTHT, XUT, XVT, XWT
+USE MODD_GRID_n, ONLY: XXHAT, XYHAT
+USE MODD_METRICS_n, ONLY: XDZZ
+USE MODD_PARAMETERS, ONLY: JPHEXT
USE MODD_RECYCL_PARAM_n
+
+USE MODE_MSG
+
USE MODI_RECYCL_FLUC
-USE MODD_LUNIT_n, ONLY : TLUOUT
-!
+
IMPLICIT NONE
!
!------------------------------------------------------------------------------
@@ -106,9 +88,6 @@ REAL, DIMENSION(:,:) ,INTENT(INOUT) :: PFLUCTUNE,PFLUCTVTE,PFLUCTVNS,PFLUCTU
!------------------------------------------------------------------------------
!
! 0.2 declaration of local variables
-INTEGER :: IIU,IJU,IKU,JCOUNT,ICOUNT,ILUOUT
-INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE,IIP
-INTEGER :: IIBG,IIEG,IJBG,IJEG,IIMAX,IJMAX
INTEGER :: PMINW,PMINE,PMINN,PMINS
INTEGER :: JIDIST,JJDIST
REAL :: Z_DELTX,Z_DELTY
@@ -116,24 +95,14 @@ REAL :: Z_DELTX,Z_DELTY
!------------------------------------------------------------------------------
!
! 0.3 allocation
-CALL GET_DIM_EXT_ll('B',IIU,IJU)
-IKU=NKMAX+2*JPVEXT
PMINW=0
PMINN=0
PMINS=0
PMINE=0
-CALL GET_OR_ll('B',IIBG,IJBG)
-IIBG = IIBG+IIB-1
-IJBG = IJBG+IJB-1
-CALL GET_GLOBALDIMS_ll( IIMAX,IJMAX)
-IIEG=IIBG+IIE-IIB
-IJEG=IJBG+IJE-IJB
Z_DELTX = XXHAT(2)-XXHAT(1)
Z_DELTY = XYHAT(2)-XYHAT(1)
-
-ILUOUT = TLUOUT%NLU
!------------------------------------------------------------------------------
!
!**** 1. Recycling distance calculation
diff --git a/src/MNH/relaxation.f90 b/src/MNH/relaxation.f90
index 7202c8ea13c44d2b708271e20e44fd158703750d..69e130288298efa89da7afaccb9fe70fd88fb11d 100644
--- a/src/MNH/relaxation.f90
+++ b/src/MNH/relaxation.f90
@@ -273,6 +273,7 @@ USE MODD_CONF, only: cconf
USE MODD_ELEC_DESCR, ONLY: LRELAX2FW_ION
USE MODD_NSV, ONLY: NSV_ELECBEG, NSV_ELECEND
USE MODD_PARAMETERS, only: jphext, jpvext
+USE MODD_PRECISION, ONLY: MNHREAL32
use mode_budget, only: Budget_store_init, Budget_store_end
USE MODE_EXTRAPOL, only: Extrapol
@@ -517,16 +518,16 @@ IF(OVE_RELAX) THEN
!
DO JK = KALBOT, IKE+1
!
- PRUS(:,:,JK) = PRUS(:,:,JK) - ZKV(JK) *(PUT(:,:,JK) -PLSUM(:,:,JK) )&
+ PRUS(:,:,JK) = PRUS(:,:,JK) - ZKV(JK) *(PUT(:,:,JK) -real(PLSUM(:,:,JK),kind=MNHREAL32) )&
* ZRHODJU(:,:,JK)
!
- PRVS(:,:,JK) = PRVS(:,:,JK) - ZKV(JK) *(PVT(:,:,JK) -PLSVM(:,:,JK) )&
+ PRVS(:,:,JK) = PRVS(:,:,JK) - ZKV(JK) *(PVT(:,:,JK) -real(PLSVM(:,:,JK),kind=MNHREAL32) )&
* ZRHODJV(:,:,JK)
!
- PRWS(:,:,JK) = PRWS(:,:,JK) - ZKVW(JK) *(PWT(:,:,JK) -PLSWM(:,:,JK) )&
+ PRWS(:,:,JK) = PRWS(:,:,JK) - ZKVW(JK) *(PWT(:,:,JK) -real(PLSWM(:,:,JK),kind=MNHREAL32) )&
* ZRHODJW(:,:,JK)
!
- PRTHS(:,:,JK) = PRTHS(:,:,JK) - ZKV(JK) *(PTHT(:,:,JK) -PLSTHM(:,:,JK) )&
+ PRTHS(:,:,JK) = PRTHS(:,:,JK) - ZKV(JK) *(PTHT(:,:,JK) -real(PLSTHM(:,:,JK),kind=MNHREAL32) )&
* PRHODJ(:,:,JK)
!
END DO
@@ -554,16 +555,16 @@ IF(OVE_RELAX_GRD) THEN
!
DO JK = 1,KALBAS
!
- PRUS(:,:,JK) = PRUS(:,:,JK) - ZKVBAS(JK) *(PUT(:,:,JK) -PLSUM(:,:,JK) )&
+ PRUS(:,:,JK) = PRUS(:,:,JK) - ZKVBAS(JK) *(PUT(:,:,JK) -real(PLSUM(:,:,JK),kind=MNHREAL32) )&
* ZRHODJU(:,:,JK)
!
- PRVS(:,:,JK) = PRVS(:,:,JK) - ZKVBAS(JK) *(PVT(:,:,JK) -PLSVM(:,:,JK) )&
+ PRVS(:,:,JK) = PRVS(:,:,JK) - ZKVBAS(JK) *(PVT(:,:,JK) -real(PLSVM(:,:,JK),kind=MNHREAL32) )&
* ZRHODJV(:,:,JK)
!
- PRWS(:,:,JK) = PRWS(:,:,JK) - ZKVWBAS(JK) *(PWT(:,:,JK) -PLSWM(:,:,JK) )&
+ PRWS(:,:,JK) = PRWS(:,:,JK) - ZKVWBAS(JK) *(PWT(:,:,JK) -real(PLSWM(:,:,JK),kind=MNHREAL32) )&
* ZRHODJW(:,:,JK)
!
- PRTHS(:,:,JK) = PRTHS(:,:,JK) - ZKVBAS(JK) *(PTHT(:,:,JK) -PLSTHM(:,:,JK) )&
+ PRTHS(:,:,JK) = PRTHS(:,:,JK) - ZKVBAS(JK) *(PTHT(:,:,JK) -real(PLSTHM(:,:,JK),kind=MNHREAL32) )&
* PRHODJ(:,:,JK)
!
END DO
diff --git a/src/MNH/subl_blowsnow.f90 b/src/MNH/subl_blowsnow.f90
index ad0cf278fa86b1eeedaf7c0a2da1344d8d86024a..96c6b9be4432ef3ee07b4f7d79368c232ca113d7 100644
--- a/src/MNH/subl_blowsnow.f90
+++ b/src/MNH/subl_blowsnow.f90
@@ -696,6 +696,7 @@ ZFCRI = MAX(ZFCRI1,ZFCRI2)
!* 2 Calculate variances of the horizontal and vertical velocity components
!
ZS0 = ZFCRI*PZZ/PVMOD
+STOP 'Bug in TURB_FLUC: ZUSTAR used but not set'
ZSIGU = 4.77 *ZUSTAR**2/ (1+33*ZS0)**0.66666
ZSIGV = 2.76 *ZUSTAR**2/ (1+9.5*ZS0)**0.66666
ZSIGW = 1.31 *ZUSTAR**2/ (1+3.12*ZS0)**0.66666
diff --git a/src/MNH/version.f90 b/src/MNH/version.f90
index 8a7e340814d6a5fc2b25c3bc61afc8b5f477c999..5c98f1946a20fa439b85ac1ef532e771ecb90a4a 100644
--- a/src/MNH/version.f90
+++ b/src/MNH/version.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2002-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -44,9 +44,9 @@ IMPLICIT NONE
!
NMNHVERSION(1)=5
NMNHVERSION(2)=5
-NMNHVERSION(3)=0
+NMNHVERSION(3)=1
NMASDEV=55
-NBUGFIX=0
+NBUGFIX=1
CBIBUSER=''
!
END SUBROUTINE VERSION
diff --git a/src/MNH/write_desfmn.f90 b/src/MNH/write_desfmn.f90
index fb24c9bae47e3b7140e3a266768cffe9a5216841..8be686e6450155250e7d03c8712d67f659048a94 100644
--- a/src/MNH/write_desfmn.f90
+++ b/src/MNH/write_desfmn.f90
@@ -434,7 +434,6 @@ IF(LBU_RRH) WRITE(UNIT=ILUSEG,NML=NAM_BU_RRH)
IF(LBU_RSV) WRITE(UNIT=ILUSEG,NML=NAM_BU_RSV)
IF(LLES_MEAN .OR. LLES_RESOLVED .OR. LLES_SUBGRID .OR. LLES_UPDRAFT &
.OR. LLES_DOWNDRAFT .OR. LLES_SPECTRA) WRITE(UNIT=ILUSEG,NML=NAM_LES)
-WRITE(UNIT=ILUSEG,NML=NAM_BLANKn)
IF(LFORCING .OR. LTRANS) WRITE(UNIT=ILUSEG,NML=NAM_FRC)
IF(CCLOUD(1:3) == 'ICE') WRITE(UNIT=ILUSEG,NML=NAM_PARAM_ICE)
IF(CCLOUD == 'C2R2' .OR. CCLOUD == 'C3R5' .OR. CCLOUD == 'KHKO') &
diff --git a/src/MNH/write_diachro.f90 b/src/MNH/write_diachro.f90
index 77a681358e063314bd31dec2ed5cf979c3fbee75..365a24ff25b0ca25b73f07f971d2876cbf0d50f5 100644
--- a/src/MNH/write_diachro.f90
+++ b/src/MNH/write_diachro.f90
@@ -89,6 +89,7 @@ subroutine Write_diachro( tpdiafile, tpbudiachro, tpfields, &
! P. Wautelet 03/03/2021: add tbudiachrometadata type (useful to pass more information to Write_diachro)
! P. Wautelet 11/03/2021: remove ptrajx/y/z optional dummy arguments of Write_diachro
! + get the trajectory data for LFI files differently
+! P. Wautelet 01/09/2021: allow NMNHDIM_LEVEL and NMNHDIM_LEVEL_W simultaneously
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -1035,7 +1036,10 @@ do jp = 2, Size( tpfields )
if ( tpfields(jp)%ndimlist(ji) /= tpfields(1)%ndimlist(ji) ) then
!For SERIES: it is possible to have NMNHDIM_NI and NMNHDIM_NI_U in the different tpfields
!For SERIES: it is possible to have NMNHDIM_SERIES_LEVEL and NMNHDIM_SERIES_LEVEL_W in the different tpfields
- if ( tpfields(jp)%ndimlist(ji) /= NMNHDIM_NI .and. tpfields(jp)%ndimlist(ji) /= NMNHDIM_NI_U .and. &
+ !For profiles: it is possible to have NMNHDIM_LEVEL and NMNHDIM_LEVEL_W in the different tpfields
+ !This check is not perfect but should catch most problems
+ if ( tpfields(jp)%ndimlist(ji) /= NMNHDIM_NI .and. tpfields(jp)%ndimlist(ji) /= NMNHDIM_NI_U .and. &
+ tpfields(jp)%ndimlist(ji) /= NMNHDIM_LEVEL .and. tpfields(jp)%ndimlist(ji) /= NMNHDIM_LEVEL_W .and. &
tpfields(jp)%ndimlist(ji) /= NMNHDIM_SERIES_LEVEL .and. tpfields(jp)%ndimlist(ji) /= NMNHDIM_SERIES_LEVEL_W ) then
call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
'some dimensions are not the same for all tpfields entries for variable '//trim(tpfields(jp)%cmnhname) )
diff --git a/src/MNH/write_lbn.f90 b/src/MNH/write_lbn.f90
index a974878fe56b1eeaead8f1bdef34dc348ff2f958..af8e011ed01410b38464704b6bb4a796128e4f5b 100644
--- a/src/MNH/write_lbn.f90
+++ b/src/MNH/write_lbn.f90
@@ -381,7 +381,7 @@ IF (NSV >=1) THEN
DO JSV = NSV_LIMA_CCN_FREE,NSV_LIMA_CCN_FREE+NMOD_CCN-1
WRITE(INDICE,'(I2.2)')(JSV - NSV_LIMA_CCN_FREE + 1)
IF(NSIZELBXSV_ll /= 0) THEN
- TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3))//INDICE)
+ TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3)))//INDICE
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBX'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'2_Y_Z_','LBXSVM',JSV
@@ -389,7 +389,7 @@ IF (NSV >=1) THEN
END IF
!
IF(NSIZELBYSV_ll /= 0) THEN
- TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3))//INDICE)
+ TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_WARM_NAMES(3)))//INDICE
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBY'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'X_2_Z_','LBYSVM',JSV
@@ -400,7 +400,7 @@ IF (NSV >=1) THEN
DO JSV = NSV_LIMA_IFN_FREE,NSV_LIMA_IFN_FREE+NMOD_IFN-1
WRITE(INDICE,'(I2.2)')(JSV - NSV_LIMA_IFN_FREE + 1)
IF(NSIZELBXSV_ll /= 0) THEN
- TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2))//INDICE)
+ TZFIELD%CMNHNAME = 'LBX_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2)))//INDICE
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBX'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'2_Y_Z_','LBXSVM',JSV
@@ -408,7 +408,7 @@ IF (NSV >=1) THEN
END IF
!
IF(NSIZELBYSV_ll /= 0) THEN
- TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2))//INDICE)
+ TZFIELD%CMNHNAME = 'LBY_'//TRIM(UPCASE(CLIMA_COLD_NAMES(2)))//INDICE
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
TZFIELD%CLBTYPE = 'LBY'
WRITE(TZFIELD%CCOMMENT,'(A6,A6,I3.3)')'X_2_Z_','LBYSVM',JSV
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index c4fc3d50bfc71bd6de4f70d3d3674e6e9da7ea86..b2b31cf6f73aaf332df655a701499dbd80df676c 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -2689,9 +2689,10 @@ IF ( LMEAN_POVO ) THEN
IWORK1(:,:)=0
ZWORK21(:,:)=0.
IF (XMEAN_POVO(1)>XMEAN_POVO(2)) THEN
- XMEAN_POVO(1) = ZX0D
- XMEAN_POVO(2) = XMEAN_POVO(1)
- ZX0D = XMEAN_POVO(2)
+ !Invert values (smallest must be first)
+ ZX0D = XMEAN_POVO(1)
+ XMEAN_POVO(1) = XMEAN_POVO(2)
+ XMEAN_POVO(2) = ZX0D
END IF
DO JK=IKB,IKE
WHERE((XPABST(:,:,JK)>XMEAN_POVO(1)).AND.(XPABST(:,:,JK)<XMEAN_POVO(2)))
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index d397bec0bf30b6aaea185e0a96fe8ae3bb4c3838..d9fb11aa284dd60a5669711c58b87c4a7cd93e21 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -1552,6 +1552,11 @@ CALL WRITE_LB_n(TPFILE)
!
!
CALL IO_Field_write(TPFILE,'DRYMASST',XDRYMASST)
+IF (CPROGRAM == 'MESONH') THEN
+ CALL IO_Field_write(TPFILE,'DRYMASSS',XDRYMASSS)
+ELSE
+ CALL IO_Field_write(TPFILE,'DRYMASSS',0.)
+END IF
!
IF( CTURB /= 'NONE' .AND. CTOM=='TM06') THEN
CALL IO_Field_write(TPFILE,'BL_DEPTH',XBL_DEPTH)
diff --git a/src/MNH/write_profilern.f90 b/src/MNH/write_profilern.f90
index f3fbf5aa0363207f897d13df23a40bf4f24ed78c..8cabe0d3b0e0570e8c28c24e775368e2ec4233f1 100644
--- a/src/MNH/write_profilern.f90
+++ b/src/MNH/write_profilern.f90
@@ -17,6 +17,8 @@
! P. Wautelet 11/03/2021: bugfix: correct name for NSV_LIMA_IMM_NUCL
! P. Wautelet 05/07/2021: reorganisation to store point values correctly (not in vertical profiles)
! M. Taufour 07/2021: modify RARE for hydrometeors containing ice and add bright band calculation for RARE
+! P. Wautelet 01/09/2021: fix: correct vertical dimension for ALT and W
+! P. Wautelet 19/11/2021: bugfix in units for LIMA variables
!-----------------------------------------------------------------
! ###########################
MODULE MODE_WRITE_PROFILER_n
@@ -85,7 +87,7 @@ USE MODD_CH_AEROSOL, ONLY: CAERONAMES, LORILAM, JPMODE
USE MODD_CH_M9_n, ONLY: CNAMES
USE MODD_CST, ONLY: XRV
USE MODD_ELEC_DESCR, ONLY: CELECNAMES
-use modd_field, only: NMNHDIM_LEVEL, NMNHDIM_PROFILER_TIME, NMNHDIM_PROFILER_PROC, NMNHDIM_UNUSED, &
+use modd_field, only: NMNHDIM_LEVEL, NMNHDIM_LEVEL_W, NMNHDIM_PROFILER_TIME, NMNHDIM_PROFILER_PROC, NMNHDIM_UNUSED, &
tfield_metadata_base, TYPEREAL
USE MODD_ICE_C1R3_DESCR, ONLY: C1R3NAMES
USE MODD_IO, ONLY: TFILEDATA
@@ -114,10 +116,12 @@ INTEGER, INTENT(IN) :: KI
character(len=2) :: yidx
character(len=100) :: ycomment
character(len=100) :: yname
+character(len=40) :: yunits
CHARACTER(LEN=:), allocatable :: YGROUP ! group title
INTEGER :: IKU
INTEGER :: IPROC ! number of variables records
INTEGER :: JPROC
+integer :: jproc_alt, jproc_w
INTEGER :: JRR ! loop counter
INTEGER :: JSV ! loop counter
integer :: ji
@@ -161,11 +165,15 @@ call Add_profile( 'RARE', 'Radar reflectivity', 'dBZ', tprofil
call Add_profile( 'RAREatt', 'Radar attenuated reflectivity', 'dBZ', tprofiler%crare_att )
call Add_profile( 'P', 'Pressure', 'Pa', tprofiler%p )
call Add_profile( 'ALT', 'Altitude', 'm', tprofiler%zz )
+!Store position of ALT in the field list. Useful because it is not computed on the same Arakawa-grid points
+jproc_alt = jproc
call Add_profile( 'ZON_WIND', 'Zonal wind', 'm s-1', tprofiler%zon )
call Add_profile( 'MER_WIND', 'Meridional wind', 'm s-1', tprofiler%mer )
call Add_profile( 'FF', 'Wind intensity', 'm s-1', tprofiler%ff )
call Add_profile( 'DD', 'Wind direction', 'degree', tprofiler%dd )
call Add_profile( 'W', 'Air vertical speed', 'm s-1', tprofiler%w )
+!Store position of W in the field list. Useful because it is not computed on the same Arakawa-grid points
+jproc_w = jproc
if ( ldiag_in_run ) &
call Add_profile( 'TKE_DISS', 'TKE dissipation rate', 'm2 s-2', tprofiler% tke_diss )
@@ -207,6 +215,7 @@ if ( Size( tprofiler%sv, 4 ) > 0 ) then
end do
! LIMA variables
do jsv = nsv_lima_beg, nsv_lima_end
+ yunits = 'kg-1'
if ( jsv == nsv_lima_nc ) then
yname = Trim( clima_warm_names(1) ) // 'T'
else if ( jsv == nsv_lima_nr ) then
@@ -219,6 +228,7 @@ if ( Size( tprofiler%sv, 4 ) > 0 ) then
yname = Trim( clima_warm_names(4) ) // yidx // 'T'
else if ( jsv == nsv_lima_scavmass ) then
yname = Trim( caero_mass(1) ) // 'T'
+ yunits = 'kg kg-1'
else if ( jsv == nsv_lima_ni ) then
yname = Trim( clima_cold_names(1) ) // 'T'
else if ( jsv >= nsv_lima_ifn_free .and. jsv < nsv_lima_ifn_free + nmod_ifn ) then
@@ -235,7 +245,7 @@ if ( Size( tprofiler%sv, 4 ) > 0 ) then
else if ( jsv == nsv_lima_spro ) then
yname = Trim( clima_warm_names(5) ) // 'T'
end if
- call Add_profile( yname, '', 'kg-1', tprofiler%sv(:,:,:,jsv) )
+ call Add_profile( yname, '', yunits, tprofiler%sv(:,:,:,jsv) )
end do
! electrical scalar variables
do jsv = nsv_elecbeg, nsv_elecend
@@ -389,6 +399,8 @@ tzfields(:)%ndims = 3
tzfields(:)%ndimlist(1) = NMNHDIM_UNUSED
tzfields(:)%ndimlist(2) = NMNHDIM_UNUSED
tzfields(:)%ndimlist(3) = NMNHDIM_LEVEL
+tzfields(jproc_alt)%ndimlist(3) = NMNHDIM_LEVEL_W
+tzfields(jproc_w)%ndimlist(3) = NMNHDIM_LEVEL_W
tzfields(:)%ndimlist(4) = NMNHDIM_PROFILER_TIME
tzfields(:)%ndimlist(5) = NMNHDIM_UNUSED
tzfields(:)%ndimlist(6) = NMNHDIM_PROFILER_PROC
diff --git a/src/Makefile.MESONH.mk b/src/Makefile.MESONH.mk
index 5d019e2fe151007cb19c03d0776d646b52d7cb72..6ad235831de52676df58510c69e06282150f02c8 100644
--- a/src/Makefile.MESONH.mk
+++ b/src/Makefile.MESONH.mk
@@ -337,7 +337,7 @@ INC_MPI = -I$(B)$(DIR_MPI)
DIR_MASTER += $(DIR_MPI)
OBJS_LISTE_MASTER += mpivide.o
INC += $(INC_MPI)
-mpivide.o : CPPFLAGS += -DFUJI -DMNH_INT=$(MNH_INT) -DMNH_REAL=$(MNH_REAL) \
+mpivide.o : CPPFLAGS += -DMNH_INT=$(MNH_INT) -DMNH_REAL=$(MNH_REAL) \
-I$(DIR_MPI)/include
VPATH += $(DIR_MPI)
endif
@@ -648,8 +648,8 @@ endif
#
ifeq "$(VER_CDF)" "CDFCTI"
CDF_PATH?=/usr
-INC_NETCDF = -I${CDF_PATH}/include
-LIB_NETCDF = -L${CDF_PATH}/lib64 -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 -lsz -lz
+INC_NETCDF ?= $(shell nf-config --fflags)
+LIB_NETCDF ?= $(shell nf-config --flibs)
INC += $(INC_NETCDF)
LIBS += $(LIB_NETCDF)
endif
diff --git a/src/SURFEX/init_megann.F90 b/src/SURFEX/init_megann.F90
index ee3b097f8a6b1a5ee6374c4458225f0cd1d16d81..cfb05aecfc8ad0e6d332681ec41d1da73c983a61 100644
--- a/src/SURFEX/init_megann.F90
+++ b/src/SURFEX/init_megann.F90
@@ -26,6 +26,7 @@ SUBROUTINE INIT_MEGAN_n(IO, S, K, NP, MSF, MGN, PLAT, HSV, PMEGAN_FIELDS)
!! Original: 25/10/14
!! Modified: 06/2017, J. Pianezze, adaptation for SurfEx v8.0
!! Modified: 06/2018, P. Tulet, add PFT and LAI
+!! Modified: 11/2021, P. Tulet, update PFT with Ecoclimap-SG
!!
!!
!! EXTERNAL
@@ -40,10 +41,7 @@ USE MODD_MEGAN_n, ONLY : MEGAN_t
USE MODD_ISBA_OPTIONS_n, ONLY : ISBA_OPTIONS_t
USE MODD_ISBA_n, ONLY : ISBA_S_t, ISBA_P_t, ISBA_K_t, ISBA_NP_t
!
-USE MODD_DATA_COVER_PAR, ONLY : NVT_C4, NVT_TRBE, NVT_TRBD, NVT_TEBE, &
- NVT_TEBD, NVT_TENE, NVT_BOBD, NVT_BONE, NVT_BOND, &
- NVT_BOGR, NVT_SHRB, NVT_GRAS, NVT_TROG, NVT_C3, &
- NVT_NO, NVT_ROCK, NVT_SNOW, NVT_IRR, NVT_PARK
+USE MODD_DATA_COVER_PAR
!
USE MODD_SURF_PAR, ONLY : XUNDEF
USE MODD_DATA_COVER, ONLY : XDATA_LAI
@@ -94,6 +92,7 @@ REAL,DIMENSION(SIZE(K%XCLAY,1)) :: ZLAI
!
ALLOCATE(MGN%XPFT (N_MGN_PFT,SIZE(K%XCLAY,1)))
ALLOCATE(MGN%XEF (N_MGN_SPC,SIZE(K%XCLAY,1)))
+ALLOCATE(MGN%XLAI (SIZE(K%XCLAY,1)))
ALLOCATE(MGN%NSLTYP (SIZE(K%XCLAY,1)))
ALLOCATE(MGN%XBIOFLX(SIZE(K%XCLAY,1)))
MGN%XBIOFLX(:) = 0.
@@ -313,134 +312,91 @@ ENDDO
!
! 1 Needleleaf evergreen temperate tree
! -------------------------------------
-! utilisation de la classe NVT_CONI pour 30 < LAT < 60
-WHERE ((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.))
- MGN%XPFT(1,:) = S%XVEGTYPE(:,NVT_TENE)
-END WHERE
-WHERE ((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.))
- MGN%XPFT(1,:) = S%XVEGTYPE(:,NVT_TENE)
-END WHERE
+! utilisation de la classe NVT_TENE
+MGN%XPFT(1,:) = S%XVEGTYPE(:,NVT_TENE)
!
! 2 Needleleaf evergreen boreal tree
! -------------------------------------
-!utilisation de la classe NVT_CONI pour LAT > 60
-WHERE ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.))
- MGN%XPFT(2,:) = S%XVEGTYPE(:,NVT_BONE)
-END WHERE
+!utilisation de la classe NVT_BONE
+MGN%XPFT(2,:) = S%XVEGTYPE(:,NVT_BONE)
!
!3 Needleleaf deciduous boreal tree
! -------------------------------------
-!utilisation de la classe NVT_TREE pour LAT > 60
-WHERE ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.))
- MGN%XPFT(3,:) = S%XVEGTYPE(:,NVT_BOND)
-END WHERE
+!utilisation de la classe NVT_BOND
+MGN%XPFT(3,:) = S%XVEGTYPE(:,NVT_BOND)
!
!4 Broadleaf evergreen tropical tree
! -------------------------------------
-!utilisation de la classe NVT_EVER pour -30 < LAT < 30
-! et une hauteur d'arbre supérieur à 3 m
-WHERE (((PLAT(:) .GE. -30.) .AND. (PLAT(:) .LE. 30.)).AND.&
- (ZH_TREE(:,IP_TRBE) .GE. 3.).AND.(ZH_TREE(:,IP_TRBE) .NE. XUNDEF))
+!utilisation de la classe NVT_TRBE
MGN%XPFT(4,:) = S%XVEGTYPE(:,NVT_TRBE)
-END WHERE
!
!5 Broadleaf evergreen temperate tree
! -------------------------------------
-! utilisation de la classe NVT_EVER pour 30 < LAT < 60
-! et une hauteur d'arbre supérieur à 3 m.
-WHERE (((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)).AND.&
- (ZH_TREE(:,IP_TEBE) .GE. 3.).AND.(ZH_TREE(:,IP_TEBE) .NE. XUNDEF))
MGN%XPFT(5,:) = S%XVEGTYPE(:,NVT_TEBE)
-END WHERE
-WHERE (((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.)).AND.&
- (ZH_TREE(:,IP_TEBE) .GE. 3.).AND.(ZH_TREE(:,IP_TEBE) .NE. XUNDEF))
-MGN%XPFT(5,:) = S%XVEGTYPE(:,NVT_TEBE)
-END WHERE
!
!6 Broadleaf deciduous tropical tree
! -------------------------------------
-!utilisation de la classe NVT_TREE pour -30 < LAT < 30
-WHERE ((PLAT(:) .GE. -30.) .AND. (PLAT(:) .LE. 30.))
MGN%XPFT(6,:) = S%XVEGTYPE(:,NVT_TRBD)
-END WHERE
!
!7 Broadleaf deciduous temperate tree
! -------------------------------------
-!utilisation de la classe NVT_TREE pour 30 < LAT < 60
-! en utilisant une hauteur d'arbre supérieur à 3 m
-WHERE (((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)).AND.&
- (ZH_TREE(:,IP_TEBD) .GE. 3.).AND.(ZH_TREE(:,IP_TEBD) .NE. XUNDEF))
MGN%XPFT(7,:) = S%XVEGTYPE(:,NVT_TEBD)
-END WHERE
-WHERE (((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.)).AND.&
- (ZH_TREE(:,IP_TEBD) .GE. 3.).AND.(ZH_TREE(:,IP_TEBD) .NE. XUNDEF))
-MGN%XPFT(7,:) = S%XVEGTYPE(:,NVT_TEBD)
-END WHERE
!
!8 Broadleaf deciduous boreal tree
! -------------------------------------
-!utilisation de la classe NVT_TREE pour LAT > 60
-WHERE (((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.)).AND.&
- (ZH_TREE(:,IP_BOBD) .GE. 3.).AND.(ZH_TREE(:,IP_BOBD) .NE. XUNDEF))
MGN%XPFT(8,:) = S%XVEGTYPE(:,NVT_BOBD)
-END WHERE
!
!9 Broadleaf evergreen shrub
! -------------------------------------
-!utilisation de la classe NVT_EVER pour une hauteur d'arbre inférieure à 3 m
-WHERE (ZH_TREE(:,IP_SHRB) .LT. 3.)
+!utilisation de la classe NVT_SHBR pour -30 < LAT < 30
+WHERE (((PLAT(:) .GE. -30.) .AND. (PLAT(:) .LE. 30.)))
MGN%XPFT(9,:) = S%XVEGTYPE(:,NVT_SHRB)
+ELSE WHERE
+MGN%XPFT(9,:) = 0.
END WHERE
!
!10 Broadleaf deciduous temperate shrub
! -------------------------------------
-!utilisation de la classe NVT_TREE pour une hauteur d'arbre inférieure à 3 m
-! et pour 30 < LAT < 60
-WHERE ((ZH_TREE(:,IP_SHRB) .LT. 3.) .AND. (ZH_TREE(:,IP_SHRB).NE. XUNDEF) .AND. &
- ((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)))
-MGN%XPFT(10,:) = S%XVEGTYPE(:,NVT_SHRB)
-END WHERE
-WHERE ((ZH_TREE(:,IP_SHRB) .LT. 3.) .AND. (ZH_TREE(:,IP_SHRB).NE. XUNDEF) .AND. &
+!utilisation de la classe NVT_SHBR pour 30 < LAT < 60
+WHERE (((PLAT(:) .GE. 30.) .AND. (PLAT(:) .LT. 60.)).OR.&
((PLAT(:) .LE. -30.) .AND. (PLAT(:) .GT. -60.)))
MGN%XPFT(10,:) = S%XVEGTYPE(:,NVT_SHRB)
+ELSE WHERE
+MGN%XPFT(10,:) = 0.
END WHERE
!
!11 Broadleaf deciduous boreal_shrub
! -------------------------------------
-!utilisation de la classe NVT_TREE pour une hauteur d'arbre inférieure à 3 m
-! et pour LAT > 60
-WHERE ((ZH_TREE(:,IP_SHRB) .LT. 3.) .AND. (ZH_TREE(:,IP_SHRB).NE. XUNDEF) .AND. &
- ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.)))
+!utilisation de la classe NVT_SHBR pour LAT > 60
+WHERE (((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.)))
MGN%XPFT(11,:) = S%XVEGTYPE(:,NVT_SHRB)
+ELSE WHERE
+MGN%XPFT(11,:) = 0.
END WHERE
!
!12 C3 arctic grass
! -------------------------------------
-!utilisation de la classe NVT_GRAS + NVT_PARK pour LAT > 60
-WHERE ((PLAT(:) .GE. 60.) .OR. (PLAT(:) .LE. -60.))
-MGN%XPFT(12,:) = S%XVEGTYPE(:,NVT_GRAS) + S%XVEGTYPE(:,NVT_PARK)
-ELSEWHERE
+MGN%XPFT(12,:) = S%XVEGTYPE(:,NVT_BOGR)
!
!13 C3 non-arctic grass
! -------------------------------------
-!utilisation de la classe NVT_GRAS + NVT_PARK ailleur
-MGN%XPFT(13,:) = S%XVEGTYPE(:,NVT_GRAS) + S%XVEGTYPE(:,NVT_PARK)
-END WHERE
+MGN%XPFT(13,:) = S%XVEGTYPE(:,NVT_GRAS)
!
!14 C4 grass
! -------------------------------------
-! utilisation de la classe NVT_TROG
MGN%XPFT(14,:) = S%XVEGTYPE(:,NVT_TROG)
!
!15 Corn
! -------------------------------------
-! utilisation de la classe NVT_C4
MGN%XPFT(15,:) = S%XVEGTYPE(:,NVT_C4)
!
!16 Wheat
! -------------------------------------
-! utilisation de la classe NVT_C3
+IF (NVT_C3W .NE. 0 ) THEN ! use ecoclimap_sg
+MGN%XPFT(16,:) = S%XVEGTYPE(:,NVT_C3W) + S%XVEGTYPE(:,NVT_C3S)
+ELSE ! use ecaclimap2.0
MGN%XPFT(16,:) = S%XVEGTYPE(:,NVT_C3)
+END IF
!
! Emission factor
MGN%XEF(:,:) = 0.
@@ -499,7 +455,6 @@ DO JSV=1, MSF%NMEGAN_NBR
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "EFBIDER") MGN%XEF(18,:) = PMEGAN_FIELDS(:,JSV)
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "EFSTRESS") MGN%XEF(19,:) = PMEGAN_FIELDS(:,JSV)
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "EFOTHER") MGN%XEF(20,:) = PMEGAN_FIELDS(:,JSV)
-! IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "LAI") PLAI(:,1) = PMEGAN_FIELDS(:,JSV)
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT1") MGN%XPFT(1,:) = PMEGAN_FIELDS(:,JSV)
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT2") MGN%XPFT(2,:) = PMEGAN_FIELDS(:,JSV)
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT3") MGN%XPFT(3,:) = PMEGAN_FIELDS(:,JSV)
@@ -516,6 +471,7 @@ DO JSV=1, MSF%NMEGAN_NBR
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT14") MGN%XPFT(14,:) = PMEGAN_FIELDS(:,JSV)
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT15") MGN%XPFT(15,:) = PMEGAN_FIELDS(:,JSV)
IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "PFT16") MGN%XPFT(16,:) = PMEGAN_FIELDS(:,JSV)
+ IF (TRIM(MSF%CMEGAN_NAME(JSV)) == "LAI") MGN%XLAI(:) = PMEGAN_FIELDS(:,JSV)
END DO
#endif
diff --git a/src/configure b/src/configure
index c0845c0d5f5e52acd0d7cca1d31c9316e03c1c79..450f9c2ca241923e3eb3b4d7960aa232182e8ce1 100755
--- a/src/configure
+++ b/src/configure
@@ -347,7 +347,7 @@ export CC=gcc
export OPTLEVEL=${OPTLEVEL:-DEBUG}
export MVWORK=${MVWORK:-NO}
export VER_CDF=${VER_CDF:-CDFCTI}
- export NEED_NCARG=${NEED_NCARG:-NO}
+ export NEED_NCARG=${NEED_NCARG:-YES}
export NEED_TOOLS=NO
;;
'Linux nuwa'*)