From 5249fdd3e1299ca2689362a3b4492a0a608f8b6b Mon Sep 17 00:00:00 2001
From: Gaelle Tanguy <gaelle.tanguy@meteo.fr>
Date: Thu, 26 Nov 2015 15:05:51 +0000
Subject: [PATCH] C. Barthe and JP Pinty 11/2015: Update electricity
---
src/MNH/elec_fieldn.f90 | 29 +-
src/MNH/elec_tridz.f90 | 812 ++++++++++++++++++++++++++
src/MNH/flash_geom_elec.f90 | 981 ++++++++++++++++++++++++--------
src/MNH/gamma.f90 | 89 ++-
src/MNH/ini_elecn.f90 | 28 +-
src/MNH/ini_flash_geom_elec.f90 | 20 +-
src/MNH/ini_param_elec.f90 | 332 ++++++++++-
src/MNH/ini_segn.f90 | 5 +
src/MNH/ion_attach_elec.f90 | 69 +--
src/MNH/ion_drift.f90 | 219 ++-----
src/MNH/lesn.f90 | 61 ++
src/MNH/modd_elec_descr.f90 | 13 +-
src/MNH/modd_elec_flash.f90 | 46 ++
src/MNH/modd_elec_param.f90 | 25 +-
src/MNH/modd_elecn.f90 | 11 +-
src/MNH/modd_lma_simulator.f90 | 64 +++
src/MNH/modn_elec.f90 | 6 +-
src/MNH/read_desfmn.f90 | 14 +-
src/MNH/read_exsegn.f90 | 35 +-
src/MNH/resolved_elecn.f90 | 279 +++++----
src/MNH/spawn_field2.f90 | 7 +
src/MNH/write_lfin.f90 | 94 ++-
22 files changed, 2551 insertions(+), 688 deletions(-)
create mode 100644 src/MNH/elec_tridz.f90
create mode 100644 src/MNH/modd_elec_flash.f90
create mode 100644 src/MNH/modd_lma_simulator.f90
diff --git a/src/MNH/elec_fieldn.f90 b/src/MNH/elec_fieldn.f90
index 77b81ed9b..6b7aaf8f3 100644
--- a/src/MNH/elec_fieldn.f90
+++ b/src/MNH/elec_fieldn.f90
@@ -1,3 +1,4 @@
+
!MNH_LIC Copyright 1994-2014 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
@@ -43,12 +44,15 @@ USE MODD_CST, ONLY : XCPD
USE MODD_DYN_n, ONLY: XTRIGSX, XTRIGSY, XDXHATM, XDYHATM, NIFAXX, NIFAXY
USE MODD_METRICS_n
USE MODD_ELEC_DESCR, ONLY : CLSOL, NLAPITR_ELEC, XEPSILON, XEPOTFW_TOP
-USE MODD_ELEC_n, ONLY : XESOURCEFW, XRHOM_E, XAF_E, XCF_E, XBFY_E
+USE MODD_ELEC_n, ONLY : XRHOM_E, XAF_E, XCF_E, XBFY_E, &
+ XBFB_E, XBF_SXP2_YP1_Z_E
!
USE MODI_FLAT_INV
+USE MODI_FLAT_INVZ
USE MODI_RICHARDSON
USE MODI_CONJGRAD
USE MODI_CONRESOL
+USE MODI_CONRESOLZ
USE MODI_GRADIENT_M
!
USE MODD_ARGSLIST_ll, ONLY : LIST_ll
@@ -164,10 +168,18 @@ CALL CLEANLIST_ll(TZFIELDS_ll)
!
IF (LFLAT .AND. LCARTESIAN) THEN
! flat cartesian LHE case -> exact solution
- CALL FLAT_INV (ZLBCX, ZLBCY, XDXHATM, XDYHATM, &
- XRHOM_E, XAF_E, XBFY_E, XCF_E, &
- XTRIGSX, XTRIGSY, NIFAXX, NIFAXY, &
- ZDV_SOURCE, ZPHIT)
+ IF ( CLSOL /= "ZRESI" ) THEN
+ CALL FLAT_INV (ZLBCX, ZLBCY, XDXHATM, XDYHATM, &
+ XRHOM_E, XAF_E, XBFY_E, XCF_E, &
+ XTRIGSX, XTRIGSY, NIFAXX, NIFAXY, &
+ ZDV_SOURCE, ZPHIT)
+ ELSE
+ CALL FLAT_INVZ (ZLBCX, ZLBCY, XDXHATM, XDYHATM, &
+ XRHOM_E, XAF_E, XBFY_E, XCF_E, &
+ XTRIGSX, XTRIGSY, NIFAXX, NIFAXY, &
+ ZDV_SOURCE, ZPHIT, &
+ XBFB_E, XBF_SXP2_YP1_Z_E)
+ END IF
ELSE
SELECT CASE(CLSOL)
CASE('RICHA') ! Richardson's method
@@ -189,6 +201,13 @@ ELSE
XRHOM_E, XAF_E, XBFY_E, XCF_E, XTRIGSX, XTRIGSY, &
NIFAXX, NIFAXY, NLAPITR_ELEC, &
ZDV_SOURCE, ZPHIT)
+ CASE('ZRESI') ! Conjugate Residual method
+ CALL CONRESOLZ (ZLBCX, ZLBCY, XDXX, XDYY, XDZX, XDZY, XDZZ,&
+ PRHODJ,ZTHETAV, XDXHATM, XDYHATM, &
+ XRHOM_E, XAF_E, XBFY_E, XCF_E, XTRIGSX, XTRIGSY, &
+ NIFAXX, NIFAXY, NLAPITR_ELEC, &
+ ZDV_SOURCE, ZPHIT, &
+ XBFB_E, XBF_SXP2_YP1_Z_E)
END SELECT
END IF
!
diff --git a/src/MNH/elec_tridz.f90 b/src/MNH/elec_tridz.f90
new file mode 100644
index 000000000..524765eae
--- /dev/null
+++ b/src/MNH/elec_tridz.f90
@@ -0,0 +1,812 @@
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$
+! MASDEV4_7 init 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+! ######################
+ MODULE MODI_ELEC_TRIDZ
+! ######################
+!
+INTERFACE
+!
+ SUBROUTINE ELEC_TRIDZ(HLUOUT,HLBCX,HLBCY, &
+ PMAP,PDXHAT,PDYHAT,PDXHATM,PDYHATM,PRHOM, &
+ PAF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ PRHODJ,PTHVREF,PZZ,PBFY,PEPOTFW_TOP, &
+ PBFB,PBF_SXP2_YP1_Z)!++cb - Z_SPLITTING
+!
+IMPLICIT NONE
+!
+CHARACTER (LEN=*), INTENT(IN) :: HLUOUT ! output listing name
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Potential
+ ! Temperature of the reference state
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! scale factor
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height z
+!
+REAL, DIMENSION(:), INTENT(IN) :: PDXHAT ! Stretching in x direction
+REAL, DIMENSION(:), INTENT(IN) :: PDYHAT ! Stretching in y direction
+!
+REAL, INTENT(OUT) :: PDXHATM ! mean grid increment in the x
+ ! direction
+REAL, INTENT(OUT) :: PDYHATM ! mean grid increment in the y
+ ! direction
+!
+REAL, DIMENSION (:), INTENT(OUT) :: PRHOM ! mean of XRHODJ on the plane
+ ! x y localized at a mass
+ ! level
+!
+REAL, DIMENSION(:), INTENT(OUT) :: PAF,PCF ! vectors giving the nonvanishing
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFY ! elements (yslice) of the tri-diag.
+ ! matrix in the pressure eq.
+!++cb - Z_SPLITTING
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements (bsplit slide) of the tri-diag.
+ ! matrix in the pressure eq.
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
+ ! matrix in the pressure eq
+!++cb
+!
+REAL, DIMENSION(:,:), INTENT(INOUT) :: PEPOTFW_TOP ! top boundary condition of
+ ! the tri-diag. system; it
+ ! has the dimension of an
+ ! electrical potential
+!
+ ! arrays of sin or cos values
+ ! for the FFT :
+REAL, DIMENSION(:), INTENT(OUT) :: PTRIGSX ! - along x
+REAL, DIMENSION(:), INTENT(OUT) :: PTRIGSY ! - along y
+!
+ ! decomposition in prime
+ ! numbers for the FFT:
+INTEGER, DIMENSION(19), INTENT(OUT) :: KIFAXX ! - along x
+INTEGER, DIMENSION(19), INTENT(OUT) :: KIFAXY ! - along y
+
+!
+END SUBROUTINE ELEC_TRIDZ
+!
+END INTERFACE
+!
+END MODULE MODI_ELEC_TRIDZ
+!
+! ###################################################################
+ SUBROUTINE ELEC_TRIDZ(HLUOUT,HLBCX,HLBCY, &
+ PMAP,PDXHAT,PDYHAT,PDXHATM,PDYHATM,PRHOM, &
+ PAF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ PRHODJ,PTHVREF,PZZ,PBFY,PEPOTFW_TOP, &
+ PBFB,PBF_SXP2_YP1_Z)
+! ####################################################################
+!
+!!**** *ELEC_TRIDZ * - Compute coefficients for the flat operator to get the
+!! electric potential from the Gauss equation
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the vertical time independent
+! coefficients a(k), b(k), c(k) required for the calculation of the "flat"
+! (i.e. neglecting the orography) operator Laplacian. RHOJ is averaged on
+! the whole horizontal domain. The form of the eigenvalues of the flat
+! operator depends on the lateral boundary conditions. Furthermore, this
+! routine initializes TRIGS and IFAX arrays required for the FFT transform
+! used in the routine PRECOND.
+! ELEC_TRIDZ (to invert the Gauss equation) differs from TRID (to solve the
+! pressure equation) by the bottom boundary condition, here Dirichlet
+! instead of Neumann, because the earth surface is conductive so it is a
+! surface with an electrical equipotential referenced to zero.
+!
+!!** METHOD
+!! ------
+!! The forms of the eigenvalues of the horizontal Laplacian are given by:
+!! Cyclic conditions:
+!! -----------------
+!! <rhoj> 2 ( pi ) <rhoj> 2 ( pi )
+!! b(m,n) = -4 ----------- sin (----- m ) -4 ----------- sin (----- n )
+!! <dxx> <dxx> ( imax ) <dyy> <dyy> ( jmax )
+!!
+!! Open conditions:
+!! -----------------
+!! <rhoj> 2 ( pi ) <rhoj> 2 ( pi )
+!! b(m,n) = -4 ----------- sin (----- m ) -4 ----------- sin (----- n )
+!! <dxx> <dxx> ( 2imax ) <dyy> <dyy> ( 2jmax )
+!!
+!! Cyclic condition along x and open condition along y:
+!! ------------------------------------------------------
+!! <rhoj> 2 ( pi ) <rhoj> 2 ( pi )
+!! b(m,n) = -4 ----------- sin (----- m ) -4 ----------- sin (----- n )
+!! <dxx> <dxx> ( imax ) <dyy> <dyy> ( 2jmax )
+!!
+!! Open condition along x and cyclic condition along y:
+!! ------------------------------------------------------
+!! <rhoj> 2 ( pi ) <rhoj> 2 ( pi )
+!! b(m,n) = -4 ----------- sin (----- m ) -4 ----------- sin (----- n )
+!! <dxx> <dxx> ( 2imax ) <dyy> <dyy> ( jmax )
+!!
+!! where m = 0,1,2....imax-1, n = 0,1,2....jmax-1
+!! Note that rhoj contains the Jacobian J = Deltax*Deltay*Deltaz = volume of
+!! an elementary mesh.
+
+!!
+!! EXTERNAL
+!! --------
+!! Function FFTFAX: initialization of TRIGSX,IFAXX,TRIGSY,IFAXY for
+!! the FFT transform
+!! GET_DIM_EXT_ll : get extended sub-domain sizes
+!! GET_INDICE_ll : get physical sub-domain bounds
+!! GET_DIM_PHYS_ll : get physical sub-domain sizes
+!! GET_GLOBALDIMS_ll : get physical global domain sizes
+!! GET_OR_ll : get origine coordonates of the physical sub-domain
+!! in global indices
+!! REDUCESUM_ll : sum into a scalar variable
+!! GET_SLICE_ll : get a slice of the global domain
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : define constants
+!! XPI : pi
+!! XCPD
+!! Module MODD_PARAMETERS: declaration of parameter variables
+!! JPHEXT, JPVEXT: define the number of marginal points out of the
+!! physical domain along horizontal and vertical directions respectively
+!! Module MODD_CONF: model configurations
+!! LCARTESIAN: logical for CARTESIAN geometry
+!! .TRUE. = Cartesian geometry used
+!! L2D: logical for 2D model version
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (routine ELEC_TRIDZ)
+!!
+!! AUTHOR
+!! ------
+!! P. HÃ…reil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/07/94
+!! 14/04/95 (J. Stein) bug in the ZDZM computation
+!! ( stretched case)
+!! 8/07/96 (P. Jabouille) change the FFT initialization
+!! which now works for odd number.
+!! 14/01/97 Durran anelastic equation (Stein,Lafore)
+!! 15/06/98 (D.Lugato, R.Guivarch) Parallelisation
+!! 10/08/98 (N. Asencio) add parallel code
+!! use PDXHAT, PDYHAT and not PXHAT,PYHAT
+!! PBFY is initialized
+!! 20/08/00 (J. Stein, J. Escobar) optimisation of the solver
+!! PBFY transposition
+!! 14/03/02 (P. Jabouille) set values for meaningless spectral coefficients
+!! (to avoid problem in bouissinesq configuration)
+!! 01/07/12 (J-P. Pinty) add a non-homogeneous fair-weather
+!! top boundary condition (Neuman)
+!! 16/10/14 (C. Bovalo) add the z-splitting
+!------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODD_CST
+USE MODD_CONF
+USE MODD_PARAMETERS
+!
+USE MODE_ll
+USE MODE_IO_ll
+USE MODE_FM
+!++cb - Z_SPLITTING
+USE MODE_SPLITTINGZ_ll , ONLY : GET_DIM_EXTZ_ll,GET_ORZ_ll,LWESTZ_ll,LSOUTHZ_ll
+!
+USE MODE_REPRO_SUM
+!++cb -
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments
+!
+!
+!
+CHARACTER (LEN=*), INTENT(IN) :: HLUOUT ! output listing name
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Potential
+ ! Temperature of the reference state
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PMAP ! scale factor
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height z
+!
+REAL, DIMENSION(:), INTENT(IN) :: PDXHAT ! Stretching in x direction
+REAL, DIMENSION(:), INTENT(IN) :: PDYHAT ! Stretching in y direction
+!
+REAL, INTENT(OUT) :: PDXHATM ! mean grid increment in the x
+ ! direction
+REAL, INTENT(OUT) :: PDYHATM ! mean grid increment in the y
+ ! direction
+!
+REAL, DIMENSION (:), INTENT(OUT) :: PRHOM ! mean of XRHODJ on the plane
+ ! x y localized at a mass
+ ! level
+!
+REAL, DIMENSION(:), INTENT(OUT) :: PAF,PCF ! vectors giving the nonvanishing
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFY ! elements (yslice) of the tri-diag.
+! matrix in the pressure eq. which is transposed. PBFY is a y-slices structure
+!++cb - Z_SPLITTING
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements (bsplit slide) of the tri-diag.
+ ! matrix in the pressure eq.
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
+ ! matrix in the pressure eq.
+!++cb
+!
+REAL, DIMENSION(:,:), INTENT(INOUT) :: PEPOTFW_TOP ! top boundary condition of
+ ! the tri-diag. system; it
+ ! has the dimension of an
+ ! electrical potential
+!
+ ! arrays of sin or cos values
+ ! for the FFT :
+REAL, DIMENSION(:), INTENT(OUT) :: PTRIGSX ! - along x
+REAL, DIMENSION(:), INTENT(OUT) :: PTRIGSY ! - along y
+!
+ ! decomposition in prime
+ ! numbers for the FFT:
+INTEGER, DIMENSION(19), INTENT(OUT) :: KIFAXX ! - along x
+INTEGER, DIMENSION(19), INTENT(OUT) :: KIFAXY ! - along y
+
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: IRESP ! FM return code
+INTEGER :: ILUOUT ! Logical unit number for
+ ! output-listing
+INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE ! indice values of the physical subdomain
+INTEGER :: IKU ! size of the arrays along z
+INTEGER :: IIB_ll,IIE_ll,IJB_ll,IJE_ll ! indice values of the physical global domain
+INTEGER :: IIMAX,IJMAX ! Number of points of the physical subdomain
+INTEGER :: IIMAX_ll,IJMAX_ll ! Number of points of Global physical domain
+!
+INTEGER :: JI,JJ,JK ! loop indexes
+!
+INTEGER :: INN ! temporary result for the computation of array TRIGS
+!
+REAL, DIMENSION (:,:), ALLOCATABLE :: ZEIGEN_ll ! eigenvalues b(m,n) in global representation
+REAL, DIMENSION (:), ALLOCATABLE :: ZEIGENX_ll ! used for the computation of ZEIGEN_ll
+!
+REAL, DIMENSION( SIZE(PDXHAT)) :: ZWORKX ! work array to compute PDXHATM
+REAL, DIMENSION( SIZE(PDYHAT)) :: ZWORKY ! work array to compute PDYHATM
+!
+REAL :: ZGWNX,ZGWNY ! greater wave numbers allowed by the model
+ ! configuration in x and y directions respectively
+!
+REAL, DIMENSION (SIZE(PZZ,3)) :: ZDZM ! mean of deltaz on the plane x y
+ ! localized at a w-level
+!
+REAL :: ZANGLE,ZDEL ! needed for the initialization of the arrays used by the FFT
+!
+REAL :: ZINVMEAN ! inverse of inner points number in an horizontal grid
+!
+INTEGER :: IINFO_ll ! return code of parallel routine
+REAL, DIMENSION (SIZE(PMAP,1)) :: ZXMAP ! extraction of PMAP array along x
+REAL, DIMENSION (SIZE(PMAP,2)) :: ZYMAP ! extraction of PMAP array along y
+INTEGER :: IORXY_ll,IORYY_ll ! origin's coordinates of the y-slices subdomain
+INTEGER :: IIUY_ll,IJUY_ll ! dimensions of the y-slices subdomain
+INTEGER :: IXMODE_ll,IYMODE_ll ! number of modes in the x and y direction for global point of view
+INTEGER :: IXMODEY_ll,IYMODEY_ll ! number of modes in the x and y direction for y_slice point of view
+!++cb - Z_SPLITTING
+INTEGER :: IORXB_ll,IORYB_ll ! origin's coordinates of the b-slices subdomain
+INTEGER :: IIUB_ll,IJUB_ll ! dimensions of the b-slices subdomain
+INTEGER :: IXMODEB_ll,IYMODEB_ll ! number of modes in the x and y direction for b_slice point of view
+!
+INTEGER :: IORX_SXP2_YP1_Z_ll,IORY_SXP2_YP1_Z_ll ! origin's coordinates of the b-slices subdomain
+INTEGER :: IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll ! dimensions of the b-slices subdomain
+INTEGER :: IXMODE_SXP2_YP1_Z_ll,IYMODE_SXP2_YP1_Z_ll ! number of modes in the x and y direction for b_slice point of view
+!++cb
+!
+!++JUAN-16
+!TYPE(DOUBLE_DOUBLE) , DIMENSION (SIZE(PZZ,3)) :: ZRHOM_ll , ZDZM_ll
+REAL, ALLOCATABLE, DIMENSION(:,:) :: ZRHOM_2D , ZDZM_2D
+!++JUAN-16
+!
+!
+!
+!
+!
+!------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATION
+! --------------
+!
+!* 1.1 retrieve a logical unit number
+! ------------------------------
+!
+CALL FMLOOK_ll(HLUOUT,HLUOUT,ILUOUT,IRESP)
+!
+!* 1.2 compute loop bounds
+! -------------------
+!
+! extended sub-domain
+CALL GET_DIM_EXT_ll ('Y',IIUY_ll,IJUY_ll)
+!++cb - Z_SPLITTING
+CALL GET_DIM_EXT_ll ('B',IIUB_ll,IJUB_ll)
+! P1/P2 splitting
+CALL GET_DIM_EXTZ_ll('SXP2_YP1_Z',IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll)
+!++cb
+IKU=SIZE(PRHODJ,3)
+!
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IKB=1 +JPVEXT
+IKE=IKU -JPVEXT
+! physical sub-domain
+CALL GET_DIM_PHYS_ll ( 'B',IIMAX,IJMAX)
+!
+! global physical domain limits
+CALL GET_GLOBALDIMS_ll ( IIMAX_ll, IJMAX_ll)
+IIB_ll = 1 + JPHEXT
+IIE_ll = IIMAX_ll + JPHEXT
+IJB_ll = 1 + JPHEXT
+IJE_ll = IJMAX_ll + JPHEXT
+!
+! the use of local array ZEIGENX and ZEIGEN would require some technical modifications
+!
+ALLOCATE (ZEIGENX_ll(IIMAX_ll + 2*JPHEXT))
+ALLOCATE (ZEIGEN_ll(IIMAX_ll + 2*JPHEXT, IJMAX_ll + 2*JPHEXT))
+ZEIGEN_ll = 0.0
+! Get the origin coordinates of the extended sub-domain in global landmarks
+CALL GET_OR_ll('Y',IORXY_ll,IORYY_ll)
+!++cb - Z_SPLITTING
+CALL GET_OR_ll('B',IORXB_ll,IORYB_ll)
+! P1/P2 Splitting
+CALL GET_ORZ_ll('SXP2_YP1_Z',IORX_SXP2_YP1_Z_ll,IORY_SXP2_YP1_Z_ll)
+!++cb
+!
+!* 1.3 allocate x-slice array
+
+!
+!* 1.4 variables for the eigenvalues computation
+!
+ZGWNX = XPI/REAL(IIMAX_ll)
+ZGWNY = XPI/REAL(IJMAX_ll)
+!
+!------------------------------------------------------------------------------
+!
+!* 2. COMPUTE THE AVERAGE OF RHOJ*CPD*THETAVREF ALONG XY
+! --------------------------------------------------
+!
+ZINVMEAN = 1./REAL(IIMAX_ll*IJMAX_ll)
+!JUAN-16
+ALLOCATE(ZRHOM_2D(IIB:IIE, IJB:IJE))
+!
+DO JK = 1,SIZE(PZZ,3)
+ IF ( CEQNSYS == 'DUR' .OR. CEQNSYS == 'MAE' ) THEN
+ DO JJ = IJB,IJE
+ DO JI = IIB,IIE
+ ZRHOM_2D(JI,JJ) = PRHODJ(JI,JJ,JK)*XCPD*PTHVREF(JI,JJ,JK)*ZINVMEAN
+ END DO
+ END DO
+ ELSEIF ( CEQNSYS == 'LHE' ) THEN
+ DO JJ = IJB,IJE
+ DO JI = IIB,IIE
+ ZRHOM_2D(JI,JJ) = PRHODJ(JI,JJ,JK)*ZINVMEAN
+ END DO
+ END DO
+ END IF
+ ! global sum
+ PRHOM(JK) = SUM_DD_R2_ll(ZRHOM_2D)
+END DO
+
+!
+! global sum
+!CALL REDUCESUM_ll(ZRHOM_ll,IINFO_ll)
+!PRHOM = ZRHOM_ll
+!JUAN-16
+!
+!------------------------------------------------------------------------------
+!
+!* 3. COMPUTE THE MEAN INCREMENT BETWEEN Z LEVELS
+! -------------------------------------------
+!
+!JUAN-16
+!ZDZM_ll = 0.
+ALLOCATE(ZDZM_2D(IIB:IIE, IJB:IJE))
+!
+DO JK = IKB-1,IKE
+ DO JJ = IJB,IJE
+ DO JI = IIB,IIE
+ ZDZM_2D(JI,JJ) = (PZZ(JI,JJ,JK+1)-PZZ(JI,JJ,JK))*ZINVMEAN
+ END DO
+ END DO
+ ZDZM(JK) = SUM_DD_R2_ll(ZDZM_2D)
+END DO
+ZDZM(IKE+1) = ZDZM(IKE)
+!
+! global sum
+!CALL REDUCESUM_ll(ZDZM_ll,IINFO_ll)
+!ZDZM = ZDZM_ll
+!JUAN-16
+!
+!
+! vertical average to arrive at a w-level
+DO JK = IKE+1,IKB,-1
+ ZDZM(JK) = (ZDZM(JK) + ZDZM(JK-1))*0.5
+END DO
+!
+ZDZM(IKB-1) = -999.
+!
+!------------------------------------------------------------------------------
+!
+!* 4. COMPUTE THE MEAN INCREMENT BETWEEN X LEVELS
+! -------------------------------------------
+!
+PDXHATM =0.
+! . local sum
+IF (LCARTESIAN) THEN
+ PDXHATM = SUM_1DFIELD_ll ( PDXHAT,'X',IIB_ll,IIE_ll,IINFO_ll)
+ELSE
+ ! Extraction of x-slice PMAP at j=(IJB_ll+IJE_ll)/2
+ CALL GET_SLICE_ll (PMAP,'X',(IJB_ll+IJE_ll)/2,ZXMAP(IIB:IIE) &
+ ,IIB,IIE,IINFO_ll)
+ ! initialize the work array = PDXHAT/ZXMAP
+ ZWORKX(IIB:IIE) = PDXHAT(IIB:IIE)/ ZXMAP (IIB:IIE)
+ PDXHATM = SUM_1DFIELD_ll ( ZWORKX,'X',IIB_ll,IIE_ll,IINFO_ll)
+END IF
+! . division to complete sum
+PDXHATM = PDXHATM / REAL(IIMAX_ll)
+!
+!------------------------------------------------------------------------------
+!
+!* 5. COMPUTE THE MEAN INCREMENT BETWEEN Y LEVELS
+! -------------------------------------------
+!
+PDYHATM = 0.
+IF (LCARTESIAN) THEN
+ PDYHATM = SUM_1DFIELD_ll ( PDYHAT,'Y',IJB_ll,IJE_ll,IINFO_ll)
+ELSE
+ ! Extraction of y-slice PMAP at i=IIB_ll+IIE_ll/2
+ CALL GET_SLICE_ll (PMAP,'Y',(IIB_ll+IIE_ll)/2,ZYMAP(IJB:IJE) &
+ ,IJB,IJE,IINFO_ll)
+ ! initialize the work array = PDYHAT / ZYMAP
+ ZWORKY(IJB:IJE) = PDYHAT(IJB:IJE) / ZYMAP (IJB:IJE)
+ PDYHATM = SUM_1DFIELD_ll ( ZWORKY,'Y',IJB_ll,IJE_ll,IINFO_ll)
+END IF
+! . division to complete sum
+PDYHATM= PDYHATM / REAL(IJMAX_ll)
+!
+!------------------------------------------------------------------------------
+!
+!* 6. COMPUTE THE OUT-DIAGONAL ELEMENTS A AND C OF THE MATRIX
+! -------------------------------------------------------
+!
+DO JK = IKB,IKE
+ PAF(JK) = 0.5 * ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2
+ PCF(JK) = 0.5 * ( PRHOM(JK) + PRHOM(JK+1) ) / ZDZM(JK+1) **2
+END DO
+!
+! at the upper level PAF and PCF are computed using the Neumann
+! condition applying on the vertical component of the gradient
+! while at the lower level, the Dirichlet condition is used
+!
+!
+! Neumann boundary condition (top of atmosphere)
+!
+PAF(IKE+1) = -0.5 * ( PRHOM(IKE) + PRHOM(IKE+1) ) / ZDZM(IKE+1) **2
+!
+! Dirichlet boundary condition (earth surface)
+!
+PCF(IKB-1) = 0.0
+!
+PAF(IKB-1) = 999.
+PCF(IKE+1) = 999.
+!
+!------------------------------------------------------------------------------
+!* 7. COMPUTE THE DIAGONAL ELEMENTS B OF THE MATRIX
+! ---------------------------------------------
+!
+!* 7.1 compute the horizontal eigenvalues
+!
+!
+!* 7.1.1 compute the eigenvalues along the x direction
+!
+SELECT CASE (HLBCX(1))
+! in the cyclic case, the eigenvalues are the same for two following JM values:
+! it corresponds to the real and complex parts of the FFT
+ CASE('CYCL') ! cyclic case
+ IXMODE_ll = IIMAX_ll+2
+ IXMODEY_ll = IIUY_ll
+ IXMODEB_ll = IIUB_ll !++cb - Z_SPLITTING
+!
+ DO JI = 1,IXMODE_ll
+ ZEIGENX_ll(JI) = - ( 2. * SIN ( (JI-1)/2*ZGWNX ) / PDXHATM )**2
+ END DO
+ CASE DEFAULT ! other cases
+ IXMODE_ll = IIMAX_ll
+!
+!
+ IF (LEAST_ll(HSPLITTING='Y')) THEN
+ IXMODEY_ll = IIUY_ll - 2
+ ELSE
+ IXMODEY_ll = IIUY_ll
+ END IF
+!++cb - Z_SPLITTING
+ IF (LEAST_ll(HSPLITTING='B')) THEN
+ IXMODEB_ll = IIUB_ll - 2
+ ELSE
+ IXMODEB_ll = IIUB_ll
+ END IF
+!++cb
+!
+!
+ DO JI = 1,IXMODE_ll
+ ZEIGENX_ll(JI) = - ( 2. *SIN (0.5*REAL(JI-1)*ZGWNX ) / PDXHATM )**2
+ END DO
+END SELECT
+!
+!* 7.1.2 compute the eventual eigenvalues along the y direction
+!
+IF (.NOT. L2D) THEN
+!
+! y lateral boundary conditions for three-dimensional cases
+!
+ SELECT CASE (HLBCY(1))
+! in the cyclic case, the eigenvalues are the same for two following JN values:
+! it corresponds to the real and complex parts of the FFT result
+!
+ CASE('CYCL') ! 3D cyclic case
+ IYMODE_ll = IJMAX_ll+2
+ IYMODEY_ll = IJUY_ll
+ IYMODEB_ll = IJUB_ll !++cb - Z_SPLITTING
+!
+ DO JJ = 1,IYMODE_ll
+ DO JI = 1,IXMODE_ll
+ ZEIGEN_ll(JI,JJ) = ZEIGENX_ll(JI) - &
+ ( 2.* SIN ( (JJ-1)/2*ZGWNY ) / PDYHATM )**2
+ END DO
+ END DO
+!
+ CASE DEFAULT ! 3D non-cyclic cases
+ IYMODE_ll = IJMAX_ll
+ IYMODEY_ll = IJUY_ll - 2
+ IYMODEB_ll = IJUB_ll - 2 !++cb - Z_SPLITTING
+!
+ DO JJ = 1,IYMODE_ll
+ DO JI = 1,IXMODE_ll
+ ZEIGEN_ll(JI,JJ) = ZEIGENX_ll(JI) - ( 2.* SIN (0.5*REAL(JJ-1)*ZGWNY ) / &
+ PDYHATM )**2
+ END DO
+ END DO
+!
+ END SELECT
+ELSE
+!
+! copy the x eigenvalue array in a 2D array for a 2D case
+!
+ IYMODE_ll = 1
+ IYMODEY_ll = 1
+ ZEIGEN_ll(1:IXMODE_ll,1)=ZEIGENX_ll(1:IXMODE_ll)
+!
+END IF
+!
+DEALLOCATE(ZEIGENX_ll)
+!
+!* 7.2 compute the matrix diagonal elements
+!
+!
+PBFY = 1.
+PBFB = 1. ! ++cb - Z_SLIDE
+PBF_SXP2_YP1_Z = 1. ! ++cb - Z_SLIDE
+!++cb
+IF (L2D) THEN
+ DO JK= IKB,IKE
+ DO JJ= 1, IYMODEY_ll
+ DO JI= 1, IXMODEY_ll
+ PBFY(JI,JJ,JK) = PRHOM(JK)* ZEIGEN_ll(JI+IORXY_ll-1,JJ+IORYY_ll-1) - 0.5 * &
+ ( ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2 &
+ +( PRHOM(JK) + PRHOM(JK+1) ) / ZDZM(JK+1)**2 )
+ END DO
+ END DO
+ END DO
+! at the upper level PBFY is computed using the Neumann condition
+! while at the lower level, the Dirichlet condition is used
+!
+! lower level: Dirichlet condition
+ PBFY(1:IXMODEY_ll,1:IYMODEY_ll,IKB) = PBFY(1:IXMODEY_ll,1:IYMODEY_ll,IKB) - &
+ PAF(IKB)
+ PAF(IKB) = 0.0
+ PBFY(1:IXMODEY_ll,1:IYMODEY_ll,IKB-1) = 1.0
+ !
+! upper level: Neumann condition
+ PBFY(1:IXMODEY_ll,1:IYMODEY_ll,IKE+1) = 0.5 * ( PRHOM(IKE) + PRHOM(IKE+1) ) / &
+ ZDZM(IKE+1) **2
+ !
+ELSE
+ DO JK= IKB,IKE
+ DO JJ= 1, IYMODEY_ll
+ DO JI= 1, IXMODEY_ll
+ PBFY(JJ,JI,JK) = PRHOM(JK)* ZEIGEN_ll(JI+IORXY_ll-1,JJ+IORYY_ll-1) - 0.5 * &
+ ( ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2 &
+ +( PRHOM(JK) + PRHOM(JK+1) ) / ZDZM(JK+1)**2 )
+ END DO
+ END DO
+ END DO
+! at the upper level PBFY is computed using the Neumann condition
+! while at the lower level, the Dirichlet condition is used
+!
+! lower level: Dirichlet
+ PBFY(1:IYMODEY_ll,1:IXMODEY_ll,IKB) = PBFY(1:IYMODEY_ll,1:IXMODEY_ll,IKB) - &
+ PAF(IKB)
+! PAF(IKB) = 0.0 ! done later
+ PBFY(1:IYMODEY_ll,1:IXMODEY_ll,IKB-1) = 1.0
+!
+! upper level: Neumann
+ PBFY(1:IYMODEY_ll,1:IXMODEY_ll,IKE+1) = 0.5 * ( PRHOM(IKE) + PRHOM(IKE+1) ) / &
+ ZDZM(IKE+1) **2
+!
+!++cb - Z_SPLITTING
+!++cb - for Z splitting we need to do the vertical substitution
+!++cb - in Bsplitting slides so need PBFB
+ DO JK= IKB,IKE
+ DO JJ= IJB,IJE
+ DO JI= IIB, IIE
+ PBFB(JI,JJ,JK) = PRHOM(JK)* ZEIGEN_ll(JI+IORXB_ll-IIB_ll,JJ+IORYB_ll-IJB_ll) - 0.5 * &
+ ( ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2 &
+ +( PRHOM(JK) + PRHOM(JK+1) ) / ZDZM(JK+1)**2 )
+ END DO
+ END DO
+ END DO
+! at the upper level PBFB is computed using the Neumann condition
+! while at the lower level, the Dirichlet is used
+!
+! lower level: Dirichlet
+ PBFB(IIB:IIE,IJB:IJE,IKB) = PBFB(IIB:IIE,IJB:IJE,IKB) - &
+ PAF(IKB)
+!
+! PAF(IKB) = 0.0
+ PBFB(IIB:IIE,IJB:IJE,IKB-1) = 1.0
+ !
+! upper level: Neumann
+ PBFB(IIB:IIE,IJB:IJE,IKE+1) = 0.5 * ( PRHOM(IKE) + PRHOM(IKE+1) ) / &
+ ZDZM(IKE+1) **2
+!
+IF (HLBCX(1) == 'CYCL' .AND. .NOT.(L2D) ) THEN
+ !++cb
+ ! fill unused 2 coef with NI+1 coef (lost in Z transposition elsewhere )
+ JI = IXMODE_ll -1
+ ZEIGEN_ll(2,:) = ZEIGEN_ll(JI,:)
+END IF
+IF (HLBCY(1) == 'CYCL' .AND. .NOT.(L2D) ) THEN
+ !++cb
+ ! fill unused (:,2,:) coef with NJ+1 coef (lost in Z transposition elsewhere )
+ JJ = IYMODE_ll -1
+ ZEIGEN_ll(:,2) = ZEIGEN_ll(:,JJ)
+END IF
+ !
+!++cb - Z_SPLITTING
+!++cb - Z_SPLITTING
+!++cb - for Z splitting we need to do the vertical substitution
+!++cb - in _SXP2_YP1_Zsplitting slides so need PBF_SXP2_YP1_Z
+ DO JK=IKB,IKE
+ DO JJ= 1, IJU_SXP2_YP1_Z_ll
+ DO JI= 1, IIU_SXP2_YP1_Z_ll
+ PBF_SXP2_YP1_Z(JI,JJ,JK) = PRHOM(JK)* ZEIGEN_ll(JI+IORX_SXP2_YP1_Z_ll-IIB_ll,JJ+IORY_SXP2_YP1_Z_ll-IJB_ll) - 0.5 * &
+ ( ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2 &
+ +( PRHOM(JK) + PRHOM(JK+1) ) / ZDZM(JK+1)**2 )
+ END DO
+ END DO
+ END DO
+! at the upper level PBF_SXP2_YP1_Z is computed using the Neumann condition
+! while at the lower level, the Dirichlet condition is used
+!
+! lower level: Dirichlet
+ PBF_SXP2_YP1_Z(1:IIU_SXP2_YP1_Z_ll,1:IJU_SXP2_YP1_Z_ll,IKB) = PBF_SXP2_YP1_Z(1:IIU_SXP2_YP1_Z_ll,1:IJU_SXP2_YP1_Z_ll,IKB) - &
+ PAF(IKB)
+ PAF(IKB) = 0.0
+ PBF_SXP2_YP1_Z(1:IIU_SXP2_YP1_Z_ll,1:IJU_SXP2_YP1_Z_ll,IKB-1) = 1.0
+ !
+! upper level: Neumann
+ PBF_SXP2_YP1_Z(1:IIU_SXP2_YP1_Z_ll,1:IJU_SXP2_YP1_Z_ll,IKE+1) = 0.5 * ( PRHOM(IKE) + PRHOM(IKE+1) ) / &
+ ZDZM(IKE+1) **2
+!++cb
+END IF
+!
+! second coefficent is meaningless in cyclic case
+IF (HLBCX(1) == 'CYCL' .AND. L2D .AND. SIZE(PBFY,1) .GE. 2 ) PBFY(2,:,:)=1.
+IF (HLBCX(1) == 'CYCL' .AND. .NOT.(L2D) .AND. LWEST_ll(HSPLITTING='Y') .AND. SIZE(PBFY,2) .GE. 2 ) &
+ PBFY(:,2,:)=1.
+IF (HLBCY(1) == 'CYCL' .AND. .NOT.(L2D) .AND. SIZE(PBFY,1) .GE. 2 ) PBFY(2,:,:)=1.
+!++cb - Z_SPLITTING
+! second coefficent is meaningless in cyclic case
+IF (HLBCX(1) == 'CYCL' .AND. L2D .AND. SIZE(PBFB,1) .GE. 2 ) PBFB(2,:,:)=1.
+IF (HLBCX(1) == 'CYCL' .AND. .NOT.(L2D) .AND. LWEST_ll(HSPLITTING='B') .AND. SIZE(PBFB,2) .GE. 2 ) &
+ PBFB(:,2,:)=1.
+IF (HLBCY(1) == 'CYCL' .AND. .NOT.(L2D) .AND. SIZE(PBFB,1) .GE. 2 ) PBFB(2,:,:)=1.
+!++cb
+!
+DEALLOCATE(ZEIGEN_ll)
+!
+!
+!------------------------------------------------------------------------------
+!* 8. INITIALIZATION OF THE TRIGS AND IFAX ARRAYS FOR THE FFT
+! -------------------------------------------------------
+!
+! 8.1 x lateral boundary conditions
+!
+CALL SET99(PTRIGSX,KIFAXX,IIMAX_ll)
+!
+! test on the value of KIMAX: KIMAX must be factorizable as a product
+! of powers of 2,3 and 5. KIFAXX(10) is equal to IIMAX if the decomposition
+! is correct, then KIFAXX(1) contains the number of decomposition factors
+! of KIMAX.
+!
+IF (KIFAXX(10) /= IIMAX_ll) THEN
+ WRITE(UNIT=ILUOUT,FMT="(' ERROR',/, &
+ &' : THE FORM OF THE FFT USED FOR THE INVERSION OF THE FLAT ',/,&
+ &' OPERATOR REQUIRES THAT KIMAX MUST BE FACTORIZABLE' ,/,&
+ & ' AS A PRODUCT OF POWERS OF 2, 3 AND 5.')")
+ !callabortstop
+ CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
+ CALL ABORT
+ STOP
+END IF
+!
+IF (HLBCX(1) /= 'CYCL') THEN
+!
+! extra trigs for shifted (co) sine transform (FFT55)
+!
+ INN=2*(IIMAX_ll)
+ ZDEL=ASIN(1.0)/REAL(IIMAX_ll)
+ DO JI=1,IIMAX_ll
+ ZANGLE=REAL(JI)*ZDEL
+ PTRIGSX(INN+JI)=SIN(ZANGLE)
+ END DO
+!
+ENDIF
+!
+! 8.2 y lateral boundary conditions
+!
+IF (.NOT. L2D) THEN
+ CALL SET99(PTRIGSY,KIFAXY,IJMAX_ll)
+ !
+ ! test on the value of KJMAX: KJMAX must be factorizable as a product
+ ! of powers of 2,3 and 5. KIFAXY(10) is equal to IJMAX_ll if the decomposition
+ ! is correct, then KIFAXX(1) contains the number of decomposition factors
+ ! of IIMAX_ll.
+ !
+ IF (KIFAXY(10) /= IJMAX_ll) THEN
+ WRITE(UNIT=ILUOUT,FMT="(' ERROR',/, &
+ &' : THE FORM OF THE FFT USED FOR THE INVERSION OF THE FLAT ',/,&
+ &' OPERATOR REQUIRES THAT KJMAX MUST BE FACTORIZABLE' ,/,&
+ & ' AS A PRODUCT OF POWERS OF 2, 3 AND 5.')")
+ !callabortstop
+ CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
+ CALL ABORT
+ STOP
+ END IF
+ !
+ ! 8.3 top boundary conditions
+ !
+ PEPOTFW_TOP(:,:) = PEPOTFW_TOP(:,:) / ZDZM(IKE+1)
+
+ !
+ ! other cases
+ !
+ IF (HLBCY(1) /= 'CYCL') THEN
+ !
+ ! extra trigs for shifted (co) sine transform
+ !
+ INN=2*(IJMAX_ll)
+ ZDEL=ASIN(1.0)/REAL(IJMAX_ll)
+ DO JJ=1,IJMAX_ll
+ ZANGLE=REAL(JJ)*ZDEL
+ PTRIGSY(INN+JJ)=SIN(ZANGLE)
+ END DO
+ !
+ ENDIF
+ !
+ENDIF
+!
+!------------------------------------------------------------------------------
+!
+END SUBROUTINE ELEC_TRIDZ
diff --git a/src/MNH/flash_geom_elec.f90 b/src/MNH/flash_geom_elec.f90
index baa20675a..c6d39b48f 100644
--- a/src/MNH/flash_geom_elec.f90
+++ b/src/MNH/flash_geom_elec.f90
@@ -7,15 +7,18 @@
! #############################
!
INTERFACE
- SUBROUTINE FLASH_GEOM_ELEC_n (KTCOUNT, KRR, PTSTEP, PRHODJ, PRHODREF, &
- PRT, PCIT, PRSVS, PRS, PTHT, PPABST, &
- PEFIELDU, PEFIELDV, PEFIELDW, &
- PZZ, PTOWN, PSEA)
+ SUBROUTINE FLASH_GEOM_ELEC_n (KTCOUNT, KMI, KRR, PTSTEP, OEXIT, &
+ PRHODJ, PRHODREF, &
+ PRT, PCIT, PRSVS, PRS, PTHT, PPABST, &
+ PEFIELDU, PEFIELDV, PEFIELDW, &
+ PZZ, PSVS_LINOX, PTOWN, PSEA )
!
INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop counter
+INTEGER, INTENT(IN) :: KMI ! current model index
INTEGER, INTENT(IN) :: KRR ! number of moist variables
REAL, INTENT(IN) :: PTSTEP ! Double time step except for
! cold start
+LOGICAL, INTENT(IN) :: OEXIT ! switch for the end of the temporal loop
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! Reference dry air density
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! Moist variables at time t
@@ -28,6 +31,7 @@ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRS ! Moist variables vol. sourc
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta (K) at time t
REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Absolute pressure at t
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSVS_LINOX ! NOx source term
REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PTOWN ! town fraction
REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA ! Land-sea mask
!
@@ -36,12 +40,13 @@ END INTERFACE
END MODULE MODI_FLASH_GEOM_ELEC_n
!
!
-! #######################################################################
- SUBROUTINE FLASH_GEOM_ELEC_n (KTCOUNT, KRR, PTSTEP, PRHODJ, PRHODREF, &
- PRT, PCIT, PRSVS, PRS, PTHT, PPABST, &
- PEFIELDU, PEFIELDV, PEFIELDW, &
- PZZ, PTOWN, PSEA)
-! #######################################################################
+! ####################################################################
+ SUBROUTINE FLASH_GEOM_ELEC_n (KTCOUNT, KMI, KRR, PTSTEP, OEXIT, &
+ PRHODJ, PRHODREF, &
+ PRT, PCIT, PRSVS, PRS, PTHT, PPABST, &
+ PEFIELDU, PEFIELDV, PEFIELDW, &
+ PZZ, PSVS_LINOX, PTOWN, PSEA )
+! ####################################################################
!
!!**** * -
!!
@@ -72,28 +77,42 @@ END MODULE MODI_FLASH_GEOM_ELEC_n
!! Original : Jan. 2010
!! Modifications:
!! M. Chong * LA * Juin 2010 : add small ions
-!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! J-P Pinty * LA * Feb. 2013 : add LMA storage
+!! J-P Pinty * LA * Nov. 2013 : add flash map storage
+!! M. Chong * LA * Juin 2010 : add LiNOx
+!! C. Barthe * LACy * Jan. 2015 : convert trig. pt into lat,lon in ascii file
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
-USE MODD_CONF, ONLY : CEXP
-USE MODD_PARAMETERS, ONLY : JPVEXT
+USE MODD_CST, ONLY : XAVOGADRO, XMD
+USE MODD_CONF, ONLY : CEXP, LCARTESIAN
+USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT
+USE MODD_GRID, ONLY : XLATORI,XLONORI
USE MODD_GRID_n, ONLY : XXHAT, XYHAT, XZHAT
USE MODD_DYN_n, ONLY : XDXHATM, XDYHATM, NSTOP
+USE MODD_METRICS_n, ONLY : XDXX, XDYY, XDZZ ! in linox_production
USE MODD_ELEC_DESCR
USE MODD_ELEC_PARAM, ONLY : XFQLIGHTR, XEXQLIGHTR, &
XFQLIGHTI, XEXQLIGHTI, &
XFQLIGHTS, XEXQLIGHTS, &
XFQLIGHTG, XEXQLIGHTG, &
+ XFQLIGHTH, XEXQLIGHTH, &
XFQLIGHTC
USE MODD_RAIN_ICE_DESCR, ONLY : XLBR, XLBEXR, XLBS, XLBEXS, &
- XLBG, XLBEXG, XLBDAS_MAX, XRTMIN, &
- XLBDAR_MAX, XLBDAG_MAX
+ XLBG, XLBEXG, XLBH, XLBEXH, &
+ XRTMIN
USE MODD_NSV, ONLY : NSV_ELECBEG, NSV_ELECEND, NSV_ELEC
USE MODD_VAR_ll, ONLY : NPROC
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+USE MODD_PRINT_ELEC, ONLY : NLU_fgeom_diag, NLU_fgeom_coord, &
+ NIOSTAT_fgeom_diag, NIOSTAT_fgeom_coord
+USE MODD_SUB_ELEC_n
+USE MODD_TIME_n
+USE MODD_LMA_SIMULATOR
+USE MODD_ELEC_FLASH
!
USE MODI_SHUMAN
USE MODI_TO_ELEC_FIELD_n
@@ -106,10 +125,7 @@ USE MODE_PACK_PGI
!
USE MODE_ll
USE MODE_ELEC_ll
-USE MODD_ARGSLIST_ll, ONLY : LIST_ll
-USE MODD_PRINT_ELEC, ONLY : NLU_fgeom_diag, NLU_fgeom_coord, &
- NIOSTAT_fgeom_diag, NIOSTAT_fgeom_coord
-USE MODD_SUB_ELEC_n
+USE MODE_GRIDPROJ
!
IMPLICIT NONE
!
@@ -117,9 +133,11 @@ IMPLICIT NONE
! 0.1 Declaration of arguments
!
INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop counter
+INTEGER, INTENT(IN) :: KMI ! current model index
INTEGER, INTENT(IN) :: KRR ! number of moist variables
REAL, INTENT(IN) :: PTSTEP ! Double time step except for
! cold start
+LOGICAL, INTENT(IN) :: OEXIT ! switch for the end of the temporal loop
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! Reference dry air density
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! Moist variables at time t
@@ -132,6 +150,7 @@ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRS ! Moist variables vol. sourc
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta (K) at time t
REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Absolute pressure at t
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSVS_LINOX ! NOx source term
REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PTOWN ! town fraction
REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA ! Land-sea mask
!
@@ -143,6 +162,7 @@ INTEGER :: IJB, IJE ! index values of the first and last inner mass points alon
INTEGER :: IKB, IKE ! index values of the first and last inner mass points along z
INTEGER :: IKU
INTEGER :: II, IJ, IK, IL, IM, IPOINT ! loop indexes
+INTEGER :: IX, IY, IZ
INTEGER :: IXOR, IYOR ! origin of the extended subdomain
INTEGER :: INB_CELL ! Number of detected electrified cells
INTEGER :: IPROC_CELL ! Proc with the center of the cell
@@ -180,14 +200,15 @@ INTEGER :: IWORK
INTEGER :: ICHOICE
INTEGER :: IIMIN, IIMAX, IJMIN, IJMAX, IKMIN, IKMAX
INTEGER :: IPOS_LEADER, INEG_LEADER
-INTEGER :: INBSEG_GLOB ! global number of segments
INTEGER :: INBLIGHT
INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ITYPE ! flash type (IC, CGN or CGP)
INTEGER, DIMENSION(:), ALLOCATABLE :: INBSEG_LEADER ! number of segments in the leader
INTEGER, DIMENSION(:), ALLOCATABLE :: ISIGNE_EZ ! sign of the vertical electric field
! component at the trig. pt
INTEGER, DIMENSION(:), ALLOCATABLE :: IPROC_TRIG ! proc that contains the triggering point
-INTEGER, DIMENSION(:), ALLOCATABLE :: INBSEG ! Number of segments per flash
+INTEGER, DIMENSION(:), ALLOCATABLE :: INBSEG ! Number of segments per flash
+INTEGER, DIMENSION(:), ALLOCATABLE :: INBSEG_ALL ! Number of segments, all processes
+INTEGER, DIMENSION(NPROC) :: INBSEG_PROC ! ------------------ per process
INTEGER, DIMENSION(:), ALLOCATABLE :: INB_FLASH ! Number of flashes per time step / cell
INTEGER, DIMENSION(:), ALLOCATABLE :: INB_FL_REAL ! Effective Number of flashes per timestep/cell
INTEGER, DIMENSION(:), ALLOCATABLE :: IHIST_LOC ! local nb of possible branches at [r,r+dr]
@@ -219,18 +240,14 @@ REAL :: ZE_TRIG_THRES ! Triggering Electric field threshold corrected for
! pressure
REAL :: ZMAXE ! Max electric field module (V/m)
REAL :: ZEMOD_BL ! E module at the tip of the last segment of the leader (V/m)
-REAL :: IICOORD_GLOB ! global i index of the cell point
-REAL :: IJCOORD_GLOB ! global j index of the cell point
REAL :: ZMEAN_GRID ! mean grid size
REAL :: ZMAX_DIST ! max distance between the triggering pt and the possible branches
REAL :: ZMIN_DIST ! min distance between the triggering pt and the possible branches
-REAL :: ZRAD_INF ! minimum radius from which to build the histograms
-REAL :: ZRAD_SUP ! maximum radius up to which the histograms are build
REAL :: ZRANDOM ! random number
REAL :: ZQNET ! net charge carried by the flash (C/kg)
REAL :: ZCLOUDLIM ! cloud limit
REAL :: ZSIGMIN ! min efficient cross section
-!
+REAL :: ZLAT, ZLON ! lat,lon coordinates of the triggering points if not lcartesian
!
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZQMT ! mass charge density (C/kg)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZCELL ! define the electrified cells
@@ -238,10 +255,10 @@ REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZSIGMA ! efficient cross section of hyd
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZDQDT ! charge to neutralize at each pt (C/kg)
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZFLASH ! = 1 if the flash leader reaches this pt
! = 2 if the flash branch is concerned
-REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZQTRANSFER ! Charge distributed on ions
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZLBDAR ! Lambda for rain
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZLBDAS ! Lambda for snow
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZLBDAG ! Lambda for graupel
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZLBDAH ! Lambda for hail
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZQMTOT ! total mass charge density (C/kg)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZCLOUD ! total mixing ratio (kg/kg)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEMODULE ! Electric field module (V/m)
@@ -249,11 +266,18 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDIST ! distance between the trig. pt
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSIGLOB ! sum of the cross sections
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZQFLASH ! total charge in excess of xqexcess (C/kg)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWORK
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOORD_TRIG ! Global coordinates of triggering point
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOORD_SEG ! Global coordinates of segments
-REAL, DIMENSION(:,:), ALLOCATABLE :: ZCELL_GLOB ! coordinates of the cell 'center' (m)
REAL, DIMENSION(:), ALLOCATABLE :: ZEM_TRIG ! Electric field module at the triggering pt
REAL, DIMENSION(:), ALLOCATABLE :: ZNEUT_POS ! Positive charge neutralized at each segment
REAL, DIMENSION(:), ALLOCATABLE :: ZNEUT_NEG ! Negative charge neutralized at each segment
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: ISEG_GLOB ! Global indexes of LMA segments
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: ILMA_SEG_ALL ! Global indexes of LMA segments
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZLMA_QMT ! Particle charge at neutralization point
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZLMA_PRT ! Particle mixing ratio at neutralization point
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZLMA_NEUT_POS
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZLMA_NEUT_NEG
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOORD_SEG_ALL
REAL, DIMENSION(:), ALLOCATABLE :: ZEMAX ! Max electric field in each cell
REAL, DIMENSION(:), ALLOCATABLE :: ZHIST_PERCENT ! percentage of possible branches at [r,r+dr] on each proc
REAL, DIMENSION(:), ALLOCATABLE :: ZMAX_BRANCH ! max nb of branches at [r,r+dr]
@@ -266,6 +290,15 @@ INTEGER, SAVE :: ISAVE_STATUS ! 0: print and save
!
TYPE(LIST_ll), POINTER :: TZFIELDS_ll=> NULL() ! list of fields to exchange
!
+! Storage for the localization of the flashes
+LOGICAL :: GFIRSTFLASH
+INTEGER,DIMENSION(SIZE(PRT,1),SIZE(PRT,2)) :: IMAP2D
+!
+! Storage for the NOx production terms
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZLNOX
+REAL :: ZLGHTLENGTH, ZCOEF
+INTEGER :: IFLASH_COUNT, IFLASH_COUNT_GLOB ! Total number of flashes within the timestep
+!
!-------------------------------------------------------------------------------
!
!* 1. INITIALIZATION
@@ -275,7 +308,10 @@ CALL MYPROC_ELEC_ll(IPROC)
!* 1.1 subdomains indexes
!
! beginning and end indexes of the physical subdomain
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IIB = 1 + JPHEXT
+IIE = SIZE(PRT,1) - JPHEXT
+IJB = 1 + JPHEXT
+IJE = SIZE(PRT,2) - JPHEXT
IKB = 1 + JPVEXT
IKE = SIZE(PRT,3) - JPVEXT
IKU = SIZE(PRT,3)
@@ -297,6 +333,16 @@ IF (GEFIRSTCALL) THEN
ALLOCATE (ZYMASS(SIZE(XYHAT)))
ALLOCATE (ZZMASS(SIZE(PZZ,1), SIZE(PZZ,2), SIZE(PZZ,3)))
ALLOCATE (ZPRES_COEF(SIZE(PZZ,1), SIZE(PZZ,2), SIZE(PZZ,3)))
+ IF(LLMA) THEN
+ ALLOCATE (ZLMA_LAT(NFLASH_WRITE, NBRANCH_MAX))
+ ALLOCATE (ZLMA_LON(NFLASH_WRITE, NBRANCH_MAX))
+ ALLOCATE (ZSLMA_NEUT_POS(NFLASH_WRITE, NBRANCH_MAX))
+ ALLOCATE (ZSLMA_NEUT_NEG(NFLASH_WRITE, NBRANCH_MAX))
+ ALLOCATE (ISLMA_SEG_GLOB(NFLASH_WRITE, NBRANCH_MAX, 3))
+ ALLOCATE (ZSLMA_QMT(NFLASH_WRITE, NBRANCH_MAX, SIZE(PRSVS,4)))
+ ALLOCATE (ZSLMA_PRT(NFLASH_WRITE, NBRANCH_MAX, SIZE(PRSVS,4)))
+ ISLMA_SEG_GLOB(:,:,:) = 0
+ END IF
ALLOCATE (ZSCOORD_SEG(NFLASH_WRITE, NBRANCH_MAX, 3)) ! NFLASH_WRITE nb of flash to be stored
! before writing in files
! NBRANCH_MAX=5000 default
@@ -325,7 +371,6 @@ ALLOCATE (ZCLOUD(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
ALLOCATE (GPOSS(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
ALLOCATE (ZEMODULE(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
ALLOCATE (ZCELL(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NMAX_CELL))
-ALLOCATE (ZQTRANSFER(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
!
ZQMT(:,:,:,:) = 0.
ZQMTOT(:,:,:) = 0.
@@ -340,6 +385,7 @@ ZCELL(:,:,:,:) = 0.
!
PRSVS(:,:,:,1) = XECHARGE * PRSVS(:,:,:,1) ! C /(m3 s)
PRSVS(:,:,:,NSV_ELEC) = -1. * XECHARGE * PRSVS(:,:,:,NSV_ELEC) ! C /(m3 s)
+!
CALL PT_DISCHARGE
!
!
@@ -350,7 +396,7 @@ DO II = 1, NSV_ELEC
ZQMT(:,:,:,II) = PRSVS(:,:,:,II) * PTSTEP / PRHODJ(:,:,:)
!
! total mass charge density (C/kg)
- ZQMTOT(:,:,:) = ZQMTOT(:,:,:) + PRSVS(:,:,:,II) * PTSTEP / PRHODJ(:,:,:)
+ ZQMTOT(:,:,:) = ZQMTOT(:,:,:) + PRSVS(:,:,:,II) * PTSTEP / PRHODJ(:,:,:)
END DO
!
! total mixing ratio (g/kg)
@@ -371,12 +417,10 @@ ZSIGMIN = 1.E-12
! ------------------------------------
!
ALLOCATE (ZEMAX(NMAX_CELL))
-ALLOCATE (ICELL_LOC(NMAX_CELL,4))
-ALLOCATE (ZCELL_GLOB(NMAX_CELL,3))
+ALLOCATE (ICELL_LOC(4,NMAX_CELL))
!
ZEMAX(:) = 0.
ICELL_LOC(:,:) = 0
-ZCELL_GLOB(:,:) = 0.
!
WHERE (ZCLOUD(IIB:IIE,IJB:IJE,IKB:IKE) .LE. ZCLOUDLIM)
GPOSS(IIB:IIE,IJB:IJE,IKB:IKE) = .FALSE.
@@ -399,7 +443,6 @@ ZEMODULE(IIB:IIE,IJB:IJE,IKB:IKE) = ZPRES_COEF(IIB:IIE,IJB:IJE,IKB:IKE)* &
PEFIELDV(IIB:IIE,IJB:IJE,IKB:IKE)**2 + &
PEFIELDW(IIB:IIE,IJB:IJE,IKB:IKE)**2)**0.5
!
-!
DO WHILE (.NOT. GEND_DOMAIN .AND. INB_CELL .LT. NMAX_CELL)
!
! find the maximum electric field on each proc
@@ -414,19 +457,17 @@ DO WHILE (.NOT. GEND_DOMAIN .AND. INB_CELL .LT. NMAX_CELL)
!
IF (ZMAXE .GT. ZE_TRIG_THRES) THEN
INB_CELL = INB_CELL + 1 ! one cell is detected
-!
ZEMAX(INB_CELL) = ZMAXE
! local coordinates of the maximum electric field
-! ICELL_LOC(INB_CELL,1:3) = MAXLOC(ZEMODULE(:,:,:),MASK=GPOSS(:,:,:))
- ICELL_LOC(INB_CELL,1:3) = MAXLOC(ZEMODULE(IIB:IIE,IJB:IJE,IKB:IKE), &
+ ICELL_LOC(1:3,INB_CELL) = MAXLOC(ZEMODULE(IIB:IIE,IJB:IJE,IKB:IKE), &
MASK=GPOSS(IIB:IIE,IJB:IJE,IKB:IKE))
- IICOORD = ICELL_LOC(INB_CELL,1)
- IJCOORD = ICELL_LOC(INB_CELL,2)
- IKCOORD = ICELL_LOC(INB_CELL,3)
- ICELL_LOC(INB_CELL,4) = IPROC_CELL
+ IICOORD = ICELL_LOC(1,INB_CELL)
+ IJCOORD = ICELL_LOC(2,INB_CELL)
+ IKCOORD = ICELL_LOC(3,INB_CELL)
+ ICELL_LOC(4,INB_CELL) = IPROC_CELL
!
! Broadcast the center of the cell to all procs
- CALL MPI_BCAST (ICELL_LOC(INB_CELL,:), 4, MPI_INTEGER, IPROC_CELL, &
+ CALL MPI_BCAST (ICELL_LOC(:,INB_CELL), 4, MPI_INTEGER, IPROC_CELL, &
MPI_COMM_WORLD, IERR)
!
!
@@ -454,27 +495,27 @@ DO WHILE (.NOT. GEND_DOMAIN .AND. INB_CELL .LT. NMAX_CELL)
!
DO II = IIB, IIE
DO IJ = IJB, IJE
- IF ((ZCELL(II,IJ,IK,INB_CELL) .EQ. 0.) .AND. &
- (GPOSS(II,IJ,IK)) .AND. &
- (ZCLOUD(II,IJ,IK) .GT. 1.E-5) .AND. &
- ((ABS(ZQMT(II,IJ,IK,2)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
- (ABS(ZQMT(II,IJ,IK,3)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
- (ABS(ZQMT(II,IJ,IK,4)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
- (ABS(ZQMT(II,IJ,IK,5)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
- (ABS(ZQMT(II,IJ,IK,6)) * PRHODREF(II,IJ,IK) .GT. XQEXCES)) )THEN
-
- IF ((ZCELL(II-1,IJ, IK,INB_CELL) .EQ. 1.) .OR. &
- (ZCELL(II+1,IJ, IK,INB_CELL) .EQ. 1.) .OR. &
- (ZCELL(II, IJ-1,IK,INB_CELL) .EQ. 1.) .OR. &
- (ZCELL(II, IJ+1,IK,INB_CELL) .EQ. 1.) .OR. &
- (ZCELL(II-1,IJ-1,IK,INB_CELL) .EQ. 1.) .OR. &
- (ZCELL(II-1,IJ+1,IK,INB_CELL) .EQ. 1.) .OR. &
- (ZCELL(II+1,IJ+1,IK,INB_CELL) .EQ. 1.) .OR. &
- (ZCELL(II+1,IJ-1,IK,INB_CELL) .EQ. 1.)) THEN
- ZCELL(II,IJ,IK,INB_CELL) = 1.
- GPOSS(II,IJ,IK) = .FALSE.
- END IF
- END IF
+ IF ((ZCELL(II,IJ,IK,INB_CELL) .EQ. 0.) .AND. &
+ (GPOSS(II,IJ,IK)) .AND. &
+ (ZCLOUD(II,IJ,IK) .GT. 1.E-5) .AND. &
+ ((ABS(ZQMT(II,IJ,IK,2)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
+ (ABS(ZQMT(II,IJ,IK,3)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
+ (ABS(ZQMT(II,IJ,IK,4)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
+ (ABS(ZQMT(II,IJ,IK,5)) * PRHODREF(II,IJ,IK) .GT. XQEXCES).OR. &
+ (ABS(ZQMT(II,IJ,IK,6)) * PRHODREF(II,IJ,IK) .GT. XQEXCES)) )THEN
+!
+ IF ((ZCELL(II-1,IJ, IK,INB_CELL) .EQ. 1.) .OR. &
+ (ZCELL(II+1,IJ, IK,INB_CELL) .EQ. 1.) .OR. &
+ (ZCELL(II, IJ-1,IK,INB_CELL) .EQ. 1.) .OR. &
+ (ZCELL(II, IJ+1,IK,INB_CELL) .EQ. 1.) .OR. &
+ (ZCELL(II-1,IJ-1,IK,INB_CELL) .EQ. 1.) .OR. &
+ (ZCELL(II-1,IJ+1,IK,INB_CELL) .EQ. 1.) .OR. &
+ (ZCELL(II+1,IJ+1,IK,INB_CELL) .EQ. 1.) .OR. &
+ (ZCELL(II+1,IJ-1,IK,INB_CELL) .EQ. 1.)) THEN
+ ZCELL(II,IJ,IK,INB_CELL) = 1.
+ GPOSS(II,IJ,IK) = .FALSE.
+ END IF
+ END IF
END DO
END DO
!
@@ -519,13 +560,37 @@ IF (INB_CELL .GE. 1) THEN
! laisse tel quel pour le moment
!
ALLOCATE (ISEG_LOC(3*SIZE(PRT,3), INB_CELL)) ! 3 coord indices of the leader
+ ALLOCATE (ZCOORD_TRIG(3, INB_CELL))
ALLOCATE (ZCOORD_SEG(NBRANCH_MAX*3, INB_CELL))
! NBRANCH_MAX=5000 default
! 3= 3 coord index
+ ALLOCATE (ZCOORD_SEG_ALL(NBRANCH_MAX*3, INB_CELL))
+ ALLOCATE (ISEG_GLOB(NBRANCH_MAX*3, INB_CELL))
+ ISEG_GLOB(:,:) = 0
+!
+ IF(LLMA) THEN
+ ALLOCATE (ILMA_SEG_ALL (NBRANCH_MAX*3, INB_CELL))
+ ALLOCATE (ZLMA_QMT(NBRANCH_MAX*NSV_ELEC, INB_CELL)) ! charge des part.
+ ! a neutraliser
+ ALLOCATE (ZLMA_PRT(NBRANCH_MAX*NSV_ELEC, INB_CELL)) ! mixing ratio
+ ALLOCATE (ZLMA_NEUT_POS(NBRANCH_MAX, INB_CELL))
+ ALLOCATE (ZLMA_NEUT_NEG(NBRANCH_MAX, INB_CELL))
+ ZLMA_QMT(:,:) = 0.
+ ZLMA_PRT(:,:) = 0.
+ ZLMA_NEUT_POS(:,:) = 0.
+ ZLMA_NEUT_NEG(:,:) = 0.
+ END IF
+!
+ IF (LLNOX_EXPLICIT) THEN
+ ALLOCATE (ZLNOX(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
+ ZLNOX(:,:,:) = 0.
+ END IF
+!
ALLOCATE (ZEM_TRIG(INB_CELL))
ALLOCATE (INB_FLASH(INB_CELL))
ALLOCATE (INB_FL_REAL(INB_CELL))
ALLOCATE (INBSEG(INB_CELL))
+ ALLOCATE (INBSEG_ALL(INB_CELL))
ALLOCATE (ITYPE(INB_CELL))
ALLOCATE (INBSEG_LEADER(INB_CELL))
ALLOCATE (ZDQDT(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),SIZE(PRT,4)+1))
@@ -533,6 +598,7 @@ IF (INB_CELL .GE. 1) THEN
ALLOCATE (ZLBDAR(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
ALLOCATE (ZLBDAS(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
ALLOCATE (ZLBDAG(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
+ IF (KRR == 7) ALLOCATE (ZLBDAH(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
ALLOCATE (ZSIGLOB(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
ALLOCATE (ZFLASH(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),INB_CELL))
ALLOCATE (ZDIST(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
@@ -540,6 +606,7 @@ IF (INB_CELL .GE. 1) THEN
ALLOCATE (GATTACH(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
!
ISEG_LOC(:,:) = 0
+ ZCOORD_TRIG(:,:) = 0.
ZCOORD_SEG(:,:) = 0.
ZDQDT(:,:,:,:) = 0.
ZSIGMA(:,:,:,:) = 0.
@@ -554,6 +621,8 @@ IF (INB_CELL .GE. 1) THEN
INB_FLASH(:) = 0
INB_FL_REAL(:) = 0
INBSEG(:) = 0
+ INBSEG_ALL(:) = 0
+ INBSEG_PROC(:) = 0
INBSEG_LEADER(:) = 0
ITYPE(:) = 1 ! default = IC
!
@@ -577,7 +646,7 @@ IF (INB_CELL .GE. 1) THEN
MAX(PRT(:,:,:,3),XRTMIN(3)))**XLBEXR
END WHERE
!
- WHERE (PRT(:,:,:,3) > ZCLOUDLIM .AND. ZLBDAR(:,:,:) < XLBDAR_MAX .AND. &
+ WHERE (PRT(:,:,:,3) > ZCLOUDLIM .AND. ZLBDAR(:,:,:) < XLBDAR_MAXE .AND. &
ZLBDAR(:,:,:) > 0.)
ZSIGMA(:,:,:,2) = XFQLIGHTR * ZLBDAR(:,:,:)**XEXQLIGHTR
END WHERE
@@ -594,12 +663,12 @@ IF (INB_CELL .GE. 1) THEN
!* 3.4 for snow
!
WHERE (PRT(:,:,:,5) > 0.0)
- ZLBDAS(:,:,:) = MIN(XLBDAS_MAX, &
+ ZLBDAS(:,:,:) = MIN(XLBDAS_MAXE, &
XLBS * (PRHODREF(:,:,:) * &
MAX(PRT(:,:,:,5),XRTMIN(5)))**XLBEXS)
END WHERE
!
- WHERE (PRT(:,:,:,5) > ZCLOUDLIM .AND. ZLBDAS(:,:,:) < XLBDAS_MAX .AND. &
+ WHERE (PRT(:,:,:,5) > ZCLOUDLIM .AND. ZLBDAS(:,:,:) < XLBDAS_MAXE .AND. &
ZLBDAS(:,:,:) > 0.)
ZSIGMA(:,:,:,4) = XFQLIGHTS * ZLBDAS(:,:,:)**XEXQLIGHTS
ENDWHERE
@@ -611,7 +680,7 @@ IF (INB_CELL .GE. 1) THEN
ZLBDAG(:,:,:) = XLBG * (PRHODREF(:,:,:) * MAX(PRT(:,:,:,6),XRTMIN(6)))**XLBEXG
END WHERE
!
- WHERE (PRT(:,:,:,6) > ZCLOUDLIM .AND. ZLBDAG(:,:,:) < XLBDAG_MAX .AND. &
+ WHERE (PRT(:,:,:,6) > ZCLOUDLIM .AND. ZLBDAG(:,:,:) < XLBDAG_MAXE .AND. &
ZLBDAG(:,:,:) > 0.)
ZSIGMA(:,:,:,5) = XFQLIGHTG * ZLBDAG(:,:,:)**XEXQLIGHTG
ENDWHERE
@@ -619,12 +688,25 @@ IF (INB_CELL .GE. 1) THEN
!
!* 3.6 for hail
!
+ IF (KRR == 7) THEN
+ WHERE (PRT(:,:,:,7) > 0.0)
+ ZLBDAH(:,:,:) = XLBH * (PRHODREF(:,:,:) * &
+ MAX(PRT(:,:,:,7), XRTMIN(7)))**XLBEXH
+ END WHERE
+!
+ WHERE (PRT(:,:,:,7) > ZCLOUDLIM .AND. ZLBDAH(:,:,:) < XLBDAH_MAXE .AND. &
+ ZLBDAH(:,:,:) > 0.)
+ ZSIGMA(:,:,:,6) = XFQLIGHTH * ZLBDAH(:,:,:)**XEXQLIGHTH
+ ENDWHERE
+ END IF
!
!
!* 3.7 sum of the efficient cross sections
!
ZSIGLOB(:,:,:) = ZSIGMA(:,:,:,1) + ZSIGMA(:,:,:,2) + ZSIGMA(:,:,:,3) + &
ZSIGMA(:,:,:,4) + ZSIGMA(:,:,:,5)
+!
+ IF (KRR == 7) ZSIGLOB(:,:,:) = ZSIGLOB(:,:,:) + ZSIGMA(:,:,:,6)
!
!
!-------------------------------------------------------------------------------
@@ -651,6 +733,9 @@ IF (INB_CELL .GE. 1) THEN
!
!* 4. FLASH TRIGGERING
! ----------------
+!
+ IFLASH_COUNT = 0
+ IFLASH_COUNT_GLOB = 0
!
DO WHILE (GNEW_FLASH_GLOB)
!
@@ -662,7 +747,7 @@ IF (INB_CELL .GE. 1) THEN
! update lightning informations
INB_FLASH(IL) = INB_FLASH(IL) + 1 ! nb of flashes / cell / time step
INB_FL_REAL(IL) = INB_FL_REAL(IL) + 1 ! nb of flashes / cell / time step
- INBSEG(IL) = 1 ! nb of segments / flash
+ INBSEG(IL) = 0 ! nb of segments / flash
ITYPE(IL) = 1
!
IF (IPROC .EQ. IPROC_TRIG(IL)) THEN
@@ -671,7 +756,8 @@ IF (INB_CELL .GE. 1) THEN
IJBL_LOC = ISEG_LOC(2,IL)
IKBL = ISEG_LOC(3,IL)
!
- ZFLASH(IIBL_LOC,IJBL_LOC,IKBL,IL) = 1.
+ INBSEG(IL) = 1 ! nb of segments / flash
+ ZFLASH(IIBL_LOC,IJBL_LOC,IKBL,IL) = 1.
ENDIF
!
GCG = .FALSE.
@@ -702,13 +788,12 @@ IF (INB_CELL .GE. 1) THEN
CALL ONE_LEADER
!
INBSEG_LEADER(IL) = INBSEG(IL)
- INEG_LEADER = INBSEG_LEADER(IL) - IPOS_LEADER -1
+ INEG_LEADER = INBSEG_LEADER(IL) - IPOS_LEADER - 1
!
! Eliminate this flash if only positive or negative leader exists
-
IF (IPROC .EQ. IPROC_TRIG(IL)) THEN
IF (IPOS_LEADER .EQ. 0 .OR. INEG_LEADER .EQ. 0) THEN
- ZFLASH(IIBL_LOC,IJBL_LOC,IKB:IKE,IL)=0.
+ ZFLASH(IIBL_LOC,IJBL_LOC,IKB:IKE,IL) = 0.
INB_FL_REAL(IL) = INB_FL_REAL(IL) - 1
GNEW_FLASH(IL) = .FALSE.
ELSE ! return to actual Triggering electrical field
@@ -775,9 +860,9 @@ IF (INB_CELL .GE. 1) THEN
DO IJ = IJB, IJE
DO IK = IKB, IKE
IF (GPROP(II,IJ,IK,IL)) THEN
- ZDIST(II,IJ,IK) = ((ZXMASS(II) - ZCOORD_SEG(1,IL))**2 + &
- (ZYMASS(IJ) - ZCOORD_SEG(2,IL))**2 + &
- (ZZMASS(II,IJ,IK) - ZCOORD_SEG(3,IL))**2)**0.5
+ ZDIST(II,IJ,IK) = ((ZXMASS(II) - ZCOORD_TRIG(1,IL))**2 + &
+ (ZYMASS(IJ) - ZCOORD_TRIG(2,IL))**2 + &
+ (ZZMASS(II,IJ,IK) - ZCOORD_TRIG(3,IL))**2)**0.5
END IF
END DO
END DO
@@ -792,7 +877,7 @@ IF (INB_CELL .GE. 1) THEN
! transform the min and max distances into min and max increments
IIND_MIN = 1
IIND_MAX = MAX(1, INT((ZMAX_DIST-ZMIN_DIST)/ZMEAN_GRID +1.))
- IDELTA_IND = IIND_MAX +1
+ IDELTA_IND = IIND_MAX + 1
!
ALLOCATE (IHIST_LOC(IDELTA_IND))
ALLOCATE (ZHIST_PERCENT(IDELTA_IND))
@@ -860,9 +945,6 @@ IF (INB_CELL .GE. 1) THEN
! => the max number of branches / proc is proportional to the percentage of
! available points / proc at this distance
!
-! this line must be commented if branch_geom is called once
-! IMAX_BRANCH(IM) = INT(ANINT(ZMAX_BRANCH(IM) * ZHIST_PERCENT(IM) / 2.))
-! this line must be commented if branch_geom is called twice
IMAX_BRANCH(IM) = INT(ANINT(ZMAX_BRANCH(IM)))
END DO
!
@@ -875,12 +957,6 @@ IF (INB_CELL .GE. 1) THEN
!* 8.3 distribute the branches
!
ALLOCATE (ZWORK(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3)))
-!
-! upper branches
-! CALL BRANCH_GEOM(IK_TRIG+1, IKE)
-!
-! lower branches
-! CALL BRANCH_GEOM(IKB, IK_TRIG-1)
!
CALL BRANCH_GEOM(IKB, IKE)
!
@@ -952,15 +1028,15 @@ IF (INB_CELL .GE. 1) THEN
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) < 0.)
! increase negative ion charge
ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,NSV_ELEC) = &
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,NSV_ELEC) + &
- ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,NSV_ELEC) + &
+ ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
ENDWHERE
!
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) > 0.)
! Increase positive ion charge
ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,1) = &
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,1) + &
- ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,1) + &
+ ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
ENDWHERE
!
!
@@ -972,21 +1048,21 @@ IF (INB_CELL .GE. 1) THEN
IF (ITYPE(IL) .EQ. 3) THEN
DO II = 1, NSV_ELEC
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) > 0.)
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,II) = &
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,II) - &
- ZQMT(IIB:IIE,IJB:IJE,IKB:IKE,II)
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,II) = &
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,II) - &
+ ZQMT(IIB:IIE,IJB:IJE,IKB:IKE,II)
END WHERE
ENDDO
!
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) > 0.)
- ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)=0.
+ ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)=0.
END WHERE
!
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) < 0.)
! Increase negative ion charge
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,NSV_ELEC) = &
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,NSV_ELEC) + &
- ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,NSV_ELEC) = &
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,NSV_ELEC) + &
+ ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
ENDWHERE
ELSE
!
@@ -994,21 +1070,21 @@ IF (INB_CELL .GE. 1) THEN
!
DO II = 1, NSV_ELEC
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) < 0.)
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,II) = &
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,II) = &
ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,II) - &
ZQMT(IIB:IIE,IJB:IJE,IKB:IKE,II)
END WHERE
ENDDO
!
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) < 0.)
- ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)=0.
+ ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)=0.
END WHERE
!
WHERE (ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE) > 0.)
! Increase positive ion charge
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,1) = &
- ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,1) + &
- ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,1) = &
+ ZDQDT(IIB:IIE,IJB:IJE,IKB:IKE,1) + &
+ ZQFLASH(IIB:IIE,IJB:IJE,IKB:IKE)
ENDWHERE
END IF ! GCG_POS
END IF ! NOT(GCG)
@@ -1020,39 +1096,77 @@ IF (INB_CELL .GE. 1) THEN
!
CALL MPI_BCAST (INB_NEUT_OK,1, MPI_INTEGER, IPROC_TRIG(IL), &
MPI_COMM_WORLD, IERR)
- END IF ! GNEUTRALIZATION
-!
-!
-!* 9.5 update the source term
!
- DO II = IIB, IIE
- DO IJ = IJB, IJE
- DO IK = IKB, IKE
- DO IM = 1, NSV_ELEC
- IF (ZDQDT(II,IJ,IK,IM) .NE. 0.) THEN
- PRSVS(II,IJ,IK,IM) = PRSVS(II,IJ,IK,IM) + &
- ZDQDT(II,IJ,IK,IM) * &
- PRHODJ(II,IJ,IK) / PTSTEP
- END IF
+!* 9.5 Gather lightning information from all processes
+!* Save the particule charge and total pos/neg charge neutralization points.
+!* the coordinates of all flash branch points
!
+ INBSEG_PROC(IPROC+1) = INBSEG(IL)
+ DO IK = 0, NPROC-1
+ CALL MPI_BCAST (INBSEG_PROC(IK+1), 1, MPI_INTEGER, IK, &
+ MPI_COMM_WORLD, IERR)
+ END DO
+
+ INBSEG_ALL(IL) = INBSEG(IL)
+ CALL SUM_ELEC_ll(INBSEG_ALL(IL))
+
+ CALL GATHER_ALL_BRANCH
!
-!* 9.6 update the positive and negative charge neutralized
+!* 9.6 update the source term
!
- IF (ZDQDT(II,IJ,IK,IM) .LT. 0.) THEN
- ZNEUT_NEG(IL) = ZNEUT_NEG(IL) + ZDQDT(II,IJ,IK,IM)* &
- PRHODJ(II,IJ,IK)
- ELSE IF (ZDQDT(II,IJ,IK,IM) .GT. 0.) THEN
- ZNEUT_POS(IL) = ZNEUT_POS(IL) + ZDQDT(II,IJ,IK,IM)* &
- PRHODJ(II,IJ,IK)
- END IF
+ DO II = IIB, IIE
+ DO IJ = IJB, IJE
+ DO IK = IKB, IKE
+ DO IM = 1, NSV_ELEC
+ IF (ZDQDT(II,IJ,IK,IM) .NE. 0.) THEN
+ PRSVS(II,IJ,IK,IM) = PRSVS(II,IJ,IK,IM) + &
+ ZDQDT(II,IJ,IK,IM) * &
+ PRHODJ(II,IJ,IK) / PTSTEP
+ END IF
+!
+!
+!* 9.7 update the positive and negative charge neutralized
+!
+ IF (ZDQDT(II,IJ,IK,IM) .LT. 0.) THEN
+ ZNEUT_NEG(IL) = ZNEUT_NEG(IL) + ZDQDT(II,IJ,IK,IM) * &
+ PRHODJ(II,IJ,IK)
+ ELSE IF (ZDQDT(II,IJ,IK,IM) .GT. 0.) THEN
+ ZNEUT_POS(IL) = ZNEUT_POS(IL) + ZDQDT(II,IJ,IK,IM) * &
+ PRHODJ(II,IJ,IK)
+ END IF
+ END DO
END DO
END DO
END DO
- END DO
!
- CALL SUM_ELEC_ll(ZNEUT_POS(IL))
- CALL SUM_ELEC_ll(ZNEUT_NEG(IL))
+ CALL SUM_ELEC_ll(ZNEUT_POS(IL))
+ CALL SUM_ELEC_ll(ZNEUT_NEG(IL))
+!
!
+!* 9.8 compute the NOx production
+!
+!! The lightning length is first computed. The number of NOx molecules per
+!! meter of lightning flash is taken from Wang et al. (1998). It is a linear
+!! function of the pressure. No distinction is made between ICs and CGs.
+
+ IF (LLNOX_EXPLICIT) THEN
+ IFLASH_COUNT_GLOB = IFLASH_COUNT_GLOB + 1
+ IF (INBSEG(IL) .NE. 0) THEN
+ DO II = 0, INBSEG(IL)-1
+ IM = 3 * II
+ IX = ISEG_GLOB(IM+1,IL) - IXOR + 1
+ IY = ISEG_GLOB(IM+2,IL) - IYOR + 1
+ IZ = ISEG_GLOB(IM+3,IL)
+ ZLGHTLENGTH = (XDXX(IX,IY,IZ) * XDYY(IX,IY,IZ) * &
+ XDZZ(IX,IY,IZ))**(1./3.)
+ ZLNOX(IX, IY, IZ) = ZLNOX(IX, IY, IZ) + &
+ (XWANG_A + XWANG_B * PPABST(IX,IY,IZ)) * &
+ ZLGHTLENGTH
+ ENDDO
+ IFLASH_COUNT = IFLASH_COUNT + 1
+ END IF
+ END IF
+ END IF ! GNEUTRALIZATION
END IF ! end if gnew_flash
END DO ! end loop il
!
@@ -1065,6 +1179,8 @@ IF (INB_CELL .GE. 1) THEN
! ---------------------------------------------------
!
! Synchronizing all processes
+! CALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! A ACTIVER SI PB.
+!
IF (IPROC .EQ. 0) THEN
INBLIGHT = COUNT(GNEW_FLASH(1:INB_CELL))
IF (INBLIGHT .NE. 0) THEN
@@ -1073,9 +1189,9 @@ IF (INB_CELL .GE. 1) THEN
DO IL = 1, INB_CELL
IF (GNEW_FLASH(IL)) THEN
NNBLIGHT = NNBLIGHT + 1
- ISFLASH_NUMBER(NNBLIGHT) = ISFLASH_NUMBER(NNBLIGHT-1) +1
+ ISFLASH_NUMBER(NNBLIGHT) = ISFLASH_NUMBER(NNBLIGHT-1) + 1
ISNB_FLASH(NNBLIGHT) = INB_FL_REAL(IL)
- ISNBSEG(NNBLIGHT) = INBSEG(IL)
+ ISNBSEG(NNBLIGHT) = INBSEG_ALL(IL)
ISCELL_NUMBER(NNBLIGHT) = IL
ISTCOUNT_NUMBER(NNBLIGHT) = KTCOUNT
ISTYPE(NNBLIGHT) = ITYPE(IL)
@@ -1083,13 +1199,24 @@ IF (INB_CELL .GE. 1) THEN
ZSNEUT_POS(NNBLIGHT) = ZNEUT_POS(IL)
ZSNEUT_NEG(NNBLIGHT) = ZNEUT_NEG(IL)
!
- DO II = 1, INBSEG(IL)
- DO IJ = 1,3
- ZSCOORD_SEG(NNBLIGHT, II, IJ) = ZCOORD_SEG(3*(II-1)+IJ, IL)
- END DO
+ DO II = 1, INBSEG_ALL(IL)
+ IM = 3 * (II - 1)
+ ZSCOORD_SEG(NNBLIGHT,II,1:3) = ZCOORD_SEG_ALL(IM+1:IM+3,IL)
ENDDO
- END IF
- ENDDO
+!
+ IF(LLMA) THEN
+ DO II = 1, INBSEG_ALL(IL)
+ IM = 3 * (II - 1)
+ ISLMA_SEG_GLOB(NNBLIGHT,II,1:3) = ILMA_SEG_ALL(IM+1:IM+3,IL)
+ IM = NSV_ELEC * (II - 1)
+ ZSLMA_QMT(NNBLIGHT,II,2:6) = ZLMA_QMT(IM+2:IM+6,IL)
+ ZSLMA_PRT(NNBLIGHT,II,2:6) = ZLMA_PRT(IM+2:IM+6,IL)
+ ZSLMA_NEUT_POS(NNBLIGHT,II) = ZLMA_NEUT_POS(II,IL)
+ ZSLMA_NEUT_NEG(NNBLIGHT,II) = ZLMA_NEUT_NEG(II,IL)
+ END DO
+ END IF ! llma
+ END IF ! gnew_flash
+ END DO ! end loop il
!
IF (NNBLIGHT .EQ. NFLASH_WRITE) ISAVE_STATUS = 0
!
@@ -1097,8 +1224,11 @@ IF (INB_CELL .GE. 1) THEN
ISAVE_STATUS = 2
END IF
!
- IF (ISAVE_STATUS.EQ. 0 .OR. ISAVE_STATUS.EQ. 2) THEN
+ IF (ISAVE_STATUS .EQ. 0 .OR. ISAVE_STATUS .EQ. 2) THEN
CALL WRITE_OUT_ASCII
+ IF(LLMA) THEN
+ CALL WRITE_OUT_LMA
+ END IF
ISFLASH_NUMBER(0) = ISFLASH_NUMBER(NNBLIGHT)
END IF
!
@@ -1106,10 +1236,10 @@ IF (INB_CELL .GE. 1) THEN
NNBLIGHT = 0
DO IL = 1, INB_CELL
IF (GNEW_FLASH(IL)) THEN
- NNBLIGHT = NNBLIGHT + 1
- ISFLASH_NUMBER(NNBLIGHT) = ISFLASH_NUMBER(NNBLIGHT-1) +1
+ NNBLIGHT = NNBLIGHT + 1
+ ISFLASH_NUMBER(NNBLIGHT) = ISFLASH_NUMBER(NNBLIGHT-1) + 1
ISNB_FLASH(NNBLIGHT) = INB_FL_REAL(IL)
- ISNBSEG(NNBLIGHT) = INBSEG(IL)
+ ISNBSEG(NNBLIGHT) = INBSEG_ALL(IL)
ISCELL_NUMBER(NNBLIGHT) = IL
ISTCOUNT_NUMBER(NNBLIGHT) = KTCOUNT
ISTYPE(NNBLIGHT) = ITYPE(IL)
@@ -1117,11 +1247,22 @@ IF (INB_CELL .GE. 1) THEN
ZSNEUT_POS(NNBLIGHT) = ZNEUT_POS(IL)
ZSNEUT_NEG(NNBLIGHT) = ZNEUT_NEG(IL)
!
- DO II = 1, INBSEG(IL)
- DO IJ = 1,3
- ZSCOORD_SEG(NNBLIGHT, II, IJ)=ZCOORD_SEG(3*(II-1)+IJ, IL)
- ENDDO
+ DO II = 1, INBSEG_ALL(IL)
+ IM = 3 * (II - 1)
+ ZSCOORD_SEG(NNBLIGHT, II, 1:3) = ZCOORD_SEG_ALL(IM+1:IM+3, IL)
ENDDO
+!
+ IF(LLMA) THEN
+ DO II = 1, INBSEG_ALL(IL)
+ IM = 3 * (II - 1)
+ ISLMA_SEG_GLOB(NNBLIGHT,II,1:3) = ILMA_SEG_ALL(IM+1:IM+3,IL)
+ IM = NSV_ELEC*(II-1)
+ ZSLMA_QMT(NNBLIGHT,II,2:6) = ZLMA_QMT(IM+2:IM+6,IL)
+ ZSLMA_PRT(NNBLIGHT,II,2:6) = ZLMA_PRT(IM+2:IM+6,IL)
+ ZSLMA_NEUT_POS(NNBLIGHT,II) = ZLMA_NEUT_POS(II,IL)
+ ZSLMA_NEUT_NEG(NNBLIGHT,II) = ZLMA_NEUT_NEG(II,IL)
+ END DO
+ END IF
END IF
ENDDO
END IF
@@ -1132,6 +1273,45 @@ IF (INB_CELL .GE. 1) THEN
END IF ! INBLIGHT
END IF ! IPROC
!
+! Save flash location statistics in all processes
+ IF (INBLIGHT .NE. 0) THEN
+ DO IL = 1, INB_CELL
+ IF (GNEW_FLASH(IL)) THEN
+ IMAP2D(:,:) = 0
+ DO IK = IKB, IKE
+ IMAP2D(:,:) = IMAP2D(:,:) + ZFLASH(:,:,IK,IL)
+ END DO
+!
+! Detect Trig/Impact X,Y location
+ IX = 0
+ IY = 0
+ GFIRSTFLASH = .FALSE.
+ DO II = IIB, IIE
+ DO IJ = IJB, IJE
+ DO IK = IKB, IKE
+ IF (GFIRSTFLASH) EXIT
+ IF (ZFLASH(II,IJ,IK,IL)==1.) THEN
+ IX = II
+ IY = IJ
+ GFIRSTFLASH = .TRUE.
+ END IF
+ END DO
+ END DO
+ END DO
+!
+! Store
+ IF (ITYPE(IL)==1) THEN ! IC
+ IF (IX*IY/=0) NMAP_TRIG_IC(IX,IY) = NMAP_TRIG_IC(IX,IY) + 1
+ NMAP_2DAREA_IC(:,:) = NMAP_2DAREA_IC(:,:) + MIN(1,IMAP2D(:,:))
+ NMAP_3DIC(:,:,:) = NMAP_3DIC(:,:,:) + ZFLASH(:,:,:,IL)
+ ELSE ! CGN & CGP
+ IF (IX*IY/=0) NMAP_IMPACT_CG(IX,IY) = NMAP_IMPACT_CG(IX,IY) + 1
+ NMAP_2DAREA_CG(:,:) = NMAP_2DAREA_CG(:,:) + MIN(1,IMAP2D(:,:))
+ NMAP_3DCG(:,:,:) = NMAP_3DCG(:,:,:) + ZFLASH(:,:,:,IL)
+ END IF
+ END IF
+ ENDDO
+ END IF ! INBLIGHT
!
!------------------------------------------------------------------------------
!
@@ -1187,6 +1367,13 @@ IF (INB_CELL .GE. 1) THEN
PEFIELDV(IIB:IIE,IJB:IJE,IKB:IKE)**2 + &
PEFIELDW(IIB:IIE,IJB:IJE,IKB:IKE)**2)**0.5
ENDIF
+!
+ ISEG_LOC(:,:) = 0
+ ZCOORD_TRIG(:,:) = 0.
+ ZCOORD_SEG(:,:) = 0.
+ IPROC_TRIG(:) = 0
+ ISIGNE_EZ(:) = 0
+!
CALL TRIG_POINT
ELSE
GNEW_FLASH_GLOB = .FALSE.
@@ -1194,7 +1381,33 @@ IF (INB_CELL .GE. 1) THEN
!
ZNEUT_POS(:) = 0.
ZNEUT_NEG(:) = 0.
+!
+ IF (LLMA) THEN
+ ZLMA_NEUT_POS(:,:) = 0.
+ ZLMA_NEUT_NEG(:,:) = 0.
+ END IF
END DO ! end loop do while
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 13. COMPUTE THE NOX SOURCE TERM
+! ---------------------------
+!
+ IF (LLNOX_EXPLICIT) THEN
+ IF (IFLASH_COUNT_GLOB .NE. 0) THEN
+ ZCOEF = XMD / XAVOGADRO
+ XLNOX_ECLAIR = 0.
+ IF (IFLASH_COUNT .NE. 0) THEN
+ XLNOX_ECLAIR = SUM(ZLNOX(:,:,:))
+ PSVS_LINOX(:,:,:) = PSVS_LINOX(:,:,:) + ZLNOX(:,:,:) * ZCOEF ! PRHODJ is
+ ! implicit
+ END IF
+ CALL SUM_ELEC_ll (XLNOX_ECLAIR)
+ XLNOX_ECLAIR = XLNOX_ECLAIR / (XAVOGADRO * FLOAT(IFLASH_COUNT_GLOB))
+ END IF
+ DEALLOCATE (ZLNOX)
+ END IF
!
DEALLOCATE (ZNEUT_POS)
DEALLOCATE (ZNEUT_NEG)
@@ -1202,6 +1415,7 @@ IF (INB_CELL .GE. 1) THEN
DEALLOCATE (ZLBDAR)
DEALLOCATE (ZLBDAS)
DEALLOCATE (ZLBDAG)
+ IF (KRR == 7) DEALLOCATE (ZLBDAH)
DEALLOCATE (ZSIGLOB)
DEALLOCATE (ZDQDT)
DEALLOCATE (ZDIST)
@@ -1211,24 +1425,47 @@ IF (INB_CELL .GE. 1) THEN
DEALLOCATE (ISIGNE_EZ)
DEALLOCATE (GNEW_FLASH)
DEALLOCATE (INBSEG)
+ DEALLOCATE (INBSEG_ALL)
DEALLOCATE (INBSEG_LEADER)
DEALLOCATE (INB_FLASH)
DEALLOCATE (INB_FL_REAL)
DEALLOCATE (ZEM_TRIG)
DEALLOCATE (ITYPE)
DEALLOCATE (ISEG_LOC)
+ DEALLOCATE (ZCOORD_TRIG)
DEALLOCATE (ZCOORD_SEG)
+ DEALLOCATE (ZCOORD_SEG_ALL)
+ DEALLOCATE (ISEG_GLOB)
DEALLOCATE (GATTACH)
+ IF(LLMA) THEN
+ DEALLOCATE (ILMA_SEG_ALL)
+ DEALLOCATE (ZLMA_QMT)
+ DEALLOCATE (ZLMA_PRT)
+ DEALLOCATE (ZLMA_NEUT_POS)
+ DEALLOCATE (ZLMA_NEUT_NEG)
+ END IF
END IF ! (inb_cell .ge. 1)
!
!
!-------------------------------------------------------------------------------
!
!* 13. PRINT LIGHTNING INFORMATIONS FOR THE LAST TIMESTEP
+! OR LMA_TIME_SAVE IS REACHED IF LLMA OPTION IS USED
! --------------------------------------------------
-
-IF (NNBLIGHT .NE. 0 .AND. IPROC .EQ. 0 .AND. KTCOUNT .EQ. NSTOP) THEN
+!
+IF (LLMA) THEN
+ IF( IPROC .EQ. 0 .AND. TDTCUR%TIME >= TDTLMA%TIME - PTSTEP ) THEN
+ CALL WRITE_OUT_ASCII
+ CALL WRITE_OUT_LMA
+ ISFLASH_NUMBER(0) = ISFLASH_NUMBER(NNBLIGHT)
+ NNBLIGHT = 0
+ END IF
+END IF
+!
+IF (NNBLIGHT .NE. 0 .AND. ((IPROC .EQ. 0 .AND. OEXIT) .OR. &
+ (KTCOUNT == NSTOP .AND. KMI==1))) THEN
CALL WRITE_OUT_ASCII
+ IF(LLMA) CALL WRITE_OUT_LMA
END IF
!
!
@@ -1237,14 +1474,12 @@ END IF
!* 14. DEALLOCATE
! ----------
!
-DEALLOCATE (ZCELL_GLOB)
DEALLOCATE (ICELL_LOC)
DEALLOCATE (ZQMT)
DEALLOCATE (ZQMTOT)
DEALLOCATE (ZCLOUD)
DEALLOCATE (ZCELL)
DEALLOCATE (ZEMODULE)
-DEALLOCATE (ZQTRANSFER)
!
!
!-------------------------------------------------------------------------------
@@ -1353,15 +1588,21 @@ DO IL = 1, INB_CELL
ISEG_LOC(1,IL) = II_TRIG_LOC
ISEG_LOC(2,IL) = IJ_TRIG_LOC
ISEG_LOC(3,IL) = IK_TRIG
- ZCOORD_SEG(1,IL) = ZXMASS(II_TRIG_LOC)
- ZCOORD_SEG(2,IL) = ZYMASS(IJ_TRIG_LOC)
- ZCOORD_SEG(3,IL) = ZZMASS(II_TRIG_LOC, IJ_TRIG_LOC, IK_TRIG)
+!
+ ISEG_GLOB(1,IL) = II_TRIG_GLOB
+ ISEG_GLOB(2,IL) = IJ_TRIG_GLOB
+ ISEG_GLOB(3,IL) = IK_TRIG
+!
+ ZCOORD_TRIG(1,IL) = ZXMASS(II_TRIG_LOC)
+ ZCOORD_TRIG(2,IL) = ZYMASS(IJ_TRIG_LOC)
+ ZCOORD_TRIG(3,IL) = ZZMASS(II_TRIG_LOC, IJ_TRIG_LOC, IK_TRIG)
+!
+ ZCOORD_SEG(1:3,IL) = ZCOORD_TRIG(1:3,IL)
!
! electric field module at the triggering point
ZEM_TRIG(IL) = ZEMODULE(II_TRIG_LOC,IJ_TRIG_LOC,IK_TRIG)
!
! sign of Ez at the triggering point
-!
ISIGNE_EZ(IL) = 0
IF (PEFIELDW(II_TRIG_LOC,IJ_TRIG_LOC,IK_TRIG) .GT. 0.) THEN
ISIGNE_EZ(IL) = 1
@@ -1373,7 +1614,7 @@ DO IL = 1, INB_CELL
!
! broadcast IFOUND and find the proc where IFOUND = 1
CALL MAX_ELEC_ll (IFOUND, IPROC_TRIG(IL))
-
+!
END DO
!
!
@@ -1381,13 +1622,13 @@ DO IL = 1, INB_CELL
!* 4. BROADCAST USEFULL PARAMETERS
! ----------------------------
!
- CALL MPI_BCAST (ZEM_TRIG(IL), 1,&
+ CALL MPI_BCAST (ZEM_TRIG(IL), 1, &
MPI_PRECISION, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
- CALL MPI_BCAST (ISEG_LOC(:,IL), 3*SIZE(PRT,3), &
+ CALL MPI_BCAST (ISEG_LOC(:,IL), 3*SIZE(PRT,3), &
MPI_INTEGER, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
- CALL MPI_BCAST (ZCOORD_SEG(:,IL), NBRANCH_MAX*3, &
+ CALL MPI_BCAST (ZCOORD_TRIG(:,IL), 3, &
MPI_PRECISION, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
- CALL MPI_BCAST (ISIGNE_EZ(IL), 1,&
+ CALL MPI_BCAST (ISIGNE_EZ(IL), 1, &
MPI_INTEGER, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
!
!
@@ -1440,18 +1681,24 @@ IF (IPROC .EQ. IPROC_TRIG(IL)) THEN
IKBL < IKE .AND. ISTOP .EQ. 0 .AND. &
INBSEG(IL) .LE. (NLEADER_MAX-1))
!
-! local coordinates of the new segment
-!
+! local coordinates of the new segment
IIBL_LOC = ISEG_LOC(1,IL)
IJBL_LOC = ISEG_LOC(2,IL)
IKBL = IKBL + IKSTEP
- IIDECAL = INBSEG(IL)*3
+ IIDECAL = INBSEG(IL) * 3
+!
ISEG_LOC(IIDECAL+1,IL) = IIBL_LOC
ISEG_LOC(IIDECAL+2,IL) = IJBL_LOC
ISEG_LOC(IIDECAL+3,IL) = IKBL
+!
+ ISEG_GLOB(IIDECAL+1,IL) = IIBL_LOC + IXOR - 1
+ ISEG_GLOB(IIDECAL+2,IL) = IJBL_LOC + IYOR - 1
+ ISEG_GLOB(IIDECAL+3,IL) = IKBL
+!
ZCOORD_SEG(IIDECAL+1,IL) = ZXMASS(IIBL_LOC)
ZCOORD_SEG(IIDECAL+2,IL) = ZYMASS(IJBL_LOC)
ZCOORD_SEG(IIDECAL+3,IL) = ZZMASS(IIBL_LOC, IJBL_LOC, IKBL)
+!
INBSEG(IL) = INBSEG(IL) + 1
!
!
@@ -1493,9 +1740,11 @@ IF (IPROC .EQ. IPROC_TRIG(IL)) THEN
PRINT*,'THE LIGHTNING FLASH HAS REACHED THE GROUND '
ISTOP = 1
GCG = .TRUE.
+ NNB_CG = NNB_CG + 1
IF (ISIGN_LEADER > 0) THEN
GCG_POS = .TRUE.
ITYPE(IL) = 3 ! CGP
+ NNB_CG_POS = NNB_CG_POS + 1
ELSE
ITYPE(IL) = 2 ! CGN
END IF
@@ -1511,10 +1760,8 @@ IF (IPROC .EQ. IPROC_TRIG(IL)) THEN
! -------------------------
!
IF (.NOT. GCG) THEN
-! IF ((0.0005*(XZHAT(IKBL)+XZHAT(IKBL+1))) <= XALT_CG .AND. &
IF ( (ZZMASS(IIBL_LOC,IJBL_LOC,IKBL)-PZZ(IIBL_LOC,IJBL_LOC,IKB)) <= &
- XALT_CG .AND. ZCLOUD(IIBL_LOC,IJBL_LOC,IKBL) <= ZCLOUDLIM .AND. &
- INBSEG(IL) .GT. 1 .AND. IKSTEP .LT. 0) THEN
+ XALT_CG .AND. INBSEG(IL) .GT. 1 .AND. IKSTEP .LT. 0) THEN
!
!
!* 2.1 the channel is prolongated to the ground if
@@ -1522,41 +1769,41 @@ IF (IPROC .EQ. IPROC_TRIG(IL)) THEN
!
DO WHILE (IKBL > IKB)
IKBL = IKBL - 1
-
+!
! local coordinates of the new segment
- IIDECAL = INBSEG(IL)*3
+ IIDECAL = INBSEG(IL) * 3
+!
ISEG_LOC(IIDECAL+1,IL) = IIBL_LOC
ISEG_LOC(IIDECAL+2,IL) = IJBL_LOC
ISEG_LOC(IIDECAL+3,IL) = IKBL
+!
+ ISEG_GLOB(IIDECAL+1:IIDECAL+2,IL) = ISEG_GLOB(IIDECAL-2:IIDECAL-1,IL)
+ ISEG_GLOB(IIDECAL+3,IL) = IKBL
+!
ZCOORD_SEG(IIDECAL+1:IIDECAL+2,IL) = ZCOORD_SEG(IIDECAL-2:IIDECAL-1,IL)
ZCOORD_SEG(IIDECAL+3,IL) = ZZMASS(IIBL_LOC, IJBL_LOC, IKBL)
-
+!
! Increment number of segments
INBSEG(IL) = INBSEG(IL) + 1 ! Nb of segments
ZFLASH(IIBL_LOC,IJBL_LOC,IKBL,IL) = 1.
ZCELL(IIBL_LOC,IJBL_LOC,IKBL,IL) = 0.
END DO
+!
+!
+!* 2.2 update the number of CG flashes
!
GCG = .TRUE.
+ NNB_CG = NNB_CG + 1
ISTOP = 1
!
IF (ISIGN_LEADER > 0) THEN
GCG_POS = .TRUE.
+ NNB_CG_POS = NNB_CG_POS + 1
ITYPE(IL) = 3
ELSE
ITYPE(IL) = 2
END IF
END IF
-!
-!
-!* 2.2 update the number of CG flashes
-!
- IF (GCG) THEN
- NNB_CG = NNB_CG + 1
- IF (GCG_POS) THEN
- NNB_CG_POS = NNB_CG_POS + 1
- END IF
- END IF
END IF
END IF ! end if ISTOP=0
END DO ! end loop leader
@@ -1566,13 +1813,9 @@ END IF ! only iproc_trig was working
!* 3. BROADCAST THE INFORMATIONS TO ALL PROCS
! ---------------------------------------
!
-CALL MPI_BCAST (INBSEG(IL), 1, &
- MPI_INTEGER, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
-CALL MPI_BCAST (ISEG_LOC(:,IL), 3*SIZE(PRT,3), &
+CALL MPI_BCAST (ISEG_LOC(:,IL), 3*SIZE(PRT,3), &
MPI_INTEGER, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
-CALL MPI_BCAST (ZCOORD_SEG(:,IL), NBRANCH_MAX*3, &
- MPI_PRECISION, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
-CALL MPI_BCAST (ITYPE(IL), 1,&
+CALL MPI_BCAST (ITYPE(IL), 1, &
MPI_INTEGER, IPROC_TRIG(IL), MPI_COMM_WORLD, IERR)
!
!
@@ -1622,7 +1865,7 @@ DO IL = 1, INB_CELL
END IF
IF (GNEW_FLASH(IL) .AND. IPROC .EQ. IPROC_TRIG(IL)) THEN
DO II = 1, INBSEG(IL)
- IIDECAL = 3*(II-1)
+ IIDECAL = 3 * (II - 1)
IIBL_LOC = ISEG_LOC(IIDECAL+1,IL)
IJBL_LOC = ISEG_LOC(IIDECAL+2,IL)
IKBL = ISEG_LOC(IIDECAL+3,IL)
@@ -1781,13 +2024,13 @@ DO WHILE (IM .LE. IDELTA_IND .AND. ISTOP .NE. 1)
CALL SUM_ELEC_ll (IPT_DIST_GLOB)
!
IF (IPT_DIST_GLOB .LE. INB_SEG_TO_BRANCH) THEN
- IF (IPT_DIST_GLOB .LE. IMAX_BRANCH(IM)) THEN
- GRANDOM = .FALSE.
- ELSE
- GRANDOM = .TRUE.
- END IF
+ IF (IPT_DIST_GLOB .LE. IMAX_BRANCH(IM)) THEN
+ GRANDOM = .FALSE.
+ ELSE
+ GRANDOM = .TRUE.
+ END IF
ELSE
- GRANDOM = .TRUE.
+ GRANDOM = .TRUE.
END IF
!
!
@@ -1796,18 +2039,16 @@ DO WHILE (IM .LE. IDELTA_IND .AND. ISTOP .NE. 1)
!
IF (IPT_DIST_GLOB .GT. 0 .AND. INB_SEG_TO_BRANCH .NE. 0) THEN
IF (.NOT. GRANDOM) THEN
-
- INB_SEG_TO_BRANCH = INB_SEG_TO_BRANCH - IPT_DIST_GLOB
+ INB_SEG_TO_BRANCH = INB_SEG_TO_BRANCH - IPT_DIST_GLOB
!
!* 2.1 all points are selected
!
- IF(IPT_DIST .GT. 0) THEN
- WHERE (IMASKQ_DIST(IIB:IIE,IJB:IJE,IKB:IKE) .EQ. IM)
- ZFLASH(IIB:IIE,IJB:IJE,IKB:IKE,IL) = 2.
- ZCELL(IIB:IIE,IJB:IJE,IKB:IKE,IL) = 0.
- END WHERE
- END IF
-
+ IF(IPT_DIST .GT. 0) THEN
+ WHERE (IMASKQ_DIST(IIB:IIE,IJB:IJE,IKB:IKE) .EQ. IM)
+ ZFLASH(IIB:IIE,IJB:IJE,IKB:IKE,IL) = 2.
+ ZCELL(IIB:IIE,IJB:IJE,IKB:IKE,IL) = 0.
+ END WHERE
+ END IF
ELSE
!
!* 2.2 the gridpoints are chosen randomly
@@ -1878,7 +2119,7 @@ DO WHILE (IM .LE. IDELTA_IND .AND. ISTOP .NE. 1)
ZWORK(:,:,:) = 0.
ZWORK(:,:,:) = UNPACK(ZAUX(:), MASK=(IMASKQ_DIST(:,:,:).EQ.IM), FIELD=0.0)
WHERE (ZWORK(IIB:IIE,IJB:IJE,IKB:IKE) .EQ. 1.)
- ZFLASH(IIB:IIE,IJB:IJE,IKB:IKE,IL) = 2.
+ ZFLASH(IIB:IIE,IJB:IJE,IKB:IKE,IL) = 2.
ZCELL(IIB:IIE,IJB:IJE,IKB:IKE,IL) = 0.
END WHERE
DEALLOCATE (ZVECT)
@@ -1894,36 +2135,220 @@ DO WHILE (IM .LE. IDELTA_IND .AND. ISTOP .NE. 1)
! next distance
IM = IM + 1
END DO ! end loop / do while / radius IM
-
+!
INB_SEG_AFT = COUNT (ZFLASH(IIB:IIE,IJB:IJE,IKB:IKE,IL) .NE. 0.)
CALL SUM_ELEC_ll(INB_SEG_AFT)
+!
IF (INB_SEG_AFT .GT. INB_SEG_BEF) THEN
- DO IPLOOP = 0, NPROC-1
- IF (IPROC .EQ. IPLOOP) THEN
- DO II = IIB, IIE
- DO IJ = IJB, IJE
- DO IK = IKB, IKE
- IF (ZFLASH(II,IJ,IK,IL) .EQ. 2.) THEN
- IIDECALB = INBSEG(IL)*3
- ZCOORD_SEG(IIDECALB+1,IL) = ZXMASS(II)
- ZCOORD_SEG(IIDECALB+2,IL) = ZYMASS(IJ)
- ZCOORD_SEG(IIDECALB+3,IL) = ZZMASS(II, IJ, IK)
- INBSEG(IL) = INBSEG(IL) + 1
- END IF
- END DO
- END DO
+ DO II = IIB, IIE
+ DO IJ = IJB, IJE
+ DO IK = IKB, IKE
+ IF (ZFLASH(II,IJ,IK,IL) .EQ. 2.) THEN
+ IIDECALB = INBSEG(IL) * 3
+!
+ ISEG_GLOB(IIDECALB+1,IL) = II + IXOR - 1
+ ISEG_GLOB(IIDECALB+2,IL) = IJ + IYOR - 1
+ ISEG_GLOB(IIDECALB+3,IL) = IK
+!
+ ZCOORD_SEG(IIDECALB+1,IL) = ZXMASS(II)
+ ZCOORD_SEG(IIDECALB+2,IL) = ZYMASS(IJ)
+ ZCOORD_SEG(IIDECALB+3,IL) = ZZMASS(II,IJ,IK)
+ INBSEG(IL) = INBSEG(IL) + 1
+ END IF
END DO
- END IF
- CALL MPI_BCAST (INBSEG(IL), 1,&
- MPI_INTEGER, IPLOOP, MPI_COMM_WORLD, IERR)
- CALL MPI_BCAST (ZCOORD_SEG(:,IL), NBRANCH_MAX*3, &
- MPI_PRECISION, IPLOOP, MPI_COMM_WORLD, IERR)
+ END DO
END DO
END IF
!
END SUBROUTINE BRANCH_GEOM
!
-!------- ------------------------------------------------------------------------
+!--------------------------------------------------------------------------------
+!
+ SUBROUTINE GATHER_ALL_BRANCH
+!
+!!
+!! Purpose:
+!!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+IMPLICIT NONE
+!
+INTEGER :: INSEGPROC, INSEGCELL ! number of segments in the process,
+ ! and number of segments in the cell
+INTEGER :: ISAVEDECAL
+INTEGER :: INSEGTRIG, IPROCTRIG
+REAL, DIMENSION(:), ALLOCATABLE :: ZLMAQMT, ZLMAPRT, ZLMAPOS, ZLMANEG
+REAL, DIMENSION(:), ALLOCATABLE :: ZSEND, ZRECV
+INTEGER, DIMENSION(:), ALLOCATABLE :: ISEND, IRECV
+INTEGER, DIMENSION(NPROC) :: IDECAL, IDECAL3, IDECALN
+INTEGER, DIMENSION(NPROC) :: INBSEG_PROC_X3, INBSEG_PROC_XNSV
+!
+!
+IPROCTRIG = IPROC_TRIG(IL)
+INSEGCELL = INBSEG_ALL(IL)
+INSEGPROC = INBSEG_PROC(IPROC+1)
+INSEGTRIG = INBSEG_PROC(IPROCTRIG+1)
+!
+IDECAL(1) = INSEGTRIG
+DO IK = 2, NPROC
+ IDECAL(IK) = IDECAL(IK-1) + INBSEG_PROC(IK-1)
+END DO
+!
+IF(IPROCTRIG .EQ. 0) ISAVEDECAL = IDECAL(IPROCTRIG+1)
+!
+IDECAL(IPROCTRIG+1) = 0
+DO IK = IPROCTRIG+2, NPROC
+ IF(IPROCTRIG .EQ. 0) THEN
+ IDECAL(IK) = IDECAL(IK) - ISAVEDECAL
+ ELSE
+ IDECAL(IK) = IDECAL(IK) - IDECAL(1)
+ END IF
+END DO
+!
+IDECAL3(:) = 3 * IDECAL(:)
+!
+!
+!* 1. BRANCH COORDINATES
+!
+ALLOCATE (ZRECV(INSEGCELL*3))
+ALLOCATE (ZSEND(INSEGPROC*3))
+!
+IF (INSEGPROC .NE. 0) THEN
+ ZSEND(1:3*INSEGPROC) = ZCOORD_SEG(1:3*INSEGPROC,IL)
+END IF
+!
+IF (IPROC .EQ. 0) THEN
+ INBSEG_PROC_X3(:) = 3 * INBSEG_PROC(:)
+END IF
+!
+CALL MPI_GATHERV (ZSEND, 3*INSEGPROC, MPI_PRECISION, ZRECV, INBSEG_PROC_X3, &
+ IDECAL3, MPI_PRECISION, 0, MPI_COMM_WORLD, IERR)
+!
+IF (IPROC .EQ. 0) THEN
+ ZCOORD_SEG_ALL(1:3*INSEGCELL,IL) = ZRECV(1:3*INSEGCELL)
+END IF
+!
+DEALLOCATE (ZRECV)
+DEALLOCATE (ZSEND)
+!
+!
+!* 2. FOR LMA-LIKE RESULTS: Charge, mixing ratio,
+!* neutralized positive/negative charge
+!* and grid index
+!
+IF (LLMA) THEN
+ ALLOCATE (ISEND(3*INSEGPROC))
+ ALLOCATE (ZLMAQMT(INSEGPROC*NSV_ELEC))
+ ALLOCATE (ZLMAPRT(INSEGPROC*NSV_ELEC))
+ ALLOCATE (ZLMAPOS(INSEGPROC))
+ ALLOCATE (ZLMANEG(INSEGPROC))
+!
+ ISEND (:) = 0
+ ZLMAPOS(:) = 0.
+ ZLMANEG(:) = 0.
+ ZLMAQMT(:) = 0.
+ ZLMAPRT(:) = 0.
+!
+ IF (INSEGPROC .NE. 0) THEN
+ DO II = 1, INSEGPROC
+ IM = 3 * (II - 1)
+ IX = ISEG_GLOB(IM+1,IL) - IXOR + 1
+ IY = ISEG_GLOB(IM+2,IL) - IYOR + 1
+ IZ = ISEG_GLOB(IM+3,IL)
+!
+ IM = NSV_ELEC * (II - 1)
+ IF (IX .LE. IIE .AND. IX .GE. IIB .AND. &
+ IY .LE. IJE .AND. IY .GE. IJB) THEN
+ ZLMAQMT(IM+2:IM+6) = ZQMT(IX,IY,IZ,2:6)
+ ZLMAPRT(IM+2:IM+6) = PRT(IX,IY,IZ,2:6)
+ DO IJ = 1, NSV_ELEC
+ IF (ZDQDT(IX,IY,IZ,IJ) .GT. 0.) THEN
+ ZLMAPOS(II) = ZLMAPOS(II) + &
+ ZDQDT(IX,IY,IZ,IJ) * PRHODJ(IX,IY,IZ)
+ ELSE IF (ZDQDT(IX,IY,IZ,IJ) .LT. 0.) THEN
+ ZLMANEG(II) = ZLMANEG(II) + &
+ ZDQDT(IX,IY,IZ,IJ) * PRHODJ(IX,IY,IZ)
+ END IF
+ END DO
+ END IF
+ END DO
+!
+ ISEND(1:3*INSEGPROC) = ISEG_GLOB(1:3*INSEGPROC, IL)
+ END IF
+!
+! Grid Indexes
+!
+ ALLOCATE (IRECV(3*INSEGCELL))
+!
+ CALL MPI_GATHERV (ISEND, 3*INSEGPROC, MPI_INTEGER, IRECV, INBSEG_PROC_X3, &
+ IDECAL3, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
+!
+ IF (IPROC .EQ. 0) THEN
+ ILMA_SEG_ALL(1:3*INSEGCELL,IL) = IRECV(1:3*INSEGCELL)
+ END IF
+!
+ DEALLOCATE (IRECV)
+ DEALLOCATE (ISEND)
+!
+! Neutralized charge at grid points
+!
+ ALLOCATE (ZRECV(INSEGCELL))
+!
+ CALL MPI_GATHERV (ZLMAPOS, INSEGPROC, MPI_PRECISION, ZRECV, INBSEG_PROC, &
+ IDECAL, MPI_PRECISION, 0, MPI_COMM_WORLD, IERR)
+!
+ IF (IPROC .EQ. 0) THEN
+ ZLMA_NEUT_POS(1:INSEGCELL,IL) = ZRECV(1:INSEGCELL)
+ END IF
+!
+ CALL MPI_GATHERV (ZLMANEG, INSEGPROC, MPI_PRECISION, ZRECV, INBSEG_PROC, &
+ IDECAL, MPI_PRECISION, 0, MPI_COMM_WORLD, IERR)
+!
+ IF (IPROC .EQ. 0) THEN
+ ZLMA_NEUT_NEG(1:INSEGCELL,IL) = ZRECV(1:INSEGCELL)
+ END IF
+!
+ DEALLOCATE (ZLMAPOS)
+ DEALLOCATE (ZLMANEG)
+ DEALLOCATE (ZRECV)
+!
+! Charge and mixing ratios at neutralized points
+!
+ ALLOCATE (ZRECV(NSV_ELEC*INSEGCELL))
+!
+ IDECALN(:) = IDECAL(:) * NSV_ELEC
+!
+ IF (IPROC .EQ. 0) THEN
+ INBSEG_PROC_XNSV(:) = NSV_ELEC * INBSEG_PROC(:)
+ END IF
+!
+ CALL MPI_GATHERV (ZLMAQMT, NSV_ELEC*INSEGPROC, MPI_PRECISION, ZRECV, &
+ INBSEG_PROC_XNSV, &
+ IDECALN, MPI_PRECISION, 0, MPI_COMM_WORLD, IERR )
+!
+ IF (IPROC .EQ. 0) THEN
+ ZLMA_QMT(1:NSV_ELEC*INSEGCELL,IL) = ZRECV(1:NSV_ELEC*INSEGCELL)
+ END IF
+!
+ CALL MPI_GATHERV (ZLMAPRT, NSV_ELEC*INSEGPROC, MPI_PRECISION, ZRECV, &
+ INBSEG_PROC_XNSV, &
+ IDECALN, MPI_PRECISION, 0, MPI_COMM_WORLD, IERR)
+!
+ IF (IPROC .EQ. 0) THEN
+ ZLMA_PRT(1:NSV_ELEC*INSEGCELL,IL) = ZRECV(1:NSV_ELEC*INSEGCELL)
+ END IF
+!
+ DEALLOCATE (ZLMAQMT)
+ DEALLOCATE (ZLMAPRT)
+ DEALLOCATE (ZRECV)
+!
+END IF
+!
+END SUBROUTINE GATHER_ALL_BRANCH
+!
+!--------------------------------------------------------------------------------
!
SUBROUTINE PT_DISCHARGE
!
@@ -1938,24 +2363,23 @@ IMPLICIT NONE
!
!
WHERE (ABS(PEFIELDW(:,:,IKB)) > XECORONA .AND. PEFIELDW(:,:,IKB) > 0.)
- PRSVS(:,:,IKB,1) = PRSVS(:,:,IKB,1) + &
+ PRSVS(:,:,IKB,1) = PRSVS(:,:,IKB,1) + &
XFCORONA * PEFIELDW(:,:,IKB) * (ABS(PEFIELDW(:,:,IKB)) - &
XECORONA)**2 / (PZZ(:,:,IKB+1) - PZZ(:,:,IKB))
ENDWHERE
-
+!
WHERE (ABS(PEFIELDW(:,:,IKB)) > XECORONA .AND. PEFIELDW(:,:,IKB) < 0.)
PRSVS(:,:,IKB,NSV_ELEC) = PRSVS(:,:,IKB,NSV_ELEC) + &
XFCORONA * PEFIELDW(:,:,IKB) * (ABS(PEFIELDW(:,:,IKB)) - &
XECORONA)**2 / (PZZ(:,:,IKB+1) - PZZ(:,:,IKB))
ENDWHERE
-
+!
END SUBROUTINE PT_DISCHARGE
!
-
-!-----------------------------------------------------------------------------------
+!----------------------------------------------------------------------------------
!
SUBROUTINE WRITE_OUT_ASCII
-
+!
!!
!! Purpose:
!!
@@ -1969,27 +2393,51 @@ INTEGER :: I1, I2
!
!
!* 1. FLASH PARAMETERS
-! ---------------
+! ----------------
+!
OPEN (UNIT=NLU_fgeom_diag, FILE=CEXP//"_fgeom_diag.asc", ACTION="WRITE", &
STATUS="OLD", FORM="FORMATTED", POSITION="APPEND")
!
! Ecriture ascii dans CEXP//'_fgeom_diag.asc" defini dans RESOLVED_ELEC
!
-DO I1 = 1, NNBLIGHT
- WRITE (NLU_fgeom_diag,FMT='(I8,F9.1,I4,I6,I4,I6,F9.3,F9.3,F9.3,F7.3,F8.2,F9.2,f9.4)') &
- ISFLASH_NUMBER(I1), &
- ISTCOUNT_NUMBER(I1) * PTSTEP, &
- ISCELL_NUMBER(I1), &
- ISNB_FLASH(I1), &
- ISTYPE(I1), &
- ISNBSEG(I1), &
- ZSEM_TRIG(I1), &
- ZSCOORD_SEG(I1,1,1)*1.E-3, &
- ZSCOORD_SEG(I1,1,2)*1.E-3, &
- ZSCOORD_SEG(I1,1,3)*1.E-3, &
- ZSNEUT_POS(I1), &
- ZSNEUT_NEG(I1), ZSNEUT_POS(I1)+ZSNEUT_NEG(I1)
-ENDDO
+IF (LCARTESIAN) THEN
+ DO I1 = 1, NNBLIGHT
+ WRITE (UNIT=NLU_fgeom_diag,FMT='(I8,F9.1,I4,I6,I4,I6,F9.3,F12.3,F12.3,F9.3,F8.2,F9.2,f9.4)') &
+ ISFLASH_NUMBER(I1), &
+ ISTCOUNT_NUMBER(I1) * PTSTEP, &
+ ISCELL_NUMBER(I1), &
+ ISNB_FLASH(I1), &
+ ISTYPE(I1), &
+ ISNBSEG(I1), &
+ ZSEM_TRIG(I1), &
+ ZSCOORD_SEG(I1,1,1)*1.E-3, &
+ ZSCOORD_SEG(I1,1,2)*1.E-3, &
+ ZSCOORD_SEG(I1,1,3)*1.E-3, &
+ ZSNEUT_POS(I1), &
+ ZSNEUT_NEG(I1), ZSNEUT_POS(I1)+ZSNEUT_NEG(I1)
+ END DO
+ELSE
+ DO I1 = 1, NNBLIGHT
+! compute latitude and longitude of the triggering point
+ CALL SM_LATLON(XLATORI,XLONORI,ZSCOORD_SEG(I1,1,1),&
+ ZSCOORD_SEG(I1,1,2),&
+ ZLAT,ZLON)
+!
+ WRITE (UNIT=NLU_fgeom_diag,FMT='(I8,F9.1,I4,I6,I4,I6,F9.3,F12.3,F12.3,F9.3,F8.2,F9.2,f9.4)') &
+ ISFLASH_NUMBER(I1), &
+ ISTCOUNT_NUMBER(I1) * PTSTEP, &
+ ISCELL_NUMBER(I1), &
+ ISNB_FLASH(I1), &
+ ISTYPE(I1), &
+ ISNBSEG(I1), &
+ ZSEM_TRIG(I1), &
+ ZLAT, &
+ ZLON, &
+ ZSCOORD_SEG(I1,1,3)*1.E-3, &
+ ZSNEUT_POS(I1), &
+ ZSNEUT_NEG(I1), ZSNEUT_POS(I1)+ZSNEUT_NEG(I1)
+ END DO
+END IF
!
CLOSE (UNIT=NLU_fgeom_diag)
!
@@ -1998,6 +2446,7 @@ CLOSE (UNIT=NLU_fgeom_diag)
! -------------------------
!
IF (LSAVE_COORD) THEN
+!
! Ecriture ascii dans CEXP//'_fgeom_coord.asc" defini dans RESOLVED_ELEC
!
OPEN (UNIT=NLU_fgeom_coord, FILE=CEXP//"_fgeom_coord.asc", ACTION="WRITE", &
@@ -2022,6 +2471,58 @@ END SUBROUTINE WRITE_OUT_ASCII
!
!-------------------------------------------------------------------------------
!
+SUBROUTINE WRITE_OUT_LMA
+!
+!!
+!! Purpose:
+!!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+IMPLICIT NONE
+!
+INTEGER :: I1, I2
+!
+!
+!* 1. LMA SIMULATOR
+! -------------
+!
+CALL SM_LATLON(XLATORI,XLONORI,ZSCOORD_SEG(:,:,1),ZSCOORD_SEG(:,:,2), &
+ ZLMA_LAT(:,:),ZLMA_LON(:,:))
+DO I1 = 1, NNBLIGHT
+ DO I2 = 1, ISNBSEG(I1)
+ WRITE (ILMA_UNIT,FMT='(I6,F12.1,I6,2(F15.6),3(F15.3),3(I6),12(E15.4))') &
+ ISFLASH_NUMBER(I1), &
+ ISTCOUNT_NUMBER(I1) * PTSTEP, &
+ ISTYPE(I1), &
+ ZLMA_LAT(I1,I2), &
+ ZLMA_LON(I1,I2), &
+ ZSCOORD_SEG(I1,I2,1)*1.E-3, &
+ ZSCOORD_SEG(I1,I2,2)*1.E-3, &
+ ZSCOORD_SEG(I1,I2,3)*1.E-3, &
+ ISLMA_SEG_GLOB(I1,I2,1), &
+ ISLMA_SEG_GLOB(I1,I2,2), &
+ ISLMA_SEG_GLOB(I1,I2,3), &
+ ZSLMA_PRT(I1,I2,2), &
+ ZSLMA_PRT(I1,I2,3), &
+ ZSLMA_PRT(I1,I2,4), &
+ ZSLMA_PRT(I1,I2,5), &
+ ZSLMA_PRT(I1,I2,6), &
+ ZSLMA_QMT(I1,I2,2), &
+ ZSLMA_QMT(I1,I2,3), &
+ ZSLMA_QMT(I1,I2,4), &
+ ZSLMA_QMT(I1,I2,5), &
+ ZSLMA_QMT(I1,I2,6), &
+ ZSLMA_NEUT_POS(I1,I2), &
+ ZSLMA_NEUT_NEG(I1,I2)
+ END DO
+END DO
+!
+END SUBROUTINE WRITE_OUT_LMA
+!
+!-------------------------------------------------------------------------------
+!
END SUBROUTINE FLASH_GEOM_ELEC_n
!
!-------------------------------------------------------------------------------
diff --git a/src/MNH/gamma.f90 b/src/MNH/gamma.f90
index 53bae0f08..bec3e8412 100644
--- a/src/MNH/gamma.f90
+++ b/src/MNH/gamma.f90
@@ -1,16 +1,6 @@
-!MNH_LIC Copyright 1994-2014 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.
-!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 operators 2006/05/18 13:07:25
-!-----------------------------------------------------------------
-! #################
+!########################
MODULE MODI_GAMMA
-! #################
+!########################
!
INTERFACE GAMMA
!
@@ -20,14 +10,19 @@ REAL :: PGAMMA
END FUNCTION GAMMA_X0D
!
FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA)
-REAL, DIMENSION(:), INTENT(IN) :: PX
-REAL, DIMENSION(SIZE(PX)) :: PGAMMA
+REAL, DIMENSION(:), INTENT(IN) :: PX
+REAL, DIMENSION(SIZE(PX)) :: PGAMMA
END FUNCTION GAMMA_X1D
!
-END INTERFACE GAMMA
-!
+END INTERFACE
END MODULE MODI_GAMMA
!
+!--------------------------------------------------------------------------
+!
+!
+!* 1. FUNCTION GAMMA FOR SCALAR VARIABLE
+!
+!
! ######################################
FUNCTION GAMMA_X0D(PX) RESULT(PGAMMA)
! ######################################
@@ -63,11 +58,9 @@ END MODULE MODI_GAMMA
!!
!! MODIFICATIONS
!! -------------
-!! Original 7/12/95
+!! Original 7/11/95
!! C. Barthe 9/11/09 add a function for 1D arguments
-!!
-!!-------------------------------------------------------------------------------
-
+!
!* 0. DECLARATIONS
! ------------
!
@@ -130,6 +123,10 @@ END FUNCTION GAMMA_X0D
!
!-------------------------------------------------------------------------------
!
+!
+!* 1. FUNCTION GAMMA FOR 1D ARRAY
+!
+!
! ######################################
FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA)
! ######################################
@@ -176,14 +173,14 @@ IMPLICIT NONE
!
!* 0.1 declarations of arguments and result
!
-REAL,DIMENSION(:), INTENT(IN) :: PX
-REAL,DIMENSION(SIZE(PX)) :: PGAMMA
+REAL, DIMENSION(:), INTENT(IN) :: PX
+REAL, DIMENSION(SIZE(PX)) :: PGAMMA
!
!* 0.2 declarations of local variables
!
INTEGER :: JJ ! Loop index
-INTEGER :: JI ! Loop index
-REAL :: ZSER,ZSTP,ZTMP,ZX,ZY,ZCOEF(6)
+REAL, DIMENSION(SIZE(PX)) :: ZSER,ZSTP,ZTMP,ZX,ZY
+REAL :: ZCOEF(6)
REAL :: ZPI
!
!-------------------------------------------------------------------------------
@@ -200,34 +197,24 @@ ZCOEF(6) = -0.5395239384953E-5
ZSTP = 2.5066282746310005
!
ZPI = 3.141592654
-!
-!-------------------------------------------------------------------------------
-!
-!* 2. COMPUTE GAMMA
-! -------------
-!
-DO JI = 1, SIZE(PX)
- IF (PX(JI) .LT. 0.) THEN
- ZX = 1. - PX(JI)
- ELSE
- ZX = PX(JI)
- END IF
- ZY = ZX
- ZTMP = ZX + 5.5
- ZTMP = (ZX + 0.5) * ALOG(ZTMP) - ZTMP
- ZSER = 1.000000000190015
-!
- DO JJ = 1, 6
- ZY = ZY + 1.0
- ZSER = ZSER + ZCOEF(JJ) / ZY
- END DO
-!
- IF (PX(JI) .LT. 0.) THEN
- PGAMMA = ZPI / SIN(ZPI*PX(JI)) / EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
- ELSE
- PGAMMA = EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
- END IF
+ZX(:) = PX(:)
+WHERE ( PX(:)<0.0 )
+ ZX(:) = 1.- PX(:)
+END WHERE
+ZY(:) = ZX(:)
+ZTMP(:) = ZX(:) + 5.5
+ZTMP(:) = (ZX(:) + 0.5)*ALOG(ZTMP(:)) - ZTMP(:)
+ZSER(:) = 1.000000000190015
+!
+DO JJ = 1 , 6
+ ZY(:) = ZY(:) + 1.0
+ ZSER(:) = ZSER(:) + ZCOEF(JJ)/ZY(:)
END DO
+!
+PGAMMA(:) = EXP( ZTMP(:) + ALOG( ZSTP*ZSER(:)/ZX(:) ) )
+WHERE ( PX(:)<0.0 )
+ PGAMMA(:) = ZPI/SIN(ZPI*PX(:))/PGAMMA(:)
+END WHERE
RETURN
!
END FUNCTION GAMMA_X1D
diff --git a/src/MNH/ini_elecn.f90 b/src/MNH/ini_elecn.f90
index db0dc1218..5ef8d1ad7 100644
--- a/src/MNH/ini_elecn.f90
+++ b/src/MNH/ini_elecn.f90
@@ -73,6 +73,7 @@ END MODULE MODI_INI_ELEC_n
!! J.-P. Pinty 13/04/12 Add elec_trid to initialise the tridiagonal syst.
!! J.-P. Pinty 01/07/12 Add a non-homogeneous Neuman fair-weather
!! boundary condition at the top
+!! J.-P. Pinty 15/11/13 Initialize the flash maps
!!
!-------------------------------------------------------------------------------
!
@@ -87,7 +88,7 @@ USE MODE_FMREAD
!
USE MODD_DIM_n, ONLY : NIMAX_ll, NJMAX_ll
USE MODD_ELEC_DESCR
-USE MODD_ELEC_n, ONLY : XRHOM_E, XAF_E, XCF_E, XBFY_E
+USE MODD_ELEC_n, ONLY : XRHOM_E, XAF_E, XCF_E, XBFY_E, XBFB_E, XBF_SXP2_YP1_Z_E
USE MODD_CONF_n, ONLY : NRR
USE MODD_PARAMETERS, ONLY : JPVEXT, JPHEXT
USE MODD_CST
@@ -95,6 +96,7 @@ USE MODD_CONF, ONLY : CEQNSYS,CCONF,CPROGRAM
USE MODD_DYN
USE MODD_REF
USE MODD_TIME
+USE MODD_ELEC_FLASH
USE MODD_GET_n, ONLY : CGETINPRC, CGETINPRR, CGETINPRS, CGETINPRG, CGETINPRH, &
CGETCLOUD, CGETSVT
USE MODD_PRECIP_n, ONLY : XINPRR, XACPRR, XINPRS, XACPRS, XINPRG, XACPRG, &
@@ -102,8 +104,8 @@ USE MODD_PRECIP_n, ONLY : XINPRR, XACPRR, XINPRS, XACPRS, XINPRG, XACPRG, &
USE MODD_CLOUDPAR_n, ONLY : NSPLITR
USE MODD_REF_n, ONLY : XRHODJ, XTHVREF
USE MODD_GRID_n, ONLY : XMAP, XDXHAT, XDYHAT
-USE MODD_DYN_n, ONLY : XRHOM, XTRIGSX, XTRIGSY, XAF, XCF, XBFY, XDXHATM, &
- XDYHATM, NIFAXX, NIFAXY
+USE MODD_DYN_n, ONLY : XRHOM, XTRIGSX, XTRIGSY, XAF, XCF, XBFY, XBFB, XDXHATM, &
+ XDYHATM, NIFAXX, NIFAXY, XBF_SXP2_YP1_Z
USE MODD_LBC_n, ONLY : CLBCX, CLBCY
USE MODI_ELEC_TRID
@@ -113,6 +115,7 @@ USE MODI_INI_RAIN_ICE_ELEC
USE MODI_INI_FIELD_ELEC
USE MODI_INI_FLASH_GEOM_ELEC
USE MODI_READ_PRECIP_FIELD
+USE MODI_ELEC_TRIDZ
!
!
IMPLICIT NONE
@@ -291,14 +294,27 @@ ALLOCATE (XRHOM_E(SIZE(XRHOM)))
ALLOCATE (XAF_E(SIZE(XAF)))
ALLOCATE (XCF_E(SIZE(XCF)))
ALLOCATE (XBFY_E(SIZE(XBFY,1),SIZE(XBFY,2),SIZE(XBFY,3)))
+ALLOCATE (XBFB_E(SIZE(XBFB,1),SIZE(XBFB,2),SIZE(XBFB,3)))
+ALLOCATE (XBF_SXP2_YP1_Z_E(SIZE(XBF_SXP2_YP1_Z,1),SIZE(XBF_SXP2_YP1_Z,2),&
+ SIZE(XBF_SXP2_YP1_Z,3)))
!
-CALL ELEC_TRID (HLUOUT,CLBCX,CLBCY, &
- XMAP,XDXHAT,XDYHAT,XDXHATM,XDYHATM,XRHOM_E,XAF_E, &
+CALL ELEC_TRIDZ (HLUOUT,CLBCX,CLBCY, &
+ XMAP,XDXHAT,XDYHAT,XDXHATM,XDYHATM,XRHOM_E,XAF_E, &
XCF_E,XTRIGSX,XTRIGSY,NIFAXX,NIFAXY, &
- XRHODJ,XTHVREF,PZZ,XBFY_E,XEPOTFW_TOP)
+ XRHODJ,XTHVREF,PZZ,XBFY_E,XEPOTFW_TOP, &
+ XBFB_E,XBF_SXP2_YP1_Z_E)
!
CEQNSYS=YEQNSYS
!
+!* 3.7 initialize the flash maps
+!
+ALLOCATE( NMAP_TRIG_IC(IIU,IJU) ); NMAP_TRIG_IC(:,:) = 0
+ALLOCATE( NMAP_IMPACT_CG(IIU,IJU) ); NMAP_IMPACT_CG(:,:) = 0
+ALLOCATE( NMAP_2DAREA_IC(IIU,IJU) ); NMAP_2DAREA_IC(:,:) = 0
+ALLOCATE( NMAP_2DAREA_CG(IIU,IJU) ); NMAP_2DAREA_CG(:,:) = 0
+ALLOCATE( NMAP_3DIC(IIU,IJU,IKU) ); NMAP_3DIC(:,:,:) = 0
+ALLOCATE( NMAP_3DCG(IIU,IJU,IKU) ); NMAP_3DCG(:,:,:) = 0
+!
!-------------------------------------------------------------------------------
!
!
diff --git a/src/MNH/ini_flash_geom_elec.f90 b/src/MNH/ini_flash_geom_elec.f90
index a4c289cf1..b3dae9562 100644
--- a/src/MNH/ini_flash_geom_elec.f90
+++ b/src/MNH/ini_flash_geom_elec.f90
@@ -45,6 +45,9 @@ END MODULE MODI_INI_FLASH_GEOM_ELEC
!! -------------
!! Original 29/11/02
!!
+!! Modifications
+!! J.-P. Pinty jan 2015 : add LMA simulator
+!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -55,6 +58,8 @@ USE MODD_RAIN_ICE_DESCR
USE MODD_ELEC_DESCR
USE MODD_ELEC_PARAM
USE MODD_DIM_n, ONLY : NKMAX
+USE MODD_TIME_n, ONLY : TDTCUR
+USE MODD_LMA_SIMULATOR, ONLY : LLMA, TDTLMA, LWRITE_LMA, XDTLMA, CLMA_FILE
!
USE MODI_MOMG
!
@@ -105,7 +110,7 @@ XE_THRESH = 35.E3 ! (V/m)
NLEADER_MAX = NKMAX
!
! the maximum number of branches is arbitriraly set to 5000
-NBRANCH_MAX = 5000
+NBRANCH_MAX = 50000
!
! the maximum number of electrified cells in the domain is arbitrarily
! set to 10
@@ -113,7 +118,7 @@ NMAX_CELL = 10
!
! the altitude for CG to be prolongated to the ground is set to 2 km
! this threshold could be modified once ions will be taken into account
-XALT_CG = 2000. ! m
+XALT_CG = 2000. ! m
!
!
!----------------------------------------------------------------------------
@@ -128,4 +133,15 @@ NNB_CG_POS = 0
!
!----------------------------------------------------------------------------
!
+!* 4. INITIALIZE LMA RECORDS
+! ----------------------
+!
+! needs LLMA = .TRUE. to operate
+XDTLMA = 600.
+TDTLMA = TDTCUR
+LWRITE_LMA = .FALSE.
+CLMA_FILE(1:5) = "BEGIN"
+!
+!----------------------------------------------------------------------------
+!
END SUBROUTINE INI_FLASH_GEOM_ELEC
diff --git a/src/MNH/ini_param_elec.f90 b/src/MNH/ini_param_elec.f90
index 59682834f..3b7808bba 100644
--- a/src/MNH/ini_param_elec.f90
+++ b/src/MNH/ini_param_elec.f90
@@ -80,6 +80,7 @@ END MODULE MODI_INI_PARAM_ELEC
!! M. Chong 26/01/10 Small ions parameters
!! +Fair weather field from Helsdon-Farley
!! (JGR, 1987, 5661-5675)
+!! J.-P. Pinty jan 2015 tabulate the equations for Saunders
!!
!-------------------------------------------------------------------------------
!
@@ -136,6 +137,9 @@ CHARACTER (LEN=100) :: YCOMMENT
CHARACTER (LEN=16) :: YRECFM
CHARACTER (LEN=2) :: YDIR
!
+INTEGER :: JLWC, JTEMP
+REAL, DIMENSION(:), ALLOCATABLE :: ZT, ZLWCC, ZEW
+!
!-------------------------------------------------------------------------------
! constants for electricity
!
@@ -383,29 +387,65 @@ END IF
!* 8.2 Saunders et al. (1991) and
!* Saunders and Peck (1998) parameterizations
!
-IF (CNI_CHARGING == 'SAUN1' .OR. CNI_CHARGING == 'SAUN2' .OR. &
- CNI_CHARGING == 'SAP98') THEN
+IF (CNI_CHARGING == 'SAUN1' .OR. CNI_CHARGING == 'SAUN2' .OR. &
+ CNI_CHARGING == 'SAP98' .OR. &
+ CNI_CHARGING == 'BSMP1' .OR. CNI_CHARGING == 'BSMP2' .OR. &
+ CNI_CHARGING == 'TEEWC' .OR. CNI_CHARGING == 'TERAR') THEN
!
! ice particle = the smallest particle (I-S and I-G collisions)
- XIMP = 3.76
+ XIMP = 3.76 ! for positive charge
XINP = 2.5
XIKP = 4.92E13
- XIMN = 2.54
+ XIKP_TAK = 6.1E12 ! for Takahashi
+ XIMN = 2.54 ! for negative charge
XINN = 2.8
XIKN = 5.25E8
+ XIKN_TAK = 4.3E7 ! for Takahashi
!
! snow = the smallest particle (S-G collisions)
- XSMP = 0.44
+ XSMP = 0.44 ! for positive charge
XSNP = 2.5
XSKP = 52.8
- XSMN = 0.5
+ XSKP_TAK = 6.5 ! for Takahashi
+ XSMN = 0.5 ! for negative charge
XSNN = 2.8
- XSKN = 24
+ XSKN = 24.
+ XSKN_TAK = 2.0 ! for Takahashi
+!
+ XFQIAGGSP = XIKP * XCS**(1. + XINP) * &
+ MOMG(XALPHAS, XNUS, 2.+XDS*(1.+XINP)) * &
+ MOMG(XALPHAI, XNUI, XIMP)
+ XFQIAGGSN = XIKN * XCS**(1. + XINN) * &
+ MOMG(XALPHAS, XNUS, 2.+XDS*(1.+XINN)) * &
+ MOMG(XALPHAI, XNUI, XIMN)
+!
+ XFQIDRYGBSP = XIKP * XCG**(1. + XINP) * &
+ MOMG(XALPHAG, XNUG, 2.+XDG*(1.+XINP)) * &
+ MOMG(XALPHAI, XNUI, XIMP)
+ XFQIDRYGBSN = XIKN * XCG**(1. + XINN) * &
+ MOMG(XALPHAG, XNUG, 2.+XDG*(1.+XINN)) * &
+ MOMG(XALPHAI, XNUI, XIMN)
+!
+ XFQIAGGSP_TAK = XFQIAGGSP * XIKP_TAK / XIKP
+ XFQIAGGSN_TAK = XFQIAGGSN * XIKN_TAK / XIKN
+ XFQIDRYGBSP_TAK = XFQIDRYGBSP * XIKP_TAK / XIKP
+ XFQIDRYGBSN_TAK = XFQIDRYGBSN * XIKN_TAK / XIKN
+!
+ XAIGAMMABI = XAI * MOMG(XALPHAI, XNUI, XBI)
+!
+ XLBQSDRYGB1SP = MOMG(XALPHAG,XNUG,2.) * MOMG(XALPHAS, XNUS, XSMP)
+ XLBQSDRYGB1SN = MOMG(XALPHAG,XNUG,2.) * MOMG(XALPHAS, XNUS, XSMN)
+ XLBQSDRYGB2SP = 2. * MOMG(XALPHAG,XNUG,1.) * MOMG(XALPHAS, XNUS, 1.+XSMP)
+ XLBQSDRYGB2SN = 2. * MOMG(XALPHAG,XNUG,1.) * MOMG(XALPHAS, XNUS, 1.+XSMN)
+ XLBQSDRYGB3SP = MOMG(XALPHAS, XNUS, 2.+XSMP)
+ XLBQSDRYGB3SN = MOMG(XALPHAS, XNUS, 2.+XSMN)
ENDIF
!
IF (CNI_CHARGING == 'SAP98' .OR. CNI_CHARGING == 'TERAR' .OR. &
- CNI_CHARGING == 'BSMP1' .OR. CNI_CHARGING == 'BSMP2') &
- XINISAP = PRHO00**XCEXVT
+ CNI_CHARGING == 'BSMP1' .OR. CNI_CHARGING == 'BSMP2') THEN
+ XVSCOEF = XCS * MOMG(XALPHAS, XNUS, XBS+XDS) / MOMG(XALPHAS, XNUS, XBS)
+ XVGCOEF = XCG * MOMG(XALPHAG, XNUG, XBG+XDG) / MOMG(XALPHAG, XNUG, XBG)
+END IF
!
!
!* 8.3 Takahashi (1978) parameterization
@@ -526,6 +566,247 @@ IF (CNI_CHARGING == 'TAKAH') THEN
XMANSELL(:,:) = XMANSELL(:,:) * 1.E-15 ! in C
END IF
!
+!
+!* 8.4 Saunders et al. (1991) parameterization
+! Idem for Brooks et al. (1997), but with EW = ZRAR/3.
+!
+!
+IF (CNI_CHARGING == 'SAUN1' .OR. CNI_CHARGING == 'SAUN2' .OR. &
+ CNI_CHARGING == 'BSMP1' .OR. CNI_CHARGING == 'BSMP2') THEN
+!
+ NIND_TEMP = 31
+ NIND_LWC = 28
+!
+ IF( .NOT.ALLOCATED(XSAUNDER)) ALLOCATE(XSAUNDER(NIND_LWC+1,NIND_TEMP+1))
+ ALLOCATE(ZT(NIND_TEMP+1)) ! Kelvin
+ ALLOCATE(ZLWCC(NIND_TEMP+1))
+ DO JTEMP = 1, NIND_TEMP+1
+ ZT(JTEMP)=1.0-FLOAT(JTEMP)+XTT
+ END DO
+ ZLWCC(:) = MIN( MAX( -0.49 + 6.64E-2*(XTT-ZT(:)),0.22 ),1.1 ) ! (g m^-3)
+ ALLOCATE(ZEW(NIND_LWC+1))
+!
+! LWC index (0.01 g.m^-3 --> 10 g.m^-3)
+! 0.01 to 0.09 every 0.01 (9 values)
+! 0.10 to 0.90 every 0.10 (9 values)
+! 1.00 to 10.0 every 1.00 (10 values)
+ DO JLWC = 1, 9
+ ZEW(JLWC)=0.01*FLOAT(JLWC)
+ END DO
+ DO JLWC = 10, 18
+ ZEW(JLWC)=0.1 + 0.1*FLOAT(JLWC-10)
+ END DO
+ DO JLWC = 19, NIND_LWC+1
+ ZEW(JLWC)=1.0 + FLOAT(JLWC-19)
+ END DO
+!
+!
+ XSAUNDER(:,:) = 0.0
+ DO JTEMP = 1, NIND_TEMP+1
+ DO JLWC = 1, NIND_LWC+1
+!
+! region S4 : positive
+ IF (ZT(JTEMP) <= (XTT-7.35) .AND. ZT(JTEMP) > (XTT-23.9458) .AND. &
+ ZEW(JLWC) > ZLWCC(JTEMP)) THEN
+ XSAUNDER(JLWC,JTEMP) = MAX( 0., &
+ 20.22*ZEW(JLWC)+1.36*(ZT(JTEMP)-XTT)+10.05 )
+ ENDIF
+!
+! region S1 : positive --> linear interpolation
+ IF (ZT(JTEMP) > (XTT-7.35) .AND. ZT(JTEMP) < XTT .AND. &
+ ZEW(JLWC) > ZLWCC(JTEMP)) THEN
+ XSAUNDER(JLWC,JTEMP) = MAX( 0.,-(2.75*ZEW(JLWC)+0.007)*(ZT(JTEMP)-XTT) )
+ ENDIF
+!
+! region S8 : positive
+ IF (ZT(JTEMP) <= (XTT-23.9458) .AND. ZT(JTEMP) > (XTT-40.0) .AND. &
+ ZEW(JLWC) > ZLWCC(JTEMP)) THEN
+ XSAUNDER(JLWC,JTEMP) = MAX( 0.,20.22*ZEW(JLWC)-22.26 )
+ ENDIF
+!
+! region S7 : negative
+ IF (ZT(JTEMP) <= (XTT-7.35) .AND. ZT(JTEMP) > (XTT-40.0) .AND. &
+ ZEW(JLWC) >= 0.104149 .AND. ZEW(JLWC) < ZLWCC(JTEMP)) THEN
+ XSAUNDER(JLWC,JTEMP) = MIN( 0.,3.02-31.76*ZEW(JLWC)+26.53*ZEW(JLWC)**2 )
+ ENDIF
+ END DO
+ END DO
+END IF
+!
+! SAUN1 doesn't take into account marginal positive and negative regions at
+! low LWC
+!
+IF (CNI_CHARGING == 'SAUN1' .OR. CNI_CHARGING == 'BSMP1') THEN
+ DO JTEMP = 1, NIND_TEMP+1
+ DO JLWC = 1, NIND_LWC+1
+!
+! region S1 : negative --> linear interpolation
+ IF (ZT(JTEMP) > (XTT-7.35) .AND. ZT(JTEMP) < XTT .AND. &
+ ZEW(JLWC) < ZLWCC(JTEMP) .AND. ZEW(JLWC) >= 0.104149) THEN
+ XSAUNDER(JLWC,JTEMP) = MIN( 0., &
+ (-0.41+4.32*ZEW(JLWC)-3.61*ZEW(JLWC)**2)*(ZT(JTEMP)-XTT) )
+ ENDIF
+ END DO
+ END DO
+!
+ XSAUNDER(:,:) = XSAUNDER(:,:) * 1.E-15 ! in C
+!
+END IF
+!
+! SAUN2 takes into account marginal positive and negative regions at low LWC
+!
+IF (CNI_CHARGING == 'SAUN2' .OR. CNI_CHARGING == 'BSMP2') THEN
+ DO JTEMP = 1, NIND_TEMP+1
+ DO JLWC = 1, NIND_LWC+1
+!
+! region S2 : negative
+ IF (ZT(JTEMP) <= (XTT-7.35) .AND. ZT(JTEMP) > (XTT-16.0) .AND. &
+ ZEW(JLWC) >= 0.026 .AND. ZEW(JLWC) < 0.14) THEN
+ XSAUNDER(JLWC,JTEMP) = MIN( 0.,-314.4*ZEW(JLWC) + 7.92 )
+ ENDIF
+!
+! region S3 : negative
+ IF (ZT(JTEMP) <= (XTT-7.35) .AND. ZT(JTEMP) > (XTT-16.0) .AND. &
+ ZEW(JLWC) >= 0.14 .AND. ZEW(JLWC) < 0.22) THEN
+ XSAUNDER(JLWC,JTEMP) = MIN( 0.,419.4 * ZEW(JLWC) - 92.64 )
+ ENDIF
+!
+! region S5 : positive
+ IF (ZT(JTEMP) < (XTT-20.0) .AND. ZT(JTEMP) > (XTT-40.0) .AND. &
+ ZEW(JLWC) >= 0.063034 .AND. ZEW(JLWC) < 0.12) THEN
+ XSAUNDER(JLWC,JTEMP) = MAX( 0.,2041.76*ZEW(JLWC) - 128.7 )
+ ENDIF
+!
+! region S6 : positive
+ IF (ZT(JTEMP) < (XTT-20.0) .AND. ZT(JTEMP) > (XTT-40.0) .AND. &
+ ZEW(JLWC) >= 0.12 .AND. ZEW(JLWC) < 0.1596) THEN
+ XSAUNDER(JLWC,JTEMP) = MAX( 0.,-2900.22*ZEW(JLWC) + 462.91 )
+ ENDIF
+!
+! region S1 : negative --> linear interpolation of S3
+ IF (ZT(JTEMP) > (XTT-7.35) .AND. ZT(JTEMP) < XTT .AND. &
+ ZEW(JLWC) >= 0.14 .AND. ZEW(JLWC) < ZLWCC(JTEMP)) THEN
+ XSAUNDER(JLWC,JTEMP) = MIN( 0.,(-57.06*ZEW(JLWC)+12.6)*(ZT(JTEMP)-XTT) )
+ ENDIF
+!
+! region S1 : negative --> linear interpolation of S2
+ IF (ZT(JTEMP) > (XTT-7.35) .AND. ZT(JTEMP) < XTT .AND. &
+ ZEW(JLWC) >= 0.026 .AND. ZEW(JLWC) < 0.14) THEN
+ XSAUNDER(JLWC,JTEMP) = MIN( 0.,(42.8*ZEW(JLWC)-1.08)*(ZT(JTEMP)-XTT) )
+ ENDIF
+ END DO
+ END DO
+!
+ XSAUNDER(:,:) = XSAUNDER(:,:) * 1.E-15 ! in C
+!
+END IF
+!
+!* 8.5 Takahashi with EW or ZRAR (Tsenova and Mitzeva, 2009, 2011)
+! here ZRAR = 9 * EW
+! Temperature index (0C --> -30C)
+! LWC index (0.01 g.m^-3 --> 10 g.m^-3)
+! 0.01 to 0.09 every 0.01 (9 values)
+! 0.10 to 0.90 every 0.10 (9 values)
+! 1.00 to 10.0 every 1.00 (10 values)
+!
+IF (CNI_CHARGING == 'TEEWC' .OR. CNI_CHARGING == 'TERAR') THEN
+!
+ NIND_TEMP = 31
+ NIND_LWC = 28
+!
+ IF( .NOT.ALLOCATED(XTAKA_TM)) ALLOCATE(XTAKA_TM(NIND_LWC+1,NIND_TEMP+1))
+ ALLOCATE(ZT(NIND_TEMP+1)) ! Kelvin
+ ALLOCATE(ZEW(NIND_LWC+1))
+ DO JTEMP = 1, NIND_TEMP+1
+ ZT(JTEMP) = 1.0 - FLOAT(JTEMP) + XTT
+ END DO
+
+ DO JLWC = 1, 9
+ ZEW(JLWC) = 0.01 * FLOAT(JLWC)
+ END DO
+ DO JLWC = 10, 18
+ ZEW(JLWC) = 0.1 + 0.1 * FLOAT(JLWC-10)
+ END DO
+ DO JLWC = 19, NIND_LWC+1
+ ZEW(JLWC) = 1.0 + FLOAT(JLWC-19)
+ END DO
+!
+ XTAKA_TM(:,:) = 0.0
+ DO JTEMP = 1, NIND_TEMP+1
+ DO JLWC = 1, NIND_LWC+1
+!
+! Eq. 1: >0
+ IF ( ZT(JTEMP) > (XTT - 10.) .AND. ZEW(JLWC) <= 1.6) THEN
+ XTAKA_TM(JLWC, JTEMP) = 146.981 * ZEW(JLWC) - 116.37 * ZEW(JLWC)**2 &
+ + 29.76 * ZEW(JLWC)**3 &
+ - 0.03 * (ZT(JTEMP) - XTT)**3 * ZEW(JLWC) &
+ - 2.58 * (ZT(JTEMP) - XTT) &
+ - 0.21 * (ZT(JTEMP) - XTT)**3 * ZEW(JLWC)**3 &
+ + 0.36 * (ZT(JTEMP) - XTT)**3 * ZEW(JLWC)**2 &
+ + 0.15 * (ZT(JTEMP) - XTT)**2 &
+ + 2.92 * (ZT(JTEMP) - XTT) * ZEW(JLWC)**3 &
+ - 4.22 * (ZT(JTEMP) - XTT) * ZEW(JLWC) - 8.506
+ END IF
+!
+! Eq. 2: >0
+ IF ( ZT(JTEMP) > (XTT - 10.) .AND. &
+ ZEW(JLWC) > 1.6 .AND. ZEW(JLWC) <= 8.) THEN
+ XTAKA_TM(JLWC, JTEMP) = 4.179 * (ZT(JTEMP) - XTT) &
+ - 0.005 * (ZT(JTEMP) - XTT)**2 * ZEW(JLWC)**2 &
+ + 0.916 * ZEW(JLWC)**2 &
+ - 1.333 * (ZT(JTEMP) - XTT) * ZEW(JLWC) &
+ - 7.465 * ZEW(JLWC) &
+ + 0.109 * (ZT(JTEMP) - XTT) * ZEW(JLWC)**2 &
+ + 0.001 * (ZT(JTEMP) - XTT)**2 * ZEW(JLWC)**3 &
+ - 0.035 * ZEW(JLWC)**3 + 50.84454
+ END IF
+!
+! Eq. 8: > 0
+ IF ( ZEW(JLWC) <= 0.4 .AND. &
+ ZT(JTEMP) <= (XTT - 10.) .AND. ZT(JTEMP) >= (XTT - 40.)) THEN
+ XTAKA_TM(JLWC, JTEMP) = - 3.3515 * (ZT(JTEMP) - XTT) &
+ + 95.957 * (ZT(JTEMP) - XTT) * ZEW(JLWC)**2 &
+ + 511.83 * ZEW(JLWC) &
+ + 17.448 * (ZT(JTEMP) - XTT)**2 * ZEW(JLWC)**3 &
+ - 0.0007 * (ZT(JTEMP) - XTT)**3 &
+ + 20.570 * (ZT(JTEMP) - XTT) * ZEW(JLWC) &
+ + 0.1656 * (ZT(JTEMP) - XTT)**2 * ZEW(JLWC) &
+ + 0.4954 * (ZT(JTEMP) - XTT)**3 * ZEW(JLWC)**3 &
+ - 0.0975 * (ZT(JTEMP) - XTT)**3 * ZEW(JLWC)**2 &
+ + 67.457 * (ZT(JTEMP) - XTT) * ZEW(JLWC)**3 &
+ - 0.1066 * (ZT(JTEMP) - XTT)**2 - 24.5715
+ END IF
+!
+! Eq. 9: < 0
+ IF ( ZT(JTEMP) <= (XTT - 10.) .AND. ZT(JTEMP) >= (XTT - 40.) .AND. &
+ ZEW(JLWC) > 0.4 .AND. ZEW(JLWC) <= 3.2) THEN
+ XTAKA_TM(JLWC, JTEMP) = - 1.5676 * (ZT(JTEMP) - XTT) * ZEW(JLWC) &
+ + 0.2484 * (ZT(JTEMP) - XTT) * ZEW(JLWC)**3 &
+ + 0.0112 * (ZT(JTEMP) - XTT)**3 &
+ + 19.199 * (ZT(JTEMP) - XTT) &
+ + 0.8051 * (ZT(JTEMP) - XTT)**2 &
+ - 83.4 * ZEW(JLWC) &
+ + 15.4 * ZEW(JLWC)**2 &
+ + 5.97 * ZEW(JLWC)**3 + 167.9278
+ END IF
+!
+! Eq. 10: > 0
+ IF ( ZT(JTEMP) <= (XTT - 10.) .AND. ZT(JTEMP) >= (XTT - 40.) .AND. &
+ ZEW(JLWC) > 3.2 .AND. ZEW(JLWC) <= 8. ) THEN
+ XTAKA_TM(JLWC, JTEMP) = 4.2127 * (ZT(JTEMP) - XTT) &
+ - 0.8311 * (ZT(JTEMP) - XTT) * ZEW(JLWC) &
+ + 0.0670 * (ZT(JTEMP) - XTT) * ZEW(JLWC) **2 &
+ + 0.0042 * (ZT(JTEMP) - XTT)**2 * ZEW(JLWC) &
+ + 40.9642
+ END IF
+ END DO
+ END DO
+!
+ XTAKA_TM(:,:) = XTAKA_TM(:,:) * 1.E-15 ! in C
+!
+END IF
+!
+!
!-------------------------------------------------------------------------------
!
!* 9. NON_INDUCTIVE PROCESS: AGGREGATION OF ICE ON SNOW
@@ -551,10 +832,10 @@ XFQIAGGSBS = (XPI / 4.0) * XCCS
!* 9.4 Takahashi (1978) parameterization
!
IF (CNI_CHARGING == 'TAKAH') THEN
- XFQIAGGSBT1 = (XPI / 4.0) * XCCS * (PRHO00**XCEXVT) * XCS
+ XFQIAGGSBT1 = (XPI / 4.0) * XCCS * XCS
XFQIAGGSBT2 = 10 * MOMG(XALPHAS,XNUS,2.+XDS)
- XFQIAGGSBT3 = 5. * XCS * (PRHO00**XCEXVT) * MOMG(XALPHAI,XNUI,2.) * &
- MOMG(XALPHAS,XNUS,2.+2*XDS) / ((1.E-4)**2 * 8. * &
+ XFQIAGGSBT3 = 5. * XCS * MOMG(XALPHAI,XNUI,2.) * &
+ MOMG(XALPHAS,XNUS,2.+2*XDS) / ((1.E-4)**2 * 8. * &
(XAI * MOMG(XALPHAI,XNUI,XBI))**(2 / XBI))
END IF
!
@@ -694,21 +975,21 @@ ENDIF
IF (CNI_CHARGING == 'TAKAH') THEN
!
! IDRYG_boun
- XFQIDRYGBT1 = (XPI / 4.0) * XCCG * (PRHO00**XCEXVT) * XCG
+ XFQIDRYGBT1 = (XPI / 4.0) * XCCG * XCG
XFQIDRYGBT2 = 10.0 * MOMG(XALPHAG,XNUG,2.+XDG)
- XFQIDRYGBT3 = 5.0 * XCG * (PRHO00**XCEXVT) * MOMG(XALPHAI,XNUI,2.) * &
- MOMG(XALPHAG,XNUG,2.+2.*XDG) / ((2.E-4)**2 * 8. * &
+ XFQIDRYGBT3 = 5.0 * XCG * MOMG(XALPHAI,XNUI,2.) * &
+ MOMG(XALPHAG,XNUG,2.+2.*XDG) / ((2.E-4)**2 * 8. * &
(XAI * MOMG(XALPHAI,XNUI,XBI))**(2 / XBI))
!
! SDRYG_boun
- XFQSDRYGBT1 = (XPI / 4.0) * XCCG * (PRHO00**XCEXVT) * XCCS
+ XFQSDRYGBT1 = (XPI / 4.0) * XCCG * XCCS
XFQSDRYGBT2 = XCG * MOMG(XALPHAG,XNUG,XDG) * MOMG(XALPHAS,XNUS,2.)
XFQSDRYGBT3 = XCS * MOMG(XALPHAS,XNUS,2.+XDS)
XFQSDRYGBT4 = XCG * MOMG(XALPHAG,XNUG,2.+XDG)
XFQSDRYGBT5 = XCS * MOMG(XALPHAG,XNUG,2.) * MOMG(XALPHAS,XNUS,XDS)
XFQSDRYGBT6 = 2. * XCG * MOMG(XALPHAG,XNUG,1.+XDG) * MOMG(XALPHAS,XNUS,1.)
XFQSDRYGBT7 = 2. * XCS * MOMG(XALPHAG,XNUG,1.) * MOMG(XALPHAS,XNUS,1.+XDS)
- XFQSDRYGBT8 = 5. * (PRHO00**XCEXVT) / ((1.E-4)**2 * 8.)
+ XFQSDRYGBT8 = 5. / ((1.E-4)**2 * 8.)
XFQSDRYGBT9 = MOMG(XALPHAG,XNUG,2.) * MOMG(XALPHAS,XNUS,2.)
XFQSDRYGBT10 = MOMG(XALPHAS,XNUS,4.)
XFQSDRYGBT11 = 2. * MOMG(XALPHAG,XNUG,1.) * MOMG(XALPHAS,XNUS,3.)
@@ -728,7 +1009,7 @@ IF (CNI_CHARGING == 'TAKAH' .OR. CNI_CHARGING == 'SAP98' .OR. &
CNI_CHARGING == 'GARDI' .OR. &
CNI_CHARGING == 'BSMP1' .OR. CNI_CHARGING == 'BSMP2' .OR. &
CNI_CHARGING == 'TEEWC' .OR. CNI_CHARGING == 'TERAR') THEN
- XAUX_LIM = (XPI / 4.0) * XCCG * XCCS * (PRHO00**(XCEXVT))
+ XAUX_LIM = (XPI / 4.0) * XCCG * XCCS
XAUX_LIM1 = MOMG(XALPHAS,XNUS,2.)
XAUX_LIM2 = 2. * MOMG(XALPHAS,XNUS,1.) * MOMG(XALPHAG,XNUG,1.)
XAUX_LIM3 = MOMG(XALPHAG,XNUG,2.)
@@ -790,7 +1071,7 @@ XALPHA_IND = 0.07 ! moderate inductive charging
XCOS_THETA = 0.2
!
XIND1 = (XPI**3 / 8.) * (15.E-6)**2 * &
- PRHO00**(XCEXVT) * XCG * 400.E6 * XCCG * &
+ XCG * 400.E6 * XCCG * &
XCOLCG_IND * XEBOUND * XALPHA_IND
XIND2 = XPI * XEPSILON * XCOS_THETA * MOMG(XALPHAG,XNUG,2.+XDG)
XIND3 = MOMG(XALPHAG,XNUG,XDG+XFG) / 3.
@@ -815,14 +1096,13 @@ XEXQLIGHTS = XCXS - 2.
XFQLIGHTG = XPI * XCCG * MOMG(XALPHAG,XNUG,2.)
XEXQLIGHTG = XCXG - 2.
!
-IF (IP == 1) THEN
- IF( .NOT.ALLOCATED(XNEUT_POS)) ALLOCATE( XNEUT_POS(NLGHTMAX) )
- IF( .NOT.ALLOCATED(XNEUT_NEG)) ALLOCATE( XNEUT_NEG(NLGHTMAX) )
- XNEUT_POS(:) = 0.
- XNEUT_NEG(:) = 0.
-ENDIF
+XFQLIGHTH = XPI * XCCH * MOMG(XALPHAH,XNUH,2.)
+XEXQLIGHTH = XCXH - 2.
!
-XALT_CG = 2000. ! m
+IF( .NOT.ALLOCATED(XNEUT_POS)) ALLOCATE( XNEUT_POS(NLGHTMAX) )
+IF( .NOT.ALLOCATED(XNEUT_NEG)) ALLOCATE( XNEUT_NEG(NLGHTMAX) )
+XNEUT_POS(:) = 0.
+XNEUT_NEG(:) = 0.
!
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/ini_segn.f90 b/src/MNH/ini_segn.f90
index 69418af98..05ab6e943 100644
--- a/src/MNH/ini_segn.f90
+++ b/src/MNH/ini_segn.f90
@@ -6,6 +6,7 @@
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$
+! masdev4_7 BUG1 2007/06/15 17:47:27
!-----------------------------------------------------------------
! ###################
MODULE MODI_INI_SEG_n
@@ -162,6 +163,7 @@ END MODULE MODI_INI_SEG_n
!! test (Escobar)
!! 10/02/15 remove ABORT in parallel case for SPAWNING
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! 01/2015 add GLNOX_EXPLICIT (C. Barthe)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -236,6 +238,7 @@ LOGICAL :: GFOREFIRE
#endif
LOGICAL :: GCONDSAMP
LOGICAL :: GCHTRANS
+LOGICAL :: GLNOX_EXPLICIT ! flag for LNOx
! These variables
! are used to locally store
INTEGER :: ISV ! the value read in DESFM
@@ -443,6 +446,7 @@ CALL READ_DESFM_n(KMI,YDESFM,HLUOUT,YCONF,GFLAT,GUSERV,GUSERC, &
#ifdef MNH_FOREFIRE
GFOREFIRE, &
#endif
+ GLNOX_EXPLICIT, &
GCONDSAMP, IRIMX,IRIMY,ISV, &
YTURB,YTOM,GRMC01,YRAD,YDCONV,YSCONV,YCLOUD,YELEC,YEQNSYS )
!
@@ -463,6 +467,7 @@ CALL READ_EXSEG_n(KMI,YEXSEG,HLUOUT,YCONF,GFLAT,GUSERV,GUSERC, &
#ifdef MNH_FOREFIRE
GFOREFIRE, &
#endif
+ GLNOX_EXPLICIT, &
GCONDSAMP, IRIMX,IRIMY,ISV, &
YTURB,YTOM,GRMC01,YRAD,YDCONV,YSCONV,YCLOUD,YELEC,YEQNSYS, &
PTSTEP_ALL,CSTORAGE_TYPE,CINIFILEPGD_n )
diff --git a/src/MNH/ion_attach_elec.f90 b/src/MNH/ion_attach_elec.f90
index 670113ae4..33963eb14 100644
--- a/src/MNH/ion_attach_elec.f90
+++ b/src/MNH/ion_attach_elec.f90
@@ -74,7 +74,6 @@ END MODULE MODI_ION_ATTACH_ELEC
!! -------------
!! Original 2010
!! Modifications:
-!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!!
!-------------------------------------------------------------------------------
!
@@ -92,7 +91,6 @@ USE MODD_RAIN_ICE_PARAM
USE MODD_NSV, ONLY : NSV_ELECBEG, NSV_ELEC
USE MODD_BUDGET, ONLY : LBU_RSV
USE MODD_REF, ONLY : XTHVREFZ
-USE MODE_ll
!
USE MODI_BUDGET
USE MODI_MOMG
@@ -151,7 +149,10 @@ REAL :: ZCOMB ! Recombination
ZCQD = 4 * XPI * XEPSILON * XBOLTZ / XECHARGE
ZCDIF = XBOLTZ /XECHARGE
!
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IIB = 1 + JPHEXT
+IIE = SIZE(PTHT,1) - JPHEXT
+IJB = 1 + JPHEXT
+IJE = SIZE(PTHT,2) - JPHEXT
IKB = 1 + JPVEXT
IKE = SIZE(PTHT,3) - JPVEXT
!
@@ -160,24 +161,24 @@ IKE = SIZE(PTHT,3) - JPVEXT
!
IVALID = 0
DO IK = IKB, IKE
- DO IJ = IJB, IJE
- DO II = IIB, IIE
- IF (GATTACH(II,IJ,IK)) THEN
+ DO IJ = IJB, IJE
+ DO II = IIB, IIE
+ IF (GATTACH(II,IJ,IK)) THEN
! Recombination
- ZCOMB = XIONCOMB * (PSVS(II,IJ,IK,1)*PTSTEP) * &
- (PSVS(II,IJ,IK,NSV_ELEC)*PTSTEP) * &
- PRHODREF(II,IJ,IK) / PRHODJ(II,IJ,IK)
- ZCOMB = MIN(ZCOMB, PSVS(II,IJ,IK,1), PSVS(II,IJ,IK,NSV_ELEC))
- PSVS(II,IJ,IK,1) = PSVS(II,IJ,IK,1) - ZCOMB
- PSVS(II,IJ,IK,NSV_ELEC) = PSVS(II,IJ,IK,NSV_ELEC) - ZCOMB
+ ZCOMB = XIONCOMB * (PSVS(II,IJ,IK,1)*PTSTEP) * &
+ (PSVS(II,IJ,IK,NSV_ELEC)*PTSTEP) * &
+ PRHODREF(II,IJ,IK) / PRHODJ(II,IJ,IK)
+ ZCOMB = MIN(ZCOMB, PSVS(II,IJ,IK,1), PSVS(II,IJ,IK,NSV_ELEC))
+ PSVS(II,IJ,IK,1) = PSVS(II,IJ,IK,1) - ZCOMB
+ PSVS(II,IJ,IK,NSV_ELEC) = PSVS(II,IJ,IK,NSV_ELEC) - ZCOMB
! Counting
- IVALID = IVALID + 1
- IGI(IVALID) = II
- IGJ(IVALID) = IJ
- IGK(IVALID) = IK
- END IF
- ENDDO
- ENDDO
+ IVALID = IVALID + 1
+ IGI(IVALID) = II
+ IGJ(IVALID) = IJ
+ IGK(IVALID) = IK
+ END IF
+ ENDDO
+ ENDDO
ENDDO
!
!* 1.2 Compute the temperature
@@ -185,8 +186,8 @@ ENDDO
IF( IVALID /= 0 ) THEN
ALLOCATE (ZT(IVALID))
DO II = 1, IVALID
- ZT(II) = PTHT(IGI(II),IGJ(II),IGK(II)) * &
- (PPABST(IGI(II),IGJ(II),IGK(II)) / XP00) ** (XRD / XCPD)
+ ZT(II) = PTHT(IGI(II),IGJ(II),IGK(II)) * &
+ (PPABST(IGI(II),IGJ(II),IGK(II)) / XP00) ** (XRD / XCPD)
ENDDO
END IF
!
@@ -221,9 +222,9 @@ IF( IVALID /= 0 ) THEN
ITYPE = 2
IF (PRESENT(PSEA)) THEN
- CALL HYDROPARAM (IGI, IGJ, IGK, ZCONC, ZVIT, ZRADIUS, ITYPE, PSEA, PTOWN)
+ CALL HYDROPARAM (IGI, IGJ, IGK, ZCONC, ZVIT, ZRADIUS, ITYPE, PSEA, PTOWN)
ELSE
- CALL HYDROPARAM (IGI, IGJ, IGK, ZCONC, ZVIT, ZRADIUS, ITYPE)
+ CALL HYDROPARAM (IGI, IGJ, IGK, ZCONC, ZVIT, ZRADIUS, ITYPE)
ENDIF
!
CALL DIFF_COND (IGI, IGJ, IGK, PSVS(:,:,:,1), PSVS(:,:,:,NSV_ELEC), &
@@ -233,10 +234,10 @@ IF( IVALID /= 0 ) THEN
! and hail (if activated)
!
DO ITYPE = 3, KRR
- CALL HYDROPARAM (IGI, IGJ, IGK, ZCONC, ZVIT, ZRADIUS, ITYPE)
+ CALL HYDROPARAM (IGI, IGJ, IGK, ZCONC, ZVIT, ZRADIUS, ITYPE)
!
- CALL DIFF_COND (IGI, IGJ, IGK, PSVS(:,:,:,1), PSVS(:,:,:,NSV_ELEC), &
- PSVS(:,:,:,ITYPE))
+ CALL DIFF_COND (IGI, IGJ, IGK, PSVS(:,:,:,1), PSVS(:,:,:,NSV_ELEC), &
+ PSVS(:,:,:,ITYPE))
END DO
!
DEALLOCATE (ZCONC, ZVIT, ZRADIUS)
@@ -363,7 +364,7 @@ SELECT CASE (ITYPE)
!* 2. PARAMETERS FOR RAIN
! -------------------
CASE (3)
- ZEXP1 = XEXSEDR
+ ZEXP1 = XEXSEDR - 1.
ZEXP2 = ZEXP1 - XCEXVT
!
DO IV = 1, IVALID
@@ -376,7 +377,7 @@ SELECT CASE (ITYPE)
ZRADIUS (IV) = 0.5 / ZLAMBDA
ZCONC (IV) = XCCR / ZLAMBDA
ZVIT (IV) = XFSEDR * PRHODREF(JI,JJ,JK)**ZEXP2 &
- * PRS(JI,JJ,JK,3)**(ZEXP1-1.)
+ * PRS(JI,JJ,JK,3)**ZEXP1
END IF
ENDDO
!
@@ -418,7 +419,7 @@ SELECT CASE (ITYPE)
!
CASE (5)
!
- ZEXP1 = XEXSEDS
+ ZEXP1 = XEXSEDS - 1.
ZEXP2 = ZEXP1 - XCEXVT
!
DO IV = 1, IVALID
@@ -431,7 +432,7 @@ SELECT CASE (ITYPE)
ZRADIUS (IV) = 0.5 / ZLAMBDA
ZCONC (IV) = XCCS * ZLAMBDA**XCXS
ZVIT (IV) = XFSEDS * PRHODREF(JI,JJ,JK)**ZEXP2 &
- * PRS(JI,JJ,JK,5)**(ZEXP1-1.)
+ * PRS(JI,JJ,JK,5)**ZEXP1
END IF
ENDDO
!
@@ -441,7 +442,7 @@ SELECT CASE (ITYPE)
!
CASE (6)
!
- ZEXP1 = XEXSEDG
+ ZEXP1 = XEXSEDG - 1.
ZEXP2 = ZEXP1 - XCEXVT
!
DO IV = 1, IVALID
@@ -454,7 +455,7 @@ SELECT CASE (ITYPE)
ZRADIUS (IV) = 0.5 / ZLAMBDA
ZCONC (IV) = XCCG * ZLAMBDA**XCXG
ZVIT (IV) = XFSEDG * PRHODREF(JI,JJ,JK)**ZEXP2 &
- * PRS(JI,JJ,JK,6)**(ZEXP1-1.)
+ * PRS(JI,JJ,JK,6)**ZEXP1
END IF
ENDDO
!
@@ -464,7 +465,7 @@ SELECT CASE (ITYPE)
!
CASE (7)
!
- ZEXP1 = XEXSEDH
+ ZEXP1 = XEXSEDH - 1.
ZEXP2 = ZEXP1-XCEXVT
ZRAY = 0.5*MOMG(XALPHAH, XNUH, 1.)
!
@@ -478,7 +479,7 @@ SELECT CASE (ITYPE)
ZRADIUS (IV) = ZRAY / ZLAMBDA
ZCONC (IV) = XCCG * ZLAMBDA**XCXG
ZVIT (IV) = XFSEDH * PRHODREF(JI,JJ,JK)**ZEXP2 &
- * PRS(JI,JJ,JK,7)**(ZEXP1-1.)
+ * PRS(JI,JJ,JK,7)**ZEXP1
END IF
ENDDO
!
diff --git a/src/MNH/ion_drift.f90 b/src/MNH/ion_drift.f90
index 438dcb8bf..14ce9fd06 100644
--- a/src/MNH/ion_drift.f90
+++ b/src/MNH/ion_drift.f90
@@ -1,3 +1,4 @@
+
!MNH_LIC Copyright 1994-2014 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
@@ -8,25 +9,19 @@
INTERFACE
!
- SUBROUTINE ION_DRIFT(PDRIFTP, PDRIFTM, PSVT, PRHODREF, PRHODJ, &
- HLBCX, HLBCY, KTCOUNT, PTSTEP, HDRIFT)
+ SUBROUTINE ION_DRIFT(PDRIFTP, PDRIFTM, PSVT, HLBCX, HLBCY)
!
-CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY
-CHARACTER(LEN=3), INTENT(IN) :: HDRIFT
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDRIFTP, PDRIFTM
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF, PRHODJ
-INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop counter
-REAL, INTENT(IN) :: PTSTEP
+CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDRIFTP, PDRIFTM
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT
!
END SUBROUTINE ION_DRIFT
END INTERFACE
END MODULE MODI_ION_DRIFT
!
-! ################################################################
- SUBROUTINE ION_DRIFT(PDRIFTP, PDRIFTM, PSVT, PRHODREF, PRHODJ, &
- HLBCX, HLBCY, KTCOUNT, PTSTEP, HDRIFT)
-! ################################################################
+! ##########################################################
+ SUBROUTINE ION_DRIFT(PDRIFTP, PDRIFTM, PSVT, HLBCX, HLBCY)
+! ##########################################################
!
!! PURPOSE
!! -------
@@ -50,19 +45,16 @@ END MODULE MODI_ION_DRIFT
USE MODD_PARAMETERS
USE MODD_CONF
USE MODD_METRICS_n, ONLY : XDXX, XDYY, XDZX, XDZY, XDZZ
-USE MODD_NSV, ONLY: NSV_ELECBEG, NSV_ELECEND
+USE MODD_NSV, ONLY: NSV_ELECBEG, NSV_ELECEND, XSVMIN
USE MODD_ELEC_n, ONLY : XCION_POS_FW, XCION_NEG_FW, &
- XMOBIL_POS, XMOBIL_NEG, &
- XEFIELDU, XEFIELDV, XEFIELDW
+ XMOBIL_POS, XMOBIL_NEG, &
+ XEFIELDU, XEFIELDV, XEFIELDW
USE MODD_ARGSLIST_ll, ONLY : LIST_ll
!
USE MODE_ll
USE MODE_ELEC_ll
!
USE MODI_SHUMAN
-USE MODI_CONTRAV
-USE MODI_PPM_SCALAR
-USE MODI_PPM_RHODJ
USE MODI_GDIV
USE MODI_ION_BOUND4DRIFT
!
@@ -71,13 +63,9 @@ IMPLICIT NONE
!
!* 0.1 declarations of arguments
!
-CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY
-CHARACTER(LEN=3), INTENT(IN) :: HDRIFT
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDRIFTP, PDRIFTM
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF, PRHODJ
-INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop counter
-REAL, INTENT(IN) :: PTSTEP
+CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDRIFTP, PDRIFTM
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT
!
!
!* 0.2 declarations of local variables
@@ -86,24 +74,10 @@ INTEGER :: IIB, IIE ! index of first and last inner mass points along x
INTEGER :: IJB, IJE ! index of first and last inner mass points along y
INTEGER :: IKB, IKE ! index of first and last inner mass points along z
INTEGER :: IKU
-INTEGER :: IXOR, IYOR ! origin of the extended subdomain
-INTEGER, DIMENSION(3) :: IM_LOC
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZDRIFTX
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZDRIFTY
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZDRIFTZ
-
-REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZADVS, ZSVT !for advection form
-REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZXCT, ZYCT, ZZCT !
-CHARACTER (LEN=6) :: HSV_ADV_SCHEME
- ! of drift source
-!
-REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZRHOX1,ZRHOX2
-REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZRHOY1,ZRHOY2
-REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZRHOZ1,ZRHOZ2
-!
-REAL :: ZMIN_DRIFT, ZMAX_DRIFT
-REAL :: ZMAX__POS, ZMAX__NEG
-INTEGER :: IPROC, IPROCMIN, ISV
+!
INTEGER :: IINFO_ll ! return code of parallel routine
!
TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
@@ -113,8 +87,6 @@ NULLIFY(TZFIELDS_ll)
!
!------------------------------------------------------------------------
!
-CALL MYPROC_ELEC_ll (IPROC)
-!
!* 1. COMPUTE DIMENSIONS OF ARRAYS AND OTHER INDICES
! ----------------------------------------------
!
@@ -127,7 +99,7 @@ IKU = SIZE(PSVT,3)
!
!-------------------------------------------------------------------------------
!
-!* 3. UPDATE BOUNDARY CONDITION FOR IONS ACCORDING TO THE DRIFT MOTION
+!* 2. UPDATE BOUNDARY CONDITION FOR IONS ACCORDING TO THE DRIFT MOTION
! ----------------------------------------------------------------
!
IF (LWEST_ll() ) THEN
@@ -160,151 +132,76 @@ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDS_ll)
!
! specify upper boundary ion mixing ratio
-WHERE (XEFIELDW(:,:,IKE) .GE. 0.) ! Out(In)flow for positive (negative) ions
- PSVT (:,:,IKE+1,NSV_ELECBEG) = 2. * PSVT (:,:,IKE,NSV_ELECBEG) - &
- PSVT (:,:,IKE-1,NSV_ELECBEG)
+WHERE (XEFIELDW(:,:,IKE+1) .GE. 0.) ! Out(In)flow for positive (negative) ions
+ PSVT (:,:,IKE+1,NSV_ELECBEG) = MAX(2. * PSVT (:,:,IKE,NSV_ELECBEG) - &
+ PSVT (:,:,IKE-1,NSV_ELECBEG),XSVMIN(NSV_ELECBEG))
PSVT (:,:,IKE+1,NSV_ELECEND) = XCION_NEG_FW(:,:,IKE+1)
ELSE WHERE ! In(Out)flow for positive (negative) ions
PSVT (:,:,IKE+1,NSV_ELECBEG) = XCION_POS_FW(:,:,IKE+1)
- PSVT (:,:,IKE+1,NSV_ELECEND) = 2.* PSVT (:,:,IKE,NSV_ELECEND) - &
- PSVT (:,:,IKE-1,NSV_ELECEND)
+ PSVT (:,:,IKE+1,NSV_ELECEND) = MAX(2.* PSVT (:,:,IKE,NSV_ELECEND) - &
+ PSVT (:,:,IKE-1,NSV_ELECEND),XSVMIN(NSV_ELECEND))
END WHERE
!
XEFIELDW(:,:,IKB-1) = XEFIELDW(:,:,IKB)
XEFIELDW(:,:,IKE+1) = XEFIELDW(:,:,IKE)
!
+!-------------------------------------------------------------------------------
!
-!* 4. positive ion source
+!* 3. DRIFT MOTION TAKEN AS THE DIVERGENCE OF THE DRIFT FLUXES
+! --------------------------------------------------------
+!
+!* 3.1 positive ion source (drifting along E)
!
-IF (HDRIFT /= 'PPM') THEN ! Divergence form
! x-component of div term
- ZDRIFTX(:,:,:) = -PSVT(:,:,:,NSV_ELECBEG) * XMOBIL_POS(:,:,:) * XEFIELDU(:,:,:) &
- * PRHODJ(:,:,:)
+ZDRIFTX(:,:,:) = PSVT(:,:,:,NSV_ELECBEG) * XMOBIL_POS(:,:,:)
+ZDRIFTX(:,:,:) = ZDRIFTX(:,:,:) * XEFIELDU(:,:,:)
+ZDRIFTX(:,:,:) = -MXM(ZDRIFTX(:,:,:)) ! Put components at flux sides
+!
! y-component of div term
- ZDRIFTY(:,:,:) = -PSVT(:,:,:,NSV_ELECBEG) * XMOBIL_POS(:,:,:) * XEFIELDV(:,:,:) &
- * PRHODJ(:,:,:)
-! z-component of div term
- ZDRIFTZ(:,:,:) = -PSVT(:,:,:,NSV_ELECBEG) * XMOBIL_POS(:,:,:) * XEFIELDW(:,:,:) &
- * PRHODJ(:,:,:)
-ELSE ! Advection form
- ZDRIFTX(:,:,:) = XMOBIL_POS(:,:,:) * XEFIELDU(:,:,:)
- ZDRIFTY(:,:,:) = XMOBIL_POS(:,:,:) * XEFIELDV(:,:,:)
- ZDRIFTZ(:,:,:) = XMOBIL_POS(:,:,:) * XEFIELDW(:,:,:)
-ENDIF
+ZDRIFTY(:,:,:) = PSVT(:,:,:,NSV_ELECBEG) * XMOBIL_POS(:,:,:)
+ZDRIFTY(:,:,:) = ZDRIFTY(:,:,:) * XEFIELDV(:,:,:)
+ZDRIFTY(:,:,:) = -MYM(ZDRIFTY(:,:,:)) ! Put components at flux sides
!
-! Put components at flux sides
-ZDRIFTX(:,:,:) = MXM(ZDRIFTX(:,:,:))
-ZDRIFTY(:,:,:) = MYM(ZDRIFTY(:,:,:))
-ZDRIFTZ(:,:,:) = MZM(1,IKU,1,ZDRIFTZ (:,:,:))
+! z-component of div term
+ZDRIFTZ(:,:,:) = PSVT(:,:,:,NSV_ELECBEG) * XMOBIL_POS(:,:,:)
+ZDRIFTZ(:,:,:) = ZDRIFTZ(:,:,:) * XEFIELDW(:,:,:)
+ZDRIFTZ(:,:,:) = -MZM(1,IKU,1,ZDRIFTZ(:,:,:)) ! Put components at flux sides
!
IF (LWEST_ll( )) ZDRIFTX(IIB-1,:,:) = ZDRIFTX(IIB,:,:)
+IF (LEAST_ll( )) ZDRIFTX(IIE+1,:,:) = ZDRIFTX(IIE,:,:)
IF (LSOUTH_ll( )) ZDRIFTY(:,IJB-1,:) = ZDRIFTY(:,IJB,:)
+IF (LNORTH_ll( )) ZDRIFTY(:,IJE+1,:) = ZDRIFTY(:,IJE,:)
ZDRIFTZ(:,:,IKB-1) = ZDRIFTZ(:,:,IKB)
+ZDRIFTZ(:,:,IKE+1) = ZDRIFTZ(:,:,IKE)
!
-IF (HDRIFT /= 'PPM') THEN
- CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZDRIFTX,ZDRIFTY,ZDRIFTZ,PDRIFTP)
-ELSE
- ISV = 1
- ALLOCATE (ZADVS(SIZE(PDRIFTP,1),SIZE(PDRIFTP,2),SIZE(PDRIFTP,3), ISV))
- ALLOCATE (ZSVT(SIZE(PDRIFTP,1),SIZE(PDRIFTP,2),SIZE(PDRIFTP,3), ISV))
- ALLOCATE (ZXCT(SIZE(PDRIFTP,1),SIZE(PDRIFTP,2),SIZE(PDRIFTP,3)))
- ALLOCATE (ZYCT(SIZE(PDRIFTP,1),SIZE(PDRIFTP,2),SIZE(PDRIFTP,3)))
- ALLOCATE (ZZCT(SIZE(PDRIFTP,1),SIZE(PDRIFTP,2),SIZE(PDRIFTP,3)))
-!
- CALL CONTRAV (HLBCX,HLBCY,ZDRIFTX,ZDRIFTY,ZDRIFTZ,XDXX,XDYY,XDZZ, &
- XDZX,XDZY,ZXCT,ZYCT,ZZCT,4)
-!
- ZXCT = ZXCT*PTSTEP
- ZYCT = ZYCT*PTSTEP
- ZZCT = ZZCT*PTSTEP
-!
- ZADVS(:,:,:,1) = 0.
- ZSVT(:,:,:,1) = PSVT(:,:,:,NSV_ELECBEG)
- HSV_ADV_SCHEME = 'PPM_01'
-!
- CALL PPM_RHODJ(HLBCX,HLBCY, ZXCT, ZYCT, ZZCT, &
- PTSTEP, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, &
- ZRHOZ1, ZRHOZ2 )
-!
- CALL PPM_SCALAR (HLBCX,HLBCY, ISV, KTCOUNT, ZXCT, ZYCT, ZZCT, &
- PTSTEP, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
- ZSVT, ZADVS, HSV_ADV_SCHEME )
-!
- PDRIFTP(:,:,:) = ZADVS(:,:,:,1)
-!
-! Additional term: -N . DIV(mu E)
- ZDRIFTX = ZDRIFTX*MXM(PRHODJ)
- ZDRIFTY = ZDRIFTY*MYM(PRHODJ)
- ZDRIFTZ = ZDRIFTZ*MZM(1,IKU,1,PRHODJ)
-!
- CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZDRIFTX,ZDRIFTY,ZDRIFTZ,ZXCT)
+CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZDRIFTX,ZDRIFTY,ZDRIFTZ,PDRIFTP)
!
- PDRIFTP(:,:,:) = PDRIFTP(:,:,:)-ZXCT(:,:,:)*PSVT(:,:,:,NSV_ELECBEG)
-ENDIF
-!
-!
-!* 4.2.2 negative ion source
+!* 3.2 negative ion source (drifting counter E)
!
-IF (HDRIFT /= 'PPM') THEN
! x-component of div term
- ZDRIFTX(:,:,:) = PSVT(:,:,:,NSV_ELECEND) * XMOBIL_NEG(:,:,:) * XEFIELDU(:,:,:) &
- *PRHODJ(:,:,:)
+ZDRIFTX(:,:,:) = PSVT(:,:,:,NSV_ELECEND) * XMOBIL_NEG(:,:,:)
+ZDRIFTX(:,:,:) = ZDRIFTX(:,:,:) * XEFIELDU(:,:,:)
+ZDRIFTX(:,:,:) = +MXM(ZDRIFTX(:,:,:)) ! Put components at flux sides
+!
! y-component of div term
- ZDRIFTY(:,:,:) = PSVT(:,:,:,NSV_ELECEND) * XMOBIL_NEG(:,:,:) * XEFIELDV(:,:,:) &
- * PRHODJ(:,:,:)
-! z-component of div term
- ZDRIFTZ(:,:,:) = PSVT(:,:,:,NSV_ELECEND) * XMOBIL_NEG(:,:,:) * XEFIELDW(:,:,:) &
- * PRHODJ(:,:,:)
-ELSE
- ZDRIFTX(:,:,:) = -XMOBIL_NEG(:,:,:) * XEFIELDU(:,:,:)
- ZDRIFTY(:,:,:) = -XMOBIL_NEG(:,:,:) * XEFIELDV(:,:,:)
- ZDRIFTZ(:,:,:) = -XMOBIL_NEG(:,:,:) * XEFIELDW(:,:,:)
-ENDIF
+ZDRIFTY(:,:,:) = PSVT(:,:,:,NSV_ELECEND) * XMOBIL_NEG(:,:,:)
+ZDRIFTY(:,:,:) = ZDRIFTY(:,:,:) * XEFIELDV(:,:,:)
+ZDRIFTY(:,:,:) = +MYM(ZDRIFTY(:,:,:)) ! Put components at flux sides
!
-! Put components at flux sides
-ZDRIFTX(:,:,:) = MXM(ZDRIFTX(:,:,:))
-ZDRIFTY(:,:,:) = MYM(ZDRIFTY(:,:,:))
-ZDRIFTZ(:,:,:) = MZM(1,IKU,1,ZDRIFTZ (:,:,:))
+! z-component of div term
+ZDRIFTZ(:,:,:) = PSVT(:,:,:,NSV_ELECEND) * XMOBIL_NEG(:,:,:)
+ZDRIFTZ(:,:,:) = ZDRIFTZ(:,:,:) * XEFIELDW(:,:,:)
+ZDRIFTZ(:,:,:) = +MZM(1,IKU,1,ZDRIFTZ(:,:,:)) ! Put components at flux sides
+
!
IF (LWEST_ll( )) ZDRIFTX(IIB-1,:,:) = ZDRIFTX(IIB,:,:)
+IF (LEAST_ll( )) ZDRIFTX(IIE+1,:,:) = ZDRIFTX(IIE,:,:)
IF (LSOUTH_ll( )) ZDRIFTY(:,IJB-1,:) = ZDRIFTY(:,IJB,:)
+IF (LNORTH_ll( )) ZDRIFTY(:,IJE+1,:) = ZDRIFTY(:,IJE,:)
ZDRIFTZ(:,:,IKB-1) = ZDRIFTZ(:,:,IKB)
+ZDRIFTZ(:,:,IKE+1) = ZDRIFTZ(:,:,IKE)
!
-IF (HDRIFT /= 'PPM') THEN
- CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZDRIFTX,ZDRIFTY,ZDRIFTZ,PDRIFTM)
-ELSE
- CALL CONTRAV (HLBCX,HLBCY,ZDRIFTX,ZDRIFTY,ZDRIFTZ,XDXX,XDYY,XDZZ, &
- XDZX,XDZY,ZXCT,ZYCT,ZZCT,4)
-!
- ZXCT = ZXCT * PTSTEP
- ZYCT = ZYCT * PTSTEP
- ZZCT = ZZCT * PTSTEP
-!
- ZADVS(:,:,:,1) = 0.
- ZSVT(:,:,:,1) = PSVT(:,:,:,NSV_ELECEND)
- HSV_ADV_SCHEME = 'PPM_01'
-!
- CALL PPM_RHODJ(HLBCX,HLBCY, ZXCT, ZYCT, ZZCT, &
- PTSTEP, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, &
- ZRHOZ1, ZRHOZ2 )
-!
- CALL PPM_SCALAR (HLBCX,HLBCY, ISV, KTCOUNT, ZXCT, ZYCT, ZZCT, &
- PTSTEP, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
- ZSVT, ZADVS, HSV_ADV_SCHEME )
-!
- PDRIFTM(:,:,:) = ZADVS(:,:,:,1)
-!
-! Additional term -N . DIV(-mu E)
- ZDRIFTX = ZDRIFTX * MXM(PRHODJ)
- ZDRIFTY = ZDRIFTY * MYM(PRHODJ)
- ZDRIFTZ = ZDRIFTZ * MZM(1,IKU,1,PRHODJ)
-!
- CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZDRIFTX,ZDRIFTY,ZDRIFTZ,ZXCT)
-!
- PDRIFTM(:,:,:) = PDRIFTM(:,:,:) - ZXCT(:,:,:) * PSVT(:,:,:,NSV_ELECEND)
-!
- DEALLOCATE (ZXCT, ZYCT, ZZCT, ZADVS, ZSVT)
-ENDIF
+CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZDRIFTX,ZDRIFTY,ZDRIFTZ,PDRIFTM)
!
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/lesn.f90 b/src/MNH/lesn.f90
index 3823fb966..897c1a861 100644
--- a/src/MNH/lesn.f90
+++ b/src/MNH/lesn.f90
@@ -110,6 +110,11 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE ::ZZZ_LES
REAL, DIMENSION(:,:,:), ALLOCATABLE ::ZINPRR3D_LES ! precipitation flux 3D
REAL, DIMENSION(:,:,:), ALLOCATABLE ::ZEVAP3D_LES !evaporation 3D
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZP_LES ! pres. on LES vertical grid
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDP_LES ! dynamical production TKE
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTP_LES ! thermal production TKE
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTR_LES ! transport production TKE
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDISS_LES ! dissipation TKE
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZLM_LES ! mixing length
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDPDZ_LES ! dp/dz on LES vertical grid
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDTHLDZ_LES ! dThl/dz on LES vertical grid
@@ -197,6 +202,7 @@ INTEGER :: IIMAX_ll, IJMAX_ll ! total physical domain I size
INTEGER :: JLOOP
!
INTEGER :: IMASK ! mask counter
+INTEGER :: IMASKUSER! mask user number
!
INTEGER :: IRESP, ILUOUT
INTEGER :: IMI ! Current model index
@@ -239,6 +245,11 @@ END IF
! -----------
!
ALLOCATE(ZP_LES (IIU,IJU,NLES_K))
+ALLOCATE(ZDP_LES (IIU,IJU,NLES_K))
+ALLOCATE(ZTP_LES (IIU,IJU,NLES_K))
+ALLOCATE(ZTR_LES (IIU,IJU,NLES_K))
+ALLOCATE(ZDISS_LES (IIU,IJU,NLES_K))
+ALLOCATE(ZLM_LES (IIU,IJU,NLES_K))
ALLOCATE(ZDTHLDZ_LES(IIU,IJU,NLES_K))
ALLOCATE(ZDTHDZ_LES(IIU,IJU,NLES_K))
ALLOCATE(ZDRTDZ_LES(IIU,IJU,NLES_K))
@@ -486,6 +497,11 @@ END IF
!
CALL LES_VER_INT( XZZ , ZZZ_LES)
CALL LES_VER_INT( XPABST, ZP_LES )
+CALL LES_VER_INT( XDYP, ZDP_LES )
+CALL LES_VER_INT( XTHP, ZTP_LES )
+CALL LES_VER_INT( XTR, ZTR_LES )
+CALL LES_VER_INT( XDISS, ZDISS_LES )
+CALL LES_VER_INT( XLEM, ZLM_LES )
CALL LES_VER_INT( GZ_M_M(1,IKU,1,XPABST,XDZZ), ZDPDZ_LES )
!
CALL LES_VER_INT( MXF(XUT) ,ZU_LES )
@@ -642,6 +658,21 @@ END IF
!
CALL LES_MEAN_ll ( ZP_LES, LLES_CURRENT_CART_MASK, &
XLES_MEAN_P(:,NLES_CURRENT_TCOUNT,1) )
+!
+ CALL LES_MEAN_ll ( ZDP_LES, LLES_CURRENT_CART_MASK, &
+ XLES_MEAN_DP(:,NLES_CURRENT_TCOUNT,1) )
+!
+ CALL LES_MEAN_ll ( ZTP_LES, LLES_CURRENT_CART_MASK, &
+ XLES_MEAN_TP(:,NLES_CURRENT_TCOUNT,1) )
+!
+ CALL LES_MEAN_ll ( ZTR_LES, LLES_CURRENT_CART_MASK, &
+ XLES_MEAN_TR(:,NLES_CURRENT_TCOUNT,1) )
+!
+ CALL LES_MEAN_ll ( ZDISS_LES, LLES_CURRENT_CART_MASK, &
+ XLES_MEAN_DISS(:,NLES_CURRENT_TCOUNT,1) )
+!
+ CALL LES_MEAN_ll ( ZLM_LES, LLES_CURRENT_CART_MASK, &
+ XLES_MEAN_LM(:,NLES_CURRENT_TCOUNT,1) )
!
CALL LES_MEAN_ll ( ZRHO_LES, LLES_CURRENT_CART_MASK, &
XLES_MEAN_RHO(:,NLES_CURRENT_TCOUNT,1) )
@@ -998,6 +1029,11 @@ DEALLOCATE(ZTHL_LES)
DEALLOCATE(ZRT_LES)
DEALLOCATE(ZSV_LES)
DEALLOCATE(ZP_LES )
+DEALLOCATE(ZDP_LES )
+DEALLOCATE(ZTP_LES )
+DEALLOCATE(ZTR_LES )
+DEALLOCATE(ZDISS_LES )
+DEALLOCATE(ZLM_LES )
DEALLOCATE(ZDPDZ_LES)
DEALLOCATE(ZLWP_ANOM)
DEALLOCATE(ZWORK2D)
@@ -1277,6 +1313,31 @@ IF (LLES_MEAN .AND. IMASK > 1) THEN
CALL LES_MEAN_ll ( ZP_LES, OMASK, &
XLES_MEAN_P(:,NLES_CURRENT_TCOUNT,IMASK) )
!
+!* dynamical production TKE
+!
+ CALL LES_MEAN_ll ( ZDP_LES, OMASK, &
+ XLES_MEAN_DP(:,NLES_CURRENT_TCOUNT,IMASK) )
+!
+!* thermal production TKE
+!
+ CALL LES_MEAN_ll ( ZTP_LES, OMASK, &
+ XLES_MEAN_TP(:,NLES_CURRENT_TCOUNT,IMASK) )
+!
+!* transport TKE
+!
+ CALL LES_MEAN_ll ( ZTR_LES, OMASK, &
+ XLES_MEAN_TR(:,NLES_CURRENT_TCOUNT,IMASK) )
+!
+!* dissipation TKE
+!
+ CALL LES_MEAN_ll ( ZDISS_LES, OMASK, &
+ XLES_MEAN_DISS(:,NLES_CURRENT_TCOUNT,IMASK) )
+!
+!* mixing length
+!
+ CALL LES_MEAN_ll ( ZLM_LES, OMASK, &
+ XLES_MEAN_LM(:,NLES_CURRENT_TCOUNT,IMASK) )
+!
!* density
!
CALL LES_MEAN_ll ( ZRHO_LES, OMASK, &
diff --git a/src/MNH/modd_elec_descr.f90 b/src/MNH/modd_elec_descr.f90
index d5ec89786..649d0dc14 100644
--- a/src/MNH/modd_elec_descr.f90
+++ b/src/MNH/modd_elec_descr.f90
@@ -34,6 +34,8 @@
!! M. Chong 26/01/10 Small ions parameters
!! +Option for Fair weather field from
!! Helsdon-Farley (JGR, 1987, 5661-5675)
+!! Add "Beard" effect via sedimentation process
+!! J.-P. Pinty 25/10/13 Add "Latham" effect via aggregation process
!!
!-------------------------------------------------------------------------------
!
@@ -56,7 +58,7 @@ LOGICAL :: LRELAX2FW_ION = .FALSE. ! .T.= Relaxation to fair weather ion
! layer
LOGICAL :: LFLASH_GEOM=.TRUE. ! .T.: the 'geometric' flash scheme is used
LOGICAL :: LSAVE_COORD=.FALSE. ! .T.: the flash coord are written in an ascii file
-INTEGER :: NFLASH_WRITE = 100 ! Number of flashes to be saved before writing
+INTEGER :: NFLASH_WRITE = 1000 ! Number of flashes to be saved before writing
! the diag and/or coordinates in ascii files
LOGICAL :: LINDUCTIVE=.FALSE. ! .T.: inductive process is taken into account
LOGICAL :: LLNOX_EXPLICIT=.FALSE. ! .T.: lnox production is computed
@@ -157,7 +159,7 @@ INTEGER :: NNBLIGHT=0 ! Nb of lightning flashes
REAL, DIMENSION(:), ALLOCATABLE :: XNEUT_POS, XNEUT_NEG
INTEGER :: NNB_CG ! Nb of CG flashes
INTEGER :: NNB_CG_POS ! Nb of positive CG flashes
-REAL :: XALT_CG ! Altitude at which CG are detected
+REAL :: XALT_CG ! Altitude (m) at which CG are detected
!
CHARACTER(LEN=10), DIMENSION(8) &
:: CELECNAMES=(/'QNIONP','QCELEC','QRELEC','QIELEC','QSELEC', &
@@ -169,4 +171,11 @@ REAL :: XLNOX_ECLAIR
!
REAL, DIMENSION(:,:), ALLOCATABLE :: XEPOTFW_TOP
!
+! Parameters relative to the "Beard" effect ELEC=>MICROPHYS
+!
+LOGICAL :: LSEDIM_BEARD=.FALSE. ! .T.: to enable ELEC=>MICROPHYS via
+! ! particule sedimentation rate
+LOGICAL :: LIAGGS_LATHAM=.FALSE. ! .T.: to enable ELEC=>MICROPHYS via
+! ! ice aggregation rate
+!
END MODULE MODD_ELEC_DESCR
diff --git a/src/MNH/modd_elec_flash.f90 b/src/MNH/modd_elec_flash.f90
new file mode 100644
index 000000000..6a95d19d3
--- /dev/null
+++ b/src/MNH/modd_elec_flash.f90
@@ -0,0 +1,46 @@
+! ######################
+ MODULE MODD_ELEC_FLASH
+! ######################
+!
+!!**** *MODD_ELEC_FLASH* - declaration of flash map arrays
+!!
+!! PURPOSE
+!! -------
+! The purpose of this declarative module is to define the storage of the
+! flash maps
+!
+!!
+!!** IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty *Laboratoire Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 15/11/13
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+IMPLICIT NONE
+!
+!* data records: accumulated flash locations
+!
+INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: NMAP_TRIG_IC ! X-Y trig point
+INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: NMAP_IMPACT_CG ! X-Y ground impact
+INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: NMAP_2DAREA_IC ! IC: 2D max extent
+INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: NMAP_2DAREA_CG ! CG: 2D max extent
+INTEGER, DIMENSION(:,:,:), SAVE, ALLOCATABLE :: NMAP_3DIC ! reached gridpoints
+INTEGER, DIMENSION(:,:,:), SAVE, ALLOCATABLE :: NMAP_3DCG ! reached gridpoints
+!
+!------------------------------------------------------------------------------
+!
+END MODULE MODD_ELEC_FLASH
diff --git a/src/MNH/modd_elec_param.f90 b/src/MNH/modd_elec_param.f90
index 66679ab3f..d05c1fac6 100644
--- a/src/MNH/modd_elec_param.f90
+++ b/src/MNH/modd_elec_param.f90
@@ -87,9 +87,6 @@ REAL, SAVE :: XFQSDRYG, XFQSDRYGB, XFQRDRYG ! Constant in SDRYG and RDRYG
REAL, DIMENSION(:,:), SAVE, ALLOCATABLE :: XKER_Q_LIMSG
REAL, DIMENSION(:,:), SAVE, ALLOCATABLE &
:: XKER_Q_SDRYGB, XKER_Q_SDRYGB1, XKER_Q_SDRYGB2
-REAL, DIMENSION(:,:), SAVE, ALLOCATABLE &
- :: XFTAKA ! F(T,LWC) for Takahashi (1978) charge
- ! generation parametrization
!
! Helsdon-Farley
!
@@ -111,7 +108,15 @@ REAL, SAVE :: XIMP, XINP, XIKP, & ! Parameters m, n and k
XIMN, XINN, XIKN, & ! for the NI processes
XSMP, XSNP, XSKP, & ! following
XSMN, XSNN, XSKN ! Saunders et al. (1991)
-REAL, SAVE :: XINISAP
+REAL, SAVE :: XFQIAGGSP, XFQIAGGSN, & ! Auxiliary parameters
+ XFQIDRYGBSP, XFQIDRYGBSN, & ! containing MOMG function
+ XLBQSDRYGB1SP, XLBQSDRYGB1SN, &
+ XLBQSDRYGB2SP, XLBQSDRYGB2SN, &
+ XLBQSDRYGB3SP, XLBQSDRYGB3SN, &
+ XAIGAMMABI
+REAL, SAVE :: XIKP_TAK, XIKN_TAK, XSKP_TAK, XSKN_TAK ! Using Takahashi charge
+REAL, SAVE :: XFQIAGGSP_TAK, XFQIAGGSN_TAK, XFQIDRYGBSP_TAK, XFQIDRYGBSN_TAK
+REAL, SAVE :: XVSCOEF, XVGCOEF
REAL, SAVE :: XFQIDRYGBS, XLBQIDRYGBS ! Constants in QIDRYGB
REAL, SAVE :: XFQSDRYGBS ! Constants in QSDRYGB
REAL, SAVE :: XLBQSDRYGB1S, XLBQSDRYGB2S !
@@ -120,7 +125,9 @@ REAL, SAVE :: XLBQSDRYGB1S, XLBQSDRYGB2S !
!
INTEGER, SAVE :: NIND_TEMP ! number of indexes for temperature
INTEGER, SAVE :: NIND_LWC ! number of indexes for liquid water content
-REAL, DIMENSION(:,:), SAVE, ALLOCATABLE :: XMANSELL ! F(T,LWC) for Takahashi(1978)
+REAL, DIMENSION(:,:), SAVE, ALLOCATABLE :: XMANSELL ! F(LWC, T) for Takahashi(1978) /Mansell
+REAL, DIMENSION(:,:), SAVE, ALLOCATABLE :: XSAUNDER ! F(LWC, T) for SAUN1/SAUN2, BSMP1/BSMP2
+REAL, DIMENSION(:,:), SAVE, ALLOCATABLE :: XTAKA_TM ! F(LWC, T) for Takahashi/Tsenova and Mitzeva
!
REAL, SAVE :: XFQIDRYGBT1, XFQIDRYGBT2, XFQIDRYGBT3, & ! IDRYGB
XFQSDRYGBT1, XFQSDRYGBT2, XFQSDRYGBT3, & ! SDRYGB
@@ -156,9 +163,9 @@ REAL, SAVE :: XIND1, XIND2, XIND3
!
! lightning
!
-REAL :: XFQLIGHTC, XFQLIGHTR, &
- XFQLIGHTI, XFQLIGHTS, XFQLIGHTG ! Constant for charge redistribution
-REAL :: XEXQLIGHTR, &
- XEXQLIGHTI, XEXQLIGHTS, XEXQLIGHTG ! Exponent for charge redistribution
+REAL :: XFQLIGHTC, XFQLIGHTR, XFQLIGHTI, &
+ XFQLIGHTS, XFQLIGHTG, XFQLIGHTH ! Constant for charge redistribution
+REAL :: XEXQLIGHTR, XEXQLIGHTI, &
+ XEXQLIGHTS, XEXQLIGHTG, XEXQLIGHTH ! Exponent for charge redistribution
!
END MODULE MODD_ELEC_PARAM
diff --git a/src/MNH/modd_elecn.f90 b/src/MNH/modd_elecn.f90
index 9d84f0304..9c50af2ba 100644
--- a/src/MNH/modd_elecn.f90
+++ b/src/MNH/modd_elecn.f90
@@ -54,7 +54,9 @@ TYPE ELEC_t
! Parameters for flat lapalcian operator to solve the electric field
! (see MODD_DYN_n)
REAL, DIMENSION(:), POINTER :: XRHOM_E =>NULL(), XAF_E =>NULL(), XCF_E =>NULL()
- REAL, DIMENSION(:,:,:), POINTER :: XBFY_E =>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XBFY_E =>NULL(), &
+ XBFB_E =>NULL(), XBF_SXP2_YP1_Z_E =>NULL()
+ ! Z_Splitting
!
END TYPE ELEC_t
@@ -66,7 +68,8 @@ REAL, DIMENSION(:,:,:), POINTER :: XNI_SDRYG=>NULL(), XNI_IDRYG=>NULL(), &
XESOURCEFW=>NULL(), XEFIELDV=>NULL(), XEFIELDW=>NULL(), &
XIND_RATE=>NULL(), XIONSOURCEFW =>NULL(), XEW=>NULL(), &
XCION_POS_FW =>NULL(), XCION_NEG_FW =>NULL(), &
- XMOBIL_POS =>NULL(), XMOBIL_NEG=>NULL(), XBFY_E =>NULL()
+ XMOBIL_POS =>NULL(), XMOBIL_NEG=>NULL(), XBFY_E =>NULL(), &
+ XBFB_E =>NULL(), XBF_SXP2_YP1_Z_E =>NULL()
REAL, DIMENSION(:), POINTER :: XRHOM_E =>NULL(), XAF_E =>NULL(), XCF_E =>NULL()
CONTAINS
@@ -90,6 +93,8 @@ ELEC_MODEL(KFROM)%XCION_NEG_FW=>XCION_NEG_FW
ELEC_MODEL(KFROM)%XMOBIL_POS=>XMOBIL_POS
ELEC_MODEL(KFROM)%XMOBIL_NEG=>XMOBIL_NEG
ELEC_MODEL(KFROM)%XBFY_E=>XBFY_E
+ELEC_MODEL(KFROM)%XBFB_E=>XBFB_E
+ELEC_MODEL(KFROM)%XBF_SXP2_YP1_Z_E=>XBF_SXP2_YP1_Z_E
ELEC_MODEL(KFROM)%XRHOM_E=>XRHOM_E
ELEC_MODEL(KFROM)%XAF_E=>XAF_E
ELEC_MODEL(KFROM)%XCF_E=>XCF_E
@@ -110,6 +115,8 @@ XCION_NEG_FW=>ELEC_MODEL(KTO)%XCION_NEG_FW
XMOBIL_POS=>ELEC_MODEL(KTO)%XMOBIL_POS
XMOBIL_NEG=>ELEC_MODEL(KTO)%XMOBIL_NEG
XBFY_E=>ELEC_MODEL(KTO)%XBFY_E
+XBFB_E=>ELEC_MODEL(KTO)%XBFB_E
+XBF_SXP2_YP1_Z_E=>ELEC_MODEL(KTO)%XBF_SXP2_YP1_Z_E
XRHOM_E=>ELEC_MODEL(KTO)%XRHOM_E
XAF_E=>ELEC_MODEL(KTO)%XAF_E
XCF_E=>ELEC_MODEL(KTO)%XCF_E
diff --git a/src/MNH/modd_lma_simulator.f90 b/src/MNH/modd_lma_simulator.f90
new file mode 100644
index 000000000..0724bd3ee
--- /dev/null
+++ b/src/MNH/modd_lma_simulator.f90
@@ -0,0 +1,64 @@
+! ############################
+ MODULE MODD_LMA_SIMULATOR
+! ############################
+!
+!!**** *MODD_LMA_SIMULATOR* - declaration of LMA network
+!!
+!! PURPOSE
+!! -------
+! The purpose of this declarative module is to define the storage of the
+! LMA simulator
+!
+!!
+!!** IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty *Laboratoire Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 15/02/13
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_TYPE_DATE
+!
+IMPLICIT NONE
+!
+!
+!* general information
+!
+LOGICAL, SAVE :: LLMA=.FALSE.! Flag to record LMA-like
+ ! data along the
+ ! simulation
+REAL, SAVE :: XDTLMA ! Time length of a LMA record
+TYPE (DATE_TIME), SAVE :: TDTLMA ! Date and Time of LMA file
+CHARACTER (LEN=31), SAVE :: CLMA_FILE ! File name
+INTEGER, SAVE :: ILMA_UNIT ! File information
+INTEGER, SAVE :: ILMA_IOSTAT ! File information
+!
+!* storage monitoring
+!
+INTEGER, DIMENSION(:,:,:), SAVE, ALLOCATABLE :: ISLMA_SEG_GLOB ! Global indexes
+ ! of the LMA segments
+!
+!* data records
+!
+LOGICAL, SAVE :: LWRITE_LMA
+REAL, DIMENSION(:,:), SAVE, ALLOCATABLE :: ZLMA_LAT, ZLMA_LON
+REAL, DIMENSION(:,:), SAVE, ALLOCATABLE :: ZSLMA_NEUT_POS, ZSLMA_NEUT_NEG
+REAL, DIMENSION(:,:,:), SAVE, ALLOCATABLE :: ZSLMA_QMT
+REAL, DIMENSION(:,:,:), SAVE, ALLOCATABLE :: ZSLMA_PRT
+!
+!------------------------------------------------------------------------------
+!
+END MODULE MODD_LMA_SIMULATOR
diff --git a/src/MNH/modn_elec.f90 b/src/MNH/modn_elec.f90
index 813fafaa1..7c9e53f3c 100644
--- a/src/MNH/modn_elec.f90
+++ b/src/MNH/modn_elec.f90
@@ -26,6 +26,8 @@
!! -------------
!! Original 09/11/09
!! M. Chong 26/01/10 Option for fair weather field from Helsdon-Farley
+!! M. Chong 26/08/13 Option for "Beard" effect and LLMA storage
+!! J.-P. Pinty 25/10/13 Option for "Latham" effect
!!
!! IMPLICIT ARGUMENTS
!! ------------------
@@ -37,6 +39,7 @@
USE MODD_ELEC_DESCR
USE MODD_ELEC_PARAM
+USE MODD_LMA_SIMULATOR
!
IMPLICIT NONE
!
@@ -50,6 +53,7 @@ NAMELIST /NAM_ELEC/ LOCG, LELEC_FIELD, LFLASH_GEOM, &
CLSOL, NLAPITR_ELEC, XRELAX_ELEC, &
XETRIG, XEBALANCE, XEPROP, XQEXCES, XQNEUT,&
XDFRAC_ECLAIR, XDFRAC_L, &
- XWANG_A, XWANG_B
+ XWANG_A, XWANG_B, &
+ LSEDIM_BEARD, LIAGGS_LATHAM, LLMA
!
END MODULE MODN_ELEC
diff --git a/src/MNH/read_desfmn.f90 b/src/MNH/read_desfmn.f90
index 960404a5e..5593a4ab1 100644
--- a/src/MNH/read_desfmn.f90
+++ b/src/MNH/read_desfmn.f90
@@ -21,6 +21,7 @@ INTERFACE
#ifdef MNH_FOREFIRE
OFOREFIRE, &
#endif
+ OLNOX_EXPLICIT, &
OCONDSAMP, &
KRIMX,KRIMY,KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC,HEQNSYS )
@@ -50,6 +51,7 @@ LOGICAL, INTENT(OUT) :: OPASPOL ! Passive pollutant flag
#ifdef MNH_FOREFIRE
LOGICAL, INTENT(OUT) :: OFOREFIRE! ForeFire flag
#endif
+LOGICAL, INTENT(OUT) :: OLNOX_EXPLICIT ! explicit LNOx flag
LOGICAL, INTENT(OUT) :: OCONDSAMP! Conditional sampling flag
LOGICAL, INTENT(OUT) :: OORILAM ! Orilam flag
LOGICAL, INTENT(OUT) :: OCHTRANS ! Deep convection on scalar
@@ -84,6 +86,7 @@ END MODULE MODI_READ_DESFM_n
#ifdef MNH_FOREFIRE
OFOREFIRE, &
#endif
+ OLNOX_EXPLICIT, &
OCONDSAMP, &
KRIMX,KRIMY,KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC,HEQNSYS )
@@ -186,6 +189,7 @@ END MODULE MODI_READ_DESFM_n
!! Modification 03/2006 (O.Geoffroy) Add KHKO scheme
!! Modification 04/2010 (M. Leriche) Add aqueous + ice chemistry
!! Modification 07/2013 (Bosseur & Filippi) Adds Forefire
+!! Modification 01/2015 (C. Barthe) Add explicit LNOx
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -279,6 +283,7 @@ LOGICAL, INTENT(OUT) :: OPASPOL ! Passive pollutant flag
#ifdef MNH_FOREFIRE
LOGICAL, INTENT(OUT) :: OFOREFIRE ! ForeFire flag
#endif
+LOGICAL, INTENT(OUT) :: OLNOX_EXPLICIT ! explicit LNOx flag
LOGICAL, INTENT(OUT) :: OCONDSAMP! Conditional sampling flag
LOGICAL, INTENT(OUT) :: ODUST ! Dust flag
LOGICAL, INTENT(OUT) :: OORILAM ! Dust flag
@@ -509,6 +514,7 @@ OPASPOL = LPASPOL
#ifdef MNH_FOREFIRE
OFOREFIRE = LFOREFIRE
#endif
+OLNOX_EXPLICIT = LLNOX_EXPLICIT
OCONDSAMP= LCONDSAMP
KRIMX = NRIMX
KRIMY = NRIMY
@@ -674,10 +680,10 @@ IF (NVERB >= 10) THEN
WRITE(UNIT=ILUOUT,NML=NAM_PARAM_C1R3)
END IF
!
- IF (CELEC /= 'NONE') THEN
- WRITE(UNIT=ILUOUT,FMT="('************ ELEC SCHEME **********************')")
- WRITE(UNIT=ILUOUT,NML=NAM_ELEC)
- END IF
+ IF (CELEC /= 'NONE') THEN
+ WRITE(UNIT=ILUOUT,FMT="('************ ELEC SCHEME **********************')")
+ WRITE(UNIT=ILUOUT,NML=NAM_ELEC)
+ END IF
!
END IF
!
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
index c2158e74a..7a5060ed0 100644
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -6,6 +6,7 @@
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$
+! masdev4_8 2008/07/09 16:40:30
!-----------------------------------------------------------------
! ######################
MODULE MODI_READ_EXSEG_n
@@ -21,6 +22,7 @@ INTERFACE
#ifdef MNH_FOREFIRE
OFOREFIRE, &
#endif
+ OLNOX_EXPLICIT, &
OCONDSAMP, &
KRIMX,KRIMY, KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC, &
@@ -51,6 +53,7 @@ LOGICAL, INTENT(IN) :: OPASPOL ! Passive pollutant FLAG in FMF
#ifdef MNH_FOREFIRE
LOGICAL, INTENT(IN) :: OFOREFIRE ! ForeFire FLAG in FMFILE
#endif
+LOGICAL, INTENT(IN) :: OLNOX_EXPLICIT ! explicit LNOx FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCONDSAMP ! Conditional sampling FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCHTRANS ! LCHTRANS FLAG in FMFILE
@@ -89,6 +92,7 @@ END MODULE MODI_READ_EXSEG_n
#ifdef MNH_FOREFIRE
OFOREFIRE, &
#endif
+ OLNOX_EXPLICIT, &
OCONDSAMP, &
KRIMX,KRIMY, KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC, &
@@ -179,7 +183,7 @@ END MODULE MODI_READ_EXSEG_n
!! Module MODN_PARAM1 : CTURB,CRAD,CDCONV,CSCONV
!!
!! Module MODN_LUNIT1 :
-!! Module MODN_LBC1 : CLBCX,CLBCY,NLBLX,NLBLY,XCPHASE
+!! Module MODN_LBC1 : CLBCX,CLBCY,NLBLX,NLBLY,XCPHASE,XPOND
!!
!! Module MODN_TURB_n : CTURBLEN,CTURBDIM
!!
@@ -273,6 +277,7 @@ END MODULE MODI_READ_EXSEG_n
!! Modification 12/2012 (S.Bielli) add NAM_NCOUT for netcdf output
!! Modification 02/2012 (Pialat/Tulet) add ForeFire
!! Modification 02/2012 (T.Lunet) add of new Runge-Kutta methods
+!! Modification 01/2015 (C. Barthe) add explicit LNOx
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!!------------------------------------------------------------------------------
!
@@ -384,6 +389,7 @@ LOGICAL, INTENT(IN) :: OPASPOL ! Passive pollutant FLAG in FMF
#ifdef MNH_FOREFIRE
LOGICAL, INTENT(IN) :: OFOREFIRE ! ForeFire FLAG in FMFILE
#endif
+LOGICAL, INTENT(IN) :: OLNOX_EXPLICIT ! explicit LNOx FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCONDSAMP ! Conditional sampling FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCHTRANS ! LCHTRANS FLAG in FMFILE
@@ -618,7 +624,8 @@ CALL TEST_NAM_VAR(ILUOUT,'CTURBLEN_CLOUD',CTURBLEN_CLOUD,'NONE','DEAR','DELT','B
!
! The test on the mass flux scheme for shallow convection
!
-CALL TEST_NAM_VAR(ILUOUT,'CMF_UPDRAFT',CMF_UPDRAFT,'NONE','EDKF','RHCJ')
+CALL TEST_NAM_VAR(ILUOUT,'CMF_UPDRAFT',CMF_UPDRAFT,'NONE','EDKF','RHCJ',&
+ 'HRIO','BOUT')
CALL TEST_NAM_VAR(ILUOUT,'CMF_CLOUD',CMF_CLOUD,'NONE','STAT','DIRE')
!
! The test on the CSOLVER name is made elsewhere
@@ -635,9 +642,6 @@ END IF
!* 2. FIRST INITIALIZATIONS
! ---------------------
!
-!!!!!!!!!!!!!!!!!!!! TEST CL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!IF (NPROC==1) JPHEXT=1
-!!!!!!!!!!!!!!!!!!!! TEST CL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!* 2.1 Time step in gridnesting case
!
IF (KMI /= 1 .AND. NDAD(KMI) /= KMI) THEN
@@ -1357,6 +1361,20 @@ IF (CELEC /= 'NONE') THEN
END IF
END IF
!
+! (explicit) LINOx SV case
+!
+IF (CELEC /= 'NONE' .AND. LLNOX_EXPLICIT) THEN
+ IF (HELEC /= 'NONE' .AND. OLNOX_EXPLICIT) THEN
+ CGETSVT(NSV_LNOXBEG:NSV_LNOXEND)='READ'
+ ELSE
+ WRITE(UNIT=ILUOUT,FMT=9001) KMI
+ WRITE(UNIT=ILUOUT,FMT='("THERE IS NO SCALAR VARIABLES FOR LINOX &
+ & IN INITIAL FMFILE",/,&
+ & "THE LINOX VARIABLES HAVE BEEN INITIALIZED TO ZERO ")')
+ CGETSVT(NSV_LNOXBEG:NSV_LNOXEND)='INIT'
+ END IF
+END IF
+!
! Chemical SV case (including aqueous chemical species)
!
IF (LUSECHEM) THEN
@@ -2493,13 +2511,6 @@ IF (XTNUDGING < 4.*XTSTEP) THEN
WRITE(ILUOUT,FMT=*) 'XTNUDGING is SET TO ',XTNUDGING
END IF
!
-IF (XBULEN < 2.*XTSTEP .AND. KMI==NBUMOD) THEN
- XBULEN = 2.*XTSTEP
- WRITE(UNIT=ILUOUT,FMT=9002) KMI
- WRITE(UNIT=ILUOUT,FMT='("DURATION FOR BUDGET AVERAGING CAN NOT BE SMALLER THAN", &
- & " TWO TIMES THE TIME STEP")')
- WRITE(ILUOUT,FMT=*) 'XBULEN is SET TO ',XBULEN
-END IF
!
IF (XWAY(KMI) == 3. ) THEN
XWAY(KMI) = 2.
diff --git a/src/MNH/resolved_elecn.f90 b/src/MNH/resolved_elecn.f90
index 40a40c65b..842ab3bd3 100644
--- a/src/MNH/resolved_elecn.f90
+++ b/src/MNH/resolved_elecn.f90
@@ -1,3 +1,4 @@
+
!MNH_LIC Copyright 1994-2014 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
@@ -8,7 +9,7 @@
!
INTERFACE
SUBROUTINE RESOLVED_ELEC_n (HCLOUD, HSCONV, HMF_CLOUD, &
- KRR, KSPLITR, KMI, KTCOUNT, &
+ KRR, KSPLITR, KMI, KTCOUNT, OEXIT, &
HLBCX, HLBCY, HFMFILE, HLUOUT, HRAD, HTURBDIM, &
OCLOSE_OUT, OSUBG_COND, OSIGMAS,PSIGQSAT, HSUBG_AUCV, &
PTSTEP, PZZ, PRHODJ, PRHODREF, PEXNREF, &
@@ -25,9 +26,10 @@ CHARACTER(LEN=4), INTENT(IN) :: HSCONV ! Shallow convection scheme
CHARACTER(LEN=4), INTENT(IN) :: HMF_CLOUD! Type of statistical cloud
INTEGER, INTENT(IN) :: KRR ! Number of moist variables
INTEGER, INTENT(IN) :: KSPLITR ! Number of small time step
- ! integrations for rain sedimendation
+ ! integrations for rain sedimendation
INTEGER, INTENT(IN) :: KMI ! Model index
INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop counter
+LOGICAL, INTENT(IN) :: OEXIT ! switch for the end of the temporal loop
CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY ! X and Y-direc. LBC type
CHARACTER(LEN=*), INTENT(IN) :: HFMFILE ! Name of the output FM-file
CHARACTER(LEN=*), INTENT(IN) :: HLUOUT ! Output-listing name for
@@ -36,16 +38,16 @@ CHARACTER*4, INTENT(IN) :: HRAD ! Radiation scheme name
CHARACTER*4, INTENT(IN) :: HTURBDIM ! Dimensionality of the
! turbulence scheme
LOGICAL, INTENT(IN) :: OCLOSE_OUT ! Conditional closure of
- ! the OUTPUT FM-file
+ ! the OUTPUT FM-file
LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for Subgrid Cond.
LOGICAL, INTENT(IN) :: OSIGMAS ! Switch for Sigma_s:
- ! use values computed in CONDENSATION
- ! or that from turbulence scheme
+ ! use values computed in CONDENSATION
+ ! or that from turbulence scheme
REAL, INTENT(IN) :: PSIGQSAT ! coeff applied to qsat variance contribution
CHARACTER(LEN=4), INTENT(IN) :: HSUBG_AUCV
- ! Kind of Subgrid autoconversion method
+ ! Kind of Subgrid autoconversion method
REAL, INTENT(IN) :: PTSTEP ! Double Time step
- ! (single if cold start)
+ ! (single if cold start)
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ !Dry density * Jacobian
@@ -99,7 +101,7 @@ END MODULE MODI_RESOLVED_ELEC_n
!
! #####################################################################################
SUBROUTINE RESOLVED_ELEC_n (HCLOUD, HSCONV, HMF_CLOUD, &
- KRR, KSPLITR, KMI, KTCOUNT, &
+ KRR, KSPLITR, KMI, KTCOUNT, OEXIT, &
HLBCX, HLBCY, HFMFILE, HLUOUT, HRAD, HTURBDIM, &
OCLOSE_OUT, OSUBG_COND, OSIGMAS,PSIGQSAT, HSUBG_AUCV, &
PTSTEP, PZZ, PRHODJ, PRHODREF, PEXNREF, &
@@ -168,6 +170,7 @@ END MODULE MODI_RESOLVED_ELEC_n
!! Original 06/11/09
!! Modifications:
!! M. Chong 26/01/10 Add Small ions parameters
+!! M. Chong 31/07/14 Add explicit LiNOx
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!!
!-------------------------------------------------------------------------------
@@ -189,7 +192,11 @@ USE MODD_ELEC_n
USE MODD_BUDGET
USE MODD_NSV
USE MODD_CH_MNHC_n, ONLY: LUSECHEM,LCH_CONV_LINOX
-USE MODD_DYN_n, ONLY: NSTOP
+USE MODD_DYN_n, ONLY: NSTOP, XTSTEP
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+
+USE MODD_TIME_n
+USE MODD_LMA_SIMULATOR
USE MODD_PRINT_ELEC
!
USE MODI_IO_ll
@@ -219,6 +226,7 @@ INTEGER, INTENT(IN) :: KSPLITR ! Number of small time step
! integrations for ice sedimendation
INTEGER, INTENT(IN) :: KMI ! Model index
INTEGER, INTENT(IN) :: KTCOUNT ! Temporal loop counter
+LOGICAL, INTENT(IN) :: OEXIT ! switch for the end of the temporal loop
CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY ! X and Y-direc. LBC type
CHARACTER(LEN=*), INTENT(IN) :: HFMFILE ! Name of the output FM-file
CHARACTER(LEN=*), INTENT(IN) :: HLUOUT ! Output-listing name for
@@ -339,7 +347,11 @@ REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZDRIFT
INTEGER :: IPROCMIN, IK
INTEGER :: IXOR, IYOR ! origin of the extended subdomain
CHARACTER (LEN=32) :: YASCFILE
-
+!
+REAL :: ZTEMP_DIST
+CHARACTER (LEN=18) :: YNAME
+LOGICAL :: GLMA_FILE
+!
NULLIFY(TZFIELDS_ll)
!
!------------------------------------------------------------------------------
@@ -618,7 +630,7 @@ PSVT(:,:,:,NSV_ELECBEG) = XECHARGE*PSVT(:,:,:,NSV_ELECBEG) ! 1/kg --> C/kg
PSVT(:,:,:,NSV_ELECEND) =-XECHARGE*PSVT(:,:,:,NSV_ELECEND)
!
CALL TO_ELEC_FIELD_n (PRT, PSVT(:,:,:,NSV_ELECBEG:NSV_ELECEND), PRHODJ, &
- KTCOUNT, KRR, &
+ KTCOUNT, KRR, &
XEFIELDU, XEFIELDV, XEFIELDW )
!
PSVT(:,:,:,NSV_ELECBEG) = PSVT(:,:,:,NSV_ELECBEG)/XECHARGE ! back to 1/kg
@@ -634,12 +646,12 @@ CALL MYPROC_ELEC_ll (IPROC)
!
! Hereafter, ZCPH and ZCOR are used temporarily to store the drift sources
! of the positive and negative ions, respectively
+!
+CALL ION_DRIFT(ZCPH, ZCOR, PSVT, HLBCX, HLBCY)
-CALL ION_DRIFT(ZCPH, ZCOR, PSVT, PRHODREF, PRHODJ, HLBCX, HLBCY, &
- KTCOUNT, PTSTEP, CDRIFT)
-PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) + ZCPH(:,:,:)/PRHODJ(:,:,:)
-PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) + ZCOR(:,:,:)/PRHODJ(:,:,:)
+PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) + ZCPH(:,:,:)
+PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) + ZCOR(:,:,:)
!
!* 4.3 Add Cosmic Ray source
!
@@ -779,7 +791,7 @@ SELECT CASE (HCLOUD)
!
END SELECT
!
-IF(KTCOUNT .EQ. 1 .AND. IPROC.EQ.0) PRINT *,'KSPLITR=', KSPLITR
+IF(KTCOUNT .EQ. 1 .AND. IPROC .EQ. 0) PRINT *,'KSPLITR=', KSPLITR
!
!-------------------------------------------------------------------------------
!
@@ -798,6 +810,10 @@ DO JSV = NSV_ELECBEG, NSV_ELECEND
PSVS(:,:,:,JSV) = PSVS(:,:,:,JSV) * PRHODJ(:,:,:)
ENDDO
!
+! Note that the LiNOx Conc. (in mol/mol) is PSVS (:,::,NSV_LNOXBEG)
+! but there is no need to *PRHODJ(:,:,:) as it is done implicitly
+! during unit conversion in flash_geom.
+!
PSVS(:,:,:,NSV_ELECBEG) = MAX(0., PSVS(:,:,:,NSV_ELECBEG))
PSVS(:,:,:,NSV_ELECEND) = MAX(0., PSVS(:,:,:,NSV_ELECEND))
!
@@ -810,62 +826,57 @@ GATTACH(:,:,:) = .FALSE.
GATTACH(IIB:IIE, IJB:IJE, IKB:IKE) = .TRUE.
!
IF (PRESENT(PSEA)) THEN
- CALL ION_ATTACH_ELEC(KTCOUNT, KRR, PTSTEP, PRHODREF, &
- PRHODJ, PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
- PRS, PTHT, PCIT, PPABST, XEFIELDU, &
- XEFIELDV, XEFIELDW, GATTACH, PTOWN, PSEA )
+ CALL ION_ATTACH_ELEC(KTCOUNT, KRR, PTSTEP, PRHODREF, &
+ PRHODJ, PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
+ PRS, PTHT, PCIT, PPABST, XEFIELDU, &
+ XEFIELDV, XEFIELDW, GATTACH, PTOWN, PSEA )
ELSE
- CALL ION_ATTACH_ELEC(KTCOUNT, KRR, PTSTEP, PRHODREF, &
- PRHODJ, PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
- PRS, PTHT, PCIT, PPABST, XEFIELDU, &
- XEFIELDV, XEFIELDW, GATTACH )
+ CALL ION_ATTACH_ELEC(KTCOUNT, KRR, PTSTEP, PRHODREF, &
+ PRHODJ, PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
+ PRS, PTHT, PCIT, PPABST, XEFIELDU, &
+ XEFIELDV, XEFIELDW, GATTACH )
ENDIF
!
!-------------------------------------------------------------------------------
!
-!* 9. LIGHTNING FLASHES
-! -----------------
-!
-ZQTOT(:,:,:) = XECHARGE * (PSVT(:,:,:,NSV_ELECBEG) - PSVT(:,:,:,NSV_ELECEND))
-DO JSV = NSV_ELECBEG+1, NSV_ELECEND-1
- ZQTOT(:,:,:) = ZQTOT(:,:,:) + PSVT(:,:,:,JSV)
-END DO
-!
-!
-!* 9.1 open the output ascii files
+!* 9. OPEN THE OUTPUT ASCII FILES
+! ---------------------------
!
IF (KTCOUNT .EQ. 1) THEN
IF (LFLASH_GEOM) THEN
- YASCFILE = CEXP//"_fgeom_diag.asc"
+ YASCFILE = CEXP//"_fgeom_diag.asc"
CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW", &
- FORM="FORMATTED", POSITION="APPEND", &
- UNIT=NLU_fgeom_diag, IOSTAT= NIOSTAT_fgeom_diag )
- IF ( IPROC .EQ. 0) THEN
- WRITE (NLU_fgeom_diag, FMT='(A)') '--------------------------------------------------------'
- WRITE (NLU_fgeom_diag, FMT='(A)') '*******FLASH CHARACTERISTICS FROM FLASH_GEOM_ELEC*******'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 1 : total flash number --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 2 : time (s) --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 3 : cell number --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 4 : flash number/cell/time step --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 5 : flash type (1=IC, 2=CGN, 3=CGP) --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 6 : number of segments --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 7 : triggering electric field (kV/m) --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 8 : x coord. of the triggering point (km) --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 9 : y coord. of the triggering point (km) --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 10 : z coord. of the triggering point (km) --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 11 : positive charge neutralized (C) --'
- WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 12 : negative charge neutralized (C) --'
+ FORM="FORMATTED", POSITION="APPEND", &
+ UNIT=NLU_fgeom_diag, IOSTAT= NIOSTAT_fgeom_diag )
+ IF ( IPROC .EQ. 0) THEN
WRITE (NLU_fgeom_diag, FMT='(A)') '--------------------------------------------------------'
- END IF
+ WRITE (NLU_fgeom_diag, FMT='(A)') '*FLASH CHARACTERISTICS FROM FLASH_GEOM_ELEC*'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 1 : total flash number --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 2 : time (s) --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 3 : cell number --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 4 : flash number/cell/time step --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 5 : flash type 1=IC, 2=CGN, 3=CGP '
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 6 : number of segments --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 7 : trig electric field (kV/m) --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 8 : x coord. trig. point --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 9 : y coord. trig. point --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- --> x,y in km if lcartesian=t, --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- deg otherwise --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 10 : z coord. trig. point (km) --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 11: neutr. positive charge (C) --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 12: neutr. negative charge (C) --'
+ WRITE (NLU_fgeom_diag, FMT='(A)') '--------------------------------------------'
+ END IF
!
- CALL CLOSE_ll (YASCFILE)
+ CALL CLOSE_ll (YASCFILE)
+ CALL MPI_BCAST (NLU_fgeom_diag,1, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
!
- IF (LSAVE_COORD) THEN
- YASCFILE = CEXP//"_fgeom_coord.asc"
- CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW",&
- FORM="FORMATTED", POSITION="APPEND", &
- UNIT=NLU_fgeom_coord, IOSTAT= NIOSTAT_fgeom_coord )
- IF ( IPROC .EQ. 0) THEN
+ IF (LSAVE_COORD) THEN
+ YASCFILE = CEXP//"_fgeom_coord.asc"
+ CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW",&
+ FORM="FORMATTED", POSITION="APPEND", &
+ UNIT=NLU_fgeom_coord, IOSTAT= NIOSTAT_fgeom_coord )
+ IF ( IPROC .EQ. 0) THEN
WRITE (NLU_fgeom_coord,FMT='(A)') '------------------------------------------'
WRITE (NLU_fgeom_coord,FMT='(A)') '*****FLASH COORD. FROM FLASH_GEOM_ELEC****'
WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 1 : flash number --'
@@ -875,20 +886,21 @@ IF (KTCOUNT .EQ. 1) THEN
WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 5 : coordinate along Y (km) --'
WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 6 : coordinate along Z (km) --'
WRITE (NLU_fgeom_coord,FMT='(A)') '------------------------------------------'
- END IF
-!
- CALL CLOSE_ll (YASCFILE)
END IF
+!
+ CALL CLOSE_ll (YASCFILE)
+ CALL MPI_BCAST (NLU_fgeom_coord,1, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
END IF
+ END IF
!
- IF (LSERIES_ELEC) THEN
- YASCFILE = CEXP//"_series_cloud_elec.asc"
- CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW", &
- FORM="FORMATTED", POSITION="APPEND", &
- UNIT=NLU_series_cloud_elec, IOSTAT= NIOSTAT_series_cloud_elec)
- IF ( IPROC .EQ. 0) THEN
+ IF (LSERIES_ELEC) THEN
+ YASCFILE = CEXP//"_series_cloud_elec.asc"
+ CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW", &
+ FORM="FORMATTED", POSITION="APPEND", &
+ UNIT=NLU_series_cloud_elec, IOSTAT= NIOSTAT_series_cloud_elec)
+ IF ( IPROC .EQ. 0) THEN
WRITE (NLU_series_cloud_elec, FMT='(A)') '----------------------------------------------------'
- WRITE (NLU_series_cloud_elec, FMT='(A)') '********* RESULTS FROM USE of LSERIES_ELEC *********'
+ WRITE (NLU_series_cloud_elec, FMT='(A)') '********* RESULTS FROM of LSERIES_ELEC *************'
WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 1 : Time (s) --'
WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 2 : Cloud top height / Z > 20 dBZ (m) --'
WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 3 : Cloud top height / m.r. > 1.e-4 (m) --'
@@ -909,58 +921,105 @@ IF (KTCOUNT .EQ. 1) THEN
WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 18 : Cloud volume (m3) --'
WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 19 : Maximum rain inst. precip. (mm/H) --'
WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 20 : Rain instant. precip. (mm/H) --'
-
-
WRITE (NLU_series_cloud_elec, FMT='(A)') '----------------------------------------------------'
- END IF
+ END IF
!
- CALL CLOSE_ll (YASCFILE)
+ CALL CLOSE_ll (YASCFILE)
+ CALL MPI_BCAST (NLU_series_cloud_elec,1, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
END IF
+END IF
!
- IF (LFLASH_GEOM) THEN
- CALL MPI_BCAST (NLU_fgeom_diag,1, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
- IF (LSAVE_COORD) THEN
- CALL MPI_BCAST (NLU_fgeom_coord,1, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
- END IF
- ELSE
- PRINT *,' LIGHTNING_ELEC NOT YET DEVELOPPED'
- STOP
+IF (LFLASH_GEOM .AND. LLMA) THEN
+!
+! test to see if a new LMA file should be created
+!
+ GLMA_FILE = .FALSE.
+!
+ IF (TDTCUR%TIME >= TDTLMA%TIME-PTSTEP .AND. CLMA_FILE(1:5) /= "BEGIN") THEN
+ LWRITE_LMA = .TRUE.
+ END IF
+ IF (TDTCUR%TIME >= TDTLMA%TIME) THEN
+ TDTLMA%TIME = TDTLMA%TIME + XDTLMA
+ GLMA_FILE = .TRUE.
+ LWRITE_LMA = .FALSE.
END IF
!
- IF (LSERIES_ELEC) THEN
- CALL MPI_BCAST (NLU_series_cloud_elec,1, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
- ENDIF
+ IF (GLMA_FILE) THEN
+ IF(CLMA_FILE(1:5) /= "BEGIN") THEN ! close preceeding file when existing
+ CALL CLOSE_ll (CLMA_FILE)
+ ENDIF
+!
+ TDTLMA%TIME = TDTLMA%TIME - XDTLMA
+ WRITE (YNAME,FMT='(3I2.2,A1,3I2.2,A1,I4.4)') &
+ ABS(TDTCUR%TDATE%YEAR-2000),TDTCUR%TDATE%MONTH,TDTCUR%TDATE%DAY,'_', &
+ INT(TDTLMA%TIME/3600.),INT(FLOAT(MOD(INT(TDTLMA%TIME),3600))/60.), &
+ MOD(INT(TDTLMA%TIME),60), '_', INT(XDTLMA)
+ TDTLMA%TIME = MOD(TDTLMA%TIME + XDTLMA,86400.)
+ CLMA_FILE = CEXP//"_SIMLMA_"//YNAME//".dat"
+!
+ CALL OPEN_ll (FILE=CLMA_FILE, ACTION="WRITE", FORM="FORMATTED", STATUS="NEW", &
+ UNIT=ILMA_UNIT, POSITION="APPEND", IOSTAT=ILMA_IOSTAT )
+ IF ( IPROC .EQ. 0 ) THEN
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '----------------------------------------'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '*** FLASH COORD. FROM LMA SIMULATOR ****'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 1 : flash number --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 2 : time (s) --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 3 : type --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 4 : coordinate along X (km)--'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 5 : coordinate along Y (km)--'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 6 : coordinate along Z (km)--'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 7 : cld drop. mixing ratio --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 8 : rain mixing ratio --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 9 : ice cryst mixing ratio --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 10 : snow mixing ratio --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 11 : graupel mixing ratio --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 12 : rain charge neut --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 13 : ice cryst. charge neut --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 14 : snow charge neut --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 15 : graupel charge neut --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 16 : positive ions neut --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 17 : negative ions neut --'
+ WRITE (UNIT=ILMA_UNIT,FMT='(A)') '----------------------------------------'
+ END IF
+ CALL CLOSE_ll (CLMA_FILE)
+ CALL MPI_BCAST (ILMA_UNIT,1, MPI_INTEGER, 0, MPI_COMM_WORLD, IERR)
+ END IF
END IF
!
!
-!* 9.2 lightning flashes
+!-------------------------------------------------------------------------------
+!
+!* 10. LIGHTNING FLASHES AND NOX PRODUCTION
+! ------------------------------------
!
! the lightning scheme is now called at each time step
! but only if there's electric charge in the domain
!
+ZQTOT(:,:,:) = XECHARGE * (PSVT(:,:,:,NSV_ELECBEG) - PSVT(:,:,:,NSV_ELECEND))
+DO JSV = NSV_ELECBEG+1, NSV_ELECEND-1
+ ZQTOT(:,:,:) = ZQTOT(:,:,:) + PSVT(:,:,:,JSV)
+END DO
+!
IF ((.NOT. LOCG) .AND. LELEC_FIELD .AND. MAX_ll(ABS(ZQTOT),IINFO_ll)>0.) THEN
IF (PRESENT(PSEA)) THEN
IF (LFLASH_GEOM) THEN
- CALL FLASH_GEOM_ELEC_n (KTCOUNT, KRR, PTSTEP, PRHODJ, PRHODREF, &
- PRT, PCIT, &
- PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
- PRS, PTHT, PPABST, &
- XEFIELDU, XEFIELDV, XEFIELDW, &
- PZZ, PTOWN, PSEA )
- ELSE
- PRINT *,' LIGHTNING_ELEC NOT YET DEVELOPPED'
- STOP
+ CALL FLASH_GEOM_ELEC_n (KTCOUNT, KMI, KRR, PTSTEP, OEXIT, &
+ PRHODJ, PRHODREF, &
+ PRT, PCIT, &
+ PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
+ PRS, PTHT, PPABST, &
+ XEFIELDU, XEFIELDV, XEFIELDW, &
+ PZZ, PSVS(:,:,:,NSV_LNOXBEG), PTOWN, PSEA)
END IF
ELSE
IF (LFLASH_GEOM) THEN
- CALL FLASH_GEOM_ELEC_n (KTCOUNT, KRR, PTSTEP, PRHODJ, PRHODREF, &
- PRT, PCIT, &
- PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
- PRS, PTHT, PPABST, &
- XEFIELDU, XEFIELDV, XEFIELDW, PZZ )
- ELSE
- PRINT *,' LIGHTNING_ELEC NOT YET DEVELOPPED'
- STOP
+ CALL FLASH_GEOM_ELEC_n (KTCOUNT, KMI, KRR, PTSTEP, OEXIT, &
+ PRHODJ, PRHODREF, &
+ PRT, PCIT, &
+ PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND), &
+ PRS, PTHT, PPABST, &
+ XEFIELDU, XEFIELDV, XEFIELDW, &
+ PZZ, PSVS(:,:,:,NSV_LNOXBEG) )
END IF
ENDIF
!
@@ -969,17 +1028,6 @@ IF ((.NOT. LOCG) .AND. LELEC_FIELD .AND. MAX_ll(ABS(ZQTOT),IINFO_ll)>0.) THEN
!
END IF
!
-!-------------------------------------------------------------------------------
-!
-!* 10. NOX PRODUCTION BY LIGHTNING FLASHES
-! -----------------------------------
-!
-IF (LLNOX_EXPLICIT) THEN
- PRINT *,' LINOX PRODUCTION FROM THE EXPLICIT LIGHTNING SCHEME NOT YET DEVELOPPED'
-!
- STOP
-END IF
-!
!
!-------------------------------------------------------------------------------
!
@@ -998,9 +1046,10 @@ END IF
!
! Close Ascii Files if KTCOUNT = NSTOP
-IF (KTCOUNT == NSTOP) THEN
+IF (OEXIT) THEN
IF (.NOT. LFLASH_GEOM) CLOSE (UNIT=NLU_light_diag)
IF (.NOT. LFLASH_GEOM .AND. LSAVE_COORD) CLOSE (UNIT=NLU_light_coord)
+ IF (LLMA) CLOSE (UNIT=ILMA_UNIT)
ENDIF
!
!
diff --git a/src/MNH/spawn_field2.f90 b/src/MNH/spawn_field2.f90
index 79d0874e4..37e47a7ce 100644
--- a/src/MNH/spawn_field2.f90
+++ b/src/MNH/spawn_field2.f90
@@ -141,6 +141,7 @@ END MODULE MODI_SPAWN_FIELD2
!! for 2D west african monsoon
!! Modification 07/13 (Bosseur & Filippi) Adds Forefire
!! Modification 2014 (M.Faivre)
+!! Modification 01/15 (C. Barthe) add LNOx
!! Modification 25/02/2015 (M.Moge) correction of the parallelization attempted by M.Faivre
!-------------------------------------------------------------------------------
!
@@ -843,6 +844,12 @@ IF (PRESENT(HSONFILE)) THEN
YCOMMENT,IRESP)
IF(IRESP==0) PSVT(KIB2:KIE2,KJB2:KJE2,:,JSV)=ZWORK3D(KIB1:KIE1,KJB1:KJE1,:)
END DO
+ DO JSV = NSV_LNOXBEG,NSV_LNOXEND ! LNOx Scalar Variables
+ YRECFM='LINOX'
+ CALL FMREAD(HSONFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ IF(IRESP==0) PSVT(KIB2:KIE2,KJB2:KJE2,:,JSV)=ZWORK3D(KIB1:KIE1,KJB1:KJE1,:)
+ END DO
DO JSV = NSV_PPBEG,NSV_PPEND ! Passive scalar variables
WRITE(YRECFM,'(A3,I3.3)')'SVT',JSV
CALL FMREAD(HSONFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH, &
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index 1eab31e1f..e72f6ab13 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -6,6 +6,7 @@
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$
+! masdev4_7 BUG1 2007/06/20 16:58:20
!-----------------------------------------------------------------
! #########################
MODULE MODI_WRITE_LFIFM_n
@@ -157,6 +158,7 @@ END MODULE MODI_WRITE_LFIFM_n
!! P. Tulet Nov 2014 accumulated moles of aqueous species that fall at the surface
!! M.Faivre 2014
!! C.Lac Dec.2014 writing past wind fields for centred advection
+!! J.-P. Pinty Jan 2015 add LNOx and flash map diagnostics
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!!
!-------------------------------------------------------------------------------
@@ -200,7 +202,7 @@ USE MODD_CH_PH_n
USE MODD_CH_M9_n
USE MODD_RAIN_C2R2_DESCR, ONLY: C2R2NAMES
USE MODD_ICE_C1R3_DESCR, ONLY: C1R3NAMES
-USE MODD_ELEC_DESCR, ONLY: CELECNAMES
+USE MODD_ELEC_DESCR, ONLY: CELECNAMES, LLNOX_EXPLICIT
USE MODD_LG, ONLY: CLGNAMES
USE MODD_NSV
USE MODD_AIRCRAFT_BALLOON
@@ -217,6 +219,7 @@ USE MODD_CONDSAMP
USE MODD_CH_AEROSOL
USE MODD_PAST_FIELD_n
USE MODD_ADV_n, ONLY: CUVW_ADV_SCHEME,XRTKEMS
+USE MODD_ELEC_FLASH
!
USE MODE_FMWRIT
USE MODE_ll
@@ -484,9 +487,6 @@ IGRID=3
ILENCH=LEN(YCOMMENT)
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,XYHAT,IGRID,ILENCH,YCOMMENT,IRESP)
!
-!print*,'XXHAT=',XXHAT
-!print*,'XYHAT=',XYHAT
-!print*,'XZHAT=',XZHAT
YRECFM='ZHAT'
YDIR='--'
YCOMMENT='height level without orography (METERS)'
@@ -647,7 +647,6 @@ YDIR='XY'
! CALL EXTRAPOL('N',XUT)
! CALL EXTRAPOL('S',XUT)
CALL MPPDB_CHECK3D(XUT,"write_lfifmn before FMWRIT::XUT",PRECISION)
-!
YRECFM='UT'
YCOMMENT='X_Y_Z_U component of wind (m/s)'
IGRID=2
@@ -721,91 +720,91 @@ END IF
IF (MEAN_COUNT /= 0) THEN
!
- YRECFM='UMMEAN'
+ YRECFM='UMME'
YCOMMENT='X_Y_Z_U component of mean wind (m/s)'
IGRID=2
ILENCH=LEN(YCOMMENT)
ZWORK3D = XUM_MEAN/MEAN_COUNT
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='U2MEAN'
+ YRECFM='U2ME'
YCOMMENT='X_Y_Z_U component of mean wind (m/s)'
IGRID=2
ILENCH=LEN(YCOMMENT)
ZWORK3D = XU2_MEAN/MEAN_COUNT-XUM_MEAN**2/MEAN_COUNT**2
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='VMMEAN'
+ YRECFM='VMME'
YCOMMENT='X_Y_Z_V component of mean wind (m/s)'
IGRID=3
ILENCH=LEN(YCOMMENT)
ZWORK3D = XVM_MEAN/MEAN_COUNT
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='V2MEAN'
+ YRECFM='V2ME'
YCOMMENT='X_Y_Z_V component of mean wind (m/s)'
IGRID=3
ILENCH=LEN(YCOMMENT)
ZWORK3D = XV2_MEAN/MEAN_COUNT-XVM_MEAN**2/MEAN_COUNT**2
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='WMMEAN'
+ YRECFM='WMME'
YCOMMENT='X_Y_Z_vertical mean wind (m/s)'
IGRID=4
ILENCH=LEN(YCOMMENT)
ZWORK3D = XWM_MEAN/MEAN_COUNT
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='W2MEAN'
+ YRECFM='W2ME'
YCOMMENT='X_Y_Z_vertical mean wind (m/s)'
IGRID=4
ILENCH=LEN(YCOMMENT)
ZWORK3D = XW2_MEAN/MEAN_COUNT-XWM_MEAN**2/MEAN_COUNT**2
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='THMMEAN'
+ YRECFM='THMME'
YCOMMENT='X_Y_Z_mean potential temperature (K)'
IGRID=1
ILENCH=LEN(YCOMMENT)
ZWORK3D = XTHM_MEAN/MEAN_COUNT
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='TH2MEAN'
+ YRECFM='TH2ME'
YCOMMENT='X_Y_Z_mean potential temperature (K)'
IGRID=1
ILENCH=LEN(YCOMMENT)
ZWORK3D = XTH2_MEAN/MEAN_COUNT-XTHM_MEAN**2/MEAN_COUNT**2
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='TEMPMMEAN'
+ YRECFM='TEMPMME'
YCOMMENT='X_Y_Z_mean potential temperature (K)'
IGRID=1
ILENCH=LEN(YCOMMENT)
ZWORK3D= XTEMPM_MEAN/MEAN_COUNT
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='TEMP2MEAN'
+ YRECFM='TEMP2ME'
YCOMMENT='X_Y_Z_mean potential temperature (K)'
IGRID=1
ILENCH=LEN(YCOMMENT)
ZWORK3D = XTEMP2_MEAN/MEAN_COUNT-XTEMPM_MEAN**2/MEAN_COUNT**2
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='PABSMMEAN'
+ YRECFM='PABSMME'
YCOMMENT='X_Y_Z_mean ABSolute Pressure (Pa)'
IGRID=1
ILENCH=LEN(YCOMMENT)
ZWORK3D= XPABSM_MEAN/MEAN_COUNT
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='PABS2MEAN'
+ YRECFM='PABS2ME'
YCOMMENT='X_Y_Z_mean ABSolute Pressure (Pa)'
IGRID=1
ILENCH=LEN(YCOMMENT)
ZWORK3D = XPABS2_MEAN/MEAN_COUNT-XPABSM_MEAN**2/MEAN_COUNT**2
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D,IGRID,ILENCH,YCOMMENT,IRESP)
!
- YRECFM='TKEMMEAN'
+ YRECFM='TKEMME'
YCOMMENT='X_Y_Z_mean ABSolute Pressure (Pa)'
IGRID=1
ILENCH=LEN(YCOMMENT)
@@ -822,11 +821,13 @@ IF (CTURB /= 'NONE') THEN
ILENCH=LEN(YCOMMENT)
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,XTKET,IGRID,ILENCH,YCOMMENT,IRESP)
!
+ IF (CPROGRAM == 'MESONH') THEN
YRECFM='TKEMS'
YCOMMENT='X_Y_Z_Turbulent Kinetic Energy adv source (M**2/S**3)'
IGRID=1
ILENCH=LEN(YCOMMENT)
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,XRTKEMS,IGRID,ILENCH,YCOMMENT,IRESP)
+ END IF
END IF
!
!
@@ -1013,6 +1014,63 @@ IF (NSV >=1) THEN
ZWORK3D(:,:,:) = XIND_RATE(:,:,:) * 1.E12
CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D(:,:,:),IGRID,ILENCH, &
YCOMMENT,IRESP)
+ !
+ ZWORK2D(:,:) = 0.
+ YRECFM='TRIG_IC'
+ YCOMMENT='X_Y_Z_FLASH_MAP_TRIG_IC (no unit)'
+ ILENCH=LEN(YCOMMENT)
+ ZWORK2D(:,:) = FLOAT(NMAP_TRIG_IC(:,:))
+ CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK2D(:,:),IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ !
+ ZWORK2D(:,:) = 0.
+ YRECFM='IMPACT_CG'
+ YCOMMENT='X_Y_Z_FLASH_MAP_IMPACT_CG (no unit)'
+ ILENCH=LEN(YCOMMENT)
+ ZWORK2D(:,:) = FLOAT(NMAP_IMPACT_CG(:,:))
+ CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK2D(:,:),IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ !
+ ZWORK2D(:,:) = 0.
+ YRECFM='AREA_CG'
+ YCOMMENT='X_Y_Z_FLASH_MAP_2DAREA_CG (no unit)'
+ ILENCH=LEN(YCOMMENT)
+ ZWORK2D(:,:) = FLOAT(NMAP_2DAREA_CG(:,:))
+ CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK2D(:,:),IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ !
+ ZWORK2D(:,:) = 0.
+ YRECFM='AREA_IC'
+ YCOMMENT='X_Y_Z_FLASH_MAP_2DAREA_IC (no unit)'
+ ILENCH=LEN(YCOMMENT)
+ ZWORK2D(:,:) = FLOAT(NMAP_2DAREA_IC(:,:))
+ CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK2D(:,:),IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ !
+ ZWORK3D(:,:,:) = 0.
+ YRECFM='FLASH_3DCG'
+ YCOMMENT='X_Y_Z_FLASH_MAP_3DCG (no unit)'
+ ILENCH=LEN(YCOMMENT)
+ ZWORK3D(:,:,:) = FLOAT(NMAP_3DCG(:,:,:))
+ CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D(:,:,:),IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ !
+ ZWORK3D(:,:,:) = 0.
+ YRECFM='FLASH_3DIC'
+ YCOMMENT='X_Y_Z_FLASH_MAP_3DIC (no unit)'
+ ILENCH=LEN(YCOMMENT)
+ ZWORK3D(:,:,:) = FLOAT(NMAP_3DIC(:,:,:))
+ CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK3D(:,:,:),IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ !
+ IF (LLNOX_EXPLICIT) THEN
+ YRECFM='LINOX'
+ WRITE(YCOMMENT,'(A6,A3,I3.3,A10)')'X_Y_Z_','SVT',JSV,' (mol/mol)'
+ ILENCH=LEN(YCOMMENT)
+ CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,XSVT(:,:,:,NSV_LNOXEND),IGRID,ILENCH, &
+ YCOMMENT,IRESP)
+ JSA=JSA+1
+ END IF
END IF
! lagrangian variables
DO JSV = NSV_LGBEG,NSV_LGEND
--
GitLab