Philippe 14/11/2016: minor: removed unused IO arguments
[MNH-git_open_source-lfs.git] / src / MNH / resolved_elecn.f90
1
2 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
3 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
4 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
5 !MNH_LIC for details. version 1.
6 !     ###########################
7       MODULE MODI_RESOLVED_ELEC_n
8 !     ###########################
9 !
10 INTERFACE
11       SUBROUTINE RESOLVED_ELEC_n (HCLOUD, HSCONV, HMF_CLOUD,                              &
12                                   KRR, KSPLITR, KMI, KTCOUNT, OEXIT,                      &
13                                   HLBCX, HLBCY, HRAD, HTURBDIM,                           &
14                                   OSUBG_COND, OSIGMAS,PSIGQSAT, HSUBG_AUCV,               &
15                                   PTSTEP, PZZ, PRHODJ, PRHODREF, PEXNREF,                 &
16                                   PPABST, PTHT, PTHS, PWT,                                & 
17                                   PRT, PRS, PSVT, PSVS, PCIT,                             & 
18                                   PSIGS, PSRCS, PCLDFR, PMFCONV, PCF_MF, PRC_MF,          &
19                                   PRI_MF, OSEDIC, OWARM,                                  &
20                                   PINPRC, PINPRR, PINPRR3D, PEVAP3D,                      &
21                                   PINPRS, PINPRG, PINPRH,                                 &
22                                   PSEA, PTOWN                                             )   
23 !
24 CHARACTER(LEN=4),         INTENT(IN)   :: HCLOUD   ! kind of cloud
25 CHARACTER(LEN=4),         INTENT(IN)   :: HSCONV   ! Shallow convection scheme
26 CHARACTER(LEN=4),         INTENT(IN)   :: HMF_CLOUD! Type of statistical cloud
27 INTEGER,                  INTENT(IN)   :: KRR      ! Number of moist variables
28 INTEGER,                  INTENT(IN)   :: KSPLITR  ! Number of small time step
29                                                    ! integrations for  rain sedimendation
30 INTEGER,                  INTENT(IN)   :: KMI      ! Model index
31 INTEGER,                  INTENT(IN)   :: KTCOUNT  ! Temporal loop counter
32 LOGICAL,                  INTENT(IN)   :: OEXIT    ! switch for the end of the temporal loop
33 CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY   ! X and Y-direc. LBC type
34 CHARACTER*4,              INTENT(IN)   :: HRAD     ! Radiation scheme name
35 CHARACTER*4,              INTENT(IN)   :: HTURBDIM ! Dimensionality of the
36                                                    ! turbulence scheme
37 LOGICAL,                  INTENT(IN)   :: OSUBG_COND ! Switch for Subgrid Cond.
38 LOGICAL,                  INTENT(IN)   :: OSIGMAS  ! Switch for Sigma_s:
39                                                    ! use values computed in CONDENSATION
40                                                    ! or that from turbulence scheme
41 REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
42 CHARACTER(LEN=4),         INTENT(IN)   :: HSUBG_AUCV
43                                                    ! Kind of Subgrid autoconversion method
44 REAL,                     INTENT(IN)   :: PTSTEP ! Double Time step
45                                                  ! (single if cold start)
46 !
47 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PZZ     ! Height (z)
48 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PRHODJ  !Dry density * Jacobian
49 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PRHODREF! Reference dry air density
50 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PEXNREF ! Reference Exner function
51 !
52 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PPABST  ! abs. pressure at time t
53 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PTHT    ! Theta at time t
54 REAL, DIMENSION(:,:,:,:), INTENT(INOUT):: PRT     ! Moist variables at time t
55 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PSIGS   ! Sigma_s at time t
56 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PMFCONV ! convective mass flux
57 !
58 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS  ! Theta source
59 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRS   ! Moist  variable sources
60 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT  ! Scalar variable at time t
61 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVS  ! Scalar variable sources
62 !
63 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS ! Second-order flux
64                                                  ! s'rc'/2Sigma_s2 at time t+1
65                                                  ! multiplied by Lambda_3
66 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PCLDFR! Cloud fraction
67 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PCIT  ! Pristine ice number
68                                                  ! concentration at time t
69 !
70 LOGICAL,                  INTENT(IN) :: OSEDIC! Switch to activate the
71                                               ! cloud droplet sedimentation
72 LOGICAL,                  INTENT(IN) :: OWARM ! Control of the rain formation
73                                               !  by slow warm microphysical
74                                               !         processes
75 !
76 REAL, DIMENSION(:,:,:), INTENT(IN) :: PCF_MF! Convective Mass Flux Cloud fraction 
77 REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC_MF! Convective Mass Flux liquid mixing ratio
78 REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI_MF! Convective Mass Flux solid mixing ratio
79 !
80 REAL, DIMENSION(:,:),   INTENT(INOUT) :: PINPRC   ! Cloud instant precip
81 REAL, DIMENSION(:,:),   INTENT(INOUT) :: PINPRR   ! Rain instant precip
82 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PINPRR3D ! sed flux of precip
83 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PEVAP3D  ! evap profile
84 REAL, DIMENSION(:,:),   INTENT(INOUT) :: PINPRS   ! Snow instant precip
85 REAL, DIMENSION(:,:),   INTENT(INOUT) :: PINPRG   ! Graupel instant precip
86 REAL, DIMENSION(:,:),   INTENT(INOUT) :: PINPRH   ! Hail instant precip
87 !
88 REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA   ! Land Sea mask
89 REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PTOWN  ! Town fraction
90 !
91 REAL, DIMENSION(:,:,:), INTENT(IN) :: PWT  ! vertical velocity at time t-dt
92 !
93 END SUBROUTINE RESOLVED_ELEC_n
94 END INTERFACE
95 END MODULE MODI_RESOLVED_ELEC_n
96 !
97 !     #####################################################################################
98       SUBROUTINE RESOLVED_ELEC_n (HCLOUD, HSCONV, HMF_CLOUD,                              &
99                                   KRR, KSPLITR, KMI, KTCOUNT, OEXIT,                      &
100                                   HLBCX, HLBCY, HRAD, HTURBDIM,                           &
101                                   OSUBG_COND, OSIGMAS,PSIGQSAT, HSUBG_AUCV,               &
102                                   PTSTEP, PZZ, PRHODJ, PRHODREF, PEXNREF,                 &
103                                   PPABST, PTHT, PTHS, PWT,                                & 
104                                   PRT, PRS, PSVT, PSVS, PCIT,                             & 
105                                   PSIGS, PSRCS, PCLDFR, PMFCONV, PCF_MF, PRC_MF,          &
106                                   PRI_MF, OSEDIC, OWARM,                                  &
107                                   PINPRC, PINPRR, PINPRR3D, PEVAP3D,                      &
108                                   PINPRS, PINPRG, PINPRH,                                 &
109                                   PSEA, PTOWN                                             )   
110 !     #####################################################################################
111 !
112 !!    PURPOSE
113 !!    -------
114 !!      The purpose of this routine is to compute the resolved clouds and 
115 !!    precipitation, the associated cloud electrification, and the charge 
116 !!    neutralization associated to lightning flashes
117 !!
118 !!
119 !!    METHOD
120 !!    ------
121 !!      The main action of this routine is to call the routines computing the
122 !!    microphysical and electrical sources. Before that:
123 !!        - it computes the real absolute pressure,
124 !!        - negative values of the current guess of all mixing ratio are removed.
125 !!          This is done by a global filling algorithm based on a multiplicative
126 !!          method (Rood, 1987), in order to conserved the total mass in the
127 !!          simulation domain.
128 !!        - Sources are transformed in physical tendencies, by removing the
129 !!          multiplicative term Rhod*J.
130 !!        - External points values are filled owing to the use of cyclic
131 !!          l.b.c., in order to performe computations on the full domain.
132 !!      After calling to microphysical and electrical routines, the physical 
133 !!    tendencies are switched back to prognostic variables.
134 !!
135 !!
136 !!    EXTERNAL
137 !!    --------
138 !!      Subroutine FMLOOK: to recover the logical unit number linked to a FMfile
139 !!      Subroutine SLOW_TERMS: Computes the explicit microphysical sources
140 !!      Subroutine FAST_TERMS: Performs the saturation adjustment for l
141 !!      Subroutine RAIN_ICE  : Computes the explicit microphysical sources for i
142 !!      Subroutine ICE_ADJUST: Performs the saturation adjustment for i+l
143 !!      MIN_ll,SUM3D_ll : distributed functions equivalent to MIN and SUM
144 !!
145 !!
146 !!    IMPLICIT ARGUMENTS
147 !!    ------------------
148 !!      Module MODD_PARAMETERS : contains declarations of parameter variables
149 !!         JPHEXT       : Horizontal external points number
150 !!         JPVEXT       : Vertical external points number
151 !!      Module MODD_CST
152 !!          XP00               ! Reference pressure
153 !!          XRD                ! Gaz  constant for dry air
154 !!          XCPD               ! Cpd (dry air)
155 !!
156 !!    REFERENCE
157 !!    ---------
158 !!
159 !!    AUTHOR
160 !!    ------
161 !!      C. Barthe       * LACy *
162 !!
163 !!    MODIFICATIONS
164 !!    -------------
165 !!      Original    06/11/09
166 !!      Modifications: 
167 !!      M. Chong      26/01/10  Add Small ions parameters
168 !!      M. Chong      31/07/14  Add explicit LiNOx
169 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
170 !!
171 !-------------------------------------------------------------------------------
172 !
173 !*       0.    DECLARATIONS
174 !              ------------
175 !
176 USE MODE_ll
177 USE MODE_FM
178 USE MODE_ELEC_ll
179 !
180 USE MODD_METRICS_n, ONLY : XDXX, XDYY, XDZX, XDZY, XDZZ 
181 USE MODD_FIELD_n, ONLY : XRSVS
182 USE MODD_CONF, ONLY : L1D, L2D, CEXP
183 USE MODD_CST
184 USE MODD_PARAMETERS, ONLY : JPVEXT
185 USE MODD_ELEC_DESCR
186 USE MODD_ELEC_n          
187 USE MODD_BUDGET
188 USE MODD_NSV
189 USE MODD_CH_MNHC_n,    ONLY: LUSECHEM,LCH_CONV_LINOX
190 USE MODD_DYN_n, ONLY: NSTOP, XTSTEP
191 USE MODD_ARGSLIST_ll, ONLY : LIST_ll
192
193 USE MODD_TIME_n
194 USE MODD_LMA_SIMULATOR
195 USE MODD_PRINT_ELEC
196 !
197 USE MODI_IO_ll
198 USE MODD_VAR_ll, ONLY : NMNH_COMM_WORLD
199 USE MODI_RAIN_ICE_ELEC
200 USE MODI_ICE_ADJUST_ELEC
201 USE MODI_TO_ELEC_FIELD_n
202 USE MODI_FLASH_GEOM_ELEC_n
203 USE MODI_SHUMAN
204 USE MODI_BUDGET
205 USE MODI_ION_ATTACH_ELEC
206 USE MODD_ARGSLIST_ll, ONLY : LIST_ll
207 !
208 USE MODI_ION_DRIFT
209 USE MODI_SERIES_CLOUD_ELEC
210 !
211 IMPLICIT NONE
212 !
213 !*       0.1   Declarations of dummy arguments :
214 !
215 CHARACTER(LEN=4),         INTENT(IN)   :: HCLOUD   ! kind of cloud
216                                                    ! paramerization
217 CHARACTER(LEN=4),         INTENT(IN)   :: HSCONV   ! Shallow convection scheme
218 CHARACTER(LEN=4),         INTENT(IN)   :: HMF_CLOUD! Type of statistical cloud
219 INTEGER,                  INTENT(IN)   :: KRR      ! Number of moist variables
220 INTEGER,                  INTENT(IN)   :: KSPLITR  ! Number of small time step
221                                        ! integrations for  rain sedimendation
222                                        ! integrations for  ice  sedimendation
223 INTEGER,                  INTENT(IN)   :: KMI      ! Model index
224 INTEGER,                  INTENT(IN)   :: KTCOUNT  ! Temporal loop counter
225 LOGICAL,                  INTENT(IN)   :: OEXIT    ! switch for the end of the temporal loop
226 CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX,HLBCY   ! X and Y-direc. LBC type
227 CHARACTER*4,              INTENT(IN)   :: HRAD     ! Radiation scheme name
228 CHARACTER*4,              INTENT(IN)   :: HTURBDIM ! Dimensionality of the
229                                                    ! turbulence scheme
230 LOGICAL,                  INTENT(IN)   :: OSUBG_COND ! Switch for Subgrid Cond.
231 LOGICAL,                  INTENT(IN)   :: OSIGMAS  ! Switch for Sigma_s:
232                                         ! use values computed in CONDENSATION
233                                         ! or that from turbulence scheme
234 REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
235 CHARACTER(LEN=4),         INTENT(IN)   :: HSUBG_AUCV
236                                         ! Kind of Subgrid autoconversion method
237 REAL,                     INTENT(IN)   :: PTSTEP   ! Double Time step
238                                                    ! (single if cold start)
239 !
240 !
241 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PZZ     ! Height (z)
242 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PRHODJ  !Dry density * Jacobian
243 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PRHODREF! Reference dry air density
244 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PEXNREF ! Reference Exner function
245 !
246 !
247 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PPABST  ! abs. pressure at time t
248 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PTHT    ! Theta at time t
249 REAL, DIMENSION(:,:,:,:), INTENT(INOUT):: PRT     ! Moist variables at time t
250 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PSIGS   ! Sigma_s at time t
251 REAL, DIMENSION(:,:,:),   INTENT(IN)   :: PMFCONV ! convective mass flux
252 !
253 !
254 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS  ! Theta source
255 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRS   ! Moist  variable sources
256 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT  ! Scalar variable at time t
257 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVS  ! Scalar variable sources
258 !
259 REAL, DIMENSION(:,:,:),   INTENT(OUT) :: PSRCS ! Second-order flux
260                                                ! s'rc'/2Sigma_s2 at time t+1
261                                                ! multiplied by Lambda_3
262 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PCLDFR! Cloud fraction
263 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PCIT  ! Pristine ice number
264                                                ! concentration at time t
265 LOGICAL,                  INTENT(IN) :: OSEDIC! Switch to activate the
266                                               ! cloud droplet sedimentation
267                                               ! for ICE3            
268 LOGICAL,                  INTENT(IN) :: OWARM ! Control of the rain formation
269                                               !  by slow warm microphysical
270                                               !         processes
271 !
272 REAL, DIMENSION(:,:,:), INTENT(IN) :: PCF_MF! Convective Mass Flux Cloud fraction 
273 REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC_MF! Convective Mass Flux liquid mixing ratio
274 REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI_MF! Convective Mass Flux solid mixing ratio
275 !
276 REAL, DIMENSION(:,:), INTENT(INOUT)   :: PINPRC ! Cloud instant precip
277 REAL, DIMENSION(:,:), INTENT(INOUT)   :: PINPRR ! Rain instant precip
278 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PINPRR3D ! sed flux of precip
279 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PEVAP3D  ! evap profile
280 REAL, DIMENSION(:,:), INTENT(INOUT)   :: PINPRS ! Snow instant precip
281 REAL, DIMENSION(:,:), INTENT(INOUT)   :: PINPRG ! Graupel instant precip
282 REAL, DIMENSION(:,:), INTENT(INOUT)   :: PINPRH ! Hail instant precip
283 !
284 REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PSEA   ! Land Sea mask
285 REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: PTOWN  ! Town fraction
286 !
287 REAL, DIMENSION(:,:,:), INTENT(IN) :: PWT  ! vertical velocity at time t-dt
288 !
289 !*       0.2   Declarations of local variables :
290 !
291 INTEGER :: JRR,JSV       ! Loop index for the moist and scalar variables
292 INTEGER :: IIB           !  Define the physical domain
293 INTEGER :: IIE           !
294 INTEGER :: IJB           !
295 INTEGER :: IJE           !
296 INTEGER :: IKB           !
297 INTEGER :: IKE           !
298 INTEGER :: IKU
299 INTEGER :: IINFO_ll      ! return code of parallel routine
300 INTEGER :: IPROC         ! my proc number
301 INTEGER :: IERR          ! error status
302 !
303 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZT,   &
304                                                        ZEXN, &
305                                                        ZLV,  &
306                                                        ZLS,  &
307                                                        ZCPH
308 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZCOR
309                                     ! for the correction of negative rv
310 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZZZ
311                                     ! model layer height
312 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZQTOT
313                                     ! total charge source term
314 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZION_NUMBER  !nearly Nb
315                          !  of elementary charge in hydrometeor charge
316 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZADD         ! ratio (0
317                          !  or 1) of ZION_NUMBER to add to positive
318                          !              or negative ion number
319 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZIONTOT
320 !
321 REAL :: ZMASSTOT         ! total mass  for one water category
322                          ! including the negative values
323 REAL :: ZMASSPOS         ! total mass  for one water category
324                          ! after removing the negative values
325 REAL :: ZRATIO           ! ZMASSTOT / ZMASSCOR
326 !
327 INTEGER, DIMENSION(3) :: IMINLOC, IMAXLOC
328 !
329 LOGICAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: GMASSCOR ! mask for
330                                                                ! mass correction
331 LOGICAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: GATTACH  ! mask for
332                                      !ion recombination and attachment
333 !
334 TYPE(LIST_ll), POINTER :: TZFIELDS_ll   ! list of fields to exchange
335
336 INTEGER, DIMENSION(3) :: IM_LOC
337 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZDRIFT
338 INTEGER :: IPROCMIN, IK
339 INTEGER :: IXOR, IYOR  ! origin of the extended subdomain
340 CHARACTER (LEN=32) :: YASCFILE
341 !
342 REAL               :: ZTEMP_DIST
343 CHARACTER (LEN=18) :: YNAME
344 LOGICAL            :: GLMA_FILE
345 !
346 NULLIFY(TZFIELDS_ll)
347 !
348 !------------------------------------------------------------------------------
349 !
350 !*       1.     PRELIMINARY COMPUTATIONS
351 !               ------------------------
352 !
353 CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
354 IKB = 1 + JPVEXT
355 IKE = SIZE(PZZ,3) - JPVEXT
356 IKU = SIZE(PZZ,3)
357 !
358 !
359 !------------------------------------------------------------------------------
360 !
361 !*       2.     MICROPHYSICS AND CLOUD ELECTRIFICATION
362 !               --------------------------------------
363 !
364 !*       2.1    Transformation into physical tendencies
365 !
366 ! X-Component per m3.s into X-Component per kg.s
367 PTHS(:,:,:) = PTHS(:,:,:) / PRHODJ(:,:,:)
368 DO JRR = 1, KRR
369   PRS(:,:,:,JRR) = PRS(:,:,:,JRR) / PRHODJ(:,:,:)
370 END DO
371 !
372 DO JSV = NSV_ELECBEG, NSV_ELECEND
373   PSVS(:,:,:,JSV) = PSVS(:,:,:,JSV) / PRHODJ(:,:,:)
374 ENDDO
375 !
376 !  complete the lateral boundaries to avoid possible problems
377 !
378 PTHS(IIB-1,:,:) = PTHS(IIB,:,:)
379 PTHS(IIE+1,:,:) = PTHS(IIE,:,:)
380 PTHS(:,IJB-1,:) = PTHS(:,IJB,:)
381 PTHS(:,IJE+1,:) = PTHS(:,IJE,:)
382 !
383 PRS(IIB-1,:,:,1) = PRS(IIB,:,:,1)
384 PRS(IIE+1,:,:,1) = PRS(IIE,:,:,1)
385 PRS(:,IJB-1,:,1) = PRS(:,IJB,:,1)
386 PRS(:,IJE+1,:,1) = PRS(:,IJE,:,1)
387 !
388 PRS(IIB-1,:,:,2:) = 0.0
389 PRS(IIE+1,:,:,2:) = 0.0
390 PRS(:,IJB-1,:,2:) = 0.0
391 PRS(:,IJE+1,:,2:) = 0.0
392 !
393 ! positive ion source
394 PSVS(IIB-1,:,:,NSV_ELECBEG) = PSVS(IIB,:,:,NSV_ELECBEG)
395 PSVS(IIE+1,:,:,NSV_ELECBEG) = PSVS(IIE,:,:,NSV_ELECBEG) 
396 PSVS(:,IJB-1,:,NSV_ELECBEG) = PSVS(:,IJB,:,NSV_ELECBEG)
397 PSVS(:,IJE+1,:,NSV_ELECBEG) = PSVS(:,IJE,:,NSV_ELECBEG)
398 ! source of hydrometeor charge
399 PSVS(IIB-1,:,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0    
400 PSVS(IIE+1,:,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0   
401 PSVS(:,IJB-1,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
402 PSVS(:,IJE+1,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
403 ! negative ion source
404 PSVS(IIB-1,:,:,NSV_ELECEND) = PSVS(IIB,:,:,NSV_ELECEND) 
405 PSVS(IIE+1,:,:,NSV_ELECEND) = PSVS(IIE,:,:,NSV_ELECEND)
406 PSVS(:,IJB-1,:,NSV_ELECEND) = PSVS(:,IJB,:,NSV_ELECEND)
407 PSVS(:,IJE+1,:,NSV_ELECEND) = PSVS(:,IJE,:,NSV_ELECEND)
408 !
409 !  complete the physical boundaries to avoid some computations
410 !
411 IF(LWEST_ll()  .AND. HLBCX(1) /= 'CYCL')  PRT(IIB-1,:,:,2:) = 0.0
412 IF(LEAST_ll()  .AND. HLBCX(2) /= 'CYCL')  PRT(IIE+1,:,:,2:) = 0.0
413 IF(LSOUTH_ll() .AND. HLBCY(1) /= 'CYCL')  PRT(:,IJB-1,:,2:) = 0.0
414 IF(LNORTH_ll() .AND. HLBCY(2) /= 'CYCL')  PRT(:,IJE+1,:,2:) = 0.0
415 !
416 IF(LWEST_ll()  .AND. HLBCX(1) /= 'CYCL')  &
417                         PSVT(IIB-1,:,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
418 IF(LEAST_ll()  .AND. HLBCX(2) /= 'CYCL')  &
419                         PSVT(IIE+1,:,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
420 IF(LSOUTH_ll() .AND. HLBCY(1) /= 'CYCL')  &
421                         PSVT(:,IJB-1,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
422 IF(LNORTH_ll() .AND. HLBCY(2) /= 'CYCL')  &
423                         PSVT(:,IJE+1,:,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
424 !
425 !  complete the vertical boundaries
426 !
427 PTHS(:,:,IKB-1) = PTHS(:,:,IKB)
428 PTHS(:,:,IKE+1) = PTHS(:,:,IKE)
429 !
430 PRS(:,:,IKB-1,1) = PRS(:,:,IKB,1)
431 PRS(:,:,IKE+1,1) = PRS(:,:,IKE,1)
432 PRS(:,:,IKB-1,2:) = 0.0
433 PRS(:,:,IKE+1,2:) = 0.0
434 !
435 PRT(:,:,IKB-1,1) = PRT(:,:,IKB,1)
436 PRT(:,:,IKE+1,1) = PRT(:,:,IKE,1)
437 PRT(:,:,IKB-1,2:) = 0.0
438 PRT(:,:,IKE+1,2:) = 0.0
439 !
440 PSVS(:,:,IKB-1,NSV_ELECBEG) = PSVS(:,:,IKB,NSV_ELECBEG)    ! Positive ion
441 PSVT(:,:,IKB-1,NSV_ELECBEG) = PSVT(:,:,IKB,NSV_ELECBEG)
442 PSVS(:,:,IKB-1,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0          ! Hydrometeor charge
443 PSVS(:,:,IKE+1,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
444 PSVT(:,:,IKB-1,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
445 PSVT(:,:,IKE+1,NSV_ELECBEG+1:NSV_ELECEND-1) = 0.0
446 PSVS(:,:,IKB-1,NSV_ELECEND) = PSVS(:,:,IKB,NSV_ELECEND)    ! Negative ion
447 PSVT(:,:,IKB-1,NSV_ELECEND) = PSVT(:,:,IKB,NSV_ELECEND)
448 !
449 ! personal comment:  tranfering these variables to the
450 !                    microphysical routines would save
451 !                    computing time
452 !
453 ZEXN(:,:,:) = (PPABST(:,:,:) / XP00)**(XRD / XCPD)
454 ZT(:,:,:)   = PTHT(:,:,:) * ZEXN(:,:,:)
455 ZLV(:,:,:)  = XLVTT + (XCPV - XCL) * (ZT(:,:,:) - XTT)
456 ZLS(:,:,:)  = XLSTT + (XCPV - XCI) * (ZT(:,:,:) - XTT)
457 ZCPH(:,:,:) = XCPD + XCPV * PTSTEP * PRS(:,:,:,1)
458 !
459 !
460 !------------------------------------------------------------------------------
461 !
462 !*       3.     REMOVE NEGATIVE VALUES
463 !               ----------------------
464 !
465 !*       3.1    Non local correction for precipitating species (Rood 87)
466 !
467 DO JRR = 3, KRR
468   SELECT CASE (JRR)
469     CASE(3,5,6,7) ! rain, snow, graupel and hail
470 !
471       IF (MIN_ll(PRS(:,:,:,JRR), IINFO_ll) < 0.0) THEN
472         GMASSCOR = PRS(:,:,:,JRR) < 0.
473 !
474 ! compute the total water mass computation
475         ZMASSTOT = MAX( 0. , SUM3D_ll( PRS(:,:,:,JRR), IINFO_ll ) )
476 !
477 ! remove the negative values
478         PRS(:,:,:,JRR) = MAX( 0., PRS(:,:,:,JRR) )
479 !
480 ! compute the new total mass
481         ZMASSPOS = MAX(XMNH_TINY,SUM3D_ll( PRS(:,:,:,JRR), IINFO_ll ) )
482 !
483 ! correct again in such a way to conserve the total mass
484         ZRATIO = ZMASSTOT / ZMASSPOS
485         PRS(:,:,:,JRR) = PRS(:,:,:,JRR) * ZRATIO
486 !
487 ! No electric charge without hydrometeors
488         WHERE( GMASSCOR )
489           PSVS(:,:,:,NSV_ELECBEG-1+JRR) = 0.
490         ENDWHERE
491       END IF
492   END SELECT
493 END DO
494 !
495 !
496 !*       3.2    Adjustement for liquid and solid cloud
497 !
498 WHERE (PRS(:,:,:,4) < 0.)                          ! ice particles
499   PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,4)
500   PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,4) * ZLS(:,:,:) /  &
501                 ZCPH(:,:,:) / ZEXN(:,:,:)
502   PRS(:,:,:,4) = 0.
503 !
504   ZION_NUMBER(:,:,:) = ABS(PSVS(:,:,:,NSV_ELECBEG+3)) / XECHARGE
505   ZADD(:,:,:) = 0.5 + SIGN(0.5, PSVS(:,:,:,NSV_ELECBEG+3))
506   PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) +  &
507                             ZADD(:,:,:) * ZION_NUMBER(:,:,:)
508   PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) +  &
509                             (1. - ZADD(:,:,:)) * ZION_NUMBER(:,:,:)
510   PSVS(:,:,:,NSV_ELECBEG+3) = 0.0
511 END WHERE
512 !
513 ! cloud
514 WHERE (PRS(:,:,:,2) < 0.)
515   PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2)
516   PTHS(:,:,:) = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) /  &
517        ZCPH(:,:,:) / ZEXN(:,:,:)
518   PRS(:,:,:,2) = 0.
519 !
520   ZION_NUMBER(:,:,:) = ABS(PSVS(:,:,:,NSV_ELECBEG+1)) / XECHARGE
521   ZADD(:,:,:) = 0.5 + SIGN(0.5, PSVS(:,:,:,NSV_ELECBEG+1))
522   PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) +  &
523                             ZADD(:,:,:) * ZION_NUMBER(:,:,:)
524   PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) +  &
525                             (1. - ZADD(:,:,:)) * ZION_NUMBER(:,:,:)
526   PSVS(:,:,:,NSV_ELECBEG+1) = 0.0
527 END WHERE
528 !
529 ! if rc or ri are positive, we can correct negative rv
530 ! cloud
531 WHERE ((PRS(:,:,:,1) < 0.) .AND. (PRS(:,:,:,2) > 0.) )
532   PRS(:,:,:,1) = PRS(:,:,:,1) + PRS(:,:,:,2)
533   PTHS(:,:,:)  = PTHS(:,:,:) - PRS(:,:,:,2) * ZLV(:,:,:) /  &
534                  ZCPH(:,:,:) / ZEXN(:,:,:)
535   PRS(:,:,:,2) = 0.
536 !
537   ZION_NUMBER(:,:,:) = ABS(PSVS(:,:,:,NSV_ELECBEG+1)) / XECHARGE
538   ZADD(:,:,:) = 0.5 + SIGN(0.5, PSVS(:,:,:,NSV_ELECBEG+1))
539   PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) +  &
540                             ZADD(:,:,:) * ZION_NUMBER(:,:,:)
541   PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) +  &
542                             (1.-ZADD(:,:,:)) * ZION_NUMBER(:,:,:)
543   PSVS(:,:,:,NSV_ELECBEG+1) = 0.0
544 END WHERE
545 !
546 ! ice
547 IF(KRR > 3) THEN
548   WHERE ((PRS(:,:,:,1) < 0.).AND.(PRS(:,:,:,4) > 0.))
549     ZCOR(:,:,:)  = MIN(-PRS(:,:,:,1),PRS(:,:,:,4))
550     PRS(:,:,:,1) = PRS(:,:,:,1) + ZCOR(:,:,:)
551     PTHS(:,:,:)  = PTHS(:,:,:) - ZCOR(:,:,:) * ZLS(:,:,:) /  &
552                    ZCPH(:,:,:) / ZEXN(:,:,:)
553     PRS(:,:,:,4) = PRS(:,:,:,4) -ZCOR(:,:,:)
554 !
555   ZION_NUMBER(:,:,:) = ABS(PSVS(:,:,:,NSV_ELECBEG+3)) / XECHARGE
556   ZADD(:,:,:) = 0.5 + SIGN(0.5, PSVS(:,:,:,NSV_ELECBEG+3))
557   PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) +  &
558                             ZADD(:,:,:) * ZION_NUMBER(:,:,:)
559   PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) +  &
560                             (1. - ZADD(:,:,:)) * ZION_NUMBER(:,:,:)
561   PSVS(:,:,:,NSV_ELECBEG+3) = 0.0
562   END WHERE
563 END IF
564 !
565 !
566 !*       3.3     cascade the electric charges in absence of hydrometeor
567 !
568 DO JRR = KRR, 5, -1
569   WHERE(PRS(:,:,:,JRR) < XRTMIN_ELEC(JRR))
570     PSVS(:,:,:,NSV_ELECBEG-2+JRR) = PSVS(:,:,:,NSV_ELECBEG-2+JRR) + &
571                                     PSVS(:,:,:,NSV_ELECBEG-1+JRR)
572     PSVS(:,:,:,NSV_ELECBEG-1+JRR) = 0.0
573   END WHERE
574 END DO
575 JRR = 3
576 WHERE(PRS(:,:,:,JRR) < XRTMIN_ELEC(JRR))
577   PSVS(:,:,:,NSV_ELECBEG-2+JRR) = PSVS(:,:,:,NSV_ELECBEG-2+JRR) + &
578                                   PSVS(:,:,:,NSV_ELECBEG-1+JRR)
579   PSVS(:,:,:,NSV_ELECBEG-1+JRR) = 0.0
580 END WHERE
581 DO JRR = 4, 2, -2
582   WHERE(PRS(:,:,:,JRR) < XRTMIN_ELEC(JRR))
583 !
584     ZION_NUMBER(:,:,:) = ABS(PSVS(:,:,:,NSV_ELECBEG-1+JRR)) / XECHARGE
585     ZADD(:,:,:) = 0.5 + SIGN(0.5, PSVS(:,:,:,NSV_ELECBEG-1+JRR))
586     PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) +  &
587                               ZADD(:,:,:) * ZION_NUMBER(:,:,:)
588     PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) +  &
589                               (1. - ZADD(:,:,:)) * ZION_NUMBER(:,:,:)
590     PSVS(:,:,:,NSV_ELECBEG-1+JRR) = 0.0
591   END WHERE
592 END DO
593 !
594 !
595 !*       3.4     store the budget terms
596 !
597 IF (LBUDGET_RV) CALL BUDGET (PRS(:,:,:,1) * PRHODJ(:,:,:), 6,'NEGA_BU_RRV')
598 IF (LBUDGET_RC) CALL BUDGET (PRS(:,:,:,2) * PRHODJ(:,:,:), 7,'NEGA_BU_RRC')
599 IF (LBUDGET_RR) CALL BUDGET (PRS(:,:,:,3) * PRHODJ(:,:,:), 8,'NEGA_BU_RRR')
600 IF (LBUDGET_RI) CALL BUDGET (PRS(:,:,:,4) * PRHODJ(:,:,:) ,9,'NEGA_BU_RRI')
601 IF (LBUDGET_RS) CALL BUDGET (PRS(:,:,:,5) * PRHODJ(:,:,:),10,'NEGA_BU_RRS')
602 IF (LBUDGET_RG) CALL BUDGET (PRS(:,:,:,6) * PRHODJ(:,:,:),11,'NEGA_BU_RRG')
603 IF (LBUDGET_RH) CALL BUDGET (PRS(:,:,:,7) * PRHODJ(:,:,:),12,'NEGA_BU_RRH')
604 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)  * PRHODJ(:,:,:), 4,'NEGA_BU_RTH')
605 !
606 IF (LBUDGET_SV) THEN
607   DO JSV = NSV_ELECBEG, NSV_ELECEND
608     CALL BUDGET (PSVS(:,:,:,JSV) * PRHODJ(:,:,:), 12+JSV, 'NEGA_BU_RSV')
609   END DO
610 END IF
611 !
612 !
613 !------------------------------------------------------------------------------
614 !
615 !*       4.     ION SOURCE FROM DRIFT MOTION AND COSMIC RAYS
616 !               ---------------------------------------------
617 !
618 !*       4.1    Compute the electric field at mass points
619 !
620 PSVT(:,:,:,NSV_ELECBEG) = XECHARGE*PSVT(:,:,:,NSV_ELECBEG)    ! 1/kg --> C/kg
621 PSVT(:,:,:,NSV_ELECEND) =-XECHARGE*PSVT(:,:,:,NSV_ELECEND)
622 !
623 CALL TO_ELEC_FIELD_n (PRT, PSVT(:,:,:,NSV_ELECBEG:NSV_ELECEND), PRHODJ, &
624                       KTCOUNT, KRR,                                     &
625                       XEFIELDU, XEFIELDV, XEFIELDW                      )
626 !
627 PSVT(:,:,:,NSV_ELECBEG) = PSVT(:,:,:,NSV_ELECBEG)/XECHARGE    ! back to 1/kg 
628 PSVT(:,:,:,NSV_ELECEND) =-PSVT(:,:,:,NSV_ELECEND)/XECHARGE
629 !
630 !
631 !*       4.2    Compute source term from -/+(Div (N.mu E)) at mass points, 
632 !               N positive or negative ion number per kg of air (= PSVT)
633 !               This is a contribution of drift motion to Source PSVS for ions
634 !               in 1/(kg.s)
635 !
636 CALL MYPROC_ELEC_ll (IPROC)
637 !
638 !     Hereafter, ZCPH and ZCOR are used temporarily to store the drift sources 
639 !          of the positive and negative ions, respectively
640 !
641 CALL ION_DRIFT(ZCPH, ZCOR, PSVT, HLBCX, HLBCY)
642
643
644 PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) + ZCPH(:,:,:)
645 PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) + ZCOR(:,:,:)
646 !
647 !*       4.3    Add Cosmic Ray source
648 !
649 PSVS(:,:,:,NSV_ELECBEG) = PSVS(:,:,:,NSV_ELECBEG) +           &
650                               XIONSOURCEFW(:,:,:) / PRHODREF(:,:,:)
651 PSVS(:,:,:,NSV_ELECEND) = PSVS(:,:,:,NSV_ELECEND) +           &
652                               XIONSOURCEFW(:,:,:) / PRHODREF(:,:,:)
653 !
654 !-------------------------------------------------------------------------------
655 !
656 SELECT CASE (HCLOUD)
657 !
658   CASE ('ICE3')
659 !
660 !*       5.     MIXED-PHASE MICROPHYSICAL SCHEME (WITH 3 ICE SPECIES)
661 !               -----------------------------------------------------
662 !
663 !*       5.1    Compute the explicit microphysical sources and
664 !*              the explicit charging rates
665 !
666     CALL RAIN_ICE_ELEC (OSEDIC, HSUBG_AUCV, OWARM,                            &
667                         KSPLITR, PTSTEP, KMI, KRR,                            &
668                         PZZ, PRHODJ, PRHODREF, PEXNREF, PPABST, PCIT, PCLDFR, &
669                         PTHT, PRT(:,:,:,1), PRT(:,:,:,2), PRT(:,:,:,3),       &
670                         PRT(:,:,:,4), PRT(:,:,:,5), PRT(:,:,:,6),             &
671                         PTHS, PRS(:,:,:,1), PRS(:,:,:,2), PRS(:,:,:,3),       &
672                         PRS(:,:,:,4), PRS(:,:,:,5), PRS(:,:,:,6),             &
673                         PINPRC, PINPRR, PINPRR3D, PEVAP3D,                    &
674                         PINPRS, PINPRG, PSIGS,                                &
675                         PSVT(:,:,:,NSV_ELECBEG),   PSVT(:,:,:,NSV_ELECBEG+1), &
676                         PSVT(:,:,:,NSV_ELECBEG+2), PSVT(:,:,:,NSV_ELECBEG+3), &
677                         PSVT(:,:,:,NSV_ELECBEG+4), PSVT(:,:,:,NSV_ELECBEG+5), &
678                         PSVT(:,:,:,NSV_ELECEND),                              &
679                         PSVS(:,:,:,NSV_ELECBEG),   PSVS(:,:,:,NSV_ELECBEG+1), &
680                         PSVS(:,:,:,NSV_ELECBEG+2), PSVS(:,:,:,NSV_ELECBEG+3), &
681                         PSVS(:,:,:,NSV_ELECBEG+4), PSVS(:,:,:,NSV_ELECBEG+5), &
682                         PSVS(:,:,:,NSV_ELECEND),                              &
683                         PSEA, PTOWN                                           )
684
685 !
686 !*       5.2    Perform the saturation adjustment over cloud ice and cloud water
687 !
688     ZZZ = MZF(1,IKU,1, PZZ )
689     CALL ICE_ADJUST_ELEC (KRR, KMI, HRAD, HTURBDIM,                           &
690                           HSCONV, HMF_CLOUD,                                  &
691                           OSUBG_COND, OSIGMAS, PTSTEP,PSIGQSAT,               &
692                           PRHODJ, PEXNREF, PSIGS, PPABST, ZZZ,                &
693                           PMFCONV, PCF_MF, PRC_MF, PRI_MF,                    &
694                           PRVT=PRT(:,:,:,1), PRCT=PRT(:,:,:,2),               &
695                           PRVS=PRS(:,:,:,1), PRCS=PRS(:,:,:,2),               &
696                           PTHS=PTHS, PSRCS=PSRCS, PCLDFR=PCLDFR,              &
697                           PRRT=PRT(:,:,:,3), PRRS=PRS(:,:,:,3),               &
698                           PRIT=PRT(:,:,:,4), PRIS=PRS(:,:,:,4),               &
699                           PRST=PRT(:,:,:,5), PRSS=PRS(:,:,:,5),               &
700                           PRGT=PRT(:,:,:,6), PRGS=PRS(:,:,:,6),               &
701                           PQPIT=PSVT(:,:,:,NSV_ELECBEG),  & !..PI.. Positive
702                           PQPIS=PSVS(:,:,:,NSV_ELECBEG),  & !  Ion Mixing Ratio
703                           PQCT=PSVT(:,:,:,NSV_ELECBEG+1), &
704                           PQCS=PSVS(:,:,:,NSV_ELECBEG+1), &
705                           PQRT=PSVT(:,:,:,NSV_ELECBEG+2), &
706                           PQRS=PSVS(:,:,:,NSV_ELECBEG+2), &
707                           PQIT=PSVT(:,:,:,NSV_ELECBEG+3), &
708                           PQIS=PSVS(:,:,:,NSV_ELECBEG+3), &
709                           PQST=PSVT(:,:,:,NSV_ELECBEG+4), &
710                           PQSS=PSVS(:,:,:,NSV_ELECBEG+4), &
711                           PQGT=PSVT(:,:,:,NSV_ELECBEG+5), &
712                           PQGS=PSVS(:,:,:,NSV_ELECBEG+5), &
713                           PQNIT=PSVT(:,:,:,NSV_ELECEND),  & !..NI.. Negative
714                           PQNIS=PSVS(:,:,:,NSV_ELECEND))    !  Ion Mixing Ratio
715 !
716 !
717 !-------------------------------------------------------------------------------
718 !
719 !*       6.     MIXED-PHASE MICROPHYSICAL SCHEME (WITH 4 ICE SPECIES)
720 !               -----------------------------------------------------
721 !
722 !*       6.1    Compute the explicit microphysical sources and
723 !*              the explicit charging rates
724 !
725   CASE ('ICE4')
726 !
727     CALL RAIN_ICE_ELEC (OSEDIC, HSUBG_AUCV, OWARM,                            &
728                         KSPLITR, PTSTEP, KMI, KRR,                            &
729                         PZZ, PRHODJ, PRHODREF, PEXNREF, PPABST, PCIT, PCLDFR, &
730                         PTHT, PRT(:,:,:,1), PRT(:,:,:,2), PRT(:,:,:,3),       &
731                         PRT(:,:,:,4), PRT(:,:,:,5), PRT(:,:,:,6),             &
732                         PTHS, PRS(:,:,:,1), PRS(:,:,:,2), PRS(:,:,:,3),       &
733                         PRS(:,:,:,4), PRS(:,:,:,5), PRS(:,:,:,6),             &
734                         PINPRC, PINPRR, PINPRR3D, PEVAP3D,                    &
735                         PINPRS, PINPRG, PSIGS,                   &
736                         PSVT(:,:,:,NSV_ELECBEG),   PSVT(:,:,:,NSV_ELECBEG+1), &
737                         PSVT(:,:,:,NSV_ELECBEG+2), PSVT(:,:,:,NSV_ELECBEG+3), &
738                         PSVT(:,:,:,NSV_ELECBEG+4), PSVT(:,:,:,NSV_ELECBEG+5), &
739                         PSVT(:,:,:,NSV_ELECEND),          &
740                         PSVS(:,:,:,NSV_ELECBEG),   PSVS(:,:,:,NSV_ELECBEG+1), &
741                         PSVS(:,:,:,NSV_ELECBEG+2), PSVS(:,:,:,NSV_ELECBEG+3), &
742                         PSVS(:,:,:,NSV_ELECBEG+4), PSVS(:,:,:,NSV_ELECBEG+5), &
743                         PSVS(:,:,:,NSV_ELECEND),          &
744                         PSEA, PTOWN,                   &
745                         PRT(:,:,:,7), PRS(:,:,:,7), PINPRH,                   &
746                         PSVT(:,:,:,NSV_ELECBEG+6), PSVS(:,:,:,NSV_ELECBEG+6)  )
747 ! Index NSV_ELECBEG: Positive ion , NSV_ELECEND: Negative ion
748 !
749 !
750 !*       6.2    Perform the saturation adjustment over cloud ice and cloud water
751 !
752     ZZZ = MZF(1,IKU,1, PZZ )
753     CALL ICE_ADJUST_ELEC (KRR, KMI, HRAD,                                     &
754                           HTURBDIM, HSCONV, HMF_CLOUD,                        &
755                           OSUBG_COND, OSIGMAS, PTSTEP,PSIGQSAT,               &
756                           PRHODJ, PEXNREF, PSIGS, PPABST, ZZZ,                &
757                           PMFCONV, PCF_MF, PRC_MF, PRI_MF,                    & 
758                           PRVT=PRT(:,:,:,1), PRCT=PRT(:,:,:,2),               &
759                           PRVS=PRS(:,:,:,1), PRCS=PRS(:,:,:,2),               &
760                           PTHS=PTHS, PSRCS=PSRCS, PCLDFR=PCLDFR,              &
761                           PRRT=PRT(:,:,:,3), PRRS=PRS(:,:,:,3),               &
762                           PRIT=PRT(:,:,:,4), PRIS=PRS(:,:,:,4),               &
763                           PRST=PRT(:,:,:,5), PRSS=PRS(:,:,:,5),               &
764                           PRGT=PRT(:,:,:,6), PRGS=PRS(:,:,:,6),               &
765                           PQPIT=PSVT(:,:,:,NSV_ELECBEG), &  !..PI.. Positive
766                           PQPIS=PSVS(:,:,:,NSV_ELECBEG),  & !  Ion Mixing Ratio
767                           PQCT=PSVT(:,:,:,NSV_ELECBEG+1), &
768                           PQCS=PSVS(:,:,:,NSV_ELECBEG+1), &
769                           PQRT=PSVT(:,:,:,NSV_ELECBEG+2), &
770                           PQRS=PSVS(:,:,:,NSV_ELECBEG+2), &
771                           PQIT=PSVT(:,:,:,NSV_ELECBEG+3), &
772                           PQIS=PSVS(:,:,:,NSV_ELECBEG+3), &
773                           PQST=PSVT(:,:,:,NSV_ELECBEG+4), &
774                           PQSS=PSVS(:,:,:,NSV_ELECBEG+4), &
775                           PQGT=PSVT(:,:,:,NSV_ELECBEG+5), &
776                           PQGS=PSVS(:,:,:,NSV_ELECBEG+5), &
777                           PQNIT=PSVT(:,:,:,NSV_ELECEND),  & !..NI.. Negative
778                           PQNIS=PSVS(:,:,:,NSV_ELECEND),  & !  Ion Mixing Ratio
779                           PRHT=PRT(:,:,:,7), PRHS=PRS(:,:,:,7),               &
780                           PQHT=PSVT(:,:,:,NSV_ELECBEG+6), &
781                           PQHS=PSVS(:,:,:,NSV_ELECBEG+6) )
782 !
783 END SELECT
784 !
785 IF(KTCOUNT .EQ. 1 .AND. IPROC .EQ. 0) PRINT *,'KSPLITR=', KSPLITR
786 !
787 !-------------------------------------------------------------------------------
788 !
789 !*      7.      SWITCH BACK TO THE PROGNOSTIC VARIABLES
790 !               ---------------------------------------
791 !
792 ! Convert source into component per m3 of air and sec., i.e. volumetric source
793 !
794 PTHS(:,:,:) = PTHS(:,:,:) * PRHODJ(:,:,:)
795 !
796 DO JRR = 1,KRR
797   PRS(:,:,:,JRR)  = PRS(:,:,:,JRR) * PRHODJ(:,:,:)
798 END DO
799 !
800 DO JSV = NSV_ELECBEG, NSV_ELECEND
801   PSVS(:,:,:,JSV) = PSVS(:,:,:,JSV) * PRHODJ(:,:,:)
802 ENDDO
803 !
804 ! Note that the LiNOx Conc. (in mol/mol) is PSVS (:,::,NSV_LNOXBEG)
805 ! but there is no need to *PRHODJ(:,:,:) as it is done implicitly
806 ! during unit conversion in flash_geom.
807 !
808 PSVS(:,:,:,NSV_ELECBEG) = MAX(0., PSVS(:,:,:,NSV_ELECBEG))
809 PSVS(:,:,:,NSV_ELECEND) = MAX(0., PSVS(:,:,:,NSV_ELECEND))
810 !
811 !-------------------------------------------------------------------------------
812 !
813 !*      8.      ION RECOMBINATION AND ATTACHMENT
814 !               --------------------------------
815 !
816 GATTACH(:,:,:) = .FALSE.
817 GATTACH(IIB:IIE, IJB:IJE, IKB:IKE) = .TRUE.
818 !
819 IF (PRESENT(PSEA)) THEN
820   CALL ION_ATTACH_ELEC(KTCOUNT, KRR, PTSTEP, PRHODREF,                   &
821                        PRHODJ, PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND),      & 
822                        PRS, PTHT, PCIT, PPABST, XEFIELDU,                &
823                        XEFIELDV, XEFIELDW, GATTACH, PTOWN, PSEA          )
824 ELSE
825   CALL ION_ATTACH_ELEC(KTCOUNT, KRR, PTSTEP, PRHODREF,                   &
826                        PRHODJ, PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND),      &
827                        PRS, PTHT, PCIT, PPABST, XEFIELDU,                &
828                        XEFIELDV, XEFIELDW, GATTACH                       )
829 ENDIF
830 !
831 !-------------------------------------------------------------------------------
832 !
833 !*      9.      OPEN THE OUTPUT ASCII FILES
834 !               ---------------------------
835 !
836 IF (KTCOUNT .EQ. 1) THEN
837   IF (LFLASH_GEOM) THEN
838     YASCFILE = CEXP//"_fgeom_diag.asc"
839     CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW",                 &
840                   FORM="FORMATTED", POSITION="APPEND",                         &
841                   UNIT=NLU_fgeom_diag, IOSTAT= NIOSTAT_fgeom_diag              )
842     IF ( IPROC .EQ. 0) THEN
843       WRITE (NLU_fgeom_diag, FMT='(A)') '--------------------------------------------------------'
844       WRITE (NLU_fgeom_diag, FMT='(A)') '*FLASH CHARACTERISTICS FROM FLASH_GEOM_ELEC*'
845       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 1 : total flash number          --'
846       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 2 : time (s)                    --'
847       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 3 : cell number                 --'
848       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 4 : flash number/cell/time step --'
849       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 5 : flash type 1=IC, 2=CGN, 3=CGP '
850       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 6 : number of segments          --'
851       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 7 : trig electric field (kV/m)  --'
852       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 8 : x coord. trig. point        --'
853       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 9 : y coord. trig. point        --'
854       WRITE (NLU_fgeom_diag, FMT='(A)') '--         --> x,y in km if lcartesian=t, --'
855       WRITE (NLU_fgeom_diag, FMT='(A)') '--                    deg otherwise       --'
856       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 10 : z coord. trig. point (km)  --'
857       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 11: neutr. positive charge (C)  --' 
858       WRITE (NLU_fgeom_diag, FMT='(A)') '-- Column 12: neutr. negative charge (C)  --'
859       WRITE (NLU_fgeom_diag, FMT='(A)') '--------------------------------------------'
860     END IF
861 !  
862     CALL CLOSE_ll (YASCFILE)
863     CALL MPI_BCAST (NLU_fgeom_diag,1, MPI_INTEGER, 0, NMNH_COMM_WORLD, IERR)
864 !
865     IF (LSAVE_COORD) THEN
866       YASCFILE = CEXP//"_fgeom_coord.asc"
867       CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW",&
868                     FORM="FORMATTED", POSITION="APPEND",                        &
869                     UNIT=NLU_fgeom_coord, IOSTAT= NIOSTAT_fgeom_coord     )
870       IF ( IPROC .EQ. 0) THEN
871         WRITE (NLU_fgeom_coord,FMT='(A)') '------------------------------------------'
872         WRITE (NLU_fgeom_coord,FMT='(A)') '*****FLASH COORD. FROM FLASH_GEOM_ELEC****'
873         WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 1 : flash number             --'
874         WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 2 : time (s)                 --'
875         WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 3 : type                     --'
876         WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 4 : coordinate along X (km)  --'
877         WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 5 : coordinate along Y (km)  --'
878         WRITE (NLU_fgeom_coord,FMT='(A)') '-- Column 6 : coordinate along Z (km)  --'
879         WRITE (NLU_fgeom_coord,FMT='(A)') '------------------------------------------'
880       END IF
881 !
882       CALL CLOSE_ll (YASCFILE)
883       CALL MPI_BCAST (NLU_fgeom_coord,1, MPI_INTEGER, 0, NMNH_COMM_WORLD, IERR)
884     END IF
885   END IF
886 !
887   IF (LSERIES_ELEC) THEN
888     YASCFILE = CEXP//"_series_cloud_elec.asc"                              
889     CALL OPEN_ll (FILE=YASCFILE, ACTION="WRITE", STATUS="NEW", &
890                   FORM="FORMATTED", POSITION="APPEND",         &
891                   UNIT=NLU_series_cloud_elec, IOSTAT= NIOSTAT_series_cloud_elec)
892     IF ( IPROC .EQ. 0) THEN
893       WRITE (NLU_series_cloud_elec, FMT='(A)') '----------------------------------------------------'
894       WRITE (NLU_series_cloud_elec, FMT='(A)') '********* RESULTS FROM of LSERIES_ELEC *************'
895       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 1 : Time (s)                            --'
896       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 2 : Cloud top height / Z > 20 dBZ (m)   --'
897       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 3 : Cloud top height / m.r. > 1.e-4 (m) --'
898       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 4 : Maximum radar reflectivity (dBZ)    --'
899       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 5 : Maximum vertical velocity (m/s)     --'
900       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 6 : Updraft volume for W > 5 m/s        --'
901       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 7 : Updraft volume for W > 10 m/s       --'
902       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 8 : Cloud water mass (kg)               --'
903       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 9 : Rain water mass (kg)                --'
904       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 10 : Ice crystal mass (kg)              --'
905       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 11 : Snow mass (kg)                     --'
906       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 12 : Graupel mass (kg)                  --'
907       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 13 : Precipitation ice mass (kg)        --'
908       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 14 : Ice mass flux product (kg2 m2/s2)  --'
909       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 15 : Precip. ice mass flux (kg m/s)     --'
910       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 16 : Non-precip. ice mass flux (kg m/s) --'
911       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 17 : Ice water path (kg/m2)             --'
912       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 18 : Cloud volume (m3)                  --'
913       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 19 : Maximum rain inst. precip. (mm/H)  --'
914       WRITE (NLU_series_cloud_elec, FMT='(A)') '-- Column 20 : Rain instant. precip. (mm/H)       --'
915       WRITE (NLU_series_cloud_elec, FMT='(A)') '----------------------------------------------------'
916     END IF
917 !
918     CALL CLOSE_ll (YASCFILE)
919     CALL MPI_BCAST (NLU_series_cloud_elec,1, MPI_INTEGER, 0, NMNH_COMM_WORLD, IERR)
920   END IF
921 END IF
922 !
923 IF (LFLASH_GEOM .AND. LLMA) THEN
924 !
925 ! test to see if a new LMA file should be created
926 !
927   GLMA_FILE = .FALSE.
928 !
929   IF (TDTCUR%TIME >= TDTLMA%TIME-PTSTEP .AND. CLMA_FILE(1:5) /= "BEGIN") THEN
930     LWRITE_LMA  = .TRUE.
931   END IF
932   IF (TDTCUR%TIME >= TDTLMA%TIME) THEN
933     TDTLMA%TIME = TDTLMA%TIME + XDTLMA
934     GLMA_FILE   = .TRUE.
935     LWRITE_LMA  = .FALSE.
936   END IF
937 !
938   IF (GLMA_FILE) THEN
939     IF(CLMA_FILE(1:5) /= "BEGIN") THEN ! close preceeding file when existing
940       CALL CLOSE_ll (CLMA_FILE)
941     ENDIF
942 !
943     TDTLMA%TIME = TDTLMA%TIME - XDTLMA
944     WRITE (YNAME,FMT='(3I2.2,A1,3I2.2,A1,I4.4)')                               &
945           ABS(TDTCUR%TDATE%YEAR-2000),TDTCUR%TDATE%MONTH,TDTCUR%TDATE%DAY,'_', &
946             INT(TDTLMA%TIME/3600.),INT(FLOAT(MOD(INT(TDTLMA%TIME),3600))/60.), &
947                                       MOD(INT(TDTLMA%TIME),60), '_', INT(XDTLMA)
948     TDTLMA%TIME = MOD(TDTLMA%TIME + XDTLMA,86400.)
949     CLMA_FILE = CEXP//"_SIMLMA_"//YNAME//".dat"
950 !
951     CALL OPEN_ll (FILE=CLMA_FILE, ACTION="WRITE",    FORM="FORMATTED", STATUS="NEW", &
952                   UNIT=ILMA_UNIT, POSITION="APPEND", IOSTAT=ILMA_IOSTAT )
953     IF ( IPROC .EQ. 0 ) THEN
954       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '----------------------------------------'
955       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '*** FLASH COORD. FROM LMA SIMULATOR ****'
956       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 1  : flash number           --'
957       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 2  : time (s)               --'
958       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 3  : type                   --'
959       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 4  : coordinate along X (km)--'
960       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 5  : coordinate along Y (km)--'
961       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 6  : coordinate along Z (km)--'
962       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 7  : cld drop. mixing ratio --'
963       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 8  : rain mixing ratio      --'
964       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 9  : ice cryst mixing ratio --'
965       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 10 : snow mixing ratio      --'
966       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 11 : graupel mixing ratio   --'
967       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 12 : rain charge neut       --'
968       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 13 : ice cryst. charge neut --'
969       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 14 : snow charge neut       --'
970       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 15 : graupel charge neut    --'
971       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 16 : positive ions neut     --'
972       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '-- Column 17 : negative ions neut     --'
973       WRITE (UNIT=ILMA_UNIT,FMT='(A)') '----------------------------------------'
974     END IF
975     CALL CLOSE_ll (CLMA_FILE)
976     CALL MPI_BCAST (ILMA_UNIT,1, MPI_INTEGER, 0, NMNH_COMM_WORLD, IERR)
977   END IF
978 END IF
979 !
980 !
981 !-------------------------------------------------------------------------------
982 !
983 !*      10.     LIGHTNING FLASHES AND NOX PRODUCTION
984 !               ------------------------------------
985 !
986 ! the lightning scheme is now called at each time step
987 ! but only if there's electric charge in the domain
988 !
989 ZQTOT(:,:,:) = XECHARGE * (PSVT(:,:,:,NSV_ELECBEG) - PSVT(:,:,:,NSV_ELECEND))
990 DO JSV = NSV_ELECBEG+1, NSV_ELECEND-1
991   ZQTOT(:,:,:) = ZQTOT(:,:,:) + PSVT(:,:,:,JSV)
992 END DO
993 !
994 IF ((.NOT. LOCG) .AND. LELEC_FIELD .AND.  MAX_ll(ABS(ZQTOT),IINFO_ll)>0.) THEN
995   IF (PRESENT(PSEA)) THEN
996     IF (LFLASH_GEOM) THEN
997       CALL FLASH_GEOM_ELEC_n (KTCOUNT, KMI, KRR, PTSTEP, OEXIT,       &
998                               PRHODJ, PRHODREF,                       &
999                               PRT, PCIT,                              &
1000                               PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND),    &
1001                               PRS, PTHT, PPABST,                      & 
1002                               XEFIELDU, XEFIELDV, XEFIELDW,           &
1003                               PZZ, PSVS(:,:,:,NSV_LNOXBEG), PTOWN, PSEA)
1004     END IF 
1005   ELSE
1006     IF (LFLASH_GEOM) THEN
1007       CALL FLASH_GEOM_ELEC_n (KTCOUNT, KMI, KRR, PTSTEP, OEXIT,       &
1008                               PRHODJ, PRHODREF,                       &
1009                               PRT, PCIT,                              &
1010                               PSVS(:,:,:,NSV_ELECBEG:NSV_ELECEND),    &
1011                               PRS, PTHT, PPABST,                      & 
1012                               XEFIELDU, XEFIELDV, XEFIELDW,           &
1013                               PZZ, PSVS(:,:,:,NSV_LNOXBEG)            )
1014     END IF
1015   ENDIF
1016 !
1017   PSVS(:,:,:,NSV_ELECBEG) = MAX(0., PSVS(:,:,:,NSV_ELECBEG))
1018   PSVS(:,:,:,NSV_ELECEND) = MAX(0., PSVS(:,:,:,NSV_ELECEND))
1019 !
1020 END IF
1021 !
1022 !
1023 !-------------------------------------------------------------------------------
1024 !
1025 !*      11.     LOOK FOR FLASH RATE PROXIES
1026 !               ---------------------------
1027 !
1028 IF (LSERIES_ELEC) THEN
1029   CALL SERIES_CLOUD_ELEC (KTCOUNT, PTSTEP,                &
1030                           PZZ, PRHODJ, PRHODREF, PEXNREF, &
1031                           PRT, PRS, PSVT,                 &
1032                           PTHT, PWT, PPABST, PCIT, PINPRR  )
1033 END IF
1034 !
1035 !
1036 !-------------------------------------------------------------------------------
1037 !
1038 !   Close Ascii Files if KTCOUNT = NSTOP
1039
1040 IF (OEXIT) THEN
1041   IF (.NOT. LFLASH_GEOM) CLOSE (UNIT=NLU_light_diag)
1042   IF (.NOT. LFLASH_GEOM .AND. LSAVE_COORD) CLOSE (UNIT=NLU_light_coord)
1043   IF (LLMA) CLOSE (UNIT=ILMA_UNIT)
1044 ENDIF
1045 !
1046 !
1047 !-------------------------------------------------------------------------------
1048 !
1049 END SUBROUTINE RESOLVED_ELEC_n