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