1module physicsFunctions_mod
2 ! MODULE physicsFunctions_mod (prefix='phf' category='8. Low-level utilities and constants')
3 !
4 !:Purpose: A collection of basic functions for various purposes
5 ! (e.g. computing saturation vapour pressure).
6 !
7 use MathPhysConstants_mod
8 use earthConstants_mod
9 use utilities_mod
10 use message_mod
11
12 implicit none
13 private
14
15 ! public procedures
16 public :: phf_FOEW8, phf_FODLE8, phf_FOQST8, phf_FODQS8, phf_FOEFQ8, phf_FOQFE8, phf_FOTVT8, phf_FOTTV8
17 public :: phf_FOHR8, phf_FOEWA8, phf_FODLA8, phf_FOQSA8, phf_FODQA8, phf_FOHRA8, phf_FOTW8, phf_FOTI8
18 public :: phf_FODTW8, phf_FODTI8, phf_FOTWI8, phf_FODTWI8, phf_FOEW8_CMAM, phf_FOEI8_CMAM, phf_FOERAT8_CMAM
19 public :: phf_FOEWI8_CMAM, phf_FODLE8_CMAM, phf_FOQST8_CMAM, phf_FOTW8_CMAM, phf_FOTI8_CMAM, phf_FODTW8_CMAM
20 public :: phf_FODTI8_CMAM, phf_FOTWI8_CMAM, phf_FODTWI8_CMAM, phf_FQBRANCH, phf_FOEFQL, phf_fotvvl, phf_FOEFQA
21 public :: phf_FOEFQPSA, phf_fottva, phf_folnqva
22 public :: phf_convert_z_to_pressure,phf_convert_z_to_gz
23 public :: phf_get_tropopause, phf_get_pbl, phf_calcDistance, phf_calcDistanceFast
24 public :: phf_height2geopotential, phf_gravityalt, phf_gravitysrf
25
26 logical :: phf_initialized = .false.
27
28 ! namelist variables:
29 character(len=20) :: saturationCurve ! saturationCurve must be one of 'Tetens_1930', 'Tetens_2018a', 'Tetens_2018'
30
31 contains
32
33 !--------------------------------------------------------------------------
34 ! tetens_coefs_switch
35 !--------------------------------------------------------------------------
36 subroutine phf_tetens_coefs_switch
37 !
38 !:Purpose: Read water saturation strategy from nml. Options are:
39 !
40 ! - 'Tetens_1930' was active before 2018.
41 ! - 'Tetens_2018a' is a partial update that missed some functions (active before IC4)
42 ! - 'Tetens_2018' completes the update to the intended 2018 specification.
43 !
44 implicit none
45
46 ! Locals:
47 integer :: NULNAM,IERR,FNOM,FCLOS
48 logical :: validOption
49 NAMELIST /NAMPHY/ saturationCurve
50
51 !$omp critical
52 if (.not.phf_initialized) then
53
54 saturationCurve = 'Tetens_2018'
55
56 if ( .not. utl_isNamelistPresent('NAMPHY','./flnml') ) then
57 call msg( 'phf_tetens_coefs_switch', &
58 'NAMPHY is missing in the namelist. Default values will be taken.', mpiAll_opt=.False.)
59 else
60 ! Reading the namelist
61 nulnam = 0
62 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
63 read(nulnam, nml=namphy, iostat=ierr)
64 if ( ierr /= 0) call utl_abort('tetens_coefs: Error reading namelist')
65 ierr = fclos(nulnam)
66 end if
67
68 validOption = (trim(saturationCurve) == 'Tetens_1930' .or. &
69 trim(saturationCurve) == 'Tetens_2018a' .or. &
70 trim(saturationCurve) == 'Tetens_2018')
71 if (.not.validOption) then
72 call utl_abort('phf_tetens: WV Saturation not in expected list')
73 end if
74
75 call msg( 'phf_tetens_coefs_switch ', saturationCurve )
76 phf_initialized = .true.
77 end if
78 !$omp end critical
79
80 end subroutine phf_tetens_coefs_switch
81
82 !--------------------------------------------------------------------------
83 ! phf_FOEW8
84 !--------------------------------------------------------------------------
85 real(8) function phf_FOEW8(TTT)
86 !
87 !:Purpose: Water vapour saturation pressure (Tetens) - EW or EI as a fct of TT
88 !
89 implicit none
90
91 ! Arguments:
92 real(8), intent(in) :: TTT
93
94 if (.not.phf_initialized) call phf_tetens_coefs_switch
95
96 if (trim(saturationCurve) == 'Tetens_2018') then
97 ! Updated coefficients 2018
98 if (TTT > MPC_K_C_DEGREE_OFFSET_R8) then
99 phf_FOEW8 = 610.94D0*EXP(17.625D0 * (TTT-MPC_TRIPLE_POINT_R8) / (TTT-30.11D0))
100 else
101 phf_FOEW8 = 610.94D0*EXP(22.587D0 * (TTT-MPC_TRIPLE_POINT_R8) / (TTT+ 0.71D0))
102 end if
103 else
104 ! Classic Tetens 1930 coefficients
105 if (TTT > MPC_TRIPLE_POINT_R8) then
106 phf_FOEW8 = 610.78D0*EXP(17.269D0 * (TTT-MPC_TRIPLE_POINT_R8) / (TTT-35.86D0))
107 else
108 phf_FOEW8 = 610.78D0*EXP(21.875D0 * (TTT-MPC_TRIPLE_POINT_R8) / (TTT- 7.66D0))
109 end if
110 end if
111 end function phf_foew8
112
113 !--------------------------------------------------------------------------
114 ! phf_FODLE8
115 !--------------------------------------------------------------------------
116 real(8) function phf_FODLE8(TTT)
117 !
118 !:Purpose: FONCTION CALCULANT LA DERIVEE SELON T DE LN EW (OU LN EI)
119 !
120 implicit none
121
122 ! Arguments:
123 real(8), intent(in) :: TTT
124
125 if (.not.phf_initialized) call phf_tetens_coefs_switch
126
127 if (trim(saturationCurve) == 'Tetens_2018') then
128 if (TTT > MPC_TRIPLE_POINT_R8) then
129 phf_FODLE8 = 4283.76D0/(TTT-30.11D0)**2
130 else
131 phf_FODLE8 = 6185.90D0/(TTT+ 0.71D0)**2
132 end if
133 else
134 if (TTT > MPC_TRIPLE_POINT_R8) then
135 phf_FODLE8 = 4097.93D0/(TTT-35.86D0)**2
136 else
137 phf_FODLE8 = 5807.81D0/(TTT- 7.66D0)**2
138 end if
139 end if
140 end function phf_FODLE8
141
142 !--------------------------------------------------------------------------
143 ! phf_FOQST8
144 !--------------------------------------------------------------------------
145 real(8) function phf_FOQST8(TTT,PRS)
146 !
147 !:Purpose: FONCTION CALCULANT L'HUMIDITE SPECIFIQUE SATURANTE (QSAT)
148 !
149 implicit none
150
151 ! Arguments:
152 real(8), intent(in) :: TTT
153 real(8), intent(in) :: PRS
154
155 phf_FOQST8=MPC_EPS1_R8/(MAX(1.D0,PRS/phf_FOEW8(TTT))-MPC_EPS2_R8)
156 end function phf_FOQST8
157
158 !--------------------------------------------------------------------------
159 ! phf_FODQS8
160 !--------------------------------------------------------------------------
161 real(8) function phf_FODQS8(QST,TTT)
162 !
163 !:Purpose: FONCTION CALCULANT LA DERIVEE DE QSAT SELON T
164 !
165 implicit none
166
167 ! Arguments:
168 real(8), intent(in) :: TTT
169 real(8), intent(in) :: QST
170
171 phf_FODQS8=QST*(1.D0+MPC_DELTA_R8*QST)*phf_FODLE8(TTT)
172 end function phf_FODQS8
173 ! QST EST LA SORTIE DE FOQST
174
175 !--------------------------------------------------------------------------
176 ! phf_FOEFQ8
177 !--------------------------------------------------------------------------
178 real(8) function phf_FOEFQ8(QQQ,PRS)
179 !
180 !:Purpose: FONCTION CALCULANT TENSION VAP (EEE) FN DE HUM SP (QQQ) ET PRS
181 !
182 implicit none
183
184 ! Arguments:
185 real(8), intent(in) :: QQQ
186 real(8), intent(in) :: PRS
187
188 phf_FOEFQ8= MIN(PRS,(QQQ*PRS) / (MPC_EPS1_R8 + MPC_EPS2_R8*QQQ))
189 end function phf_FOEFQ8
190
191 !--------------------------------------------------------------------------
192 ! FOQFE8
193 !--------------------------------------------------------------------------
194 real(8) function phf_FOQFE8(EEE,PRS)
195 !
196 !:Purpose: FONCTION CALCULANT HUM SP (QQQ) DE TENS. VAP (EEE) ET PRES (PRS)
197 !
198 implicit none
199
200 ! Arguments:
201 real(8), intent(in) :: EEE
202 real(8), intent(in) :: PRS
203
204 phf_FOQFE8= MIN(1.D0,MPC_EPS1_R8*EEE / (PRS-MPC_EPS2_R8*EEE))
205 end function phf_FOQFE8
206
207 !--------------------------------------------------------------------------
208 ! phf_FOTVT8
209 !--------------------------------------------------------------------------
210 real(8) function phf_FOTVT8(TTT,QQQ)
211 !
212 !:Purpose: FONCTION CALCULANT TEMP VIRT. (TVI) DE TEMP (TTT) ET HUM SP (QQQ)
213 !
214 implicit none
215
216 ! Arguments:
217 real(8), intent(in) :: TTT
218 real(8), intent(in) :: QQQ
219
220 phf_FOTVT8= TTT * (1.0D0 + MPC_DELTA_R8*QQQ)
221 end function phf_FOTVT8
222
223 !--------------------------------------------------------------------------
224 ! phf_FOTTV8
225 !--------------------------------------------------------------------------
226 real(8) function phf_FOTTV8(TVI,QQQ)
227 !
228 !:Purpose: FONCTION CALCULANT TTT DE TEMP VIRT. (TVI) ET HUM SP (QQQ)
229 !
230 implicit none
231
232 ! Arguments:
233 real(8), intent(in) :: TVI
234 real(8), intent(in) :: QQQ
235
236 phf_FOTTV8= TVI / (1.0D0 + MPC_DELTA_R8*QQQ)
237 end function phf_FOTTV8
238
239 !--------------------------------------------------------------------------
240 ! phf_FOHR8
241 !--------------------------------------------------------------------------
242 real(8) function phf_FOHR8(QQQ,TTT,PRS)
243 !
244 !:Purpose: FONCTION CALCULANT HUM REL DE HUM SP (QQQ), TEMP (TTT) ET PRES (PRS).
245 ! Where HR = E/ESAT
246 !
247 implicit none
248
249 ! Arguments:
250 real(8), intent(in) :: QQQ
251 real(8), intent(in) :: TTT
252 real(8), intent(in) :: PRS
253
254 phf_FOHR8 = MIN(PRS,phf_FOEFQ8(QQQ,PRS)) / phf_FOEW8(TTT)
255 end function phf_FOHR8
256
257 ! LES 5 FONCTIONS SUIVANTES SONT VALIDES DANS LE CONTEXTE OU ON
258 ! NE DESIRE PAS TENIR COMPTE DE LA PHASE GLACE DANS LES CALCULS
259 ! DE SATURATION.
260
261 !--------------------------------------------------------------------------
262 ! phf_FOEWA8
263 !--------------------------------------------------------------------------
264 real(8) function phf_FOEWA8(TTT)
265 !
266 !:Purpose: FONCTION DE VAPEUR SATURANTE (TETENS)
267 !
268 implicit none
269
270 ! Arguments:
271 real(8), intent(in) :: TTT
272
273 if (.not.phf_initialized) call phf_tetens_coefs_switch
274
275 if (trim(saturationCurve) == 'Tetens_2018') then
276 phf_FOEWA8=610.94D0*EXP(17.625D0*(TTT-MPC_TRIPLE_POINT_R8)/(TTT-30.11D0))
277 else
278 phf_FOEWA8=610.78D0*EXP(17.269D0*(TTT-MPC_TRIPLE_POINT_R8)/(TTT-35.86D0))
279 end if
280 end function phf_FOEWA8
281
282 !--------------------------------------------------------------------------
283 ! phf_FODLA8
284 !--------------------------------------------------------------------------
285 real(8) function phf_FODLA8(TTT)
286 !
287 !:Purpose: FONCTION CALCULANT LA DERIVEE SELON T DE LN EW
288 !
289 implicit none
290
291 ! Arguments:
292 real(8), intent(in) :: TTT
293
294 if (.not.phf_initialized) call phf_tetens_coefs_switch
295
296 if (trim(saturationCurve) == 'Tetens_2018') then
297 phf_FODLA8=17.625D0*(MPC_TRIPLE_POINT_R8-30.11D0)/(TTT-30.11D0)**2
298 else
299 phf_FODLA8=17.269D0*(MPC_TRIPLE_POINT_R8-35.86D0)/(TTT-35.86D0)**2
300 end if
301 end function phf_FODLA8
302
303 !--------------------------------------------------------------------------
304 ! FOQSA8
305 !--------------------------------------------------------------------------
306 real(8) function phf_FOQSA8(TTT,PRS)
307 !
308 !:Purpose: FONCTION CALCULANT L'HUMIDITE SPECIFIQUE SATURANTE
309 !
310 implicit none
311
312 ! Arguments:
313 real(8), intent(in) :: TTT
314 real(8), intent(in) :: PRS
315
316 phf_FOQSA8=MPC_EPS1_R8/(MAX(1.D0,PRS/phf_FOEWA8(TTT))-MPC_EPS2_R8)
317 end function phf_FOQSA8
318
319 !--------------------------------------------------------------------------
320 ! phf_FODQA8
321 !--------------------------------------------------------------------------
322 real(8) function phf_FODQA8(QST,TTT)
323 !
324 !:Purpose: FONCTION CALCULANT LA DERIVEE DE QSAT SELON T
325 !
326 implicit none
327
328 ! Arguments:
329 real(8), intent(in) :: QST
330 real(8), intent(in) :: TTT
331
332 phf_FODQA8=QST*(1.D0+MPC_DELTA_R8*QST)*phf_FODLA8(TTT)
333 end function phf_FODQA8
334
335 !--------------------------------------------------------------------------
336 ! phf_FOHRA8
337 !--------------------------------------------------------------------------
338 real(8) function phf_FOHRA8(QQQ,TTT,PRS)
339 !
340 !:Purpose: FONCTION CALCULANT L'HUMIDITE RELATIVE
341 !
342 implicit none
343
344 ! Arguments:
345 real(8), intent(in) :: QQQ
346 real(8), intent(in) :: TTT
347 real(8), intent(in) :: PRS
348
349 phf_FOHRA8=MIN(PRS,phf_FOEFQ8(QQQ,PRS))/phf_FOEWA8(TTT)
350 end function phf_FOHRA8
351
352 ! LES 6 FONCTIONS SUIVANTES SONT REQUISES POUR LA TEMPERATURE
353 ! EN FONCTION DE LA TENSION DE VAPEUR SATURANTE.
354 ! (AJOUTE PAR YVES J. ROCHON, ARQX/SMC, JUIN 2004)
355
356 !--------------------------------------------------------------------------
357 ! phf_FOTW8
358 !--------------------------------------------------------------------------
359 real(8) function phf_FOTW8(EEE)
360 !
361 !:Purpose: FONCTION DE LA TEMPERATURE EN FONCTION DE LA TENSION DE VAPEUR
362 ! SATURANTE PAR RAPPORT A EW.
363 !
364 implicit none
365
366 ! Arguments:
367 real(8), intent(in) :: EEE
368
369 if (.not.phf_initialized) call phf_tetens_coefs_switch
370
371 if (trim(saturationCurve) == 'Tetens_2018a' .or. trim(saturationCurve) == 'Tetens_2018') then
372 phf_FOTW8=(30.11D0*LOG(EEE/610.94D0)-17.625D0*MPC_TRIPLE_POINT_R8)/ &
373 (LOG(EEE/610.94D0)-17.625D0)
374 else
375 phf_FOTW8=(35.86D0*LOG(EEE/610.78D0)-17.269D0*MPC_TRIPLE_POINT_R8)/ &
376 (LOG(EEE/610.78D0)-17.269D0)
377 end if
378 end function phf_FOTW8
379
380 !--------------------------------------------------------------------------
381 ! phf_FOTI8
382 !--------------------------------------------------------------------------
383 real(8) function phf_FOTI8(EEE)
384 !
385 !:Purpose: FONCTION DE LA TEMPERATURE EN FONCTION DE LA TENSION DE VAPEUR
386 ! SATURANTE PAR RAPPORT A EI.
387 !
388 implicit none
389
390 ! Arguments:
391 real(8), intent(in) :: EEE
392
393 if (.not.phf_initialized) call phf_tetens_coefs_switch
394
395 if (trim(saturationCurve) == 'Tetens_2018') then
396 phf_FOTI8=(-0.71D0*LOG(EEE/610.94D0)-22.587D0*MPC_TRIPLE_POINT_R8)/ &
397 (LOG(EEE/610.94D0)-22.587D0)
398 else
399 phf_FOTI8=( 7.66D0*LOG(EEE/610.78D0)-21.875D0*MPC_TRIPLE_POINT_R8)/ &
400 (LOG(EEE/610.78D0)-21.875D0)
401 end if
402 end function phf_FOTI8
403
404 !--------------------------------------------------------------------------
405 ! phf_FODTH8
406 !--------------------------------------------------------------------------
407 real(8) function phf_FODTW8(TTT,EEE)
408 !
409 !:Purpose: FONCTION DE LA DERIVE DE LA TEMPERATURE EN FONCTION DE LA TENSION DE
410 ! VAPEUR SATURANTE (EW).
411 !
412 implicit none
413
414 ! Arguments:
415 real(8), intent(in) :: TTT
416 real(8), intent(in) :: EEE
417
418 if (.not.phf_initialized) call phf_tetens_coefs_switch
419
420 if (trim(saturationCurve) == 'Tetens_2018a' .or. trim(saturationCurve) == 'Tetens_2018') then
421 phf_FODTW8=(30.11D0-TTT)/EEE/(LOG(EEE/610.94D0)-17.625D0)
422 else
423 phf_FODTW8=(35.86D0-TTT)/EEE/(LOG(EEE/610.78D0)-17.269D0)
424 endif
425 end function phf_FODTW8
426
427 !--------------------------------------------------------------------------
428 ! phf_FODTI8
429 !--------------------------------------------------------------------------
430 real(8) function phf_FODTI8(TTT,EEE)
431 !
432 !:Purpose: FONCTION DE LA DERIVE DE LA TEMPERATURE EN FONCTION DE LA TENSION DE
433 ! VAPEUR SATURANTE (EI).
434 !
435 implicit none
436
437 ! Arguments:
438 real(8), intent(in) :: TTT
439 real(8), intent(in) :: EEE
440
441 if (.not.phf_initialized) call phf_tetens_coefs_switch
442
443 if (trim(saturationCurve) == 'Tetens_2018') then
444 phf_FODTI8=(-0.71D0-TTT)/EEE / (LOG(EEE/610.94D0)-22.587D0)
445 else
446 phf_FODTI8=( 7.66D0-TTT)/EEE / (LOG(EEE/610.78D0)-21.875D0)
447 end if
448 end function phf_FODTI8
449
450 !--------------------------------------------------------------------------
451 ! phf_FOTWI8
452 !--------------------------------------------------------------------------
453 real(8) function phf_FOTWI8(TTT,EEE)
454 !
455 !:Purpose: FONCTION DE L'AJUSTEMENT DE LA TEMPERATURE.
456 !
457 implicit none
458
459 ! Arguments:
460 real(8), intent(in) :: TTT
461 real(8), intent(in) :: EEE
462
463 phf_FOTWI8=MAX(0.0D0,SIGN(1.0D0,TTT-MPC_TRIPLE_POINT_R8))*phf_FOTW8(EEE) &
464 -MIN(0.0D0,SIGN(1.0D0,TTT-MPC_TRIPLE_POINT_R8))*phf_FOTI8(EEE)
465 end function phf_FOTWI8
466
467 !--------------------------------------------------------------------------
468 ! phf_FODWTI8
469 !--------------------------------------------------------------------------
470 real(8) function phf_FODTWI8(TTT,EEE)
471 !
472 !:Purpose: FONCTION DE L'AJUSTEMENT DE LA DERIVEE DE LA TEMPERATURE.
473 !
474 implicit none
475
476 ! Arguments:
477 real(8), intent(in) :: TTT
478 real(8), intent(in) :: EEE
479
480 phf_FODTWI8=MAX(0.0D0,SIGN(1.0D0,TTT-MPC_TRIPLE_POINT_R8))*phf_FODTW8(TTT,EEE) &
481 -MIN(0.0D0,SIGN(1.0D0,TTT-MPC_TRIPLE_POINT_R8))*phf_FODTI8(TTT,EEE)
482 end function phf_FODTWI8
483
484 ! LES 7 FONCTIONS SUIVANTES POUR EW-EI SONT REQUISES POUR LES MODELES
485 ! CMAM ET CGCM. (AJOUTE PAR YVES J. ROCHON, ARQX/SMC, JUIN 2004)
486
487 !--------------------------------------------------------------------------
488 ! phf_FOEW8_CMAM
489 !--------------------------------------------------------------------------
490 real(8) function phf_FOEW8_CMAM(TTT)
491 !
492 !:Purpose: FONCTION DE TENSION DE VAPEUR SATURANTE - EW
493 !
494 implicit none
495
496 ! Arguments:
497 real(8), intent(in) :: TTT
498
499 phf_FOEW8_CMAM= 100.D0*EXP(21.656D0-5418.D0/TTT)
500 end function phf_FOEW8_CMAM
501
502 !--------------------------------------------------------------------------
503 ! phf_FOEI8_CMAM
504 !--------------------------------------------------------------------------
505 real(8) function phf_FOEI8_CMAM(TTT)
506 !
507 !:Purpose: FONCTION DE TENSION DE VAPEUR SATURANTE - EI
508 !
509 implicit none
510
511 ! Arguments:
512 real(8), intent(in) :: TTT
513
514 phf_FOEI8_CMAM= 100.D0*EXP(24.292D0-6141.D0/TTT)
515 end function phf_FOEI8_CMAM
516
517 !--------------------------------------------------------------------------
518 ! phf_FOERAT8_CMAM
519 !--------------------------------------------------------------------------
520 real(8) function phf_FOERAT8_CMAM(TTT)
521 !
522 !:Purpose: FONCTION DE LA PROPORTION DE LA CONTRIBUTION DE EW VS EI
523 !
524 implicit none
525
526 ! Arguments:
527 real(8), intent(in) :: TTT
528
529 phf_FOERAT8_CMAM=MIN(1.0D0,0.0059D0+0.9941D0*EXP(-0.003102D0* &
530 MIN(0.0D0,TTT-MPC_TRIPLE_POINT_R8)**2))
531 end function phf_FOERAT8_CMAM
532
533 !--------------------------------------------------------------------------
534 ! phf_FOEWI8_CMAM
535 !--------------------------------------------------------------------------
536 real(8) function phf_FOEWI8_CMAM(TTT)
537 !
538 !:Purpose: FONCTION DE TENSION DE VAPEUR SATURANTE RESULTANTE - EW et EI
539 !
540 implicit none
541
542 ! Arguments:
543 real(8), intent(in) :: TTT
544
545 phf_FOEWI8_CMAM = phf_FOEW8_CMAM(TTT)*phf_FOERAT8_CMAM(TTT) &
546 +(1.0D0-phf_FOERAT8_CMAM(TTT))*phf_FOEI8_CMAM(TTT)
547 end function phf_FOEWI8_CMAM
548
549 !--------------------------------------------------------------------------
550 ! phf_FODLE8_CMAM
551 !--------------------------------------------------------------------------
552 real(8) function phf_FODLE8_CMAM(TTT)
553 !
554 !:Purpose: FONCTION DE LA DERIVE DE LN(E) PAR RAPPORT A LA TEMPERATURE
555 !
556 implicit none
557
558 ! Arguments:
559 real(8), intent(in) :: TTT
560
561 phf_FODLE8_CMAM= phf_FOERAT8_CMAM(TTT)*5418.D0/TTT/TTT &
562 +(1.0D0-phf_FOERAT8_CMAM(TTT))*6141.D0/TTT/TTT
563 end function phf_FODLE8_CMAM
564
565 !--------------------------------------------------------------------------
566 ! phf_FOQST8_CMAM
567 !--------------------------------------------------------------------------
568 real(8) function phf_FOQST8_CMAM(TTT,PRS)
569 !
570 !:Purpose: FONCTION CALCULANT L'HUMIDITE SPECIFIQUE SATURANTE (QSAT).
571 ! PREND EN COMPTE LES PHASES GLACE ET EAU.
572 !
573 implicit none
574
575 ! Arguments:
576 real(8), intent(in) :: TTT
577 real(8), intent(in) :: PRS
578
579 phf_FOQST8_CMAM=MPC_EPS1_R8/(MAX(1.0D0,PRS/ &
580 phf_FOEWI8_CMAM(TTT))-MPC_EPS2_R8)
581 end function phf_FOQST8_CMAM
582
583 ! LES 6 FONCTIONS SUIVANTES SONT REQUISES POUR LA TEMPERATURE
584 ! EN FONCTION DE LA TENSION DE VAPEUR SATURANTE POUR CMAM/CGCM
585 ! (AJOUTE PAR YVES J. ROCHON, ARQX/SMC, JUIN 2004)
586
587 !--------------------------------------------------------------------------
588 ! phf_FOTW8_CMAM
589 !--------------------------------------------------------------------------
590 real(8) function phf_FOTW8_CMAM(EEE)
591 !
592 !:Purpose: FONCTION DE LA TEMPERATURE EN FONCTION DE LA TENSION DE VAPEUR
593 ! SATURANTE PAR RAPPORT A EW.
594 !
595 implicit none
596
597 ! Arguments:
598 real(8) :: EEE
599
600 phf_FOTW8_CMAM=5418.D0/(21.656D0-LOG(EEE/100.0D0))
601 end function phf_FOTW8_CMAM
602
603 !--------------------------------------------------------------------------
604 ! phf_FOTI8_CMAM
605 !--------------------------------------------------------------------------
606 real(8) function phf_FOTI8_CMAM(EEE)
607 !
608 !:Purpose: FONCTION DE LA TEMPERATURE EN FONCTION DE LA TENSION DE VAPEUR
609 ! SATURANTE PAR RAPPORT A EI.
610 !
611 implicit none
612
613 ! Arguments:
614 real(8), intent(in) :: EEE
615
616 phf_FOTI8_CMAM=6141.D0/(24.292D0-LOG(EEE/100.0D0))
617 end function phf_FOTI8_CMAM
618
619 !--------------------------------------------------------------------------
620 ! phf_FODTW8_CMAM
621 !--------------------------------------------------------------------------
622 real(8) function phf_FODTW8_CMAM(TTT,EEE)
623 !
624 !:Purpose: FONCTION DE LA DERIVE DE LA TEMPERATURE EN FONCTION DE LA TENSION DE
625 ! VAPEUR SATURANTE (EW).
626 !
627 implicit none
628
629 ! Arguments:
630 real(8), intent(in) :: TTT
631 real(8), intent(in) :: EEE
632
633 phf_FODTW8_CMAM=TTT/EEE/(21.656D0-LOG(EEE/100.0D0))
634 end function phf_FODTW8_CMAM
635
636 !--------------------------------------------------------------------------
637 ! phf_FODTI8_CMAM
638 !--------------------------------------------------------------------------
639 real(8) function phf_FODTI8_CMAM(TTT,EEE)
640 !
641 !:Purpose: FONCTION DE LA DERIVE DE LA TEMPERATURE EN FONCTION DE LA TENSION DE
642 ! VAPEUR SATURANTE (EI).
643 !
644 implicit none
645
646 ! Arguments:
647 real(8), intent(in) :: TTT
648 real(8), intent(in) :: EEE
649
650 phf_FODTI8_CMAM=TTT/EEE/(24.292D0-LOG(EEE/100.0D0))
651 end function phf_FODTI8_CMAM
652
653 !--------------------------------------------------------------------------
654 ! phf_FOTWI8_CMAM
655 !--------------------------------------------------------------------------
656 real(8) function phf_FOTWI8_CMAM(TTT,EEE)
657 !
658 !:Purpose: FONCTION DE L'AJUSTEMENT DE LA TEMPERATURE.
659 !
660 implicit none
661
662 ! Arguments:
663 real(8), intent(in) :: TTT
664 real(8), intent(in) :: EEE
665
666 phf_FOTWI8_CMAM= phf_FOERAT8_CMAM(TTT)*phf_FOTW8_CMAM(EEE)+ &
667 (1.0D0-phf_FOERAT8_CMAM(TTT))*phf_FOTI8_CMAM(EEE)
668 end function phf_FOTWI8_CMAM
669
670 !--------------------------------------------------------------------------
671 ! phf_FODTWI8_CMAM
672 !--------------------------------------------------------------------------
673 real(8) function phf_FODTWI8_CMAM(TTT,EEE)
674 !
675 !:Purpose: FONCTION DE L'AJUSTEMENT DE LA DERIVEE DE LA TEMPERATURE.
676 !
677 implicit none
678
679 ! Arguments:
680 real(8), intent(in) :: TTT
681 real(8), intent(in) :: EEE
682
683 phf_FODTWI8_CMAM= phf_FOERAT8_CMAM(TTT)*phf_FODTW8_CMAM(TTT,EEE)+ &
684 (1.0D0-phf_FOERAT8_CMAM(TTT))*phf_FODTI8_CMAM(TTT,EEE)
685 end function phf_FODTWI8_CMAM
686
687
688 !--------------------------------------------------------------------------
689 ! phf_FOBRANCH
690 !--------------------------------------------------------------------------
691 real(8) function phf_FQBRANCH(QQQ)
692 !
693 !:Purpose: function returning 0/1 depending on the minimum q branch condition
694 ! as discussed by Brunet (1996) to prevent getting a vapour pressure that exceeds
695 ! the total pressure p when q exceeds 1.
696 !
697 implicit none
698
699 ! Arguments:
700 real(8), intent(in) :: QQQ
701
702 phf_FQBRANCH= 0.5D0+SIGN(0.5D0,1.D0-(QQQ))
703 end function phf_FQBRANCH
704
705 !=============================================================================
706 !
707 ! TLM of THERMODYNAMIC FUNCTIONS USED IN 3DVAR
708 ! CONSTANTS FROM COMMON /CTESDYN/
709 ! NOTE: ALL FUNCTIONS WORK IN S.I. UNITS
710 ! I.E PRS IN PA, QQQ IN KG/KG
711
712 !--------------------------------------------------------------------------
713 ! phf_FOEFQL
714 !--------------------------------------------------------------------------
715 real(8) function phf_FOEFQL(QQL,PRSL,QQQ,PRS,PNETA)
716 !
717 !:Purpose: TLM OF FUNCTION CALCULATING VAPOUR PRESSURE
718 !
719 ! INPUTS:
720 !
721 ! * QQL , PERTURBATION OF LN SPECIFIC HUM
722 ! * PRSL , PERTURBATION OF SURFACE PRESSURE
723 ! * QQQ , SPECIFIC HUMIDITY
724 ! * PRS , PRESSURE
725 ! * PNETA , VALUE OF ETA LEVEL
726 !
727 ! OUTPUT:
728 !
729 ! * FOEFQL, PERTURBATION OF VAPOUR PRESSURE
730 !
731 implicit none
732
733 ! Arguments:
734 real(8), intent(in) :: QQL
735 real(8), intent(in) :: PRSL
736 real(8), intent(in) :: QQQ
737 real(8), intent(in) :: PRS
738 real(8), intent(in) :: PNETA
739
740 phf_FOEFQL= phf_FQBRANCH(QQQ) &
741 * ((QQL*MPC_EPS1_R8*PRS*QQQ/(MPC_EPS1_R8+MPC_EPS2_R8*QQQ) &
742 + PRSL*PNETA*QQQ)/(MPC_EPS1_R8+MPC_EPS2_R8*QQQ)) &
743 + (1.0D0 - phf_FQBRANCH(QQQ))*PRSL*PNETA
744 end function phf_FOEFQL
745
746 !--------------------------------------------------------------------------
747 ! phf_FOTVVL
748 !--------------------------------------------------------------------------
749 real(8) function phf_fotvvl(qqq,ttt,ttl,plnql)
750 !
751 !:Purpose: Tangent-linear operator of virtual temperature
752 !
753 implicit none
754
755 ! Arguments:
756 real(8), intent(in) :: qqq ! backgroud specific humidity
757 real(8), intent(in) :: ttt ! backgroud temperature
758 real(8), intent(in) :: ttl ! temperature increment
759 real(8), intent(in) :: plnql ! increment of logarithm specific humidity (del(ln q))
760
761 phf_fotvvl=(1 + MPC_DELTA_R8*qqq)*ttl + MPC_DELTA_R8*qqq*ttt*plnql
762 end function phf_fotvvl
763
764 !=============================================================================
765 !
766 ! DEFINITION OF ADJOINTS OF THERMODYNAMIC FUNCTIONS
767 ! CONSTANTS AS IN COMMON /CTESDYN/
768 ! NOTE: ALL UNITS S.I.
769 ! I.E. PADES IN DEG K, PRS EN PA, QQQ EN KG/KG
770 !
771
772 !--------------------------------------------------------------------------
773 ! phf_FOEFQA
774 !--------------------------------------------------------------------------
775 real(8) function phf_FOEFQA(PADES,PGAMMA,QQQ,PRS)
776 !
777 !:Purpose: ADJOINT OF LN SPECIFIC HUM (QQQ) DUE TO DEWPOINT DEPRESSION CORRECTIONS
778 !
779 implicit none
780
781 ! Arguments:
782 real(8), intent(in) :: PADES ! ADJOINT OF DEWPOINT DEPRESSION
783 real(8), intent(in) :: PGAMMA ! ADOINT OF VAPOUR PRESSURE RELATIONSHIP
784 real(8), intent(in) :: QQQ ! SPECIFIC HUMIDITY
785 real(8), intent(in) :: PRS ! PRESSURE
786
787 phf_FOEFQA= PADES*PGAMMA*MPC_EPS1_R8*PRS*QQQ/ &
788 ((MPC_EPS1_R8+MPC_EPS2_R8*QQQ)*(MPC_EPS1_R8+MPC_EPS2_R8*QQQ))
789 end function phf_FOEFQA
790
791 !--------------------------------------------------------------------------
792 ! phf_FOEFQPSA
793 !--------------------------------------------------------------------------
794 real(8) function phf_FOEFQPSA(PADES,PGAMMA,QQQ,PNETA)
795 !
796 !:Purpose: ADJOINT OF SURFACE PRESSURE DUE TO DEWPOINT DEPRESSION CORRECTIONS
797 !
798 implicit none
799
800 ! Arguments:
801 real(8), intent(in) :: PADES ! ADJOINT OF DEWPOINT DEPRESSION
802 real(8), intent(in) :: PGAMMA ! ADOINT OF VAPOUR PRESSURE RELATIONSHIP
803 real(8), intent(in) :: QQQ ! SPECIFIC HUMIDITY
804 real(8), intent(in) :: PNETA ! VALUE OF NETA
805
806 phf_FOEFQPSA = PADES*PGAMMA*QQQ*PNETA/ &
807 (MPC_EPS1_R8+MPC_EPS2_R8*QQQ)
808 end function phf_FOEFQPSA
809
810 !--------------------------------------------------------------------------
811 ! phf_FOTTVa
812 !--------------------------------------------------------------------------
813 real(8) function phf_fottva(qqq,tva)
814 !
815 !:Purpose: Adjoint of temperature due to virtual temperature correction
816 !
817 implicit none
818
819 ! Arguments:
820 real(8), intent(in) :: qqq ! background specific humidity
821 real(8), intent(in) :: tva ! adjoint variable of virtual temperature
822
823 phf_fottva= (1D0 + MPC_DELTA_R8*qqq)*tva
824 end function phf_fottva
825
826 !--------------------------------------------------------------------------
827 ! phf_FOLNQVA
828 !--------------------------------------------------------------------------
829 real(8) function phf_folnqva(qqq,ttt,tva)
830 !
831 !:Purpose: Adjoint of logarithm of specific humidity due to virtual temperature correction
832 !
833 implicit none
834
835 ! Arguments:
836 real(8), intent(in) :: qqq ! background specific humidity
837 real(8), intent(in) :: ttt ! background temperature
838 real(8), intent(in) :: tva ! adjoint variable of virtual temperature
839
840 phf_folnqva = MPC_DELTA_R8*qqq*ttt*tva
841 end function phf_folnqva
842
843 !-------------------------------------------------------------------------------------------
844
845 !--------------------------------------------------------------------------
846 ! PHF_CONVERT_Z_TO_PRESSURE
847 !--------------------------------------------------------------------------
848 function phf_convert_z_to_pressure(altitude,rgz_mod,press_mod,nlev,nlev_mod,lat,success) result(press)
849 !
850 !:Purpose: Converts an array of (geometric) altitudes to pressures. Uses linear interpolation
851 ! in log(p).
852 implicit none
853
854 ! Arguments:
855 real(8), intent(in) :: altitude(nlev) ! altitudes to convert to pressures (m)
856 real(8), intent(in) :: rgz_mod(nlev_mod) ! geopotential heights on model levels (m), assumed to be in decending order
857 real(8), intent(in) :: press_mod(nlev_mod) ! pressure on model levels, assumed to be in ascending order
858 real(8), intent(in) :: lat ! latitude (rad)
859 integer, intent(in) :: nlev ! length of altitude array
860 integer, intent(in) :: nlev_mod ! number of model levels
861 logical, intent(inout) :: success(nlev) ! flag indicating success
862 ! Result:
863 real(8) :: press(nlev) ! converted pressures
864
865 ! Locals:
866 real(8) :: rgz(nlev)
867 integer :: ilev,ilev_mod
868
869 ! Check model pressures
870 if (any(press_mod.lt.0.0D0)) call utl_abort("phf_compute_z_to_pressure: Invalid model pressure.")
871
872 ! Convert altitudes to geopotential heights
873 rgz = phf_convert_z_to_gz(altitude,lat,nlev)
874
875 do ilev=1,nlev
876
877 ! Check if height is above or below model boundaries
878 if ( rgz(ilev).gt.rgz_mod(1) .or. rgz(ilev).lt.rgz_mod(nlev_mod) ) then
879 success(ilev)=.false.
880 end if
881
882 if (success(ilev)) then
883
884 ! Find model layers directly above and below rgz(ilev).
885 ! After exit of loop we will have
886 ! rgz_mod(ilev_mod) >= rgz(ilev) > rgz_mod(ilev_mod+1)
887 do ilev_mod=1,nlev_mod-1
888 if ( rgz(ilev).le.rgz_mod(ilev_mod) .and. &
889 rgz(ilev).gt.rgz_mod(ilev_mod+1) ) exit
890 end do
891
892 ! Linear interpolation in gz,log(p)
893 press(ilev) = press_mod(ilev_mod+1) * (press_mod(ilev_mod)/press_mod(ilev_mod+1))**( &
894 (rgz(ilev)-rgz_mod(ilev_mod+1))/(rgz_mod(ilev_mod)-rgz_mod(ilev_mod+1)) )
895
896 else
897 press(ilev) = 0.0
898 end if
899
900 end do
901
902 end function phf_convert_z_to_pressure
903
904 !--------------------------------------------------------------------------
905 ! PHF_CONVERT_Z_TO_GZ
906 !--------------------------------------------------------------------------
907 function phf_convert_z_to_gz(altitude,lat,nlev) result(rgz)
908 !
909 !:Purpose: Converts altitudes to geopotential heights. Uses the Helmert formula to
910 ! parameterize the latitude dependence and uses analytical result of the
911 ! integral of \int g(z)dz for the altitude dependence (see J.A. Dutton 1976,
912 ! p.65). At an altitude of 50 km, the altitude and geopotential height
913 ! differ by around 0.2-0.5 km, depending on the latitude.
914 !
915 implicit none
916
917 ! Arguments:
918 real(8), intent(in) :: altitude(nlev) ! altitudes (m)
919 real(8), intent(in) :: lat ! latitude (rad)
920 integer, intent(in) :: nlev ! number of levels
921 ! Result:
922 real(8) :: rgz(nlev) ! geopotential heights (m)
923
924 rgz = (ec_rg/9.8) * (1.-2.64D-03*cos(2.*lat)+5.9D-6*cos(2.*lat)**2) * ec_ra*altitude/(ec_ra+altitude)
925
926 end function phf_convert_z_to_gz
927
928 !--------------------------------------------------------------------------
929 ! PHF_GET_TROPOPAUSE
930 !--------------------------------------------------------------------------
931 function phf_get_tropopause(nmodlev,pressmod,tt,height,hu_opt) result(tropo_press)
932 !
933 !:Purpose: Determines pressure level of tropopause.
934 ! Final tropopause is taken as max pressure (lowest altitude) from the
935 ! water vapour and temperature based tropopauses.
936 !
937 implicit none
938
939 ! Arguments:
940 integer, intent(in) :: nmodlev ! Number of model levels
941 real(8), intent(in) :: pressmod(nmodlev) ! Model pressure array (Pa)
942 real(8), intent(in) :: tt(nmodlev) ! Model temperature (Kelvin)
943 real(8), intent(in) :: height(nmodlev) ! Model height (m)
944 real(8), intent(in), optional :: hu_opt(nmodlev) ! Model specific humidity
945 ! Result:
946 real(8) :: tropo_press ! Tropopause level in Pa (output)
947
948 ! Locals:
949 integer :: itop,i,k,ilaps
950 real(8) :: hu_ppmv1,hu_ppmv2,hu_ppmv3,xlaps,tropo_press_hu
951 real(8), parameter :: press_min=6000. ! Min tropoause pressure 60 hPa.; equivalent to ~ 20km
952 real(8), parameter :: height_min=6000.0 ! Min tropopause level in meters.
953 real(8), parameter :: ppmv_threshold=10.0
954 real(8), parameter :: tgrad_threshold=0.002 ! degrees/m (2 degrees/km)
955 real(8), parameter :: consth=0.160754938e+07 ! conversion from mass mixing ratio to ppmv; 1.0e+06 / (18.015/28.96)
956
957 tropo_press=-1.0
958 if (all(height.lt.0.0)) call utl_abort('phf_get_tropopause: Missing height for determining tropopause pressure')
959
960 ! Initialize tropopause pressure level using temperature gradient.
961 ! Thermal tropopause is defined as the lowest level (above height_min) at which (1) the lapse rate decreases
962 ! to <= 2 C/km and (2) the average lapse rate between this level and all higher levels within 2 km are <= 2 C/km.
963 ! Ref: International Meteorological Vocabulary (2nd ed.). Geneva: Secretariat of the World Meteorological
964 ! Organization. 1992. p. 636. ISBN 92-63-02182-1.
965 ! The second requirement, based on hu, may give levels that are to high (pressure too low) in the winter hemisphere.
966
967 do itop=3,nmodlev
968 if (pressmod(itop).ge.press_min) exit
969 end do
970 itop=itop-1
971
972 do i=nmodlev,itop+1,-1
973 if (height(i)-height(nmodlev).lt.height_min) cycle
974 xlaps=-(tt(i)-tt(i-1))/(height(i)-height(i-1))
975 if (xlaps.le.tgrad_threshold) then
976 ilaps=1
977 do k=i-1,itop,-1
978 if (height(k)-height(i).gt.2000.0) exit
979 xlaps=xlaps-(tt(k)-tt(k-1))/(height(k)-height(k-1))
980 ilaps=ilaps+1
981 end do
982 if (xlaps/ilaps.le.tgrad_threshold) exit
983 end if
984 enddo
985 tropo_press=pressmod(i)
986
987 ! Improve on tropopause pressure levels using specific humidity if available,
988
989 if (present(hu_opt)) then
990
991 ! Use water vapour
992
993 hu_ppmv1=0.0
994 do i=itop,nmodlev
995
996 ! Convert specific humidity to ppmv mixing ratio.
997 ! First apply r=q/(1-q) to convert to mass mixing ratio.
998
999 if (hu_opt(i).le.0.8.and.hu_opt(i).ge.0) then
1000 hu_ppmv2 = consth*hu_opt(i)/(1.0-hu_opt(i))
1001 else if (hu_opt(i).gt.0.8) then
1002 hu_ppmv2 = consth*0.8/(1.0-0.8)
1003 else if (hu_opt(i).lt.0.0) then
1004 hu_ppmv2 = 0.0
1005 end if
1006
1007 ! Check if transition point reached.
1008 ! Added requirement that levels below also satisfy this condition.
1009
1010 if (hu_ppmv2.ge.ppmv_threshold) then
1011 ilaps=1
1012 do k=i+1,nmodlev
1013 if (height(i)-height(k).gt.5000.0) exit
1014 if (hu_opt(k).le.0.8.and.hu_opt(k).ge.0) then
1015 hu_ppmv3 = consth*hu_opt(k)/(1.0-hu_opt(k))
1016 else if (hu_opt(k).gt.0.8) then
1017 hu_ppmv3 = consth*0.8/(1.0-0.8)
1018 else
1019 hu_ppmv3=0.0
1020 end if
1021 if (hu_ppmv3.lt.ppmv_threshold) ilaps=0
1022 end do
1023 if (ilaps.eq.1) exit
1024 end if
1025 hu_ppmv1=hu_ppmv2
1026 end do
1027
1028 if (hu_ppmv2.ge.ppmv_threshold.and.ilaps.eq.1) then
1029
1030 ! Interpolate between levels
1031
1032 if (abs(hu_ppmv2-hu_ppmv1).lt.0.1) hu_ppmv1=hu_ppmv2-0.1
1033! tropo_press_hu=(log(pressmod(i))*(ppmv_threshold-hu_ppmv1)+ &
1034! log(pressmod(i-1))*(hu_ppmv2-ppmv_threshold)) &
1035! /(hu_ppmv2-hu_ppmv1)
1036! tropo_press_hu=exp(tropo_press_hu)
1037 tropo_press_hu=pressmod(i)
1038
1039 tropo_press=min(tropo_press,tropo_press_hu)
1040 else
1041 write(*,*) 'phf_get_tropopause: Level and specific humidity: ',itop,hu_ppmv2
1042 call utl_abort('phf_get_tropopause: Specific humidity too small.')
1043 end if
1044
1045 end if
1046
1047 end function phf_get_tropopause
1048
1049 !--------------------------------------------------------------------------
1050 ! PHF_GET_PBL
1051 !--------------------------------------------------------------------------
1052 function phf_get_pbl(nmodlev,pressmod,tt,height,hu_opt,uu_opt,vv_opt) result(pbl_press)
1053 !
1054 !:Purpose: Determines pressure level of planetary boundary layer using
1055 ! a first threshold of 0.5 for the bulk Richadson number (after Mahrt, 1981;
1056 ! requires availability of uu and vv). Threshold reduced to largest value
1057 ! between 0.25 and 0.5 if first not satisfied.
1058 !
1059 ! If not found with this approach, applies a variant of the Heffter approach
1060 ! described in Aliabadi et al (2016), with some local variation.
1061 !
1062 ! References:
1063 !
1064 ! * Aliabadi A.A., R.M. Staebler, J. de Grandpre, A. Zadra, and P.A.
1065 ! Vaillancourt, 2016: Comparison of Estimated Atmospheric
1066 ! Boundary Layer Mixing Height in teh Arctic and Southern Great
1067 ! Plains under Statisticallt Stable Conditions: Experimental
1068 ! and Numerical Aspects, Submitted to Atmosphere-Ocean (2015).
1069 !
1070 ! * Mahrt, L. 1981: Modelling depth of the stable boundary-layer,
1071 ! Bound-Lay. Meteorol., 21, 3-19
1072 !
1073 ! * Heffter, J.L.,1980: Transport layer depth calculations, Second
1074 ! Joint Conference on Applications of Air Pollution Meteorology,
1075 ! New Orleans, LA, 24-27 March 1980. American Meteorological
1076 ! Society, Boston, MA.
1077 !
1078 ! Comments: Currently assumes (uu,vv) midlayer levels approximately
1079 ! at tt, height, and hu levels when size(uu).ne.nmodlev.
1080 !
1081 implicit none
1082
1083 ! Arguments:
1084 integer, intent(in) :: nmodlev ! Number of model levels for variables other than uu and vv
1085 real(8), intent(in) :: pressmod(nmodlev) ! Model pressure array (Pa)
1086 real(8), intent(in) :: tt(nmodlev) ! Model temperature (Kelvin)
1087 real(8), intent(in) :: height(nmodlev) ! Model height (meters)
1088 real(8), optional, intent(in) :: uu_opt(:) ! Model zonal wind component (m/s)
1089 real(8), optional, intent(in) :: vv_opt(:) ! Model meridional wind component (m/s)
1090 real(8), optional, intent(in) :: hu_opt(nmodlev) ! Specific humidity
1091 ! Result:
1092 real(8) :: pbl_press ! PBL level in Pa (output)
1093
1094 ! Locals:
1095 integer :: itop,i,id,igradmax,inv,iRiBmax
1096 real(8) :: RiB1,RiB2,RiBmax,zs,thetavs,thetavh(nmodlev),us,vs,uv,hus,huh,gradmax,grad
1097 real(8), parameter :: kappa = 287.04/1004.67 ! R/Cp
1098 real(8), parameter :: RiB_threshold=0.5, reduced=0.5
1099 ! Imposed min presssure of PBL height of 200 hPa (extreme; PBL height should normally be under 3km)
1100 real(8), parameter :: press_min=20000.
1101 real(8) :: huw(nmodlev)
1102
1103 pbl_press=-1.0
1104
1105 ! Set values for lowest prognostic level
1106
1107 i = nmodlev
1108
1109 if (all(height.lt.0.0)) call utl_abort('phf_get_pbl: Missing height for determining PBL pressure')
1110
1111 ! Convert hu to mass mixing ratio
1112
1113 if (present(hu_opt)) then
1114 huw(:)=hu_opt(:)
1115 else
1116 huw(:)=0.0
1117 end if
1118
1119 if (huw(i).le.0.8.and.huw(i).ge.0) then
1120 hus = huw(i)/(1.0-huw(i))
1121 else if (huw(i).gt.0.8) then
1122 hus = 0.8/(1.0-0.8)
1123 else if (huw(i).lt.0.0) then
1124 hus = 0.0
1125 end if
1126 zs = height(i)*0.001
1127
1128 ! Potential virtual temperature at lowest prognostic level
1129
1130 thetavs = tt(i)*(1.D5/pressmod(i))**kappa* (1.0 + 0.61*hus )
1131 thetavh(nmodlev)=thetavs
1132
1133 ! Set max vertical level
1134
1135 do itop=2,nmodlev-1
1136 if (pressmod(itop).ge.press_min) exit
1137 end do
1138
1139 RiB1=0.0
1140 RiB2=0.0
1141 RiBmax=0.0
1142 iRiBmax=0
1143 if (present(uu_opt).and.present(vv_opt)) then
1144 id=nmodlev-size(uu_opt)
1145 if (id.gt.1.or.id.lt.0) then
1146 call utl_abort('phf_get_pbl: Unexpected number of UV levels, nmodlev = ' // trim(utl_str(nmodlev)) // ' , size(uu) = ' // trim(utl_str(size(uu_opt))) )
1147 end if
1148 us = uu_opt(size(uu_opt))
1149 vs = vv_opt(size(vv_opt))
1150
1151! us,vs set to 0.0
1152! us=0.0
1153! vs=0.0
1154
1155 ! Calc RiB from near-surface to level attaining RiB_threshold
1156
1157 do i=nmodlev-1,itop,-1
1158
1159 if (huw(i).le.0.8.and.huw(i).ge.0) then
1160 huh = huw(i)/(1.0-huw(i))
1161 else if (huw(i).gt.0.8) then
1162 huh = 0.8/(1.0-0.8)
1163 else if (huw(i).lt.0.0) then
1164 huh = 0.0
1165 end if
1166 thetavh(i) = tt(i)*(1.D5/pressmod(i))**kappa* ( 1.0 + 0.61*huh )
1167
1168 if (id.eq.0) then
1169 uv = max( (uu_opt(i)-us)**2 + (vv_opt(i)-vs)**2, 1.D-8 )
1170 else
1171 ! Take layer midpoint values
1172 uv = max( ((uu_opt(i)+uu_opt(i-1))/2.0-us)**2 + ((vv_opt(i)+vv_opt(i-1))/2.0-vs)**2, 1.D-8 )
1173 end if
1174
1175 RiB2 = ec_rg * (thetavh(i)-thetavs) * (height(i)*0.001-zs) / (thetavs*uv)
1176 if (RiBmax.lt.RiB2.and.RiB2.ge.reduced*RiB_threshold) then
1177 RiBmax=RiB2
1178 iRiBmax=i
1179 end if
1180 if (RiB2.ge.RiB_threshold) exit
1181 RiB1=RiB2
1182 end do
1183 else
1184
1185 ! Calc only theta
1186
1187 do i=nmodlev-1,itop,-1
1188
1189 if (huw(i).le.0.8.and.huw(i).ge.0) then
1190 huh = huw(i)/(1.0-huw(i))
1191 else if (huw(i).gt.0.8) then
1192 huh = 0.8/(1.0-0.8)
1193 else if (huw(i).lt.0.0) then
1194 huh = 0.0
1195 end if
1196 thetavh(i) = tt(i)*(1.D5/pressmod(i))**kappa* ( 1.0 + 0.61*huh )
1197 end do
1198 end if
1199
1200 if (RiB2.ge.RiB_threshold) then
1201
1202 ! Interpolate between levels
1203
1204 pbl_press=(log(pressmod(i))*(RiB_threshold-RiB1)+ &
1205 log(pressmod(i+1))*(RiB2-RiB_threshold)) &
1206 /(RiB2-RiB1)
1207 pbl_press=exp(pbl_press)
1208 else if (RiBmax.ge.reduced*RiB_threshold) then
1209 ! Apply to level with largest RiB between reduced*RiB_threshold and RiB_threshold
1210 pbl_press=pressmod(iRiBmax)
1211 else
1212
1213 ! Estimate PBL level using the Heffter conditions:
1214 ! First find lowest inversion layer where dtheta>2K.
1215 ! If found, assign mid of layer as PBL level.
1216 ! Otherwise, assign PBL level as that with largest
1217 ! theta gradient.
1218
1219 i=nmodlev-1
1220 do while (i.gt.itop)
1221 !if (thetavh(i)-thetavh(i+1).gt.0.0) then
1222 if ((thetavh(i)-thetavh(i+1))/(height(i)-height(i+1)).ge.0.005) then
1223 ! Near bottom of inversion layer found
1224 inv=i+1
1225 i=i-1
1226 !do while (thetavh(i)-thetavh(i+1).gt.0.0.and.i.gt.itop)
1227 do while ((thetavh(i)-thetavh(i+1))/(height(i)-height(i+1)).ge.0.005.and.i.gt.itop)
1228 i=i-1
1229 end do
1230 if ((thetavh(i+1)-thetavh(inv)).gt.2.0) then
1231 ! Apply midlayer as PBL
1232 pbl_press=sqrt(pressmod(i+1)*pressmod(inv))
1233 exit
1234 end if
1235 else
1236 i=i-1
1237 end if
1238 end do
1239 if (pbl_press.le.0.0) then
1240 gradmax=-1.D30
1241 igradmax=nmodlev-1
1242 do i=nmodlev-1,itop,-1
1243 grad=(thetavh(i)-thetavh(i+1))/(height(i)-height(i+1))
1244 if (gradmax.lt.grad) then
1245 gradmax=(thetavh(i)-thetavh(i+1))/(height(i)-height(i+1))
1246 igradmax=i
1247 if (grad.ge.0.005) then ! Check next layer as well
1248 if ((thetavh(i-1)-thetavh(i))/(height(i-1)-height(i)).ge.0.005) exit
1249 end if
1250 end if
1251 end do
1252 pbl_press=pressmod(igradmax)
1253 ! write(*,*) 'phf_get_pbl: Warning2 - Max allowed altitude reached for. ',pbl_press,igradmax,gradmax,RiB2,iRiBmax,RiBmax
1254 !else
1255 ! write(*,*) 'phf_get_pbl: Warning1 - Max allowed altitude reached for. ',pbl_press,i,RiB2,iRiBmax,RiBmax
1256 end if
1257 end if
1258
1259 end function phf_get_pbl
1260
1261 !--------------------------------------------------------------------------
1262 ! phf_calcDistance
1263 !--------------------------------------------------------------------------
1264 function phf_calcDistance(lat2, lon2, lat1, lon1) result(distanceInM)
1265 !
1266 !:Purpose: Compute the distance between two point on earth: (lat1,lon1)
1267 ! and (lat2,lon2). Calcul utilisant la Formule d'Haversine.
1268 !
1269 ! Reference: R.W. Sinnott,'Virtues of Haversine',Sky and
1270 ! Telescope, vol.68, no.2, 1984, p.159)
1271 !
1272 implicit none
1273
1274 ! Arguments:
1275 real(8), intent(in) :: lat1
1276 real(8), intent(in) :: lon1
1277 real(8), intent(in) :: lat2
1278 real(8), intent(in) :: lon2
1279 ! Result:
1280 real(8) :: distanceInM
1281
1282 ! Locals:
1283 real(8) :: dlat, dlon, a, c
1284
1285 dlon = lon2 - lon1
1286 dlat = lat2 - lat1
1287
1288 a = (sin(dlat/2.d0))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2.d0))**2
1289 c = 2.d0 * atan2(sqrt(a),sqrt(1.d0-a))
1290 distanceInM = ec_ra * c
1291
1292 end function phf_calcDistance
1293
1294 !--------------------------------------------------------------------------
1295 ! phf_calcDistanceFast
1296 !--------------------------------------------------------------------------
1297 function phf_calcDistanceFast(lat2, lon2, lat1, lon1) result(distanceInM)
1298 !
1299 !:Purpose: Compute the distance between two point on earth: (lat1,lon1)
1300 ! and (lat2,lon2). Using a quick and dirty formula good for
1301 ! short distances not close to the pole.
1302 !
1303 implicit none
1304
1305 ! Arguments:
1306 real(8), intent(in) :: lat1
1307 real(8), intent(in) :: lon1
1308 real(8), intent(in) :: lat2
1309 real(8), intent(in) :: lon2
1310 ! Result:
1311 real(8) :: distanceInM
1312
1313 ! Locals:
1314 real(8) :: dlat, dlon
1315
1316 dlon = (lon2 - lon1)*cos(lat1)
1317 dlat = lat2 - lat1
1318
1319 distanceInM = ec_ra * sqrt(dlon*dlon + dlat*dlat)
1320
1321 end function phf_calcDistanceFast
1322
1323 !--------------------------------------------------------------------------
1324 ! PHF_GRAVITYSRF
1325 !--------------------------------------------------------------------------
1326 function phf_gravitysrf(sLat)
1327 !
1328 !:Purpose: Normal gravity on ellipsoidal surface
1329 !
1330 implicit none
1331
1332 ! Arguments:
1333 real(8), intent(in) :: sLat ! sin of latitude
1334 ! Result:
1335 real(8) :: phf_gravitysrf ! normal gravity (m/s2)
1336
1337 ! Locals:
1338 real(8) :: ks2
1339 real(8) :: e2s
1340
1341 ks2 = ec_wgs_TNGk * sLat*sLat
1342 e2s = 1.D0 - ec_wgs_e2 * sLat*sLat
1343 phf_gravitysrf = ec_wgs_GammaE * (1.D0 + ks2) / sqrt(e2s)
1344 end function phf_gravitysrf
1345
1346 !--------------------------------------------------------------------------
1347 ! PHF_GRAVITYALT
1348 !--------------------------------------------------------------------------
1349 function phf_gravityalt(sLat, Altitude)
1350 !
1351 !:Purpose: Normal gravity above the ellipsoidal surface
1352 !
1353 implicit none
1354
1355 ! Arguments:
1356 real(8), intent(in) :: sLat ! sin of latitude
1357 real(8), intent(in) :: Altitude ! altitude (m)
1358 ! Result:
1359 real(8) :: phf_gravityalt ! Normal gravity (m/s2)
1360
1361 ! Locals:
1362 real(8) :: C1
1363 real(8) :: C2
1364
1365 C1 =-2.D0/ec_wgs_a*(1.D0+ec_wgs_f+ec_wgs_m-2*ec_wgs_f*sLat*sLat)
1366 C2 = 3.D0/ec_wgs_a**2
1367 phf_gravityalt = phf_gravitysrf(sLat)* &
1368 (1.D0 + C1 * Altitude + C2 * Altitude**2)
1369 end function phf_gravityalt
1370
1371 !--------------------------------------------------------------------------
1372 ! PHF_HEIGHT2GEOPOTENTIAL
1373 !--------------------------------------------------------------------------
1374 subroutine phf_height2geopotential(altitude, latitude, geopotential, &
1375 printHeight_opt)
1376 !
1377 !:Purpose: Geopotential energy at a given point.
1378 ! Result is based on the WGS84 approximate expression for the
1379 ! gravity acceleration as a function of latitude and altitude,
1380 ! integrated with the trapezoidal rule.
1381 !
1382 implicit none
1383
1384 ! Arguments:
1385 real(8), intent(in) :: altitude(:) ! altitude (m)
1386 real(8), intent(in) :: latitude ! latitude (rad)
1387 real(8), intent(inout) :: geopotential(:) ! geopotential (m2/s2)
1388 logical, optional, intent(inout) :: printHeight_opt
1389
1390 ! Locals:
1391 integer :: nlev, nlev500m
1392 real(8), allocatable :: alt500m(:), gravity500m(:)
1393 real(8) :: delAlt, aveGravity, sLat, gravity, gravityM1
1394 integer :: i, ilev
1395
1396 nlev = size(altitude)
1397 sLat = sin(latitude)
1398
1399 nlev500m = int(altitude(nlev) / 500.D0)
1400 if ( nlev500m >= 1) then
1401 allocate(alt500m(0:nlev500m))
1402 allocate(gravity500m(0:nlev500m))
1403
1404 ! Calculate gravity and height of levels for 500m-layers
1405 do i = 0, nlev500m
1406 alt500m(i) = i * 500.0D0
1407 gravity500m(i) = phf_gravityalt(sLat, alt500m(i))
1408 enddo
1409
1410 geopotential(nlev) = 0.0D0
1411 ! integrate from surface on the 500m-layers untill below the desired altitude
1412 do i = 1, nlev500m
1413 delAlt = alt500m(i) - alt500m(i-1)
1414 aveGravity = 0.5D0 * (gravity500m(i) + gravity500m(i-1))
1415 geopotential(nlev) = geopotential(nlev) + delAlt * aveGravity
1416 enddo
1417
1418 ! Add the contribution from top of the last 500m-layer to the altitude
1419 delAlt = altitude(nlev) - alt500m(nlev500m)
1420 aveGravity = 0.5D0 * (phf_gravityalt(sLat, altitude(nlev)) + gravity500m(nlev500m))
1421 geopotential(nlev) = geopotential(nlev) + delAlt * aveGravity
1422 gravityM1 = phf_gravityalt(sLat, altitude(nlev))
1423
1424 deallocate(alt500m)
1425 deallocate(gravity500m)
1426
1427 else
1428 ! At surface, use local gravity to get height
1429 gravity = phf_gravityalt(sLat,altitude(nlev))
1430 geopotential(nlev) = altitude(nlev) * gravity
1431 gravityM1 = gravity
1432
1433 endif
1434
1435 ! At upper-levels, integrate on model levels to get height
1436 do ilev = nlev-1, 1, -1
1437 gravity = phf_gravityalt(sLat,altitude(ilev))
1438 aveGravity = 0.5D0 * (gravity + gravityM1)
1439 delAlt = altitude(ilev) - altitude(ilev+1)
1440 geopotential(ilev) = geopotential(ilev+1) + delAlt * aveGravity
1441 gravityM1 = gravity
1442 enddo
1443
1444 if ( present(printHeight_opt) ) then
1445 if ( printHeight_opt ) then
1446 write(*,*) 'phf_height2geopotential, Z_T:'
1447 write(*,*) geopotential(:)
1448
1449 printHeight_opt = .false.
1450 endif
1451 endif
1452
1453 end subroutine phf_height2geopotential
1454
1455end module physicsFunctions_mod
