source: LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.f90 @ 5523

Last change on this file since 5523 was 5512, checked in by yann meurdesoif, 8 days ago

Implement GPU automatic port for :

  • Thermics
  • acama_gwd_rando
  • flott_gwd_rando

YM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory.
  • Property svn:keywords set to Id
File size: 15.6 KB
Line 
1!
2! $Id: flott_gwd_rando_m.f90 5512 2025-01-28 18:07:51Z fhourdin $
3!
4module FLOTT_GWD_rando_m
5
6  USE clesphys_mod_h
7      implicit none
8    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
9    INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed
10    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
11    LOGICAL, SAVE :: firstcall = .TRUE.
12   !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
13
14contains
15
16  SUBROUTINE FLOTT_GWD_rando_first
17  use dimphy, only: klev
18  USE ioipsl_getin_p_mod, ONLY : getin_p
19  IMPLICIT NONE
20    CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m'
21    CHARACTER (LEN=80) :: abort_message
22 
23    IF (firstcall) THEN
24    ! Cle introduite pour resoudre un probleme de non reproductibilite
25    ! Le but est de pouvoir tester de revenir a la version precedenete
26    ! A eliminer rapidement
27    CALL getin_p('gwd_reproductibilite_mpiomp',gwd_reproductibilite_mpiomp)
28    IF (NW+3*NA>=KLEV) THEN
29       abort_message = 'NW+3*NA>=KLEV Probleme pour generation des ondes'
30       CALL abort_physic (modname,abort_message,1)
31    ENDIF
32    firstcall=.false.
33  ENDIF
34  END SUBROUTINE FLOTT_GWD_rando_first
35
36
37  SUBROUTINE FLOTT_GWD_rando(DTIME, PP, tt, uu, vv, prec, zustr, zvstr, d_u, &
38       d_v,east_gwstress,west_gwstress)
39
40    ! Parametrization of the momentum flux deposition due to a discrete
41    ! number of gravity waves.
42    ! Author: F. Lott
43    ! July, 12th, 2012
44    ! Gaussian distribution of the source, source is precipitation
45    ! Reference: Lott (JGR, vol 118, page 8897, 2013)
46
47    !ONLINE:
48      USE yomcst_mod_h
49use dimphy, only: klon, klev
50      use assert_m, only: assert
51      USE ioipsl_getin_p_mod, ONLY : getin_p
52      USE vertical_layers_mod, ONLY : presnivs
53      USE yoegwd_mod_h
54      USE lmdz_fake_call, ONLY : fake_call
55
56      CHARACTER (LEN=20),PARAMETER :: modname='flott_gwd_rando'
57      CHARACTER (LEN=80) :: abort_message
58
59    ! OFFLINE:
60    ! include "dimensions_mod.f90"
61    ! include "dimphy.h"
62    ! END OF DIFFERENCE ONLINE-OFFLINE
63
64    ! 0. DECLARATIONS:
65
66    ! 0.1 INPUTS
67    REAL, intent(in)::DTIME ! Time step of the Physics
68    REAL, intent(in):: pp(KLON, KLEV) ! (KLON, KLEV) Pressure at full levels
69    REAL, intent(in):: prec(KLON) ! (klon) Precipitation (kg/m^2/s)
70    REAL, intent(in):: TT(KLON, KLEV) ! (KLON, KLEV) Temp at full levels
71    REAL, intent(in):: UU(KLON, KLEV) ! (KLON, KLEV) Zonal wind at full levels
72    REAL, intent(in):: VV(KLON, KLEV) ! (KLON, KLEV) Merid wind at full levels
73
74    ! 0.2 OUTPUTS
75    REAL, intent(out):: zustr(KLON), zvstr(KLON) ! (KLON) Surface Stresses
76
77    REAL, intent(inout):: d_u(KLON, KLEV), d_v(KLON, KLEV)
78    REAL, intent(inout):: east_gwstress(KLON, KLEV) !  Profile of eastward stress
79    REAL, intent(inout):: west_gwstress(KLON, KLEV) !  Profile of westward stress
80
81    ! (KLON, KLEV) tendencies on winds
82
83    ! O.3 INTERNAL ARRAYS
84    REAL BVLOW(klon)
85    REAL DZ   !  Characteristic depth of the Source
86
87    INTEGER II, JJ, LL
88
89    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
90
91    REAL DELTAT
92
93    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
94
95    INTEGER JK, JP, JO, JW
96    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
97    REAL CMAX ! standard deviation of the phase speed distribution
98    REAL RUWMAX,SAT  ! ONLINE SPECIFIED IN run.def
99    REAL CPHA ! absolute PHASE VELOCITY frequency
100    REAL ZK(KLON, NW) ! Horizontal wavenumber amplitude
101    REAL ZP(KLON, NW) ! Horizontal wavenumber angle
102    REAL ZO(KLON, NW) ! Absolute frequency !
103
104    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
105    REAL ZOM(KLON, NW), ZOP(KLON, NW)
106
107    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
108    REAL WWM(KLON, NW), WWP(KLON, NW)
109
110    REAL RUW0(KLON, NW) ! Fluxes at launching level
111
112    REAL RUWP(KLON, NW), RVWP(KLON, NW)
113    ! Fluxes X and Y for each waves at 1/2 Levels
114
115    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude
116
117    REAL XLAUNCH ! Controle the launching altitude
118    REAL XTROP ! SORT of Tropopause altitude
119    REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
120    REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
121
122    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
123    ! for GWs parameterisation apply
124
125    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
126
127    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
128
129    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
130
131    REAL H0 ! Characteristic Height of the atmosphere
132    REAL PR, TR ! Reference Pressure and Temperature
133
134    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
135
136    REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
137    REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
138    REAL PSEC ! Security to avoid division by 0 pressure
139    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
140    REAL BVSEC ! Security to avoid negative BVF
141    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
142
143    REAL, DIMENSION(klev+1) ::HREF
144
145
146    !-----------------------------------------------------------------
147
148    ! 1. INITIALISATIONS
149
150    ! 1.1 Basic parameter
151
152    ! Are provided from elsewhere (latent heat of vaporization, dry
153    ! gaz constant for air, gravity constant, heat capacity of dry air
154    ! at constant pressure, earth rotation rate, pi).
155
156    ! 1.2 Tuning parameters of V14
157
158   
159    RDISS = 0.5 ! Diffusion parameter
160    ! ONLINE
161      RUWMAX=GWD_RANDO_RUWMAX
162      SAT=gwd_rando_sat
163    !END ONLINE
164    ! OFFLINE
165    ! RUWMAX= 1.75    ! Launched flux
166    ! SAT=0.25     ! Saturation parameter
167    ! END OFFLINE
168
169    PRMAX = 20. / 24. /3600.
170    ! maximum of rain for which our theory applies (in kg/m^2/s)
171
172 ! Characteristic depth of the source
173    DZ = 1000.
174    XLAUNCH=0.5 ! Parameter that control launching altitude
175    XTROP=0.2 ! Parameter that control tropopause altitude
176    DELTAT=24.*3600. ! Time scale of the waves (first introduced in 9b)
177    !  OFFLINE
178    !  DELTAT=DTIME
179    !  END OFFLINE
180
181    KMIN = 2.E-5
182    ! minimum horizontal wavenumber (inverse of the subgrid scale resolution)
183
184    KMAX = 1.E-3 ! Max horizontal wavenumber
185    CMAX = 30. ! Max phase speed velocity
186
187    TR = 240. ! Reference Temperature
188    PR = 101300. ! Reference pressure
189    H0 = RD * TR / RG ! Characteristic vertical scale height
190
191    BVSEC = 5.E-3 ! Security to avoid negative BVF
192    PSEC = 1.E-6 ! Security to avoid division by 0 pressure
193    ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
194
195IF (1==0) THEN
196    !ONLINE
197        call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
198         size(vv, 1), size(zustr), size(zvstr), size(d_u, 1), &
199         size(d_v, 1), &
200         size(east_gwstress, 1), size(west_gwstress, 1) /), &
201         "FLOTT_GWD_RANDO klon")
202     call assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &
203          size(vv, 2), size(d_u, 2), size(d_v, 2), &
204          size(east_gwstress,2), size(west_gwstress,2) /), &
205          "FLOTT_GWD_RANDO klev")
206    !END ONLINE
207ENDIF
208   
209    IF(DELTAT < DTIME)THEN
210       abort_message='flott_gwd_rando: deltat < dtime!'
211       CALL abort_physic(modname,abort_message,1)
212    ENDIF
213
214    IF (KLEV < NW) THEN
215       abort_message='flott_gwd_rando: you will have problem with random numbers'
216       CALL abort_physic(modname,abort_message,1)
217    ENDIF
218   
219    CALL FAKE_CALL(BVLOW)  ! to be suppress in future
220   
221    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
222
223    ! Pressure and Inv of pressure
224    DO LL = 2, KLEV
225       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
226    end DO
227    PH(:, KLEV + 1) = 0.
228    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
229
230    ! Launching altitude
231
232    !Pour revenir a la version non reproductible en changeant le nombre de process
233    IF (gwd_reproductibilite_mpiomp) THEN
234       ! Reprend la formule qui calcule PH en fonction de PP=play
235       DO LL = 2, KLEV
236          HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1))) / 2.)
237       end DO
238       HREF(KLEV + 1) = 0.
239       HREF(1) = 2. * presnivs(1) - HREF(2)
240    ELSE
241       HREF(1:KLEV)=PH(KLON/2,1:KLEV)
242    ENDIF
243
244    LAUNCH=0
245    LTROP =0
246    DO LL = 1, KLEV
247       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
248    ENDDO
249    DO LL = 1, KLEV
250       IF (HREF(LL) / HREF(1) > XTROP) LTROP = LL
251    ENDDO
252    !LAUNCH=22 ; LTROP=33
253!   print*,'LAUNCH=',LAUNCH,'LTROP=',LTROP
254
255    ! Log pressure vert. coordinate
256    DO LL = 1, KLEV + 1
257       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
258    end DO
259
260    ! BV frequency
261    DO LL = 2, KLEV
262       ! BVSEC: BV Frequency (UH USED IS AS A TEMPORARY ARRAY DOWN TO WINDS)
263       UH(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
264            * RD**2 / RCPD / H0**2 + (TT(:, LL) &
265            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0
266    end DO
267    BVLOW(:) = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &
268         * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) &
269         - TT(:, LAUNCH))/(ZH(:, LTROP )- ZH(:, LAUNCH)) * RD / H0
270
271    UH(:, 1) = UH(:, 2)
272    UH(:, KLEV + 1) = UH(:, KLEV)
273    BV(:, 1) = UH(:, 2)
274    BV(:, KLEV + 1) = UH(:, KLEV)
275    ! SMOOTHING THE BV HELPS
276    DO LL = 2, KLEV
277       BV(:, LL)=(UH(:, LL+1)+2.*UH(:, LL)+UH(:, LL-1))/4.
278    end DO
279
280    BV=MAX(SQRT(MAX(BV, 0.)), BVSEC)
281    BVLOW=MAX(SQRT(MAX(BVLOW, 0.)), BVSEC)
282
283
284    ! WINDS
285    DO LL = 2, KLEV
286       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
287       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
288    end DO
289    UH(:, 1) = 0.
290    VH(:, 1) = 0.
291    UH(:, KLEV + 1) = UU(:, KLEV)
292    VH(:, KLEV + 1) = VV(:, KLEV)
293
294    ! 3 WAVES CHARACTERISTICS CHOSEN RANDOMLY AT THE LAUNCH ALTITUDE
295
296    ! The mod functions of weird arguments are used to produce the
297    ! waves characteristics in an almost stochastic way
298
299    DO JW = 1, NW
300             ! Angle
301             DO II = 1, KLON
302                ! Angle (0 or PI so far)
303                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
304                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
305                ZP(II, JW) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
306                     * RPI / 2.
307                ! Horizontal wavenumber amplitude
308                ZK(II, JW) = KMIN + (KMAX - KMIN) *RAN_NUM_2
309                ! Horizontal phase speed
310                CPHA = 0.
311                DO JJ = 1, NA
312                    RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.)
313                    CPHA = CPHA + &
314                    CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.)
315                END DO
316                IF (CPHA.LT.0.)  THEN
317                   CPHA = -1.*CPHA
318                   ZP(II, JW) = ZP(II, JW) + RPI
319                ENDIF
320                ! Absolute frequency is imposed
321                ZO(II, JW) = CPHA * ZK(II, JW)
322                ! Intrinsic frequency is imposed
323                ZO(II, JW) = ZO(II, JW) &
324                     + ZK(II, JW) * COS(ZP(II, JW)) * UH(II, LAUNCH) &
325                     + ZK(II, JW) * SIN(ZP(II, JW)) * VH(II, LAUNCH)
326                ! Momentum flux at launch lev
327                RUW0(II, JW) = RUWMAX
328             ENDDO
329    ENDDO
330
331    ! 4. COMPUTE THE FLUXES
332
333    ! 4.1 Vertical velocity at launching altitude to ensure
334    ! the correct value to the imposed fluxes.
335
336    DO JW = 1, NW
337
338       ! Evaluate intrinsic frequency at launching altitude:
339       ZOP(:,JW) = ZO(:, JW) &
340            - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LAUNCH) &
341            - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LAUNCH)
342
343       ! VERSION WITH CONVECTIVE SOURCE
344
345       ! Vertical velocity at launch level, value to ensure the
346       ! imposed factor related to the convective forcing:
347       ! precipitations.
348
349       ! tanh limitation to values above prmax:
350       WWP(:, JW) = RUW0(:, JW) &
351            * (RD / RCPD / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
352
353       ! Factor related to the characteristics of the waves:
354       WWP(:, JW) = WWP(:, JW) * ZK(:, JW)**3 / KMIN / BVLOW(:)  &
355            / MAX(ABS(ZOP(:, JW)), ZOISEC)**3
356
357       ! Moderation by the depth of the source (dz here):
358       WWP(:, JW) = WWP(:, JW) &
359            * EXP(- BVLOW(:)**2 / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 * ZK(:, JW)**2 &
360            * DZ**2)
361
362       ! Put the stress in the right direction:
363       RUWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
364            * BV(:, LAUNCH) * COS(ZP(:, JW)) * WWP(:, JW)**2
365       RVWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 &
366            * BV(:, LAUNCH) * SIN(ZP(:, JW)) * WWP(:, JW)**2
367    end DO
368
369
370    ! 4.2 Uniform values below the launching altitude
371
372    DO LL = 1, LAUNCH
373       RUW(:, LL) = 0
374       RVW(:, LL) = 0
375       DO JW = 1, NW
376          RUW(:, LL) = RUW(:, LL) + RUWP(:, JW)
377          RVW(:, LL) = RVW(:, LL) + RVWP(:, JW)
378       end DO
379    end DO
380
381    ! 4.3 Loop over altitudes, with passage from one level to the next
382    ! done by i) conserving the EP flux, ii) dissipating a little,
383    ! iii) testing critical levels, and vi) testing the breaking.
384
385    DO LL = LAUNCH, KLEV - 1
386       ! Warning: all the physics is here (passage from one level
387       ! to the next)
388       DO JW = 1, NW
389          ZOM(:, JW) = ZOP(:,JW)
390          WWM(:, JW) = WWP(:, JW)
391          ! Intrinsic Frequency
392          ZOP(:, JW) = ZO(:, JW) - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LL + 1) &
393               - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LL + 1)
394
395          ! No breaking (Eq.6)
396          ! Dissipation (Eq. 8)
397          WWP(:, JW) = WWM(:, JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &
398               + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
399               / MAX(ABS(ZOP(:, JW) + ZOM(:, JW)) / 2., ZOISEC)**4 &
400               * ZK(:, JW)**3 * (ZH(:, LL + 1) - ZH(:, LL)))
401
402          ! Critical levels (forced to zero if intrinsic frequency changes sign)
403          ! Saturation (Eq. 12)
404          WWP(:, JW) = min(WWP(:, JW), MAX(0., &
405               SIGN(1., ZOP(:, JW) * ZOM(:, JW))) * ABS(ZOP(:, JW))**3 &
406               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2  &
407               * SAT**2 / ZK(:, JW)**4)
408       end DO
409
410       ! Evaluate EP-flux from Eq. 7 and give the right orientation to
411       ! the stress
412
413       DO JW = 1, NW
414          RUWP(:, JW) = SIGN(1., ZOP(:, JW))*COS(ZP(:, JW)) * WWP(:, JW)
415          RVWP(:, JW) = SIGN(1., ZOP(:, JW))*SIN(ZP(:, JW)) * WWP(:, JW)
416       end DO
417
418       RUW(:, LL + 1) = 0.
419       RVW(:, LL + 1) = 0.
420
421       DO JW = 1, NW
422          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:, JW)
423          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:, JW)
424          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:, JW))/REAL(NW)
425          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:, JW))/REAL(NW)
426       end DO
427    end DO
428! OFFLINE ONLY
429!   PRINT *,'SAT PROFILE:'
430!   DO LL=1,KLEV
431!   PRINT *,ZH(KLON/2,LL)/1000.,SAT*(2.+TANH(ZH(KLON/2,LL)/H0-8.))
432!   ENDDO
433
434    ! 5 CALCUL DES TENDANCES:
435
436    ! 5.1 Rectification des flux au sommet et dans les basses couches
437
438    RUW(:, KLEV + 1) = 0.
439    RVW(:, KLEV + 1) = 0.
440    RUW(:, 1) = RUW(:, LAUNCH)
441    RVW(:, 1) = RVW(:, LAUNCH)
442    DO LL = 1, LAUNCH
443       RUW(:, LL) = RUW(:, LAUNCH+1)
444       RVW(:, LL) = RVW(:, LAUNCH+1)
445       EAST_GWSTRESS(:, LL)  = EAST_GWSTRESS(:, LAUNCH)
446       WEST_GWSTRESS(:, LL)  = WEST_GWSTRESS(:, LAUNCH)
447    end DO
448
449    ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
450    DO LL = 1, KLEV
451       D_U(:, LL) = (1.-DTIME/DELTAT) * D_U(:, LL) + DTIME/DELTAT/REAL(NW) * &
452            RG * (RUW(:, LL + 1) - RUW(:, LL)) &
453            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
454       ! NO AR-1 FOR MERIDIONAL TENDENCIES
455       D_V(:, LL) =                                            1./REAL(NW) * &
456            RG * (RVW(:, LL + 1) - RVW(:, LL)) &
457            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
458    ENDDO
459
460    ! Cosmetic: evaluation of the cumulated stress
461    ZUSTR = 0.
462    ZVSTR = 0.
463    DO LL = 1, KLEV
464       ZUSTR = ZUSTR + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME
465       ZVSTR = ZVSTR + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME
466    ENDDO
467
468
469  END SUBROUTINE FLOTT_GWD_RANDO
470
471end module FLOTT_GWD_rando_m
Note: See TracBrowser for help on using the repository browser.