source: codes/icosagcm/trunk/src/dcmip/dcmip2016_cyclone.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: 10.8 KB
Line 
1MODULE dcmip2016_cyclone_mod
2
3!=======================================================================
4!
5!  Date:  July 29, 2015
6!
7!  Function for setting up idealized tropical cyclone initial conditions
8!
9!  SUBROUTINE tropical_cyclone_sample(
10!    lon,lat,p,z,zcoords,u,v,t,thetav,phis,ps,rho,q)
11!
12!  Given a point specified by:
13!      lon    longitude (radians)
14!      lat    latitude (radians)
15!      p/z    pressure (Pa) / height (m)
16!  zcoords    1 if z is specified, 0 if p is specified
17!
18!  the functions will return:
19!        p    pressure if z is specified (Pa)
20!        z    geopotential height if p is specified (m)
21!        u    zonal wind (m s^-1)
22!        v    meridional wind (m s^-1)
23!        t    temperature (K)
24!   thetav    virtual potential temperature (K)
25!     phis    surface geopotential (m^2 s^-2)
26!       ps    surface pressure (Pa)
27!      rho    density (kj m^-3)
28!        q    specific humidity (kg/kg)
29!
30!  Initial data are currently identical to:
31!
32!       Reed, K. A., and C. Jablonowski, 2011: An analytic
33!       vortex initialization technique for idealized tropical
34!       cyclone studies in AGCMs. Mon. Wea. Rev., 139, 689-710.
35!
36!  Author: Kevin A. Reed
37!          Stony Brook University
38!          Email: kevin.a.reed@stonybrook.edu
39!
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       rp         = 282000.d0,  & ! Radius for calculation of PS
67       dp         = 1115.d0,    & ! Delta P for calculation of PS
68       zp         = 7000.d0,    & ! Height for calculation of P
69       q0         = 0.021d0,    & ! q at surface from Jordan
70       gamma      = 0.007d0,    & ! lapse rate
71       Ts0        = 302.15d0,   & ! Surface temperature (SST)
72       p00        = 101500.d0,  & ! global mean surface pressure
73       cen_lat    = 10.d0,      & ! Center latitude of initial vortex
74       cen_lon    = 180.d0,     & ! Center longitufe of initial vortex
75       zq1        = 3000.d0,    & ! Height 1 for q calculation
76       zq2        = 8000.d0,    & ! Height 2 for q calculation
77       exppr      = 1.5d0,      & ! Exponent for r dependence of p
78       exppz      = 2.d0,       & ! Exponent for z dependence of p
79       ztrop      = 15000.d0,   & ! Tropopause Height
80       qtrop      = 1.d-11,     & ! Tropopause specific humidity
81       rfpi       = 1000000.d0, & ! Radius within which to use fixed-point iter.
82       constTv    = 0.608d0,    & ! Constant for Virtual Temp Conversion
83       deltaz     = 2.d-13,     & ! Small number to ensure convergence in FPI
84       epsilon    = 1.d-25,     & ! Small number to aviod dividing by zero in wind calc
85       exponent = Rd*gamma/g,   & ! exponent
86       T0    = Ts0*(1.d0+constTv*q0),             & ! Surface temp
87       Ttrop = T0 - gamma*ztrop,                  & ! Tropopause temp
88       ptrop = p00*(Ttrop/T0)**(1.d0/exponent)      ! Tropopause pressure
89
90CONTAINS 
91
92!=======================================================================
93!    Evaluate the tropical cyclone initial conditions
94!=======================================================================
95  SUBROUTINE tropical_cyclone_test(lon,lat,p,z,zcoords,u,v,t,thetav,phis,ps,rho,q) &
96    BIND(c, name = "tropical_cyclone_test")
97    use iso_c_binding
98    IMPLICIT NONE
99
100    !------------------------------------------------
101    !   Input / output parameters
102    !------------------------------------------------
103
104    REAL(KIND=C_DOUBLE), INTENT(IN) ::     &
105              lon,             &     ! Longitude (radians)
106              lat                    ! Latitude (radians)
107
108    REAL(KIND=C_DOUBLE), INTENT(INOUT) ::  &
109              p,               &     ! Pressure (Pa)
110              z                      ! Height (m)
111
112    INTEGER(KIND=C_INT32_T), INTENT(IN) :: zcoords     ! 1 if z coordinates are specified
113                                     ! 0 if p coordinates are specified
114
115    REAL(KIND=C_DOUBLE), INTENT(OUT) ::    &
116              u,               &     ! Zonal wind (m s^-1)
117              v,               &     ! Meridional wind (m s^-1)
118              t,               &     ! Temperature (K)
119              thetav,          &     ! Virtual potential temperature (K)
120              phis,            &     ! Surface Geopotential (m^2 s^-2)
121              ps,              &     ! Surface Pressure (Pa)
122              rho,             &     ! Density (kg m^-3)
123              q                      ! Specific Humidity (kg/kg)
124
125    !------------------------------------------------
126    !   Local variables
127    !------------------------------------------------
128    real(8)  :: d1, d2, d, vfac, ufac, height, zhere, gr, f, zn
129
130    integer  n
131
132    !------------------------------------------------
133    !   Define Great circle distance (gr) and
134    !   Coriolis parameter (f)
135    !------------------------------------------------
136    f  = 2.d0*omega*sin(cen_lat*deg2rad)           ! Coriolis parameter
137    gr = a*acos(sin(cen_lat*deg2rad)*sin(lat) + &  ! Great circle radius
138         (cos(cen_lat*deg2rad)*cos(lat)*cos(lon-cen_lon*deg2rad)))
139
140    !------------------------------------------------
141    !   Initialize PS (surface pressure)
142    !------------------------------------------------
143    ps = p00-dp*exp(-(gr/rp)**exppr) 
144
145    !------------------------------------------------
146    !   Initialize altitude (z) if pressure provided
147    !   or pressure if altitude (z) is provided
148    !------------------------------------------------
149    if (zcoords .eq. 1) then
150
151       height = z
152 
153       if (height > ztrop) then
154          p = ptrop*exp(-(g*(height-ztrop))/(Rd*Ttrop))
155       else
156          p = (p00-dp*exp(-(gr/rp)**exppr)*exp(-(height/zp)**exppz)) &
157              * ((T0-gamma*height)/T0)**(1/exponent)
158       end if
159 
160    else
161
162       height = (T0/gamma)*(1.d0-(p/ps)**exponent)
163
164       ! If inside a certain distance of the center of the storm
165       ! perform a Fixed-point iteration to calculate the height
166       ! more accurately
167
168       if (gr < rfpi ) then
169          zhere = height 
170          n = 1
171          20 continue
172          n = n+1
173          zn = zhere - fpiF(p,gr,zhere)/fpidFdz(gr,zhere)
174          if (n.gt.20) then
175              PRINT *,'FPI did not converge after 20 interations in q & T!!!'
176          else if ( abs(zn-zhere)/abs(zn) > deltaz) then
177              zhere = zn
178              goto 20
179          end if
180          height = zn
181       end if
182    end if
183
184    !------------------------------------------------
185    !   Initialize U and V (wind components)
186    !------------------------------------------------
187    d1 = sin(cen_lat*deg2rad)*cos(lat) - &
188         cos(cen_lat*deg2rad)*sin(lat)*cos(lon-cen_lon*deg2rad)
189    d2 = cos(cen_lat*deg2rad)*sin(lon-cen_lon*deg2rad)
190    d  = max(epsilon, sqrt(d1**2.d0 + d2**2.d0))
191    ufac = d1/d
192    vfac = d2/d
193   
194    if (height > ztrop) then
195        u = 0.d0
196        v = 0.d0
197    else
198        v = vfac*(-f*gr/2.d0+sqrt((f*gr/2.d0)**(2.d0) &
199            - exppr*(gr/rp)**exppr*Rd*(T0-gamma*height) &
200            /(exppz*height*Rd*(T0-gamma*height)/(g*zp**exppz) &
201            +(1.d0-p00/dp*exp((gr/rp)**exppr)*exp((height/zp)**exppz)))))
202        u = ufac*(-f*gr/2.d0+sqrt((f*gr/2.d0)**(2.d0) &
203            - exppr*(gr/rp)**exppr*Rd*(T0-gamma*height) &
204            /(exppz*height*Rd*(T0-gamma*height)/(g*zp**exppz) &
205            +(1.d0-p00/dp*exp((gr/rp)**exppr)*exp((height/zp)**exppz)))))
206    end if
207
208    !------------------------------------------------
209    !   Initialize water vapor mixing ratio (q)
210    !------------------------------------------------
211    if (height > ztrop) then
212        q = qtrop
213    else
214        q = q0*exp(-height/zq1)*exp(-(height/zq2)**exppz)
215    end if
216
217    !------------------------------------------------
218    !   Initialize temperature (T)
219    !------------------------------------------------
220    if (height > ztrop) then
221        t = Ttrop
222    else
223        t = (T0-gamma*height)/(1.d0+constTv*q)/(1.d0+exppz*Rd*(T0-gamma*height)*height &
224            /(g*zp**exppz*(1.d0-p00/dp*exp((gr/rp)**exppr)*exp((height/zp)**exppz))))
225    end if
226
227    !-----------------------------------------------------
228    !   Initialize virtual potential temperature (thetav)
229    !-----------------------------------------------------
230    thetav = t * (1.d0+constTv*q) * (p0/p)**(Rd/cp)
231
232    !-----------------------------------------------------
233    !   Initialize surface geopotential (PHIS)
234    !-----------------------------------------------------
235    phis = 0.d0  ! constant
236
237    !-----------------------------------------------------
238    !   Initialize density (rho)
239    !-----------------------------------------------------
240    rho = p/(Rd*t*(1.d0+constTv*q))
241
242  END SUBROUTINE tropical_cyclone_test
243
244!-----------------------------------------------------------------------
245!    First function for fixed point iterations
246!-----------------------------------------------------------------------
247  REAL(8) FUNCTION fpiF(phere, gr, zhere)
248    IMPLICIT NONE
249    REAL(8), INTENT(IN) :: phere, gr, zhere
250
251      fpiF = phere-(p00-dp*exp(-(gr/rp)**exppr)*exp(-(zhere/zp)**exppz)) &
252             *((T0-gamma*zhere)/T0)**(g/(Rd*gamma))
253
254  END FUNCTION fpiF
255
256!-----------------------------------------------------------------------
257!    Second function for fixed point iterations
258!-----------------------------------------------------------------------
259  REAL(8) FUNCTION fpidFdz(gr, zhere) 
260    IMPLICIT NONE
261    REAL(8), INTENT(IN) :: gr, zhere
262
263      fpidFdz =-exppz*zhere*dp*exp(-(gr/rp)**exppr)*exp(-(zhere/zp)**exppz)/(zp*zp)*((T0-gamma*zhere)/T0)**(g/(Rd*gamma)) &
264               +g/(Rd*T0)*(p00-dp*exp(-(gr/rp)**exppr)*exp(-(zhere/zp)**exppz))*((T0-gamma*zhere)/T0)**(g/(Rd*gamma)-1.d0)
265
266  END FUNCTION fpidFdz
267
268END MODULE dcmip2016_cyclone_mod
269
Note: See TracBrowser for help on using the repository browser.