source: lmdz_wrf/trunk/tools/module_ForDiagnostics.f90 @ 2073

Last change on this file since 2073 was 1909, checked in by lfita, 7 years ago

Adding:

  • `fog_FRAML50': fog and visibility following Gultepe and Milbrandt, (2010)
  • `var_hus': relative humidity using August-Roche-Magnus approximation [1]

Fixing: fog_K84' and fog_RUC'

File size: 33.1 KB
Line 
1!! Fortran version of different diagnostics
2! L. Fita. LMD May 2016
3! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90
4!
5! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90
6
7MODULE module_ForDiagnostics
8
9  USE module_definitions
10  USE module_generic
11  USE module_ForDiagnosticsVars
12
13  CONTAINS
14
15!!!!!!! Calculations
16! compute_cape_afwa4D: Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute
17!   CAPE, CIN, ZLFC, PLFC, LI
18! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being
19!   3rd dimension the z-dim
20! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure being
21!   3rd dimension the z-dim
22! compute_cllmh: Computation of low, medium and high cloudiness
23! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
24! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
25! compute_clt: Computation of total cloudiness
26! compute_fog_K84: Computation of fog and visibility following Kunkel, (1984)
27! compute_fog_RUC: Computation of fog and visibility following RUC method Smirnova, (2000)
28! compute_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
29! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
30! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates
31! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
32! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates
33! compute_zmla_generic4D: Subroutine to compute pbl-height following a generic method
34! compute_zwind4D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
35! compute_zwind_log4D: Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
36! compute_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog
37
38!!!
39! Calculations
40!!!
41
42  SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
43! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure
44!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
45! It should be properly done via an 'INTERFACE', but...
46
47    IMPLICIT NONE
48
49    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
50    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D, pres4D
51    REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out)        :: cllmh4D2
52
53! Local
54    INTEGER                                              :: i,j,k, zdim, Ndim
55
56!!!!!!! Variables
57! cldfra4D: 4D cloud fraction values [1]
58! pres4D: 4D pressure values [Pa]
59! Ndim: number of dimensions of the input data
60! d[1-4]: dimensions of 'cldfra4D'
61! zdim: number of the vertical-dimension within the matrix
62! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'
63
64    fname = 'compute_cllmh4D2'
65    zdim = 2
66    Ndim = 4
67
68    DO i=1, d1
69      DO j=1, d3
70        DO k=1, d4
71          cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
72        END DO
73      END DO
74    END DO
75   
76    RETURN
77
78  END SUBROUTINE compute_cllmh4D2
79
80  SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
81! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure
82!   where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
83! It should be properly done via an 'INTERFACE', but...
84
85    IMPLICIT NONE
86
87    INTEGER, INTENT(in)                                  :: d1, d2, d3
88    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D, pres3D
89    REAL(r_k), DIMENSION(3,d2,d3), INTENT(out)           :: cllmh3D1
90
91! Local
92    INTEGER                                              :: i,j,k, zdim, Ndim
93
94!!!!!!! Variables
95! cldfra3D: 3D cloud fraction values [1]
96! pres3D: 3D pressure values [Pa]
97! Ndim: number of dimensions of the input data
98! d[1-3]: dimensions of 'cldfra3D'
99! zdim: number of the vertical-dimension within the matrix
100! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'
101
102    fname = 'compute_cllmh3D1'
103    zdim = 1
104    Ndim = 3
105
106    DO i=1, d1
107      DO j=1, d2
108        cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
109      END DO
110    END DO
111   
112    RETURN
113
114  END SUBROUTINE compute_cllmh3D1
115
116  SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D,    &
117    Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2,        &
118    cllmh4D3, cllmh4D4, d1, d2, d3, d4)
119! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ
120
121    IMPLICIT NONE
122
123    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
124    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D, pres1D
125    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D, pres2D
126    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
127    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
128      INTENT(in)                                         :: cldfra4D, pres4D
129    REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out)       :: cllmh1D
130    REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out)    :: cllmh2D1
131    REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out)    :: cllmh2D2
132    REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
133    REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
134    REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
135    REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL,                                                       &
136      INTENT(out)                                        :: cllmh4D1
137    REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL,                                                       &
138      INTENT(out)                                        :: cllmh4D2
139    REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL,                                                       &
140      INTENT(out)                                        :: cllmh4D3
141    REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL,                                                       &
142      INTENT(out)                                        :: cllmh4D4
143
144! Local
145    INTEGER                                              :: i,j,k
146
147!!!!!!! Variables
148! cldfra[1-4]D: cloud fraction values [1]
149! pres[1-4]D: pressure values [Pa]
150! Ndim: number of dimensions of the input data
151! d[1-4]: dimensions of 'cldfra'
152! zdim: number of the vertical-dimension within the matrix
153! cllmh1D: low, medium and high cloudiness for the 1D cldfra
154! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
155! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
156! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
157! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
158! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
159! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
160! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
161! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
162! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'
163
164    fname = 'compute_cllmh'
165
166    SELECT CASE (Ndim)
167      CASE (1)
168        cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
169      CASE (2)
170        IF (zdim == 1) THEN
171          DO i=1, d2
172            cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
173          END DO
174        ELSE IF (zdim == 2) THEN
175          DO i=1, d1
176            cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
177          END DO
178        ELSE
179          PRINT *,TRIM(ErrWarnMsg('err'))
180          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
181          PRINT *,'    accepted values: 1,2'
182          STOP
183        END IF
184      CASE (3)
185        IF (zdim == 1) THEN
186          DO i=1, d2
187            DO j=1, d3
188              cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
189            END DO
190          END DO
191        ELSE IF (zdim == 2) THEN
192          DO i=1, d1
193            DO j=1, d3
194              cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
195            END DO
196          END DO
197        ELSE IF (zdim == 3) THEN
198          DO i=1, d1
199            DO j=1, d2
200              cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
201            END DO
202          END DO
203        ELSE
204          PRINT *,TRIM(ErrWarnMsg('err'))
205          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
206          PRINT *,'    accepted values: 1,2,3'
207          STOP
208        END IF
209      CASE (4)
210        IF (zdim == 1) THEN
211          DO i=1, d2
212            DO j=1, d3
213              DO k=1, d4
214                cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
215              END DO
216            END DO
217          END DO
218        ELSE IF (zdim == 2) THEN
219          DO i=1, d1
220            DO j=1, d3
221              DO k=1, d4
222                cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
223              END DO
224            END DO
225          END DO
226        ELSE IF (zdim == 3) THEN
227          DO i=1, d2
228            DO j=1, d3
229              DO k=1, d4
230                cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
231              END DO
232            END DO
233          END DO
234        ELSE IF (zdim == 4) THEN
235          DO i=1, d1
236            DO j=1, d2
237              DO k=1, d3
238                cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
239              END DO
240            END DO
241          END DO
242        ELSE
243          PRINT *,TRIM(ErrWarnMsg('err'))
244          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
245          PRINT *,'    accepted values: 1,2,3,4'
246          STOP
247        END IF
248      CASE DEFAULT
249        PRINT *,TRIM(ErrWarnMsg('err'))
250        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
251        STOP
252      END SELECT
253
254    RETURN
255
256  END SUBROUTINE compute_cllmh
257
258  SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
259! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
260!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
261! It should be properly done via an 'INTERFACE', but...
262
263    IMPLICIT NONE
264
265    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
266    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D
267    REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out)          :: clt4D2
268
269! Local
270    INTEGER                                              :: i,j,k, zdim, Ndim
271
272!!!!!!! Variables
273! cldfra4D: 4D cloud fraction values [1]
274! Ndim: number of dimensions of the input data
275! d[1-4]: dimensions of 'cldfra4D'
276! zdim: number of the vertical-dimension within the matrix
277! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
278
279    fname = 'compute_clt4D2'
280    zdim = 2
281    Ndim = 4
282
283    DO i=1, d1
284      DO j=1, d3
285        DO k=1, d4
286          clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
287        END DO
288      END DO
289    END DO
290   
291    RETURN
292
293  END SUBROUTINE compute_clt4D2
294
295  SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
296! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
297!   where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
298! It should be properly done via an 'INTERFACE', but...
299
300    IMPLICIT NONE
301
302    INTEGER, INTENT(in)                                  :: d1, d2, d3
303    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D
304    REAL(r_k), DIMENSION(d2,d3), INTENT(out)             :: clt3D1
305
306! Local
307    INTEGER                                              :: i,j,k, zdim, Ndim
308
309!!!!!!! Variables
310! cldfra3D: 3D cloud fraction values [1]
311! Ndim: number of dimensions of the input data
312! d[1-3]: dimensions of 'cldfra3D'
313! zdim: number of the vertical-dimension within the matrix
314! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
315
316    fname = 'compute_clt3D1'
317    zdim = 1
318    Ndim = 3
319
320    DO i=1, d2
321      DO j=1, d3
322        clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
323      END DO
324    END DO
325   
326    RETURN
327
328  END SUBROUTINE compute_clt3D1
329
330  SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2,   &
331    clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
332! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ
333
334    IMPLICIT NONE
335
336    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
337    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D
338    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D
339    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
340    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
341      INTENT(in)                                         :: cldfra4D
342    REAL(r_k), OPTIONAL, INTENT(out)                     :: clt1D
343    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out)      :: clt2D1
344    REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out)      :: clt2D2
345    REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out)   :: clt3D1
346    REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out)   :: clt3D2
347    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out)   :: clt3D3
348    REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
349    REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
350    REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
351    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4
352
353! Local
354    INTEGER                                              :: i,j,k
355
356!!!!!!! Variables
357! cldfra[1-4]D: cloud fraction values [1]
358! Ndim: number of dimensions of the input data
359! d[1-4]: dimensions of 'cldfra'
360! zdim: number of the vertical-dimension within the matrix
361! clt1D: total cloudiness for the 1D cldfra
362! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
363! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
364! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
365! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
366! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
367! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
368! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
369! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
370! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'
371
372    fname = 'compute_clt'
373
374    SELECT CASE (Ndim)
375      CASE (1)
376        clt1D = var_clt(cldfra1D, d1)
377      CASE (2)
378        IF (zdim == 1) THEN
379          DO i=1, d2
380            clt2D1(i) = var_clt(cldfra2D(:,i), d1)
381          END DO
382        ELSE IF (zdim == 2) THEN
383          DO i=1, d1
384            clt2D2(i) = var_clt(cldfra2D(:,i), d2)
385          END DO
386        ELSE
387          PRINT *,TRIM(ErrWarnMsg('err'))
388          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
389          PRINT *,'    accepted values: 1,2'
390          STOP
391        END IF
392      CASE (3)
393        IF (zdim == 1) THEN
394          DO i=1, d2
395            DO j=1, d3
396              clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
397            END DO
398          END DO
399        ELSE IF (zdim == 2) THEN
400          DO i=1, d1
401            DO j=1, d3
402              clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
403            END DO
404          END DO
405        ELSE IF (zdim == 3) THEN
406          DO i=1, d1
407            DO j=1, d2
408              clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
409            END DO
410          END DO
411        ELSE
412          PRINT *,TRIM(ErrWarnMsg('err'))
413          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
414          PRINT *,'    accepted values: 1,2,3'
415          STOP
416        END IF
417      CASE (4)
418        IF (zdim == 1) THEN
419          DO i=1, d2
420            DO j=1, d3
421              DO k=1, d4
422                clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
423              END DO
424            END DO
425          END DO
426        ELSE IF (zdim == 2) THEN
427          DO i=1, d1
428            DO j=1, d3
429              DO k=1, d4
430                clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
431              END DO
432            END DO
433          END DO
434        ELSE IF (zdim == 3) THEN
435          DO i=1, d2
436            DO j=1, d3
437              DO k=1, d4
438                clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
439              END DO
440            END DO
441          END DO
442        ELSE IF (zdim == 4) THEN
443          DO i=1, d1
444            DO j=1, d2
445              DO k=1, d3
446                clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
447              END DO
448            END DO
449          END DO
450        ELSE
451          PRINT *,TRIM(ErrWarnMsg('err'))
452          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
453          PRINT *,'    accepted values: 1,2,3,4'
454          STOP
455        END IF
456      CASE DEFAULT
457        PRINT *,TRIM(ErrWarnMsg('err'))
458        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
459        STOP
460      END SELECT
461
462    RETURN
463
464  END SUBROUTINE compute_clt
465
466  SUBROUTINE compute_massvertint1D(var, mutot, dz, deta, integral)
467    ! Subroutine to vertically integrate a 1D variable in eta vertical coordinates
468
469    IMPLICIT NONE
470
471    INTEGER, INTENT(in)                                  :: dz
472    REAL(r_k), INTENT(in)                                :: mutot
473    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: var, deta
474    REAL(r_k), INTENT(out)                               :: integral
475
476! Local
477    INTEGER                                              :: k
478
479!!!!!!! Variables
480! var: vertical variable to integrate (assuming kgkg-1)
481! mutot: total dry-air mass in column
482! dz: vertical dimension
483! deta: eta-levels difference between full eta-layers
484
485    fname = 'compute_massvertint1D'
486
487!    integral=0.
488!    DO k=1,dz
489!      integral = integral + var(k)*deta(k)
490!    END DO
491     integral = SUM(var*deta)
492
493    integral=integral*mutot/g
494
495    RETURN
496
497  END SUBROUTINE compute_massvertint1D
498
499  SUBROUTINE compute_zint4D(var4D, dlev, zweight, d1, d2, d3, d4, int3D)
500    ! Subroutine to vertically integrate a 4D variable in any vertical coordinates
501
502    IMPLICIT NONE
503
504    INTEGER, INTENT(in)                                  :: d1,d2,d3,d4
505    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: var4D, dlev, zweight
506    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: int3D
507
508! Local
509    INTEGER                                              :: i,j,l
510
511!!!!!!! Variables
512! var4D: vertical variable to integrate
513! dlev: height of layers
514! zweight: weight for each level to be applied (=1. for no effect)
515
516    fname = 'compute_zint4D'
517
518    DO i=1,d1
519      DO j=1,d2
520        DO l=1,d4
521          CALL compute_vertint1D(var4D(i,j,:,l),d3, dlev(i,j,:,l), zweight(i,j,:,l), &
522            int3D(i,j,l))
523        END DO
524      END DO
525    END DO
526
527    RETURN
528
529  END SUBROUTINE compute_zint4D
530
531  SUBROUTINE compute_vertint1D(var, dz, deta, zweight, integral)
532    ! Subroutine to vertically integrate a 1D variable in any vertical coordinates
533
534    IMPLICIT NONE
535
536    INTEGER, INTENT(in)                                  :: dz
537    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: var, deta, zweight
538    REAL(r_k), INTENT(out)                               :: integral
539
540! Local
541    INTEGER                                              :: k
542
543!!!!!!! Variables
544! var: vertical variable to integrate
545! dz: vertical dimension
546! deta: eta-levels difference between layers
547! zweight: weight for each level to be applied (=1. for no effect)
548
549    fname = 'compute_vertint1D'
550
551!    integral=0.
552!    DO k=1,dz
553!      integral = integral + var(k)*deta(k)
554!    END DO
555    integral = SUM(var*deta*zweight)
556
557    RETURN
558
559  END SUBROUTINE compute_vertint1D
560
561  SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod,    &
562    d1, d2, d3, d4)
563! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI
564
565    IMPLICIT NONE
566
567    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4, parcelmethod
568    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ta, hur, press, zg
569    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
570    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: cape, cin, zlfc, plfc, li
571 
572! Local
573    INTEGER                                              :: i, j, it
574    INTEGER                                              :: ofunc
575
576!!!!!!! Variables
577! ta: air temperature [K]
578! hur: relative humidity [%]
579! press: air pressure [Pa]
580! zg: geopotential height [gpm]
581! hgt: topographical height [m]
582! cape: Convective available potential energy [Jkg-1]
583! cin: Convective inhibition [Jkg-1]
584! zlfc: height at the Level of free convection [m]
585! plfc: pressure at the Level of free convection [Pa]
586! li: lifted index [1]
587! parcelmethod:
588!   Most Unstable = 1 (default)
589!   Mean layer = 2
590!   Surface based = 3
591
592    fname = 'compute_cape_afwa4D'
593
594    DO i=1, d1
595      DO j=1, d2
596        DO it=1, d4
597          ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it),     &
598            1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod)
599          IF (zlfc(i,j,it) /= -1.) zlfc(i,j,it) = zlfc(i,j,it) - hgt(i,j)
600        END DO
601      END DO
602    END DO
603
604    RETURN
605
606  END SUBROUTINE compute_cape_afwa4D
607
608  SUBROUTINE compute_psl_ecmwf(ps, hgt, T, press, unpress, psl, d1, d2, d4)
609! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
610
611    IMPLICIT NONE
612
613    INTEGER, INTENT(in)                                  :: d1, d2, d4
614    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: ps, T, press, unpress
615    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
616    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: psl
617 
618! Local
619    INTEGER                                              :: i, j, it
620
621!!!!!!! Variables
622! ps: surface pressure [Pa]
623! hgt: terrain height [m]
624! T: temperature at first half-mass level [K]
625! press: pressure at first full levels [Pa]
626! unpress: pressure at first mass (half) levels [Pa]
627! psl: sea-level pressure [Pa]
628
629    fname = 'compute_psl_ecmwf'
630
631    DO i=1, d1
632      DO j=1, d2
633        DO it=1, d4
634          CALL var_psl_ecmwf(ps(i,j,it), hgt(i,j), T(i,j,it), unpress(i,j,it), press(i,j,it),         &
635            psl(i,j,it))
636        END DO
637      END DO
638    END DO
639
640    RETURN
641
642  END SUBROUTINE compute_psl_ecmwf
643
644  SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4)
645! Subroutine to compute pbl-height following a generic method
646!    from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
647!    applied also in Garcia-Diez et al., 2013, QJRMS
648!   where
649!     "The technique identifies the ML height as a threshold increase of potential temperature from
650!       its minimum value within the boundary layer."
651!   here applied similarly to Garcia-Diez et al. where
652!      zmla = "...first level where potential temperature exceeds the minimum potential temperature
653!        reached in the mixed layer by more than 1.5 K"
654
655    IMPLICIT NONE
656
657    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
658    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: tpot, qratio, z
659    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
660    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: zmla3D
661 
662! Local
663    INTEGER                                              :: i, j, it
664
665!!!!!!! Variables
666! tpot: potential air temperature [K]
667! qratio: water vapour mixing ratio [kgkg-1]
668! z: height above sea level [m]
669! hgt: terrain height [m]
670! zmla3D: boundary layer height from surface [m]
671
672    fname = 'compute_zmla_generic4D'
673
674    DO i=1, d1
675      DO j=1, d2
676        DO it=1, d4
677          CALL var_zmla_generic(d3, qratio(i,j,:,it), tpot(i,j,:,it), z(i,j,:,it), hgt(i,j),          &
678            zmla3D(i,j,it))
679        END DO
680      END DO
681    END DO
682
683    RETURN
684
685  END SUBROUTINE compute_zmla_generic4D
686
687  SUBROUTINE compute_zwind4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
688! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
689
690    IMPLICIT NONE
691
692    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
693    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ua, va, z
694    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: uas, vas
695    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
696    REAL(r_k), INTENT(in)                                :: zextrap
697    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: uaz, vaz
698 
699! Local
700    INTEGER                                              :: i, j, it
701
702!!!!!!! Variables
703! tpot: potential air temperature [K]
704! qratio: water vapour mixing ratio [kgkg-1]
705! z: height above surface [m]
706! sina, cosa: local sine and cosine of map rotation [1.]
707! zmla3D: boundary layer height from surface [m]
708
709    fname = 'compute_zwind4D'
710
711    DO i=1, d1
712      DO j=1, d2
713        DO it=1, d4
714          CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it),       &
715            sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
716        END DO
717      END DO
718    END DO
719
720    RETURN
721
722  END SUBROUTINE compute_zwind4D
723
724  SUBROUTINE compute_zwind_log4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
725! Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
726
727    IMPLICIT NONE
728
729    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
730    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ua, va, z
731    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: uas, vas
732    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
733    REAL(r_k), INTENT(in)                                :: zextrap
734    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: uaz, vaz
735 
736! Local
737    INTEGER                                              :: i, j, it
738
739!!!!!!! Variables
740! tpot: potential air temperature [K]
741! qratio: water vapour mixing ratio [kgkg-1]
742! z: height above surface [m]
743! sina, cosa: local sine and cosine of map rotation [1.]
744! zmla3D: boundary layer height from surface [m]
745
746    fname = 'compute_zwind_log4D'
747
748    DO i=1, d1
749      DO j=1, d2
750        DO it=1, d4
751          CALL var_zwind_log(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it),   &
752            sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
753        END DO
754      END DO
755    END DO
756
757    RETURN
758
759  END SUBROUTINE compute_zwind_log4D
760
761  SUBROUTINE compute_zwindMO3D(d1, d2, d3, ust, znt, rmol, uas, vas, sina, cosa, newz, uznew, vznew)
762! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
763!   NOTE: only usefull for newz < 80. m
764
765    IMPLICIT NONE
766
767    INTEGER, INTENT(in)                                  :: d1, d2, d3
768    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: ust, znt, rmol
769    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: uas, vas
770    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
771    REAL(r_k), INTENT(in)                                :: newz
772    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: uznew, vznew
773 
774! Local
775    INTEGER                                              :: i, j, it
776
777!!!!!!! Variables
778! ust: u* in similarity theory [ms-1]
779! znt: thermal time-varying roughness length [m]
780! rmol: Inverse of the Obukhov length [m-1]
781! uas: x-component 10-m wind speed [ms-1]
782! vas: y-component 10-m wind speed [ms-1]
783! sina, cosa: local sine and cosine of map rotation [1.]
784
785    fname = 'compute_zwindMO3D'
786
787    DO i=1, d1
788      DO j=1, d2
789        DO it=1, d3
790          CALL var_zwind_MOtheor(ust(i,j,it), znt(i,j,it), rmol(i,j,it), uas(i,j,it), vas(i,j,it),    &
791            sina(i,j), cosa(i,j), newz, uznew(i,j,it), vznew(i,j,it))
792        END DO
793      END DO
794    END DO
795
796    RETURN
797
798  END SUBROUTINE compute_zwindMO3D
799
800  SUBROUTINE compute_potevap_orPM3D(d1, d2, d3, rho1, ust, uas, vas, tas, ps, qv1, potevap)
801! Subroutine to compute potential evapotranspiration Penman-Monteith formulation implemented in
802!   ORCHIDEE in src_sechiba/enerbil.f90
803
804    IMPLICIT NONE
805
806    INTEGER, INTENT(in)                                  :: d1, d2, d3
807    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: rho1, ust, uas, vas, tas, ps, qv1
808    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: potevap
809 
810! Local
811    INTEGER                                              :: i, j, it
812
813!!!!!!! Variables
814! rho1: atsmophere density at the first layer [kgm-3]
815! ust: u* in similarity theory [ms-1]
816! uas: x-component 10-m wind speed [ms-1]
817! vas: y-component 10-m wind speed [ms-1]
818! tas: 2-m atmosphere temperature [K]
819! ps: surface pressure [Pa]
820! qv1: 1st layer atmospheric mixing ratio [kgkg-1]
821! potevap: potential evapo transpiration [kgm-2s-1]
822
823    fname = 'compute_potevap_orPM3D'
824
825    DO i=1, d1
826      DO j=1, d2
827        DO it=1, d3
828          CALL var_potevap_orPM(rho1(i,j,it), ust(i,j,it), uas(i,j,it), vas(i,j,it), tas(i,j,it),     &
829            ps(i,j,it), qv1(i,j,it), potevap(i,j,it))
830        END DO
831      END DO
832    END DO
833
834    RETURN
835
836  END SUBROUTINE compute_potevap_orPM3D
837
838  SUBROUTINE compute_fog_K84(d1, d2, d3, qc, qi, fog, vis)
839! Subroutine to compute fog: qcloud + qice /= 0.
840! And visibility following Kunkel, B. A., (1984): Parameterization of droplet terminal velocity and
841!   extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34–41.
842
843    IMPLICIT NONE
844
845    INTEGER, INTENT(in)                                  :: d1, d2, d3
846    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: qc, qi
847    INTEGER, DIMENSION(d1,d2,d3), INTENT(out)            :: fog
848    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: vis
849 
850! Local
851    INTEGER                                              :: i, j, it
852
853!!!!!!! Variables
854! qc: cloud mixing ratio [kgkg-1]
855! qi, ice mixing ratio [kgkg-1]
856! fog: presence of fog (1: yes, 0: no)
857! vis: visibility within fog [km]
858
859    fname = 'compute_fog_K84'
860
861    DO i=1, d1
862      DO j=1, d2
863        DO it=1, d3
864          CALL var_fog_K84(qc(i,j,it), qi(i,j,it), fog(i,j,it), vis(i,j,it))
865        END DO
866      END DO
867    END DO
868
869    RETURN
870
871  END SUBROUTINE compute_fog_K84
872
873  SUBROUTINE compute_fog_RUC(d1, d2, d3, qv, ta, pres, fog, vis)
874! Subroutine to compute fog: qcloud + qice /= 0.
875! And visibility following RUC method Smirnova, T. G., S. G. Benjamin, and J. M. Brown, 2000: Case
876!   study verification of RUC/MAPS fog and visibility forecasts. Preprints, 9 th Conference on
877!   Aviation, Range, and Aerospace Meteorlogy, AMS, Orlando, FL, Sep. 2000. Paper#2.3, 6 pp.
878
879    IMPLICIT NONE
880
881    INTEGER, INTENT(in)                                  :: d1, d2, d3
882    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: qv, ta, pres
883    INTEGER, DIMENSION(d1,d2,d3), INTENT(out)            :: fog
884    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: vis
885 
886! Local
887    INTEGER                                              :: i, j, it
888
889!!!!!!! Variables
890! qv: water vapor mixing ratio [kgkg-1]
891! ta: temperature [K]
892! pres: pressure [Pa]
893! fog: presence of fog (1: yes, 0: no)
894! vis: visibility within fog [km]
895
896    fname = 'compute_fog_RUC'
897
898    DO i=1, d1
899      DO j=1, d2
900        DO it=1, d3
901          CALL var_fog_RUC(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it))
902        END DO
903      END DO
904    END DO
905
906    RETURN
907
908  END SUBROUTINE compute_fog_RUC
909
910  SUBROUTINE compute_fog_FRAML50(d1, d2, d3, qv, ta, pres, fog, vis)
911! Subroutine to compute fog (vis < 1 km) and visibility following
912!   Gultepe, I. and J.A. Milbrandt, 2010: Probabilistic Parameterizations of Visibility Using
913!     Observations of Rain Precipitation Rate, Relative Humidity, and Visibility. J. Appl. Meteor.
914!     Climatol., 49, 36-46, https://doi.org/10.1175/2009JAMC1927.1
915!   Interest is focused on a 'general' fog/visibilty approach, thus the fit at 50 % of probability is
916!     chosen
917!   Effects from precipitation are not considered
918
919    IMPLICIT NONE
920
921    INTEGER, INTENT(in)                                  :: d1, d2, d3
922    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: qv, ta, pres
923    INTEGER, DIMENSION(d1,d2,d3), INTENT(out)            :: fog
924    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: vis
925 
926! Local
927    INTEGER                                              :: i, j, it
928
929!!!!!!! Variables
930! qv: mixing ratio in [kgkg-1]
931! ta: temperature [K]
932! pres: pressure field [Pa]
933! fog: presence of fog (1: yes, 0: no)
934! vis: visibility within fog [km]
935
936    fname = 'compute_fog_FRAML50'
937
938    DO i=1, d1
939      DO j=1, d2
940        DO it=1, d3
941          CALL var_fog_FRAML50(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it))
942        END DO
943      END DO
944    END DO
945
946    RETURN
947
948  END SUBROUTINE compute_fog_FRAML50
949
950END MODULE module_ForDiagnostics
Note: See TracBrowser for help on using the repository browser.