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