M.Mazoyer+O.Thouron 26/04/2016 Arrangements
[MNH-git_open_source-lfs.git] / src / MNH / prognos.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       MODULE MODI_PROGNOS
7 !     #######################
8 !
9 INTERFACE
10 !
11 SUBROUTINE PROGNOS(HLUOUT,PDT,PDZ,PLV,PCPH,PPRES,PRHOD,PRR,PTT,PRV,PRC,PS0,PCN,PCL)
12 !
13 CHARACTER(LEN=*),         INTENT(IN)    :: HLUOUT   ! Output-listing name for
14                                                     ! model n
15 REAL,                     INTENT(IN)    :: PDT
16 REAL, DIMENSION(:),       INTENT(IN)    :: PPRES
17 REAL, DIMENSION(:),       INTENT(IN)    :: PDZ
18 REAL, DIMENSION(:),       INTENT(IN)    :: PLV
19 REAL, DIMENSION(:),       INTENT(IN)    :: PCPH
20 REAL, DIMENSION(:),       INTENT(IN)    :: PRHOD
21 REAL, DIMENSION(:),       INTENT(IN)    :: PRR
22 REAL, DIMENSION(:),       INTENT(INOUT) :: PTT
23 REAL, DIMENSION(:),       INTENT(INOUT) :: PRV
24 REAL, DIMENSION(:),       INTENT(INOUT) :: PRC
25 REAL, DIMENSION(:),       INTENT(INOUT) :: PS0
26 REAL, DIMENSION(:),       INTENT(INOUT) :: PCN
27 REAL, DIMENSION(:),       INTENT(INOUT) :: PCL
28 !
29 END SUBROUTINE PROGNOS
30 !
31 END INTERFACE
32 !
33 END MODULE MODI_PROGNOS
34 END MODULE MODI_KHKO_NOTADJUST
35 !
36 !     ###################################################################################
37       SUBROUTINE PROGNOS(HLUOUT,PDT,PDZ,PLV,PCPH,PPRES,PRHOD,PRR,PTT,PRV,PRC,PS0,PCN,PCL)
38 !     ###################################################################################
39 !
40 !!****  * -  compute pseudo-prognostic of supersaturation according to Thouron
41 !                                                                     et al. 2012
42 !!    PURPOSE
43 !!    -------
44 !!
45 !!**  METHOD
46 !!
47 !!    REFERENCE
48 !!    ---------
49 !!
50 !!      Thouron, O., J.-L. Brenguier, and F. Burnet, Supersaturation calculation
51 !!      in large eddy simulation models for prediction of the droplet number
52 !!      concentration, Geosci. Model Dev., 5, 761-772, 2012.
53 !!
54 !!    AUTHOR
55 !!    ------
56 !!
57 !!      O. Thouron   * CNRM Meteo-France* : 
58 !!
59 !!    MODIFICATIONS
60 !!    -------------
61 !!     2014 G.Delautier : remplace MODD_RAIN_C2R2_PARAM par MODD_RAIN_C2R2_KHKO_PARAM
62 !!     2015 M.Mazoyer and O.Thouron : Physical tunings
63 !!
64 !-------------------------------------------------------------------------------
65 !
66 !*       0.    DECLARATIONS
67 !
68 USE MODD_CST
69 USE MODD_RAIN_C2R2_KHKO_PARAM
70 USE MODD_PARAM_C2R2
71 USE MODI_GAMMA
72 USE MODE_IO_ll
73 !
74 IMPLICIT NONE
75 !
76 !*       0.1   Declarations of dummy arguments :
77 !
78 !
79 !
80 CHARACTER(LEN=*),         INTENT(IN)    :: HLUOUT   ! Output-listing name for
81                                                     ! model n
82 REAL,                     INTENT(IN)    :: PDT
83 REAL, DIMENSION(:),       INTENT(IN)    :: PPRES
84 REAL, DIMENSION(:),       INTENT(IN)    :: PDZ
85 REAL, DIMENSION(:),       INTENT(IN)    :: PLV
86 REAL, DIMENSION(:),       INTENT(IN)    :: PCPH
87 REAL, DIMENSION(:),       INTENT(IN)    :: PRHOD
88 REAL, DIMENSION(:),       INTENT(IN)    :: PRR
89 REAL, DIMENSION(:),       INTENT(INOUT) :: PTT
90 REAL, DIMENSION(:),       INTENT(INOUT) :: PRV
91 REAL, DIMENSION(:),       INTENT(INOUT) :: PRC
92 REAL, DIMENSION(:),       INTENT(INOUT) :: PS0
93 REAL, DIMENSION(:),       INTENT(INOUT) :: PCN
94 REAL, DIMENSION(:),       INTENT(INOUT) :: PCL
95 !
96 !
97 !*       0.2   Declarations of local variables :
98 !
99 !
100 REAL, DIMENSION(SIZE(PRHOD,1))   ::  ZZW1,ZZW2,ZDZRC2,ZDZRC,ZCPH
101 REAL, DIMENSION(SIZE(PRHOD,1))   :: ZA1,ZA2,ZB,ZC,ZG
102 REAL, DIMENSION(SIZE(PRHOD,1))   :: ZLV,ZTT1,ZRT,ZTL,ZTT1_TEMP,ZTT2_TEMP
103 REAL, DIMENSION(SIZE(PRHOD,1))   :: ZRMOY,ZRVSAT1,ZRVSAT2
104 REAL, DIMENSION(SIZE(PRHOD,1))   :: ZVEC2  ! Work vectors forinterpolations
105 INTEGER, DIMENSION(SIZE(PRHOD,1)):: IVEC2   ! Vectors of indices for interpolations
106 INTEGER :: J1,J2
107 REAL,DIMENSION(SIZE(PS0,1))      ::MEM_PS0,ADJU2
108 REAL::AER_RAD
109 REAL, DIMENSION(SIZE(PRHOD,1))   :: ZFLAG_ACT   !Flag for activation
110 !
111 INTEGER                          :: IRESP      ! Return code of FM routines
112 INTEGER                          :: ILUOUT     ! Logical unit of output listing
113
114 !minimum radius of cloud droplet
115 AER_RAD=1.0E-6
116 !
117 ZFLAG_ACT(:)=0.0
118 !ACTIVATION
119 ZVEC2(:) =0.0
120 IVEC2(:) =0.0
121 !
122 DO J1 = 1,4
123   WHERE (PS0(:).GT.0.0)
124    ZVEC2(:) = MAX( 1.00001, MIN( FLOAT(NHYP)-0.00001,      &
125                    XHYPINTP1*LOG(PS0(:))+XHYPINTP2 ) )
126    IVEC2(:) = INT( ZVEC2(:) )
127    ZVEC2(:) = ZVEC2(:) - FLOAT( IVEC2(:) )
128   END WHERE
129 END DO
130 ZZW1(:) =0.0
131 WHERE (PS0(:).GT.0.0)
132 ZZW1(:) =   XHYPF12( IVEC2(:)+1 )* ZVEC2(:)      &
133              - XHYPF12( IVEC2(:)   )*(ZVEC2(:) - 1.0)
134 END WHERE
135 !
136 ! the CCN spectra formula uses ZSMAX in percent
137 !
138 ZZW2(:)=0.0
139 IF (XCONC_CCN > 0) THEN
140 WHERE (PS0(:).GT.0.0)
141    ZZW2(:) = MIN( XCONC_CCN,XCHEN * (100.0*PS0(:))**XKHEN * ZZW1(:) ) 
142 END WHERE
143 ELSE
144 WHERE (PS0(:).GT.0.0)
145    ZZW2(:) = XCHEN * (100.0*PS0(:))**XKHEN * ZZW1(:) 
146 END WHERE
147 ENDIF
148 ZZW2(:)=MAX( (ZZW2(:)-PCN(:)),0.0 )
149 !
150 WHERE (ZZW2(:).LT.1.0)   !Non physique d'activer moins d'une particule
151 ZZW2(:)=0
152 END WHERE
153 !
154 !
155 !FLAG ACTIVE A TRUE (1.0) si on active pas
156 DO J2=1,SIZE(PRC,1)
157  IF (ZZW2(J2).EQ.0.0) THEN
158  ZFLAG_ACT(J2)=1.0
159  ENDIF
160 ENDDO
161 !
162 ! Mean radius
163 ZRMOY(:)=0.0
164 DO J2=1,SIZE(PRC,1)
165  IF (PRC(J2).NE.0.0) THEN
166  ZRMOY(J2)=(MOMG(XALPHAC,XNUC,3.0)*4.0*XPI*PCL(J2)*XRHOLW/&
167         (3.0*PRC(J2)*PRHOD(J2)))**(1.0/3.0)
168  ZRMOY(J2)=(PCL(J2)*MOMG(XALPHAC,XNUC,1.0)/ZRMOY(J2))
169  ENDIF
170  ZRMOY(J2)=ZRMOY(J2)+(ZZW2(J2)*AER_RAD)
171 ENDDO
172 !
173 PCL(:)= PCL(:)+ZZW2(:)
174 PCN(:)= PCN(:)+ZZW2(:)
175 !
176 !CALCUL DE A1 => Estimation de (drs/dt)f
177 !T(=à determiner) avant forcage; T'(=PTT) apres forcage
178 !Calcul de ZTT1: calculé en inversant S0(T)jusqu'à T:
179 ! l'erreur faite sur cette inversion est supérieur à la précision
180 ! recherchée, on applique à rs(T') pour cxalculer le DT=T'-T qui
181 ! correspond à la variation rs(T')-rs(T). Permet de recuperer une valeur
182 ! correcte de DT et donc de determiner T comme T=T'-DT         
183 !ZRVSAT1=rs(T)
184 !
185 ZRVSAT1(:)=PRV(:)/(PS0(:)+1.0)
186 !ZTT1<--es(T) de rs(T)
187 ZTT1_TEMP(:)=PPRES(:)*((((XMV / XMD)/ZRVSAT1(:))+1.0)**(-1D0))
188 !ZTT1<--T de es(T)
189 ZTT1_TEMP(:)=LOG(ZTT1_TEMP(:)/610.8)
190 ZTT1_TEMP(:)=(31.25*ZTT1_TEMP(:) -17.5688*273.15)/(ZTT1_TEMP(:) - 17.5688)
191 !es(T')
192 ZZW1(:)=EXP(XALPW-XBETAW/PTT(:)-XGAMW*LOG(PTT(:)))
193 !ZRVSAT2=rs(T')
194 ZRVSAT2(:)=(XMV / XMD)*ZZW1(:)/(PPRES(:)-ZZW1(:))
195 !ZTT2<--es(T') de rs(T')
196 ZTT2_TEMP(:)=PPRES(:)*((((XMV / XMD)/ZRVSAT2(:))+1.0)**(-1D0))
197 !ZTT2<--T' de es(T')
198 IF (MINVAL(ZTT2_TEMP).LT.0.0) THEN
199   PRINT*,'ZTT2_TEMP',MINVAL(ZTT2_TEMP),MINLOC(ZTT2_TEMP)
200   CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
201   CALL ABORT
202   STOP
203 ENDIF
204 !
205 ZTT2_TEMP(:)=LOG(ZZW1(:)/610.8)
206 ZTT2_TEMP(:)=(31.25*ZTT2_TEMP(:) -17.5688*273.15)/(ZTT2_TEMP(:) - 17.5688)
207 !ZTT1=T'-DT
208 ZTT1(:)=PTT(:)-(ZTT2_TEMP(:)-ZTT1_TEMP(:))
209 !Lv(T)
210 ZLV(:) = XLVTT+(XCPV-XCL)*(ZTT1(:)-XTT)                
211 !
212 ZA1(:)=-(((PS0(:)+1.0)**2.0)/PRV(:))*(ZRVSAT2(:)-(PRV(:)/(PS0(:)+1.0)))/PDT
213 !G
214 ZG(:)= 1.0/(XRHOLW*((XRV*ZTT1(:)/(XDIVA*EXP(XALPW-(XBETAW/ZTT1(:))-(XGAMW*LOG(ZTT1(:))))))     &
215 +((ZLV(:)/(XTHCO*ZTT1(:)))*((ZLV(:)/(ZTT1(:)*XRV))-1.0))))
216 !
217 ZC(:)=4.0*XPI*(XRHOLW/PRHOD(:))*ZG(:)
218 ZDZRC(:)=0.0
219 ZDZRC(:)=ZC(:)*PS0(:)*ZRMOY(:)
220 MEM_PS0(:)=PS0(:)
221 !CALCUL DE B => Estimation de (drs/dT)ce
222 !T(=PTT) avant condensation; T'(=à determiner) apres condensation
223 !Lv(T),Cph(T)
224 ZLV(:) = XLVTT+(XCPV-XCL)*(PTT(:)-XTT)                
225 ZCPH(:)= XCPD+XCPV*PRV(:)+XCL*(PRC(:)+PRR(:))
226 !T'=T+(DT)ce
227 ZTT1(:)=PTT(:)+(ZDZRC(:)*PDT*ZLV(:)/ZCPH(:))
228 !es(T')
229 ZZW1(:)=EXP(XALPW-XBETAW/PTT(:)-XGAMW*LOG(PTT(:)))
230 !rs(T')
231 ZZW1(:)=(XMV / XMD)*ZZW1(:)/(PPRES(:)-ZZW1(:))
232 !es(Tcond)
233 ZZW2(:)=EXP(XALPW-XBETAW/ZTT1(:)-XGAMW*LOG(ZTT1(:)))
234 !rs(Tcond)
235 ZZW2(:)=(XMV / XMD)*ZZW2(:)/(PPRES(:)-ZZW2(:))
236 !
237 WHERE (ZTT1(:).NE.PTT(:))
238  ZB(:)=(ZLV(:)/ZCPH(:))*((ZZW2(:)-ZZW1(:))/(ZTT1(:)-PTT(:)))
239 ELSEWHERE
240  ZB(:)=0.0
241  ZDZRC(:)=0.0
242 ENDWHERE
243 !Calcul de S+dS
244 PS0(:)=PS0(:)+((ZA1(:)-(((ZB(:)*(PS0(:)+1.0)+1.0)*ZDZRC(:))/ZRVSAT1(:)))*PDT)
245 !
246 !Ajustement tel que rv=(s+1)*rvs
247 ZTL(:)=PTT(:)-(PLV(:)/PCPH(:))*PRC(:)
248 ZRT(:)=PRC(:)+PRV(:)
249 ZDZRC2(:)=PRC(:)
250 DO J2=1,SIZE(ZDZRC,1)
251   IF ((ZDZRC(J2).NE.0.0).OR.(ZDZRC2(J2).NE.0.0)) THEN 
252     DO J1=1,5
253      ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
254      ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
255      ZZW1(J2)=EXP(XALPW-XBETAW/PTT(J2)-XGAMW*LOG(PTT(J2)))
256      ZRVSAT1(J2)=(XMV / XMD)*ZZW1(J2)/(PPRES(J2)-ZZW1(J2))
257      PRV(J2)=MIN(ZRT(J2),(PS0(J2)+1.0)*ZRVSAT1(J2))
258      PRC(J2)=MAX(ZRT(J2)-PRV(J2),0.0)
259      PTT(J2)=0.5*PTT(J2)+0.5*(ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2)))
260     ENDDO
261     ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
262     ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
263     PTT(J2)=ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2))
264  ENDIF
265 ENDDO
266 ADJU2(:)=0.0
267 !
268 !Correction dans les mailles où ds a été surestimée
269 ZDZRC2(:)=PRC(:)-ZDZRC2(:)
270 WHERE ((MEM_PS0(:).LE.0.0).AND.(PS0(:).GT.0.0).AND.(ZDZRC2(:).LT.0.0))
271   PS0(:)=0.0
272   ADJU2(:)=1.0
273 ENDWHERE
274 !
275 WHERE ((MEM_PS0(:).GE.0.0).AND.(PS0(:).LT.0.0).AND.(ZDZRC2(:).GT.0.0))
276   PS0(:)=0.0
277   ADJU2(:)=1.0
278 ENDWHERE
279 !
280 DO J2=1,SIZE(ADJU2,1)
281   IF (ADJU2(J2)==1) THEN 
282    DO J1=1,5
283     ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
284     ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
285     ZZW1(J2)=EXP(XALPW-XBETAW/PTT(J2)-XGAMW*LOG(PTT(J2)))
286     ZRVSAT1(J2)=(XMV / XMD)*ZZW1(J2)/(PPRES(J2)-ZZW1(J2))
287     PRV(J2)=MIN(ZRT(J2),(PS0(J2)+1.0)*ZRVSAT1(J2))
288     PRC(J2)=MAX(ZRT(J2)-PRV(J2),0.0)
289     PTT(J2)=0.5*PTT(J2)+0.5*(ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2)))
290    ENDDO
291    ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
292    ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
293    PTT(J2)=ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2))
294   ENDIF
295 ENDDO
296 !
297 !Elimination de l'eau liquide dans les mailles où le rayon des gouttelettes est
298 !inférieur à AER_RAD
299 ZRMOY(:)=0.0
300 DO J2=1,SIZE(PRC,1)
301  IF (PRC(J2).NE.0.0) THEN
302   ZRMOY(J2)=(MOMG(XALPHAC,XNUC,3.0)*4.0*XPI*PCL(J2)*XRHOLW/&
303         (3.0*PRC(J2)*PRHOD(J2)))**(1.0/3.0)
304   ZRMOY(J2)=MOMG(XALPHAC,XNUC,1.0)/ZRMOY(J2)
305 IF ((ZFLAG_ACT(J2).EQ.1.0).AND.(MEM_PS0(J2).LT.0.0).AND.(ZRMOY(J2).LT.AER_RAD)) THEN
306    PTT(J2)=ZTL(J2)
307    PRV(J2)=ZRT(J2)
308    PRC(J2)=0.0
309   ENDIF
310  ENDIF
311 ENDDO
312 !
313 !Calcul de S au regard de T et rv en fin de pas de temps
314 ZZW1=EXP(XALPW-XBETAW/PTT(:)-XGAMW*LOG(PTT(:)))
315  !rvsat
316 ZRVSAT1(:)=(XMV / XMD)*ZZW1(:)/(PPRES-ZZW1(:))
317 !
318 WHERE (PRC(:)==0.0D0)
319  PS0(:)=(PRV(:)/ZRVSAT1(:))-1D0
320 ENDWHERE
321 !
322 CONTAINS
323 !
324 FUNCTION MOMG (PALPHA,PNU,PP) RESULT (PMOMG)
325 USE MODI_GAMMA
326 IMPLICIT NONE
327 REAL     :: PALPHA ! first shape parameter of the DIMENSIONnal distribution
328 REAL     :: PNU    ! second shape parameter of the DIMENSIONnal distribution
329 REAL     :: PP     ! order of the moment
330 REAL     :: PMOMG  ! result: moment of order ZP
331 PMOMG = GAMMA(PNU+PP/PALPHA)/GAMMA(PNU)
332 !
333 END FUNCTION MOMG
334 !
335 END SUBROUTINE PROGNOS