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