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