Philippe 14/11/2016: minor: removed unused IO arguments
[MNH-git_open_source-lfs.git] / src / MNH / fast_terms.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 !     ######################
12       MODULE MODI_FAST_TERMS
13 !     ######################
14 !
15 INTERFACE
16 !
17       SUBROUTINE FAST_TERMS( KRR, KMI, HLUOUT, HRAD,                    &
18                              HTURBDIM, HSCONV, HMF_CLOUD,               &
19                              OSUBG_COND, PTSTEP,                        &
20                              PRHODJ, PSIGS, PPABST,                     &
21                              PCF_MF,PRC_MF,                             &
22                              PRVT, PRCT, PRVS, PRCS, PRRS,              &
23                              PTHS, PSRCS, PCLDFR )
24          !
25 INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
26 INTEGER,                  INTENT(IN)    :: KMI      ! Model index 
27 CHARACTER(LEN=*),         INTENT(IN)    :: HLUOUT   ! Output-listing name for
28                                                     ! model n
29 CHARACTER*4,              INTENT(IN)    :: HTURBDIM ! Dimensionality of the
30                                                     ! turbulence scheme
31 CHARACTER(LEN=4),         INTENT(IN)    :: HSCONV   ! Shallow convection scheme
32 CHARACTER(LEN=4),         INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
33 CHARACTER*4,              INTENT(IN)    :: HRAD     ! Radiation scheme name
34 LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
35                                                     ! Condensation
36 REAL,                     INTENT(IN)    :: PTSTEP   ! Time step          
37 !
38 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODJ  ! Dry density * Jacobian
39 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PSIGS   ! Sigma_s at time t
40 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PPABST  ! Absolute Pressure at t     
41 !
42 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRVT    ! Water vapor m.r. at t
43 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRCT    ! Cloud water m.r. at t
44 !
45 REAL, DIMENSION(:,:,:),   INTENT(INOUT)       :: PRVS  ! Water vapor m.r. source
46 REAL, DIMENSION(:,:,:),   INTENT(INOUT)       :: PRCS  ! Cloud water m.r. source
47 REAL, DIMENSION(:,:,:), OPTIONAL,  INTENT(IN) :: PRRS  ! Rain  water m.r. source
48 !
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 !
53 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
54 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS   ! Second-order flux
55                                                    ! s'rc'/2Sigma_s2 at time t+1
56                                                    ! multiplied by Lambda_3
57 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PCLDFR  ! Cloud fraction          
58 !
59 END SUBROUTINE FAST_TERMS
60 !
61 END INTERFACE
62 !
63 END MODULE MODI_FAST_TERMS
64 !
65 !     ##########################################################################
66       SUBROUTINE FAST_TERMS( KRR, KMI, HLUOUT, HRAD,                    &
67                              HTURBDIM, HSCONV, HMF_CLOUD,               &
68                              OSUBG_COND, PTSTEP,                        &
69                              PRHODJ, PSIGS, PPABST,                     &
70                              PCF_MF,PRC_MF,                             &                   
71                              PRVT, PRCT, PRVS, PRCS, PRRS,              &
72                              PTHS, PSRCS, PCLDFR )
73 !     ##########################################################################
74 !
75 !!****  *FAST_TERMS* -  compute the fast  microphysical sources 
76 !!
77 !!    PURPOSE
78 !!    -------
79 !!    The purpose of this routine is to compute the fast microphysical sources
80 !!    through a saturation ajustement procedure.
81 !!
82 !!
83 !!**  METHOD
84 !!    ------
85 !!    Langlois, Tellus, 1973 for the cloudless version.
86 !!    When cloud water is taken into account, refer to book 1 of the
87 !!    documentation.
88 !!
89 !!     
90 !!
91 !!    EXTERNAL
92 !!    --------
93 !!      None
94 !!     
95 !!
96 !!    IMPLICIT ARGUMENTS
97 !!    ------------------
98 !!      Module MODD_CST
99 !!         XP00               ! Reference pressure
100 !!         XMD,XMV            ! Molar mass of dry air and molar mass of vapor
101 !!         XRD,XRV            ! Gaz constant for dry air, gaz constant for vapor
102 !!         XCPD,XCPV          ! Cpd (dry air), Cpv (vapor)
103 !!         XCL                ! Cl (liquid)
104 !!         XTT                ! Triple point temperature
105 !!         XLVTT              ! Vaporization heat constant
106 !!         XALPW,XBETAW,XGAMW ! Constants for saturation vapor 
107 !!                            !  pressure  function 
108 !!      Module  MODD_CONF 
109 !!         CCONF
110 !!      Module MODD_BUDGET:
111 !!         NBUMOD 
112 !!         CBUTYPE
113 !!         NBUPROCCTR 
114 !!         LBU_RTH    
115 !!         LBU_RRV  
116 !!         LBU_RRC  
117 !!
118 !!    REFERENCE
119 !!    ---------
120 !!
121 !!      Book 1 and Book2 of documentation ( routine FAST_TERMS )
122 !!      Langlois, Tellus, 1973
123 !!    AUTHOR
124 !!    ------
125 !!      E. Richard       * Laboratoire d'Aerologie*
126 !!   
127 !!
128 !!    MODIFICATIONS
129 !!    -------------
130 !!      Original    20/12/94 
131 !!      Modifications: March 1, 1995 (J.M. Carriere) 
132 !!                                  Introduction of cloud water with order 1
133 !!                                  formulation
134 !!      Modifications: June 8, 1995 ( J.Stein )
135 !!                                  Cleaning 
136 !!      Modifications: August 30, 1995 ( J.Stein )
137 !!                                  add Lambda3 for the subgrid condensation
138 !!   
139 !!                     October 16, 1995 (J. Stein)     change the budget calls 
140 !!                     March   16, 1996 (J. Stein)     store the cloud fraction
141 !!                     April   03, 1996 (J. Stein)     displace the nebulosity
142 !!                                      computation in the all and nothing case
143 !!                     April   15, 1996 (J. Stein)     displace the lambda 3 
144 !!                         multiplication and change the nebulosity threshold
145 !!                     September 16, 1996  (J. Stein)  bug in the SG cond for
146 !!                                                     the M computation
147 !!                     October 10, 1996 (J. Stein)     reformulate the Subgrid
148 !!                                                     condensation scheme
149 !!                     October 8,  1996 (Cuxart,Sanchez) Cloud frac. LES diag (XNA)
150 !!                     December 6, 1996 (J.-P. Pinty)  correction of Delta_2
151 !!                     November 5, 1996 (J. Stein) remove Rnp<0 values
152 !!                     November 13 1996 (V. Masson) add prints in test above
153 !!                     November 6, 2002 (V. Masson) update the budget calls
154 !!                     November 6, 2002 (S. Malardel,J.Pergaud) Cloud Fract + Rc of 
155 !!                                                              Mass flux convection
156 !!                                                              Scheme
157 !!                     J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
158 !!                     J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
159 !!                     June 17, 2016 (P. Wautelet) removed unused variables
160 !-------------------------------------------------------------------------------
161 !
162 !*       0.    DECLARATIONS
163 !              ------------
164 !
165 USE MODD_PARAMETERS
166 USE MODD_CST
167 USE MODD_CONF
168 USE MODD_BUDGET
169 !
170 USE MODI_CONDENS
171 USE MODI_BUDGET
172 USE MODE_FM
173 USE MODE_FMWRIT
174 USE MODI_GET_HALO
175 !
176 IMPLICIT NONE
177 !
178 !*       0.1   Declarations of dummy arguments :
179 !
180 !
181 INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
182 INTEGER,                  INTENT(IN)    :: KMI      ! Model index 
183 CHARACTER(LEN=*),         INTENT(IN)    :: HLUOUT   ! Output-listing name for
184                                                     ! model n
185 CHARACTER*4,              INTENT(IN)    :: HTURBDIM ! Dimensionality of the
186                                                     ! turbulence scheme
187 CHARACTER(LEN=4),         INTENT(IN)    :: HSCONV   ! Shallow convection scheme
188 CHARACTER(LEN=4),         INTENT(IN)    :: HMF_CLOUD! Type of statistical cloud
189 CHARACTER*4,              INTENT(IN)    :: HRAD     ! Radiation scheme name
190 LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
191                                                     ! Condensation
192 REAL,                     INTENT(IN)    :: PTSTEP   ! Time step          
193 !
194 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODJ  ! Dry density * Jacobian
195 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PSIGS   ! Sigma_s at time t
196 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PPABST  ! Absolute Pressure at t     
197 !
198 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRVT    ! Water vapor m.r. at t
199 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRCT    ! Cloud water m.r. at t
200 !
201 REAL, DIMENSION(:,:,:),   INTENT(INOUT)       :: PRVS  ! Water vapor m.r. source
202 REAL, DIMENSION(:,:,:),   INTENT(INOUT)       :: PRCS  ! Cloud water m.r. source
203 REAL, DIMENSION(:,:,:), OPTIONAL,  INTENT(IN) :: PRRS  ! Rain  water m.r. source
204 !!
205 REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PCF_MF! Convective Mass Flux Cloud fraction 
206 REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRC_MF! Convective Mass Flux liquid mixing ratio
207 !
208 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
209 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS   ! Second-order flux
210                                                    ! s'rc'/2Sigma_s2 at time t+1
211                                                    ! multiplied by Lambda_3
212 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PCLDFR  ! Cloud fraction          
213 !
214 !
215 !
216 !*       0.2   Declarations of local variables :
217 !
218 !
219 REAL  :: ZEPS  ! Mv/Md
220 REAL  :: ZTL,ZTHETA,ZRVS
221 REAL, DIMENSION(SIZE(PRHODJ,1),SIZE(PRHODJ,2),SIZE(PRHODJ,3)) &
222                          :: ZEXNS,&      ! guess of the Exner function at t+1
223                             ZT,   &      ! guess of the temperature at t+1
224                             ZCPH, &      ! guess of the CPh for the mixing
225                             ZLV,  &      ! guess of the Lv at t+1
226                             ZW1,ZW2,ZW3  ! Work arrays for intermediate
227                                          ! fields
228 !
229 INTEGER             :: IRESP      ! Return code of FM routines
230 INTEGER             :: ILENG      ! Length of comment string in LFIFM file
231 INTEGER             :: IGRID      ! C-grid indicator in LFIFM file
232 INTEGER             :: ILENCH     ! Length of comment string in LFIFM file
233 INTEGER             :: JK         ! Var for vertical DO loops
234 INTEGER             :: JITER,ITERMAX  ! iterative loop for first order adjustment
235 INTEGER             :: ILUOUT     ! Logical unit of output listing 
236 !-------------------------------------------------------------------------------
237 !
238 !*       1.     PRELIMINARIES
239 !               -------------
240 !
241 CALL FMLOOK_ll(HLUOUT,HLUOUT,ILUOUT,IRESP)
242 ZEPS= XMV / XMD
243 !
244 IF (OSUBG_COND) THEN
245   ITERMAX=2
246 ELSE
247   ITERMAX=1
248 END IF
249 !
250 !
251 !-------------------------------------------------------------------------------
252 !
253 !*       2.     COMPUTE QUANTITIES WITH THE GUESS OF THE FUTURE INSTANT
254 !               -------------------------------------------------------
255 !
256 !*       2.1    remove negative non-precipitating negative water
257 !               ------------------------------------------------
258 !
259 IF( NVERB>5) THEN
260   IF (ANY(PRCS(:,:,:)+PRVS(:,:,:) < 0.) ) THEN
261     WRITE(ILUOUT,*) 'FAST_TERMS:  negative values of total water (reset to zero)'
262     WRITE(ILUOUT,*) '  location of minimum:', MINLOC(PRCS(:,:,:)+PRVS(:,:,:))
263     WRITE(ILUOUT,*) '  value of minimum   :', MINVAL(PRCS(:,:,:)+PRVS(:,:,:))
264   END IF
265 END IF
266 !
267 WHERE ( PRCS(:,:,:)+PRVS(:,:,:) < 0.)
268   PRVS(:,:,:) = -  PRCS(:,:,:)
269 END WHERE
270 !
271 !*       2.2    estimate the Exner function at t+1
272 !
273 ZEXNS(:,:,:) = ( PPABST(:,:,:)  / XP00 ) ** (XRD/XCPD)  
274 !
275 !    beginning of the iterative loop
276 !
277 DO JITER =1,ITERMAX
278 !
279 !*       2.3    compute the intermediate temperature at t+1, T*
280 !  
281   ZT(:,:,:) = ( PTHS(:,:,:) * PTSTEP ) * ZEXNS(:,:,:)
282 !
283 !*       2.4    compute the latent heat of vaporization Lv(T*) at t+1
284 !
285   ZLV(:,:,:) = XLVTT + ( XCPV - XCL ) * ( ZT(:,:,:) -XTT )
286 !
287 !*       2.5    compute the specific heat for moist air (Cph) at t+1
288 !
289   IF     ( KRR == 3 ) THEN 
290     ZCPH(:,:,:) =  XCPD + XCPV *PTSTEP*   PRVS(:,:,:)       &
291                         + XCL  *PTSTEP* ( PRCS(:,:,:) + PRRS(:,:,:) )
292   ELSE IF( KRR == 2 ) THEN
293     ZCPH(:,:,:) =  XCPD + XCPV *PTSTEP*   PRVS(:,:,:)       &
294                         + XCL  *PTSTEP*   PRCS(:,:,:)  
295   END IF
296 !
297 !*       2.6    compute the saturation vapor pressure at t+1
298 !
299   CALL GET_HALO(ZT)
300   ZW1(:,:,:) = EXP( XALPW - XBETAW/ZT(:,:,:) - XGAMW*ALOG( ZT(:,:,:) ) )
301 !
302 !*       2.7    compute the saturation mixing ratio at t+1
303 !
304   ZW2(:,:,:) =  ZW1(:,:,:) * ZEPS /              &
305                 (  PPABST(:,:,:) - ZW1(:,:,:) )   
306 !
307 !*       2.8    compute the saturation mixing ratio derivative (rvs')
308 !
309   ZW1(:,:,:) = ( XBETAW / ZT(:,:,:)  - XGAMW ) / ZT(:,:,:)   & 
310                * ZW2(:,:,:) * ( 1. + ZW2(:,:,:) / ZEPS )
311 !
312 !*       2.9    compute  Cph + Lv * rvs'
313 !
314   ZW3(:,:,:) = ZCPH(:,:,:) + ZLV(:,:,:) * ZW1(:,:,:)
315 !
316 !
317 !-------------------------------------------------------------------------------
318 !
319 !*       3.     FIRST ORDER SUBGRID CONDENSATION SCHEME
320 !               ---------------------------------------
321 !
322   IF ( OSUBG_COND ) THEN
323 !
324 !*       3.1    compute  Q1              
325 !
326     ZW2(:,:,:) = ( ( PRVS(:,:,:)*PTSTEP - ZW2(:,:,:) ) * ZCPH(:,:,:) / ZW3(:,:,:) &
327                   + PRCS(:,:,:)*PTSTEP                                            &
328                  ) / ( 2. * PSIGS(:,:,:) )
329
330 !
331 !*       3.2    compute s'rc'/2Sigma_s2, Rc and the nebolisity
332 !
333     CALL CONDENS(HTURBDIM, ZW2,ZW1,ZW3,PSRCS) ! ZW1 = Cloud fraction 
334                                               ! PSRC = s'rc'/(2 Sigma_s**2) 
335                                               ! ZW3 = Rc / (2 Sigma_s)
336     ZW3(:,:,:) = 2. * PSIGS(:,:,:) * ZW3(:,:,:)    ! Rc 
337 !
338 !    multiply PSRCS by the lambda3 coefficient
339 !
340     IF ( HTURBDIM == '1DIM'.AND.JITER==ITERMAX ) THEN
341       PSRCS(:,:,:) = PSRCS(:,:,:) * MIN( 3. , MAX(1.,1.-ZW2(:,:,:)) )
342     END IF       ! in the 3D case lamda_3 = 1.
343 !
344 !*       3.3    compute the variation of mixing ratio
345 !
346                                                       !         Rc - Rc*
347     ZW3(:,:,:) = (ZW3(:,:,:)/PTSTEP) - PRCS(:,:,:)       ! Pcon = ----------
348                                                       !         2 Delta t
349   ELSE
350 !
351 !
352 !*       4.     SECOND ORDER ALL OR NOTHING CONDENSATION SCHEME
353 !               -----------------------------------------------
354 !
355 !*       4.1    compute Delta 2
356 !
357     ZW1(:,:,:) = (ZW1(:,:,:) * ZLV(:,:,:) / ZW3(:,:,:) ) *                     &
358           ( ((-2.*XBETAW+XGAMW*ZT(:,:,:)) / (XBETAW-XGAMW*ZT(:,:,:))           &
359                + (XBETAW/ZT(:,:,:)-XGAMW)*(1.0+2.0*ZW2(:,:,:)/ZEPS))/ZT(:,:,:) )
360 !
361 !*       4.2    compute Delta 1
362 !
363     ZW2(:,:,:) = ZLV(:,:,:) / ZW3(:,:,:) * ( ZW2(:,:,:) - PRVS(:,:,:)*PTSTEP )
364 !
365 !*       4.3    compute the variation of mixing ratio
366 !
367     ZW3(:,:,:) = - ZW2(:,:,:) * ( 1 + 0.5 * ZW2(:,:,:) * ZW1(:,:,:) )  &
368                  * ZCPH(:,:,:) / ZLV(:,:,:) / PTSTEP
369 !
370 !  end of the IF structure on the all or nothing or statistical condensation
371 !  scheme
372 !
373   END IF
374 !
375 !
376 !*       5.     COMPUTE THE SOURCES AND STORES THE CLOUD FRACTION
377 !               -------------------------------------------------
378 !
379 !
380 !*       5.1    compute the sources 
381
382   ZW3(:,:,:) = MAX ( ZW3(:,:,:), -PRCS(:,:,:) )  
383   WHERE (ZW3(:,:,:) > 0.0)
384     ZW3(:,:,:) = MIN ( ZW3(:,:,:),  PRVS(:,:,:) )
385   END WHERE
386   PRCS(:,:,:) = PRCS(:,:,:) + ZW3(:,:,:)
387   PRVS(:,:,:) = PRVS(:,:,:) - ZW3(:,:,:) 
388   PTHS(:,:,:) = PTHS(:,:,:) + ZW3(:,:,:) * ZLV(:,:,:) / ZCPH(:,:,:)     &
389                                          / ZEXNS(:,:,:)
390 !
391 !  end of the iterative loop
392 !
393 END DO
394 !
395 !*       5.2    compute the cloud fraction PCLDFR
396 !
397 IF ( .NOT. OSUBG_COND ) THEN
398   WHERE (PRCS(:,:,:) > 1.E-12 / PTSTEP)
399     ZW1(:,:,:)  = 1.
400   ELSEWHERE
401     ZW1(:,:,:)  = 0. 
402   ENDWHERE 
403   IF ( SIZE(PSRCS,3) /= 0 ) THEN
404     PSRCS(:,:,:) = ZW1(:,:,:) 
405   END IF
406 ELSE 
407   IF ( HRAD /= 'NONE' ) THEN
408     PCLDFR(:,:,:) = ZW1(:,:,:)
409   END IF
410 !
411   IF (HSCONV == 'EDKF' .AND. HMF_CLOUD == 'DIRE') THEN
412     PCLDFR(:,:,:) = MIN(1.,(ZW1(:,:,:)+ PCF_MF(:,:,:)))
413     PRCS(:,:,:) = PRCS(:,:,:)+PRC_MF(:,:,:)/PTSTEP
414     PRVS(:,:,:) = PRVS(:,:,:)-PRC_MF(:,:,:)/PTSTEP
415     PTHS(:,:,:) = PTHS(:,:,:) + PRC_MF * ZLV(:,:,:) / ZCPH(:,:,:)     &
416                                        / ZEXNS(:,:,:) /PTSTEP
417   END IF
418 ENDIF
419 !
420 !
421 !
422 !        6.  Horizontal mean Cloud fraction (for LES uses)
423 !            ---------------------------------------------
424 !
425
426 !
427 !*       7.  STORE THE BUDGET TERMS
428 !            ----------------------
429 !
430 !
431 IF (LBUDGET_RV) CALL BUDGET (PRVS(:,:,:) * PRHODJ(:,:,:),6,'COND_BU_RRV')
432 IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:) * PRHODJ(:,:,:),7,'COND_BU_RRC')
433 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:) * PRHODJ(:,:,:),4,'COND_BU_RTH')
434 !
435 !------------------------------------------------------------------------------
436 !
437 !
438 END SUBROUTINE FAST_TERMS