Philippe 14/11/2016: minor: removed unused IO arguments
[MNH-git_open_source-lfs.git] / src / MNH / ice_adjust_elec.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 !-----------------------------------------------------------------
6 !--------------- special set of characters for RCS information
7 !-----------------------------------------------------------------
8 ! $Source$ $Revision$ $Date$
9 !-----------------------------------------------------------------
10 !     ###########################
11       MODULE MODI_ICE_ADJUST_ELEC
12 !     ###########################
13 !
14 INTERFACE
15 !
16       SUBROUTINE ICE_ADJUST_ELEC (KRR, KMI, HRAD, HTURBDIM, HSCONV, HMF_CLOUD,    &
17                                   OSUBG_COND, OSIGMAS, PTSTEP,PSIGQSAT,           &
18                                   PRHODJ, PEXNREF, PSIGS, PPABST, PZZ,            &
19                                   PMFCONV, PCF_MF, PRC_MF, PRI_MF,                &
20                                   PRVT, PRCT, PRVS, PRCS, PTHS, PSRCS, PCLDFR ,   &
21                                   PRRT, PRRS, PRIT, PRIS, PRST, PRSS, PRGT, PRGS, &
22                                   PQPIT, PQPIS, PQCT, PQCS,                       &
23                                   PQRT, PQRS, PQIT, PQIS, PQST, PQSS, PQGT, PQGS, &
24                                   PQNIT, PQNIS, PRHT, PRHS, PQHT, PQHS            )
25 !
26 INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
27 INTEGER,                  INTENT(IN)    :: KMI      ! Model index 
28 CHARACTER*4,              INTENT(IN)    :: HTURBDIM ! Dimensionality of the
29                                                     ! turbulence scheme
30 CHARACTER(LEN=4),         INTENT(IN)    :: HSCONV   ! Shallow convection scheme
31 CHARACTER(LEN=4),         INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
32 CHARACTER*4,              INTENT(IN)    :: HRAD     ! Radiation scheme name
33 LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
34                                                     ! Condensation
35 LOGICAL                                 :: OSIGMAS  ! Switch for Sigma_s: 
36                                                     ! use values computed in CONDENSATION
37                                                     ! or that from turbulence scheme
38 REAL,                     INTENT(IN)   :: PTSTEP    ! Double Time step
39                                                     ! (single if cold start)
40 REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
41 !
42 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODJ  ! Dry density * Jacobian
43 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PEXNREF ! Reference Exner function
44 !
45 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PSIGS   ! Sigma_s at time t
46 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PMFCONV ! convective mass flux
47 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PPABST  ! Absolute Pressure at t        
48 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PZZ     ! height of model layer
49 !
50 REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PCF_MF! Convective Mass Flux Cloud fraction 
51 REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRC_MF! Convective Mass Flux liquid mixing ratio
52 REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRI_MF! Convective Mass Flux solid mixing ratio
53 !
54 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRVT    ! Water vapor m.r. at t
55 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRCT    ! Cloud water m.r. at t
56 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
57 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRCS    ! Cloud water m.r. source
58 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
59 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS   ! Second-order flux
60                                                    ! s'rc'/2Sigma_s2 at time t+1
61                                                    ! multiplied by Lambda_3
62 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PCLDFR  ! Cloud fraction          
63 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRRS ! Rain water m.r. at t+1
64 REAL, DIMENSION(:,:,:),   INTENT(INOUT)::  PRIS ! Cloud ice  m.r. at t+1
65 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRSS ! Aggregate  m.r. at t+1
66 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRGS ! Graupel    m.r. at t+1
67 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRRT ! Rain water m.r. at t
68 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRIT ! Cloud ice  m.r. at t
69 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRST ! Aggregate  m.r. at t
70 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRGT ! Graupel    m.r. at t
71 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) ::  PRHT ! Hail m.r. at t
72 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) ::  PRHS ! Hail m.r. at t+1
73 !
74 REAL, DIMENSION(:,:,:),   INTENT(IN)      ::  PQPIT  ! positive ion m.r. at t
75 REAL, DIMENSION(:,:,:),   INTENT(INOUT)   ::  PQPIS  !           source
76 REAL, DIMENSION(:,:,:),   INTENT(IN)      ::  PQNIT  ! negative ion m.r. at t
77 REAL, DIMENSION(:,:,:),   INTENT(INOUT)   ::  PQNIS  ! source
78 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQCT    ! Cloud water m.r. at t
79 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PQCS    ! Cloud water m.r. source
80 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQRS ! Rain water m.r. at t+1
81 REAL, DIMENSION(:,:,:),   INTENT(INOUT)::  PQIS ! Cloud ice  m.r. at t+1
82 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQSS ! Aggregate  m.r. at t+1
83 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQGS ! Graupel    m.r. at t+1
84 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQRT ! Rain water m.r. at t
85 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQIT ! Cloud ice  m.r. at t
86 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQST ! Aggregate  m.r. at t
87 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQGT ! Graupel    m.r. at t
88 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHT ! Hail  m.r. at t
89 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHS ! Hail  m.r. at t+1
90 !
91 END SUBROUTINE ICE_ADJUST_ELEC
92 END INTERFACE
93 END MODULE MODI_ICE_ADJUST_ELEC
94 !
95 !     ########################################################################
96       SUBROUTINE ICE_ADJUST_ELEC (KRR, KMI, HRAD, HTURBDIM, HSCONV,          &
97                              HMF_CLOUD, OSUBG_COND, OSIGMAS, PTSTEP,PSIGQSAT,&
98                              PRHODJ, PEXNREF, PSIGS, PPABST, PZZ,            &
99                              PMFCONV, PCF_MF, PRC_MF, PRI_MF,                &  
100                              PRVT, PRCT, PRVS, PRCS, PTHS, PSRCS, PCLDFR ,   &
101                              PRRT, PRRS, PRIT, PRIS, PRST, PRSS, PRGT, PRGS, &
102                              PQPIT, PQPIS, PQCT, PQCS,                       &
103                              PQRT, PQRS, PQIT, PQIS, PQST, PQSS, PQGT, PQGS, &
104                              PQNIT, PQNIS, PRHT, PRHS, PQHT, PQHS            )
105 !     ########################################################################
106 !
107 !!****  *ICE_ADJUST_ELEC* -  compute the ajustment of water vapor in mixed-phase 
108 !!                           clouds
109 !!
110 !!    PURPOSE
111 !!    -------
112 !!    The purpose of this routine is to compute the fast microphysical sources
113 !!    through a saturation ajustement procedure in case of mixed-phase clouds.
114 !!
115 !!
116 !!**  METHOD
117 !!    ------
118 !!    Langlois, Tellus, 1973 for the cloudless version.
119 !!    When cloud water is taken into account, refer to book 1 of the
120 !!    documentation.
121 !!
122 !!
123 !!    EXTERNAL
124 !!    --------
125 !!      None
126 !!     
127 !!
128 !!    IMPLICIT ARGUMENTS
129 !!    ------------------
130 !!      Module MODD_CST
131 !!         XP00               ! Reference pressure
132 !!         XMD,XMV            ! Molar mass of dry air and molar mass of vapor
133 !!         XRD,XRV            ! Gaz constant for dry air, gaz constant for vapor
134 !!         XCPD,XCPV          ! Cpd (dry air), Cpv (vapor)
135 !!         XCL                ! Cl (liquid)
136 !!         XCI                ! Ci (ice)
137 !!         XTT                ! Triple point temperature
138 !!         XLVTT              ! Vaporization heat constant
139 !!         XLSTT              ! Sublimation  heat constant
140 !!         XALPW,XBETAW,XGAMW ! Constants for saturation vapor over liquid
141 !!                            !  pressure  function 
142 !!         XALPI,XBETAI,XGAMI ! Constants for saturation vapor over ice
143 !!                            !  pressure  function 
144 !!      Module  MODD_CONF 
145 !!         CCONF
146 !!      Module MODD_BUDGET:
147 !!         NBUMOD 
148 !!         CBUTYPE
149 !!         NBUPROCCTR 
150 !!         LBU_RTH    
151 !!         LBU_RRV  
152 !!         LBU_RRC  
153 !!         LBU_RRI  
154 !!
155 !!
156 !!    REFERENCE
157 !!    ---------
158 !!      Book 1 and Book2 of documentation ( routine ICE_ADJUST )
159 !!      Langlois, Tellus, 1973
160 !!
161 !!    AUTHOR
162 !!    ------
163 !!      J.-P. Pinty    * Laboratoire d'Aerologie*
164 !!   
165 !!
166 !!    MODIFICATIONS
167 !!    -------------
168 !!      Original    2002 
169 !!      C. Barthe   19/11/09   update to version 4.8.1
170 !!      M. Chong    Mar. 2010  Add small ions
171 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
172 !!
173 !-------------------------------------------------------------------------------
174 !
175 !*       0.    DECLARATIONS
176 !              ------------
177 !
178 USE MODD_PARAMETERS
179 USE MODD_CST
180 USE MODD_CONF
181 USE MODD_BUDGET
182 USE MODD_NSV, ONLY : NSV_ELECBEG, NSV_ELECEND
183 USE MODD_ELEC_DESCR, ONLY : XRTMIN_ELEC, XQTMIN, XFC, XFI, XECHARGE
184 USE MODD_RAIN_ICE_DESCR, ONLY : XRTMIN, XBI
185 !
186 USE MODI_CONDENSATION
187 USE MODI_BUDGET
188 USE MODE_FMWRIT
189 USE MODI_GET_HALO
190 !
191 IMPLICIT NONE
192 !
193 !
194 !*       0.1   Declarations of dummy arguments 
195 !
196 INTEGER,                INTENT(IN)    :: KRR      ! Number of moist variables
197 INTEGER,                INTENT(IN)    :: KMI      ! Model index 
198 CHARACTER*4,            INTENT(IN)    :: HTURBDIM ! Dimensionality of the
199                                                   ! turbulence scheme
200 CHARACTER(LEN=4),       INTENT(IN)    :: HSCONV   ! Shallow convection scheme
201 CHARACTER(LEN=4),       INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
202 CHARACTER*4,            INTENT(IN)    :: HRAD     ! Radiation scheme name
203 LOGICAL,                INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
204                                                     ! Condensation
205 LOGICAL                               :: OSIGMAS  ! Switch for Sigma_s: 
206                                                   ! use values computed in CONDENSATION
207                                                   ! or that from turbulence scheme
208 REAL,                   INTENT(IN)    :: PTSTEP   ! Double Time step
209                                                   ! (single if cold start)
210 REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
211 !
212 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ  ! Dry density * Jacobian
213 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PEXNREF ! Reference Exner function
214 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PSIGS   ! Sigma_s at time t
215 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PMFCONV ! convective mass flux
216 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PPABST  ! Absolute Pressure at t        
217 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PZZ     ! height of model layer
218 !
219 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PCF_MF! Convective Mass Flux Cloud fraction 
220 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRC_MF! Convective Mass Flux liquid mixing ratio
221 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRI_MF! Convective Mass Flux solid mixing ratio
222 !
223 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRVT    ! Water vapor m.r. at t
224 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRCT    ! Cloud water m.r. at t
225 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
226 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRCS    ! Cloud water m.r. source
227 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS    ! Theta source
228 REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PSRCS   ! Second-order flux
229                                                  ! s'rc'/2Sigma_s2 at time t+1
230                                                  ! multiplied by Lambda_3
231 REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PCLDFR  ! Cloud fraction          
232 !
233 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRRS ! Rain water m.r. at t+1
234 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRIS ! Cloud ice  m.r. at t+1
235 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRSS ! Aggregate  m.r. at t+1
236 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRGS ! Graupel    m.r. at t+1
237 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRRT ! Rain water m.r. at t
238 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRIT ! Cloud ice  m.r. at t
239 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRST ! Aggregate  m.r. at t
240 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PRGT ! Graupel    m.r. at t
241 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
242 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHS ! Hail m.r. at t+1
243 !
244 REAL, DIMENSION(:,:,:),   INTENT(INOUT)   ::  PQPIT  ! positive ion m.r. at t
245 REAL, DIMENSION(:,:,:),   INTENT(INOUT)   ::  PQPIS  !           source
246 REAL, DIMENSION(:,:,:),   INTENT(INOUT)   ::  PQNIT  ! negative ion m.r. at t
247 REAL, DIMENSION(:,:,:),   INTENT(INOUT)   ::  PQNIS  ! source
248 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PQCT    ! Cloud water m.r. at t
249 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQCS ! Cloud water m.r. source
250 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PQRS ! Rain water m.r. at t+1
251 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQIS ! Cloud ice  m.r. at t+1
252 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PQSS ! Aggregate  m.r. at t+1
253 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PQGS ! Graupel    m.r. at t+1
254 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PQRT ! Rain water m.r. at t
255 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PQIT ! Cloud ice  m.r. at t
256 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PQST ! Aggregate  m.r. at t
257 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PQGT ! Graupel    m.r. at t
258 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHT ! Hail m.r. at t
259 REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHS ! Hail m.r. at t+1
260 !
261 !
262 !*       0.2   Declarations of local variables :
263 !
264 REAL  :: ZEPS  ! Mv/Md
265 REAL  :: ZT00,ZT0   ! Min and max temperature for the mixed phase liquid and solid water
266                     ! for the coeff CND of the barycentric mixing ratio
267 REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) &
268                          :: ZEXNS,&  ! guess of the Exner functon at t+1
269                             ZT,   &  ! guess of the temperature at t+1
270                             ZCPH, &  ! guess of the CPh for the mixing
271                             ZLV,  &  ! guess of the Lv at t+1
272                             ZLS,  &  ! guess of the Ls at t+1
273       ZW1,ZW2,ZW3,ZW4,ZW5,ZW6,ZW7,&  ! Work arrays for intermediate fields
274                             ZCND     ! CND=(T-T00)/(T0-T00) cf sc doc and TAO etal (89)
275 REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) &
276                          :: ZWE1, &
277                             ZWE2
278 REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) &
279                          :: ZION_NUMBER,  &  !nearly Nb of elementary charge
280                                              ! in hydrometeor charge  
281                             ZADD             ! ratio (0 or 1) of ZION_NUMBER 
282                                              ! to add to positive
283                                              ! or negative ion number
284 !
285 INTEGER             :: IRESP      ! Return code of FM routines
286 INTEGER             :: ILENG      ! Length of comment string in LFIFM file
287 INTEGER             :: IGRID      ! C-grid indicator in LFIFM file
288 INTEGER             :: ILENCH     ! Length of comment string in LFIFM file
289 CHARACTER (LEN=100) :: YCOMMENT   ! Comment string in LFIFM file
290 CHARACTER (LEN=16)  :: YRECFM     ! Name of the desired field in LFIFM file
291 !
292 INTEGER             :: IIU,IJU,IKU! dimensions of dummy arrays
293 INTEGER             :: IIB,IJB    ! Horz index values of the first inner mass points
294 INTEGER             :: IIE,IJE    ! Horz index values of the last inner mass points
295 INTEGER             :: IKB        ! K index value of the first inner mass point
296 INTEGER             :: IKE        ! K index value of the last inner mass point
297 INTEGER             :: JITER,ITERMAX ! iterative loop for first order adjustment
298 !
299 LOGICAL             :: LPRETREATMENT, LNEW_ADJUST
300 !
301 !-------------------------------------------------------------------------------
302 !
303 !*       1.     PRELIMINARIES
304 !               -------------
305 !
306 IIU = SIZE(PEXNREF,1)
307 IJU = SIZE(PEXNREF,2)
308 IKU = SIZE(PEXNREF,3)
309 CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
310 IKB = 1 + JPVEXT
311 IKE = IKU - JPVEXT
312 !
313 ZEPS = XMV / XMD
314 !
315 ITERMAX = 1
316 !
317 LPRETREATMENT=.TRUE.     ! FALSE to retreive the previous MASDEV4_1 version
318 LNEW_ADJUST  =.TRUE.     ! FALSE to retreive the previous MASDEV4_1 version
319 ZT0  = XTT               ! Usefull if LPRETREATMENT=T or LNEW_ADJUST=T
320 ZT00 = XTT-40.           ! Usefull if LPRETREATMENT=T or LNEW_ADJUST=T
321 !
322 !-------------------------------------------------------------------------------
323 !
324 !*       2.     COMPUTE QUANTITIES WITH THE GUESS OF THE FUTURE INSTANT
325 !               -------------------------------------------------------
326 !
327 !*       2.1    estimate the pressure at t+1
328 !
329 ZEXNS(:,:,:) = ( PPABST(:,:,:)  / XP00)**(XRD/XCPD)
330 !
331 !    beginning of the iterative loop
332 !
333 DO JITER = 1, ITERMAX
334 !
335 !*       2.2    compute the intermediate temperature at t+1, T*
336 !  
337   ZT(:,:,:) = (PTHS(:,:,:) * PTSTEP) * ZEXNS(:,:,:)
338 !
339 !*       2.3    compute the latent heat of vaporization Lv(T*) at t+1
340 !                   and the latent heat of sublimation  Ls(T*) at t+1
341 !
342   ZLV(:,:,:) = XLVTT + (XCPV - XCL) * (ZT(:,:,:) - XTT)
343   ZLS(:,:,:) = XLSTT + (XCPV - XCI) * (ZT(:,:,:) - XTT)
344 !
345 !*       2.4    compute the specific heat for moist air (Cph) at t+1
346 !
347   IF     ( KRR == 7 ) THEN
348     ZCPH(:,:,:) = XCPD + XCPV *PTSTEP*  PRVS(:,:,:)                             &
349                        + XCL  *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:))              &
350                        + XCI  *PTSTEP* (PRIS(:,:,:) + PRSS(:,:,:) + PRGS(:,:,:) &
351                                                                   + PRHS(:,:,:))
352   ELSE IF( KRR == 6 ) THEN
353     ZCPH(:,:,:) = XCPD + XCPV *PTSTEP*  PRVS(:,:,:)                            &
354                        + XCL  *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:))             &
355                        + XCI  *PTSTEP* (PRIS(:,:,:) + PRSS(:,:,:) + PRGS(:,:,:))
356   ELSE IF( KRR == 5 ) THEN
357     ZCPH(:,:,:) = XCPD + XCPV *PTSTEP*  PRVS(:,:,:)                            &
358                        + XCL  *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:))             &
359                        + XCI  *PTSTEP* (PRIS(:,:,:) + PRSS(:,:,:))
360   ELSE IF( KRR == 3 ) THEN
361     ZCPH(:,:,:) = XCPD + XCPV *PTSTEP*  PRVS(:,:,:)              &
362                        + XCL  *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:))
363   ELSE IF( KRR == 2 ) THEN
364     ZCPH(:,:,:) = XCPD + XCPV *PTSTEP*  PRVS(:,:,:) &
365                        + XCL  *PTSTEP*  PRCS(:,:,:)
366   END IF
367 !
368 !
369 !*       3.     FIRST ORDER SUBGRID CONDENSATION SCHEME
370 !               ---------------------------------------
371 !
372   IF (OSUBG_COND) THEN
373 !
374 !*       3.1    compute condensate, cloud fraction
375 !
376   !   ZW3=water vapor    ZW1=rc (INOUT)  ZW2=ri (INOUT)   PSRC= s'rci'/Sigma_s^2
377     ZW3 = PRVS * PTSTEP;     ZW1 = PRCS * PTSTEP;  ZW2 = PRIS * PTSTEP
378
379     CALL CONDENSATION( IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE,1,         &
380        PPABST, PZZ, ZT, ZW3, ZW1, ZW2, PSIGS, PMFCONV, PCLDFR, PSRCS, .TRUE., &
381        OSIGMAS, PSIGQSAT )
382 !
383 !*       3.2    compute the variation of mixing ratio
384 !
385                                                         !         Rc - Rc*
386     ZW1(:,:,:) = (ZW1(:,:,:) / PTSTEP) - PRCS(:,:,:)    ! Pcon = ----------
387                                                         !         2 Delta t
388
389     ZW2(:,:,:) = (ZW2(:,:,:) / PTSTEP) - PRIS(:,:,:)    ! idem ZW1 but for Ri
390
391   ELSE
392 !
393 !
394 !*       4.     SECOND ORDER ALL OR NOTHING CONDENSATION SCHEME
395 !                            FOR MIXED-PHASE CLOUD
396 !               -----------------------------------------------
397 !
398 !
399 !*       4.1    Eventually pretreatment
400 !
401     IF (LPRETREATMENT) THEN
402 !
403 !     compute the saturation vapor pressures at t+1
404 !
405       CALL GET_HALO(ZT)
406       ZW1(:,:,:) = EXP(XALPW - XBETAW/ZT(:,:,:) - XGAMW*ALOG(ZT(:,:,:))) ! e_sw
407       ZW2(:,:,:) = EXP(XALPI - XBETAI/ZT(:,:,:) - XGAMI*ALOG(ZT(:,:,:))) ! e_si
408       ZW1(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW1(:,:,:))   ! safety limitation
409       ZW2(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW2(:,:,:))   ! safety limitation
410 !
411 !     compute the saturation mixing ratios at t+1
412 !
413       ZW3(:,:,:) = ZW1(:,:,:) * ZEPS / &
414                   ( PPABST(:,:,:) - ZW1(:,:,:))  ! r_sw
415       ZW4(:,:,:) = ZW2(:,:,:) * ZEPS / &
416                   ( PPABST(:,:,:) - ZW2(:,:,:))  ! r_si
417 !
418       WHERE(PRVS(:,:,:)*PTSTEP .LT. ZW4(:,:,:) .AND. &
419             PRCS(:,:,:) .GT. 0. .AND. ZT(:,:,:) .LT. XTT)
420 !
421 !       Subsaturation case with respect to rsi(T,P) (and case rv<0):
422 !       Evaporation of rc>0 (while enough) to decrease the lack of vapor 
423 !
424         ZW5 (:,:,:)= MIN( PRCS , ZW4(:,:,:)/PTSTEP - PRVS(:,:,:) )                ! RVCNDC
425         PRVS(:,:,:)= PRVS(:,:,:) + ZW5(:,:,:) 
426         PRCS(:,:,:)= PRCS(:,:,:) - ZW5(:,:,:) 
427         PTHS(:,:,:)= PTHS(:,:,:) - ZW5(:,:,:) * ZLV(:,:,:) /(ZCPH(:,:,:)*PEXNREF(:,:,:))
428 !
429       END WHERE
430 !  
431       WHERE (PRVS(:,:,:)*PTSTEP .GT. ZW3(:,:,:))
432 !
433 !       Supersaturation case with respect to rsw(T,P):
434 !       Condensation of the vapor that is left 
435 !
436         ZW5 (:,:,:)= PRVS(:,:,:) - ZW3(:,:,:)/PTSTEP 
437         PRVS(:,:,:)= PRVS(:,:,:) - ZW5(:,:,:)                                  ! RVCNDC 
438         PRCS(:,:,:)= PRCS(:,:,:) + ZW5(:,:,:) 
439         PTHS(:,:,:)= PTHS(:,:,:) + ZW5(:,:,:) * ZLV(:,:,:) /(ZCPH(:,:,:)*PEXNREF(:,:,:))
440 !
441       END WHERE
442 !
443       WHERE (PRCS(:,:,:) .GT. 0. .AND. ZT(:,:,:) .LT. ZT00)
444 !
445 !       Treatment of rc>0 if T<T00:
446 !
447         PRIS(:,:,:)= PRIS(:,:,:) + PRCS(:,:,:) 
448         PTHS(:,:,:)= PTHS(:,:,:) + PRCS(:,:,:) * &
449                     (ZLS(:,:,:) - ZLV(:,:,:)) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
450         PRCS(:,:,:)= 0. 
451 !
452       END WHERE
453 !
454 !*       4.2    compute the intermediate temperature at t+1, T*
455 !  
456       ZT(:,:,:) = (PTHS(:,:,:) * PTSTEP) * ZEXNS(:,:,:)
457 !
458     END IF  !end PRETREATMENT 
459 !
460 !*       4.3    compute the saturation vapor pressures at t+1
461 !
462     ZW1(:,:,:) = EXP(XALPW - XBETAW / ZT(:,:,:) - XGAMW * ALOG(ZT(:,:,:))) ! e_sw
463     ZW2(:,:,:) = EXP(XALPI - XBETAI / ZT(:,:,:) - XGAMI * ALOG(ZT(:,:,:))) ! e_si
464     ZW1(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW1(:,:,:))   ! safety limitation
465     ZW2(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW2(:,:,:))   ! safety limitation
466 !
467 !*       4.4    compute the saturation mixing ratios at t+1
468 !
469     ZW3(:,:,:) =  ZW1(:,:,:) * ZEPS     &
470                / ( PPABST(:,:,:)  - ZW1(:,:,:) ) ! r_sw
471     ZW4(:,:,:) =  ZW2(:,:,:) * ZEPS     &
472                / ( PPABST(:,:,:)  - ZW2(:,:,:) ) ! r_si
473 !
474 !*       4.5    compute the saturation mixing ratio derivatives (r'_vs)
475 !
476     ZW1(:,:,:) = (( XBETAW/ZT(:,:,:) - XGAMW ) / ZT(:,:,:)) & ! r'_sw
477                     * ZW3(:,:,:) * ( 1. + ZW3(:,:,:)/ZEPS )
478     ZW2(:,:,:) = (( XBETAI/ZT(:,:,:) - XGAMI ) / ZT(:,:,:)) & ! r'_si
479                     * ZW4(:,:,:) * ( 1. + ZW4(:,:,:)/ZEPS )
480 !
481     IF (LNEW_ADJUST) THEN
482       ZCND(:,:,:)= (ZT(:,:,:) - ZT00) / (ZT0 - ZT00)       ! Like Tao et al 89
483       ZCND(:,:,:)= MAX ( MIN(ZCND(:,:,:),1.) , 0. )
484     ELSE
485       WHERE ((PRCS(:,:,:)+PRIS(:,:,:)) .GT. 1.0E-20)      &
486       ZCND(:,:,:)= PRCS(:,:,:) / (PRCS(:,:,:) + PRIS(:,:,:)) ! Like the original version
487     END IF
488 !
489 !*       4.5    compute L_v CND + L_s DEP and F'(T)
490 !
491     WHERE ((PRCS(:,:,:)+PRIS(:,:,:)) .GT. 1.0E-20)
492 !
493       ZW5(:,:,:) = ZLS(:,:,:) + (ZLV(:,:,:) - ZLS(:,:,:)) * ZCND(:,:,:)
494       ZW6(:,:,:) = ZCPH(:,:,:) * (PRCS(:,:,:) + PRIS(:,:,:)) +          &
495                    ZW5(:,:,:)  * (PRCS(:,:,:) * ZW1(:,:,:)              &
496                                               + PRIS(:,:,:) * ZW2(:,:,:))
497 !
498 !*       4.6    compute Delta 2
499 !
500       ZW7(:,:,:) = (ZW5(:,:,:) / (ZW6(:,:,:) * ZT(:,:,:))) *                                 &
501                    (PRCS(:,:,:) * ZW1(:,:,:) *                                               &
502         ((-2. * XBETAW + XGAMW * ZT(:,:,:)) / (XBETAW - XGAMW * ZT(:,:,:)) +                 &
503                (XBETAW - XGAMW * ZT(:,:,:)) * (1.0 + 2.0 * ZW3(:,:,:) / ZEPS) / ZT(:,:,:)) + &
504                     PRIS(:,:,:) * ZW2(:,:,:) *                                               &
505         ((-2. * XBETAI + XGAMI * ZT(:,:,:)) / (XBETAI - XGAMI * ZT(:,:,:)) +                 &
506                (XBETAI - XGAMI * ZT(:,:,:)) * (1.0 + 2.0 * ZW4(:,:,:) / ZEPS) / ZT(:,:,:)))
507 !
508 !*       4.7    compute Delta 1
509 !
510       ZW6(:,:,:) = ZW5(:,:,:) * (PRCS(:,:,:) * ZW3(:,:,:) + PRIS(:,:,:) * ZW4(:,:,:) - &
511                                  PRVS(:,:,:) * PTSTEP * (PRCS(:,:,:) + PRIS(:,:,:))) / &
512                                  ZW6(:,:,:)
513 !
514 !*       4.8    compute the sources
515 !
516       ZW3(:,:,:) = (ZCPH(:,:,:) / ZW5(:,:,:)) *  &
517                    (-ZW6(:,:,:) * (1.0 + 0.5 * ZW6(:,:,:) * ZW7(:,:,:))) / PTSTEP
518       ZW1(:,:,:) = ZW3(:,:,:) * ZCND(:,:,:)                               ! RVCNDC
519       ZW2(:,:,:) = ZW3(:,:,:) * (1.0 - ZCND(:,:,:))                       ! RVDEPI
520 !
521     ELSEWHERE
522 !
523 !*       4.9    special case when both r_c and r_i are zero
524 !
525       ZW6(:,:,:) = ZCPH(:,:,:) + ZLV(:,:,:) * ZW1(:,:,:)               ! F'(T)
526       ZW7(:,:,:) = (ZLV(:,:,:) / (ZW6(:,:,:) * ZT(:,:,:))) *                & ! Delta 2
527                    (ZW1(:,:,:) *                                            &
528        ((-2. * XBETAW + XGAMW * ZT(:,:,:)) / (XBETAW - XGAMW * ZT(:,:,:)) + &
529               (XBETAW - XGAMW * ZT(:,:,:)) * (1.0 + 2.0 * ZW3(:,:,:) / ZEPS) / ZT(:,:,:)))
530       ZW6(:,:,:) = ZLV(:,:,:) * (ZW3(:,:,:) - PRVS(:,:,:) * PTSTEP) / ZW6(:,:,:) ! Delta 1
531       ZW1(:,:,:) = (ZCPH(:,:,:) / ZLV(:,:,:)) *                          & ! RVCNDC
532                    (-ZW6(:,:,:) * ( 1.0 + 0.5 * ZW6(:,:,:) * ZW7(:,:,:))) / PTSTEP
533       ZW2(:,:,:) = 0.0                                                     ! RVDEPI
534 !
535     END WHERE
536   END IF 
537 !
538 !*       5.     COMPUTE THE SOURCES AND STORES THE CLOUD FRACTION
539 !               -------------------------------------------------
540 !
541 !*       5.1    compute the sources 
542 !
543 !*       5.1.1  microphysics
544 !
545   WHERE (ZW1(:,:,:) < 0.0)
546     ZW1(:,:,:) = MAX (ZW1(:,:,:), -PRCS(:,:,:))        ! Evaporation rate
547   ELSEWHERE
548     ZW1(:,:,:) = MIN (ZW1(:,:,:),  PRVS(:,:,:))        ! Condensation rate
549   END WHERE
550 !
551   PRVS(:,:,:) = PRVS(:,:,:) - ZW1(:,:,:)
552   PRCS(:,:,:) = PRCS(:,:,:) + ZW1(:,:,:)  
553   PTHS(:,:,:) = PTHS(:,:,:) +        &
554                  ZW1(:,:,:) * ZLV(:,:,:) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
555 !
556   WHERE (ZW2(:,:,:) < 0.0)
557     ZW2(:,:,:) = MAX (ZW2(:,:,:), -PRIS(:,:,:))        ! Sublimation rate
558   ELSEWHERE
559     ZW2(:,:,:) = MIN (ZW2(:,:,:),  PRVS(:,:,:))        ! Deposition rate
560   END WHERE
561 !
562   PRVS(:,:,:) = PRVS(:,:,:) - ZW2(:,:,:)
563   PRIS(:,:,:) = PRIS(:,:,:) + ZW2(:,:,:)
564   PTHS(:,:,:) = PTHS(:,:,:) +        &
565                  ZW2(:,:,:) * ZLS(:,:,:) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
566 !
567 !*       5.1.2  electricity
568 !
569   ZWE1(:,:,:) = 0.
570   ZWE2(:,:,:) = 0.
571 !
572 ! the electrical process due to condensation is removed and
573 ! capture of ions by cloud droplets is done in ion_attach_elec routine
574 !
575 ! evaporation
576   WHERE (ABS(PRCT(:,:,:)) > XRTMIN_ELEC(2) .AND. &
577          ABS(PQCT(:,:,:)) > XQTMIN(2)      .AND. &
578               ZW1(:,:,:)  < -XRTMIN(1))
579     ZWE1(:,:,:) = (XFC / 3.) * (PQCT(:,:,:) / PRCT(:,:,:)) * (-ZW1(:,:,:))
580     ZION_NUMBER(:,:,:) = ABS(ZWE1(:,:,:)) / XECHARGE
581     ZADD(:,:,:) = 0.5 + SIGN(0.5, ZWE1(:,:,:))
582     PQPIS(:,:,:) = PQPIS(:,:,:) +  &
583                             ZADD(:,:,:) *ZION_NUMBER(:,:,:)
584     PQNIS(:,:,:) = PQNIS(:,:,:) +  &
585                             (1.-ZADD(:,:,:)) *ZION_NUMBER(:,:,:)
586     PQCS(:,:,:) = PQCS(:,:,:) - ZWE1(:,:,:)
587   END WHERE
588 !
589 ! the electrical process due to deposition is removed and
590 ! capture of ions by raindropsp is done in ion_attach_elec routine
591 !
592 ! sublimation
593   WHERE (ABS(PRIT(:,:,:)) > XRTMIN_ELEC(4) .AND. &
594          ABS(PQIT(:,:,:)) > XQTMIN(4)      .AND. &
595               ZW2(:,:,:)  < -XRTMIN(1))
596     ZWE2(:,:,:) = (XFI / XBI) * (PQIT(:,:,:) / PRIT(:,:,:)) * (-ZW2(:,:,:))
597     ZION_NUMBER(:,:,:) = ABS(ZWE2(:,:,:)) / XECHARGE
598     ZADD(:,:,:) = 0.5 + SIGN(0.5, ZWE2(:,:,:))
599     PQPIS(:,:,:) = PQPIS(:,:,:) +  &
600                             ZADD(:,:,:) *ZION_NUMBER(:,:,:)
601     PQNIS(:,:,:) = PQNIS(:,:,:) +  &
602                             (1.-ZADD(:,:,:)) *ZION_NUMBER(:,:,:)
603     PQIS(:,:,:) = PQIS(:,:,:) - ZWE2(:,:,:)
604   END WHERE
605 END DO          ! end of the iterative loop
606 !
607 !
608 !*       5.2    compute the cloud fraction PCLDFR
609 !
610 IF (.NOT. OSUBG_COND) THEN
611   WHERE (PRCS(:,:,:) + PRIS(:,:,:) > 1.E-12 / PTSTEP)
612     PCLDFR(:,:,:) = 1.
613   ELSEWHERE
614     PCLDFR(:,:,:) = 0. 
615   ENDWHERE 
616   IF (SIZE(PSRCS,3) /= 0) THEN
617     PSRCS(:,:,:) = PCLDFR(:,:,:) 
618   END IF
619 ELSE
620   IF (HSCONV == 'EDKF' .AND. HMF_CLOUD == 'DIRE') THEN
621     PCLDFR(:,:,:) = MIN(1.,PCLDFR(:,:,:)+PCF_MF(:,:,:))
622     PRCS(:,:,:)   = PRCS(:,:,:) + PRC_MF(:,:,:) / PTSTEP
623     PRIS(:,:,:)   = PRIS(:,:,:)+PRI_MF(:,:,:)/PTSTEP
624     PRVS(:,:,:)   = PRVS(:,:,:)- ( PRC_MF(:,:,:) + PRI_MF(:,:,:)) /PTSTEP
625     PTHS(:,:,:)   = PTHS(:,:,:) +  ( PRC_MF(:,:,:) * ZLV(:,:,:) +    &
626                     PRI_MF(:,:,:) * ZLS(:,:,:) ) / ZCPH(:,:,:) /     &
627                     PEXNREF(:,:,:) / PTSTEP
628   END IF
629 ENDIF
630 !
631 !
632 !
633 !*       6.  STORE THE BUDGET TERMS
634 !            ----------------------
635 !
636 IF (LBUDGET_RV) CALL BUDGET (PRVS(:,:,:) * PRHODJ(:,:,:),6,'DEPI_BU_RRV')
637 IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:) * PRHODJ(:,:,:),7,'DEPI_BU_RRC')
638 IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:) * PRHODJ(:,:,:),9,'DEPI_BU_RRI')
639 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:) * PRHODJ(:,:,:),4,'DEPI_BU_RTH')
640 IF (LBUDGET_SV) THEN
641   CALL BUDGET(PQPIS(:,:,:) * PRHODJ(:,:,:), 12+NSV_ELECBEG, 'DEPI_BU_RSV')
642   CALL BUDGET(PQNIS(:,:,:) * PRHODJ(:,:,:), 12+NSV_ELECEND, 'DEPI_BU_RSV')
643   CALL BUDGET(PQCS(:,:,:) * PRHODJ(:,:,:), 12+NSV_ELECBEG+1, 'DEPI_BU_RSV')
644   CALL BUDGET(PQIS(:,:,:) * PRHODJ(:,:,:), 12+NSV_ELECBEG+3, 'DEPI_BU_RSV')
645 END IF
646 !
647 !------------------------------------------------------------------------------
648 !
649 END SUBROUTINE ICE_ADJUST_ELEC