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