source: codes/icosagcm/trunk/src/dcmip/dcmip2016_baroclinic_wave.f90

Last change on this file was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

Fixed iso c bindings
fixed warnings with -Wall -Wno-aliasing -Wno-unused -Wno-unused-dummy-argument -Wno-maybe-uninitialized -Wno-tabs warnings
Removed all unused variables (-Wunused-variable)
vector%dot_product is now dot_product_3d to avoid compilation warning "dot_product shadows intrinsic" with GCC
File size: 15.6 KB
Line 
1MODULE dcmip2016_baroclinic_wave_mod
2
3!=======================================================================
4!
5!  Date:  July 29, 2015
6!
7!  Functions for setting up idealized initial conditions for the
8!  Ullrich, Melvin, Staniforth and Jablonowski baroclinic instability.
9!
10!  SUBROUTINE baroclinic_wave_sample(
11!    deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q)
12!
13!  Options:
14!     deep    deep atmosphere (1 = yes or 0 = no)
15!    moist    include moisture (1 = yes or 0 = no)
16!    pertt    type of perturbation (0 = exponential, 1 = stream function)
17!        X    Earth scaling factor
18!
19!  Given a point specified by:
20!      lon    longitude (radians)
21!      lat    latitude (radians)
22!      p/z    pressure (Pa) / height (m)
23!  zcoords    1 if z is specified, 0 if p is specified
24!
25!  the functions will return:
26!        p    pressure if z is specified and zcoords = 1 (Pa)
27!        u    zonal wind (m s^-1)
28!        v    meridional wind (m s^-1)
29!        t    temperature (K)
30!   thetav    virtual potential temperature (K)
31!     phis    surface geopotential (m^2 s^-2)
32!       ps    surface pressure (Pa)
33!      rho    density (kj m^-3)
34!        q    water vapor mixing ratio (kg/kg)
35!
36!
37!  Author: Paul Ullrich
38!          University of California, Davis
39!          Email: paullrich@ucdavis.edu
40!
41!=======================================================================
42  IMPLICIT NONE
43 
44!=======================================================================
45!    Physical constants
46!=======================================================================
47
48  REAL(8), PARAMETER ::               &
49       a     = 6371220.0d0,           & ! Reference Earth's Radius (m)
50       Rd    = 287.0d0,               & ! Ideal gas const dry air (J kg^-1 K^1)
51       g     = 9.80616d0,             & ! Gravity (m s^2)
52       cp    = 1004.5d0,              & ! Specific heat capacity (J kg^-1 K^1)
53       Lvap  = 2.5d6,                 & ! Latent heat of vaporization of water
54       Rvap  = 461.5d0,               & ! Ideal gas constnat for water vapor
55       Mvap  = 0.608d0,               & ! Ratio of molar mass of dry air/water
56       pi    = 3.14159265358979d0,    & ! pi
57       p0    = 100000.0d0,            & ! surface pressure (Pa)
58       kappa = 2.d0/7.d0,             & ! Ratio of Rd to cp
59       omega = 7.29212d-5,            & ! Reference rotation rate of the Earth (s^-1)
60       deg2rad  = pi/180.d0             ! Conversion factor of degrees to radians
61
62!=======================================================================
63!    Test case parameters
64!=======================================================================
65  REAL(8), PARAMETER ::               &
66       T0E        = 310.d0     ,      & ! temperature at equatorial surface (K)
67       T0P        = 240.d0     ,      & ! temperature at polar surface (K)
68       B          = 2.d0       ,      & ! jet half-width parameter
69       K          = 3.d0       ,      & ! jet width parameter
70       lapse      = 0.005d0             ! lapse rate parameter
71
72  REAL(8), PARAMETER ::               &
73       pertu0     = 0.5d0      ,      & ! SF Perturbation wind velocity (m/s)
74       pertr      = 1.d0/6.d0  ,      & ! SF Perturbation radius (Earth radii)
75       pertup     = 1.0d0      ,      & ! Exp. perturbation wind velocity (m/s)
76       pertexpr   = 0.1d0      ,      & ! Exp. perturbation radius (Earth radii)
77       pertlon    = pi/9.d0    ,      & ! Perturbation longitude
78       pertlat    = 2.d0*pi/9.d0,     & ! Perturbation latitude
79       pertz      = 15000.d0   ,      & ! Perturbation height cap
80       dxepsilon  = 1.d-5               ! Small value for numerical derivatives
81 
82  REAL(8), PARAMETER ::               &
83       moistqlat  = 2.d0*pi/9.d0,     & ! Humidity latitudinal width
84       moistqp    = 34000.d0,         & ! Humidity vertical pressure width
85       moisttr    = 0.1d0,            & ! Vertical cut-off pressure for humidity
86       moistqs    = 1.d-12,           & ! Humidity above cut-off
87       moistq0    = 0.018d0,          & ! Maximum specific humidity
88       moistqr    = 0.9d0,            & ! Maximum saturation ratio
89       moisteps   = 0.622d0,          & ! Ratio of gas constants
90       moistT0    = 273.16d0,         & ! Reference temperature (K)
91       moistE0Ast = 610.78d0            ! Saturation vapor pressure at T0 (Pa)
92
93CONTAINS
94
95!=======================================================================
96!    Generate the baroclinic instability initial conditions
97!=======================================================================
98  SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,thetav,phis,ps,rho,q) &
99    BIND(c, name = "baroclinic_wave_test")
100    use iso_c_binding
101    IMPLICIT NONE
102
103!-----------------------------------------------------------------------
104!     input/output params parameters at given location
105!-----------------------------------------------------------------------
106    INTEGER(KIND=C_INT32_T), INTENT(IN)  :: &
107                deep,       & ! Deep (1) or Shallow (0) test case
108                moist,      & ! Moist (1) or Dry (0) test case
109                pertt         ! Perturbation type
110
111    REAL(KIND=C_DOUBLE), INTENT(IN)  :: &
112                lon,        & ! Longitude (radians)
113                lat,        & ! Latitude (radians)
114                X             ! Earth scaling parameter
115
116    REAL(KIND=C_DOUBLE), INTENT(INOUT) :: &
117                p,            & ! Pressure (Pa)
118                z               ! Altitude (m)
119
120    INTEGER(KIND=C_INT32_T), INTENT(IN) :: zcoords     ! 1 if z coordinates are specified
121                                       ! 0 if p coordinates are specified
122
123    REAL(KIND=C_DOUBLE), INTENT(OUT) :: &
124                u,          & ! Zonal wind (m s^-1)
125                v,          & ! Meridional wind (m s^-1)
126                t,          & ! Temperature (K)
127                thetav,     & ! Virtual potential temperature (K)
128                phis,       & ! Surface Geopotential (m^2 s^-2)
129                ps,         & ! Surface Pressure (Pa)
130                rho,        & ! density (kg m^-3)
131                q             ! water vapor mixing ratio (kg/kg)
132
133    !------------------------------------------------
134    !   Local variables
135    !------------------------------------------------
136    REAL(8) :: aref, omegaref
137    REAL(8) :: T0, constH, constC, scaledZ, inttau2, rratio
138    REAL(8) :: inttermU, bigU, rcoslat, omegarcoslat
139    REAL(8) :: eta
140
141    !------------------------------------------------
142    !   Pressure and temperature
143    !------------------------------------------------
144    if (zcoords .eq. 1) then
145      CALL evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
146    else
147      CALL evaluate_z_temperature(deep, X, lon, lat, p, z, t)
148    end if
149
150    !------------------------------------------------
151    !   Compute test case constants
152    !------------------------------------------------
153    aref = a / X
154    omegaref = omega * X
155
156    T0 = 0.5d0 * (T0E + T0P)
157
158    constH = Rd * T0 / g
159
160    constC = 0.5d0 * (K + 2.d0) * (T0E - T0P) / (T0E * T0P)
161
162    scaledZ = z / (B * constH)
163
164    inttau2 = constC * z * exp(- scaledZ**2)
165
166    ! radius ratio
167    if (deep .eq. 0) then
168      rratio = 1.d0
169    else
170      rratio = (z + aref) / aref;
171    end if
172
173    !-----------------------------------------------------
174    !   Initialize surface pressure
175    !-----------------------------------------------------
176    ps = p0
177
178    !-----------------------------------------------------
179    !   Initialize velocity field
180    !-----------------------------------------------------
181    inttermU = (rratio * cos(lat))**(K - 1.d0) - (rratio * cos(lat))**(K + 1.d0)
182    bigU = g / aref * K * inttau2 * inttermU * t
183    if (deep .eq. 0) then
184      rcoslat = aref * cos(lat)
185    else
186      rcoslat = (z + aref) * cos(lat)
187    end if
188
189    omegarcoslat = omegaref * rcoslat
190   
191    u = - omegarcoslat + sqrt(omegarcoslat**2 + rcoslat * bigU)
192    v = 0.d0
193
194    !-----------------------------------------------------
195    !   Add perturbation to the velocity field
196    !-----------------------------------------------------
197
198    ! Exponential type
199    if (pertt .eq. 0) then
200      u = u + evaluate_exponential(lon, lat, z)
201
202    ! Stream function type
203    elseif (pertt .eq. 1) then
204      u = u - 1.d0 / (2.d0 * dxepsilon) *                       &
205          ( evaluate_streamfunction(lon, lat + dxepsilon, z)    &
206          - evaluate_streamfunction(lon, lat - dxepsilon, z))
207
208      v = v + 1.d0 / (2.d0 * dxepsilon * cos(lat)) *            &
209          ( evaluate_streamfunction(lon + dxepsilon, lat, z)    &
210          - evaluate_streamfunction(lon - dxepsilon, lat, z))
211    end if
212
213    !-----------------------------------------------------
214    !   Initialize surface geopotential
215    !-----------------------------------------------------
216    phis = 0.d0
217
218    !-----------------------------------------------------
219    !   Initialize density
220    !-----------------------------------------------------
221    rho = p / (Rd * t)
222
223    !-----------------------------------------------------
224    !   Initialize specific humidity
225    !-----------------------------------------------------
226    if (moist .eq. 1) then
227      eta = p/p0
228
229      if (eta .gt. moisttr) then
230        q = moistq0 * exp(- (lat/moistqlat)**4)          &
231                    * exp(- ((eta-1.d0)*p0/moistqp)**2)
232      else
233        q = moistqs
234      end if
235
236      ! Convert virtual temperature to temperature
237      t = t / (1.d0 + Mvap * q)
238
239    else
240      q = 0.d0
241    end if
242
243    !-----------------------------------------------------
244    !   Initialize virtual potential temperature
245    !-----------------------------------------------------
246    thetav = t * (1.d0 + 0.61d0 * q) * (p0 / p)**(Rd / cp)
247
248  END SUBROUTINE baroclinic_wave_test
249
250!-----------------------------------------------------------------------
251!    Calculate pointwise pressure and temperature
252!-----------------------------------------------------------------------
253  SUBROUTINE evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
254
255    INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
256
257    REAL(8), INTENT(IN)  :: &
258                X,          & ! Earth scaling ratio
259                lon,        & ! Longitude (radians)
260                lat,        & ! Latitude (radians)
261                z             ! Altitude (m)
262
263    REAL(8), INTENT(OUT) :: &
264                p,          & ! Pressure (Pa)
265                t             ! Temperature (K)
266
267    REAL(8) :: aref, omegaref
268    REAL(8) :: T0, constA, constB, constC, constH, scaledZ
269    REAL(8) :: tau1, tau2, inttau1, inttau2
270    REAL(8) :: rratio, inttermT
271
272    !--------------------------------------------
273    ! Constants
274    !--------------------------------------------
275    aref = a / X
276    omegaref = omega * X
277
278    T0 = 0.5d0 * (T0E + T0P)
279    constA = 1.d0 / lapse
280    constB = (T0 - T0P) / (T0 * T0P)
281    constC = 0.5d0 * (K + 2.d0) * (T0E - T0P) / (T0E * T0P)
282    constH = Rd * T0 / g
283
284    scaledZ = z / (B * constH)
285
286    !--------------------------------------------
287    !    tau values
288    !--------------------------------------------
289    tau1 = constA * lapse / T0 * exp(lapse * z / T0) &
290         + constB * (1.d0 - 2.d0 * scaledZ**2) * exp(- scaledZ**2)
291    tau2 = constC * (1.d0 - 2.d0 * scaledZ**2) * exp(- scaledZ**2)
292
293    inttau1 = constA * (exp(lapse * z / T0) - 1.d0) &
294            + constB * z * exp(- scaledZ**2)
295    inttau2 = constC * z * exp(- scaledZ**2)
296
297    !--------------------------------------------
298    !    radius ratio
299    !--------------------------------------------
300    if (deep .eq. 0) then
301      rratio = 1.d0
302    else
303      rratio = (z + aref) / aref;
304    end if
305
306    !--------------------------------------------
307    !    interior term on temperature expression
308    !--------------------------------------------
309    inttermT = (rratio * cos(lat))**K &
310             - K / (K + 2.d0) * (rratio * cos(lat))**(K + 2.d0)
311
312    !--------------------------------------------
313    !    temperature
314    !--------------------------------------------
315    t = 1.d0 / (rratio**2 * (tau1 - tau2 * inttermT))
316
317    !--------------------------------------------
318    !    hydrostatic pressure
319    !--------------------------------------------
320    p = p0 * exp(- g / Rd * (inttau1 - inttau2 * inttermT))
321
322  END SUBROUTINE evaluate_pressure_temperature
323
324!-----------------------------------------------------------------------
325!    Calculate pointwise z and temperature given pressure
326!-----------------------------------------------------------------------
327  SUBROUTINE evaluate_z_temperature(deep, X, lon, lat, p, z, t)
328   
329    INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
330
331    REAL(8), INTENT(IN)  :: &
332                X,          & ! Earth scaling ratio
333                lon,        & ! Longitude (radians)
334                lat,        & ! Latitude (radians)
335                p             ! Pressure (Pa)
336
337    REAL(8), INTENT(OUT) :: &
338                z,          & ! Altitude (m)
339                t             ! Temperature (K)
340
341    INTEGER :: ix
342
343    REAL(8) :: z0, z1, z2
344    REAL(8) :: p0, p1, p2
345
346    z0 = 0.d0
347    z1 = 10000.d0
348
349    CALL evaluate_pressure_temperature(deep, X, lon, lat, z0, p0, t)
350    CALL evaluate_pressure_temperature(deep, X, lon, lat, z1, p1, t)
351
352    DO ix = 1, 100
353      z2 = z1 - (p1 - p) * (z1 - z0) / (p1 - p0)
354
355      CALL evaluate_pressure_temperature(deep, X, lon, lat, z2, p2, t)
356
357      IF (ABS((p2 - p)/p) .lt. 1.0d-13) THEN
358        EXIT
359      END IF
360
361      z0 = z1
362      p0 = p1
363
364      z1 = z2
365      p1 = p2
366    END DO
367
368    z = z2
369
370    CALL evaluate_pressure_temperature(deep, X, lon, lat, z, p0, t)
371
372  END SUBROUTINE evaluate_z_temperature
373
374!-----------------------------------------------------------------------
375!    Exponential perturbation function
376!-----------------------------------------------------------------------
377  REAL(8) FUNCTION evaluate_exponential(lon, lat, z)
378
379    REAL(8), INTENT(IN)  :: &
380                lon,        & ! Longitude (radians)
381                lat,        & ! Latitude (radians)
382                z             ! Altitude (meters)
383
384    REAL(8) :: greatcircler, perttaper
385
386    ! Great circle distance
387    greatcircler = 1.d0 / pertexpr &
388      * acos(sin(pertlat) * sin(lat) + cos(pertlat) * cos(lat) * cos(lon - pertlon))
389
390    ! Vertical tapering of stream function
391    if (z < pertz) then
392      perttaper = 1.d0 - 3.d0 * z**2 / pertz**2 + 2.d0 * z**3 / pertz**3
393    else
394      perttaper = 0.d0
395    end if
396
397    ! Zonal velocity perturbation
398    if (greatcircler < 1.d0) then
399      evaluate_exponential = pertup * perttaper * exp(- greatcircler**2)
400    else
401      evaluate_exponential = 0.d0
402    end if
403
404  END FUNCTION evaluate_exponential
405
406!-----------------------------------------------------------------------
407!    Stream function perturbation function
408!-----------------------------------------------------------------------
409  REAL(8) FUNCTION evaluate_streamfunction(lon, lat, z)
410
411    REAL(8), INTENT(IN)  :: &
412                lon,        & ! Longitude (radians)
413                lat,        & ! Latitude (radians)
414                z             ! Altitude (meters)
415
416    REAL(8) :: greatcircler, perttaper, cospert
417
418    ! Great circle distance
419    greatcircler = 1.d0 / pertr &
420      * acos(sin(pertlat) * sin(lat) + cos(pertlat) * cos(lat) * cos(lon - pertlon))
421
422    ! Vertical tapering of stream function
423    if (z < pertz) then
424      perttaper = 1.d0 - 3.d0 * z**2 / pertz**2 + 2.d0 * z**3 / pertz**3
425    else
426      perttaper = 0.d0
427    end if
428
429    ! Horizontal tapering of stream function
430    if (greatcircler .lt. 1.d0) then
431      cospert = cos(0.5d0 * pi * greatcircler)
432    else
433      cospert = 0.d0
434    end if
435
436    evaluate_streamfunction = &
437        (- pertu0 * pertr * perttaper * cospert**4)
438
439  END FUNCTION evaluate_streamfunction
440
441END MODULE dcmip2016_baroclinic_wave_mod
442
Note: See TracBrowser for help on using the repository browser.