1 | MODULE physics_dcmip_mod |
---|
2 | USE ICOSA |
---|
3 | PRIVATE |
---|
4 | |
---|
5 | INTEGER,SAVE :: testcase |
---|
6 | PUBLIC init_physics, physics |
---|
7 | TYPE(t_field),POINTER :: f_out_i(:) |
---|
8 | TYPE(t_field),POINTER :: f_precl(:) |
---|
9 | REAL(rstd),POINTER :: out_i(:,:) |
---|
10 | |
---|
11 | CONTAINS |
---|
12 | |
---|
13 | SUBROUTINE init_physics |
---|
14 | IMPLICIT NONE |
---|
15 | |
---|
16 | testcase=1 |
---|
17 | CALL getin("dcmip_physics",testcase) |
---|
18 | CALL allocate_field(f_out_i,field_t,type_real,llm) |
---|
19 | CALL allocate_field(f_precl,field_t,type_real) |
---|
20 | |
---|
21 | END SUBROUTINE init_physics |
---|
22 | |
---|
23 | |
---|
24 | SUBROUTINE physics( it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) |
---|
25 | USE icosa |
---|
26 | IMPLICIT NONE |
---|
27 | INTEGER,INTENT(IN) :: it |
---|
28 | TYPE(t_field),POINTER :: f_phis(:) |
---|
29 | TYPE(t_field),POINTER :: f_ps(:) |
---|
30 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
31 | TYPE(t_field),POINTER :: f_ue(:) |
---|
32 | TYPE(t_field),POINTER :: f_q(:) |
---|
33 | |
---|
34 | REAL(rstd),POINTER :: phis(:) |
---|
35 | REAL(rstd),POINTER :: ps(:) |
---|
36 | REAL(rstd),POINTER :: theta_rhodz(:,:) |
---|
37 | REAL(rstd),POINTER :: ue(:,:) |
---|
38 | REAL(rstd),POINTER :: q(:,:,:) |
---|
39 | REAL(rstd),POINTER :: precl(:) |
---|
40 | INTEGER :: ind |
---|
41 | |
---|
42 | CALL transfert_request(f_ue,req_e1_vect) |
---|
43 | DO ind=1,ndomain |
---|
44 | CALL swap_dimensions(ind) |
---|
45 | CALL swap_geometry(ind) |
---|
46 | |
---|
47 | phis=f_phis(ind) |
---|
48 | ps=f_ps(ind) |
---|
49 | theta_rhodz=f_theta_rhodz(ind) |
---|
50 | ue=f_ue(ind) |
---|
51 | q=f_q(ind) |
---|
52 | out_i=f_out_i(ind) |
---|
53 | precl=f_precl(ind) |
---|
54 | |
---|
55 | CALL compute_physics( phis, ps, theta_rhodz, ue, q(:,:,1), precl) |
---|
56 | |
---|
57 | ENDDO |
---|
58 | |
---|
59 | ! CALL writefield("out_i",f_out_i) |
---|
60 | |
---|
61 | IF (mod(it,itau_out)==0 ) THEN |
---|
62 | CALL writefield("precl",f_precl) |
---|
63 | ENDIF |
---|
64 | |
---|
65 | END SUBROUTINE physics |
---|
66 | |
---|
67 | SUBROUTINE compute_physics(phis, ps, theta_rhodz, ue, q, precl) |
---|
68 | USE icosa |
---|
69 | USE pression_mod |
---|
70 | USE exner_mod |
---|
71 | USE theta2theta_rhodz_mod |
---|
72 | USE geopotential_mod |
---|
73 | USE wind_mod |
---|
74 | IMPLICIT NONE |
---|
75 | REAL(rstd) :: phis(iim*jjm) |
---|
76 | REAL(rstd) :: ps(iim*jjm) |
---|
77 | REAL(rstd) :: theta_rhodz(iim*jjm,llm) |
---|
78 | REAL(rstd) :: ue(3*iim*jjm,llm) |
---|
79 | REAL(rstd) :: q(iim*jjm,llm) |
---|
80 | REAL(rstd) :: precl(iim*jjm) |
---|
81 | |
---|
82 | REAL(rstd) :: p(iim*jjm,llm+1) |
---|
83 | REAL(rstd) :: pks(iim*jjm) |
---|
84 | REAL(rstd) :: pk(iim*jjm,llm) |
---|
85 | REAL(rstd) :: phi(iim*jjm,llm) |
---|
86 | REAL(rstd) :: T(iim*jjm,llm) |
---|
87 | REAL(rstd) :: Tfi(iim*jjm,llm) |
---|
88 | REAL(rstd) :: theta(iim*jjm,llm) |
---|
89 | |
---|
90 | REAL(rstd) :: uc(iim*jjm,3,llm) |
---|
91 | REAL(rstd) :: u(iim*jjm,llm) |
---|
92 | REAL(rstd) :: v(iim*jjm,llm) |
---|
93 | REAL(rstd) :: ufi(iim*jjm,llm) |
---|
94 | REAL(rstd) :: vfi(iim*jjm,llm) |
---|
95 | REAL(rstd) :: qfi(iim*jjm,llm) |
---|
96 | REAL(rstd) :: utemp(iim*jjm,llm) |
---|
97 | REAL(rstd) :: vtemp(iim*jjm,llm) |
---|
98 | REAL(rstd) :: lat(iim*jjm) |
---|
99 | REAL(rstd) :: lon |
---|
100 | REAL(rstd) :: pmid(iim*jjm,llm) |
---|
101 | REAL(rstd) :: pint(iim*jjm,llm+1) |
---|
102 | REAL(rstd) :: pdel(iim*jjm,llm) |
---|
103 | INTEGER :: i,j,l,ij |
---|
104 | |
---|
105 | PRINT *,'Entering in DCMIP physics' |
---|
106 | CALL compute_pression(ps,p,0) |
---|
107 | CALL compute_exner(ps,p,pks,pk,0) |
---|
108 | CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,0) |
---|
109 | CALL compute_geopotential(phis,pks,pk,theta,phi,0) |
---|
110 | CALL compute_theta_rhodz2temperature(ps,theta_rhodz,T,0) |
---|
111 | CALL compute_wind_centered(ue,uc) |
---|
112 | CALL compute_wind_centered_lonlat_compound(uc, u, v) |
---|
113 | |
---|
114 | DO j=jj_begin,jj_end |
---|
115 | DO i=ii_begin,ii_end |
---|
116 | ij=(j-1)*iim+i |
---|
117 | CALL xyz2lonlat(xyz_i(ij,:),lon,lat(ij)) |
---|
118 | ENDDO |
---|
119 | ENDDO |
---|
120 | |
---|
121 | DO l=1,llm+1 |
---|
122 | DO j=jj_begin,jj_end |
---|
123 | DO i=ii_begin,ii_end |
---|
124 | ij=(j-1)*iim+i |
---|
125 | pint(ij,l)=p(ij,llm+2-l) |
---|
126 | ENDDO |
---|
127 | ENDDO |
---|
128 | ENDDO |
---|
129 | |
---|
130 | DO l=1,llm |
---|
131 | DO j=jj_begin,jj_end |
---|
132 | DO i=ii_begin,ii_end |
---|
133 | ij=(j-1)*iim+i |
---|
134 | pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) |
---|
135 | ENDDO |
---|
136 | ENDDO |
---|
137 | ENDDO |
---|
138 | |
---|
139 | DO l=1,llm |
---|
140 | DO j=jj_begin,jj_end |
---|
141 | DO i=ii_begin,ii_end |
---|
142 | ij=(j-1)*iim+i |
---|
143 | pdel(ij,l)=pint(ij,l+1)-pint(ij,l) |
---|
144 | ENDDO |
---|
145 | ENDDO |
---|
146 | ENDDO |
---|
147 | |
---|
148 | |
---|
149 | ! ufi=u |
---|
150 | ! vfi=v |
---|
151 | |
---|
152 | DO l=1,llm |
---|
153 | DO j=jj_begin,jj_end |
---|
154 | DO i=ii_begin,ii_end |
---|
155 | ij=(j-1)*iim+i |
---|
156 | T(ij,l)=T(ij,l)/(1+0.608*q(ij,l)) |
---|
157 | ENDDO |
---|
158 | ENDDO |
---|
159 | ENDDO |
---|
160 | |
---|
161 | DO l=1,llm |
---|
162 | DO j=jj_begin,jj_end |
---|
163 | DO i=ii_begin,ii_end |
---|
164 | ij=(j-1)*iim+i |
---|
165 | Tfi(ij,l)=T(ij,llm+1-l) |
---|
166 | ufi(ij,l)=u(ij,llm+1-l) |
---|
167 | vfi(ij,l)=v(ij,llm+1-l) |
---|
168 | qfi(ij,l)=q(ij,llm+1-l) |
---|
169 | ENDDO |
---|
170 | ENDDO |
---|
171 | ENDDO |
---|
172 | |
---|
173 | ! q=0 |
---|
174 | ! out_i=T |
---|
175 | |
---|
176 | CALL simple_physics(iim*jjm, llm, dt, lat, tfi, qfi , ufi, vfi, pmid, pint, pdel, 1/pdel, ps, precl, testcase) |
---|
177 | |
---|
178 | DO l=1,llm |
---|
179 | DO j=jj_begin,jj_end |
---|
180 | DO i=ii_begin,ii_end |
---|
181 | ij=(j-1)*iim+i |
---|
182 | T(ij,l)=Tfi(ij,llm+1-l) |
---|
183 | utemp(ij,l)=ufi(ij,llm+1-l) |
---|
184 | vtemp(ij,l)=vfi(ij,llm+1-l) |
---|
185 | q(ij,l)=qfi(ij,llm+1-l) |
---|
186 | ENDDO |
---|
187 | ENDDO |
---|
188 | ENDDO |
---|
189 | |
---|
190 | |
---|
191 | DO l=1,llm |
---|
192 | DO j=jj_begin,jj_end |
---|
193 | DO i=ii_begin,ii_end |
---|
194 | ij=(j-1)*iim+i |
---|
195 | T(ij,l)=T(ij,l)*(1+0.608*q(ij,l)) |
---|
196 | ENDDO |
---|
197 | ENDDO |
---|
198 | ENDDO |
---|
199 | |
---|
200 | ! out_i=q |
---|
201 | |
---|
202 | utemp=utemp-u |
---|
203 | vtemp=vtemp-v |
---|
204 | |
---|
205 | DO l=1,llm |
---|
206 | DO j=jj_begin,jj_end |
---|
207 | DO i=ii_begin,ii_end |
---|
208 | ij=(j-1)*iim+i |
---|
209 | uc(ij,:,l)=utemp(ij,l)*elon_i(ij,:)+vtemp(ij,l)*elat_i(ij,:) |
---|
210 | ENDDO |
---|
211 | ENDDO |
---|
212 | ENDDO |
---|
213 | |
---|
214 | ! out_i=ufi |
---|
215 | |
---|
216 | DO l=1,llm |
---|
217 | DO j=jj_begin,jj_end |
---|
218 | DO i=ii_begin,ii_end |
---|
219 | ij=(j-1)*iim+i |
---|
220 | ue(ij+u_right,l)=ue(ij+u_right,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) |
---|
221 | ue(ij+u_lup,l)=ue(ij+u_lup,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) ) |
---|
222 | ue(ij+u_ldown,l)=ue(ij+u_ldown,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) ) |
---|
223 | ENDDO |
---|
224 | ENDDO |
---|
225 | ENDDO |
---|
226 | |
---|
227 | CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) |
---|
228 | |
---|
229 | |
---|
230 | END SUBROUTINE compute_physics |
---|
231 | |
---|
232 | |
---|
233 | subroutine simple_physics (pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, rpdel, ps, precl, test) |
---|
234 | !----------------------------------------------------------------------- |
---|
235 | ! |
---|
236 | ! Purpose: Simple Physics Package |
---|
237 | ! |
---|
238 | ! Author: K. A. Reed (University of Michigan, kareed@umich.edu) |
---|
239 | ! version 5 |
---|
240 | ! July/8/2012 |
---|
241 | ! |
---|
242 | ! Change log: |
---|
243 | ! v2: removal of some NCAR CAM-specific 'use' associations |
---|
244 | ! v3: corrected precl(i) computation, the precipitation rate is now computed via a vertical integral, the previous single-level computation in v2 was a bug |
---|
245 | ! v3: corrected dtdt(i,1) computation, the term '-(i,1)' was missing the temperature variable: '-t(i,1)' |
---|
246 | ! v4: modified and enhanced parameter list to make the routine truly standalone, the number of columns and vertical levels have been added: pcols, pver |
---|
247 | ! v4: 'ncol' has been removed, 'pcols' is used instead |
---|
248 | ! v5: the sea surface temperature (SST) field Tsurf is now an array, the SST now depends on the latitude |
---|
249 | ! v5: addition of the latitude array 'lat' and the flag 'test' in the parameter list |
---|
250 | ! if test = 0: constant SST is used, correct setting for the tropical cyclone test case 5-1 |
---|
251 | ! if test = 1: newly added latitude-dependent SST is used, correct setting for the moist baroclinic wave test with simple-physics (test 4-3) |
---|
252 | ! |
---|
253 | ! Description: Includes large-scale precipitation, surface fluxes and |
---|
254 | ! boundary-leyer mixing. The processes are time-split |
---|
255 | ! in that order. A partially implicit formulation is |
---|
256 | ! used to foster numerical stability. |
---|
257 | ! The routine assumes that the model levels are ordered |
---|
258 | ! in a top-down approach, e.g. level 1 denotes the uppermost |
---|
259 | ! full model level. |
---|
260 | ! |
---|
261 | ! This routine is based on an implementation which was |
---|
262 | ! developed for the NCAR Community Atmosphere Model (CAM). |
---|
263 | ! Adjustments for other models will be necessary. |
---|
264 | ! |
---|
265 | ! The routine provides both updates of the state variables |
---|
266 | ! u, v, T, q (these are local copies of u,v,T,q within this physics |
---|
267 | ! routine) and also collects their time tendencies. |
---|
268 | ! The latter might be used to couple the physics and dynamics |
---|
269 | ! in a process-split way. For a time-split coupling, the final |
---|
270 | ! state should be given to the dynamical core for the next time step. |
---|
271 | ! Test: 0 = Reed and Jablonowski (2011) tropical cyclone test case (test 5-1) |
---|
272 | ! 1 = Moist baroclinic instability test (test 4-3) |
---|
273 | ! |
---|
274 | ! |
---|
275 | ! Reference: Reed, K. A. and C. Jablonowski (2012), Idealized tropical cyclone |
---|
276 | ! simulations of intermediate complexity: A test case for AGCMs, |
---|
277 | ! J. Adv. Model. Earth Syst., Vol. 4, M04001, doi:10.1029/2011MS000099 |
---|
278 | !----------------------------------------------------------------------- |
---|
279 | ! use physics_types , only: physics_dme_adjust ! This is for CESM/CAM |
---|
280 | ! use cam_diagnostics, only: diag_phys_writeout ! This is for CESM/CAM |
---|
281 | |
---|
282 | implicit none |
---|
283 | |
---|
284 | integer, parameter :: r8 = selected_real_kind(12) |
---|
285 | |
---|
286 | ! |
---|
287 | ! Input arguments - MODEL DEPENDENT |
---|
288 | ! |
---|
289 | integer, intent(in) :: pcols ! Set number of atmospheric columns |
---|
290 | integer, intent(in) :: pver ! Set number of model levels |
---|
291 | real(r8), intent(in) :: dtime ! Set model physics timestep |
---|
292 | real(r8), intent(in) :: lat(pcols) ! Latitude |
---|
293 | integer, intent(in) :: test ! Test number |
---|
294 | |
---|
295 | ! |
---|
296 | ! Input/Output arguments |
---|
297 | ! |
---|
298 | ! pcols is the maximum number of vertical columns per 'chunk' of atmosphere |
---|
299 | ! |
---|
300 | real(r8), intent(inout) :: t(pcols,pver) ! Temperature at full-model level (K) |
---|
301 | real(r8), intent(inout) :: q(pcols,pver) ! Specific Humidity at full-model level (kg/kg) |
---|
302 | real(r8), intent(inout) :: u(pcols,pver) ! Zonal wind at full-model level (m/s) |
---|
303 | real(r8), intent(inout) :: v(pcols,pver) ! Meridional wind at full-model level (m/s) |
---|
304 | real(r8), intent(inout) :: pmid(pcols,pver) ! Pressure is full-model level (Pa) |
---|
305 | real(r8), intent(inout) :: pint(pcols,pver+1) ! Pressure at model interfaces (Pa) |
---|
306 | real(r8), intent(inout) :: pdel(pcols,pver) ! Layer thickness (Pa) |
---|
307 | real(r8), intent(in) :: rpdel(pcols,pver) ! Reciprocal of layer thickness (1/Pa) |
---|
308 | real(r8), intent(inout) :: ps(pcols) ! Surface Pressue (Pa) |
---|
309 | |
---|
310 | ! |
---|
311 | ! Output arguments |
---|
312 | ! |
---|
313 | real(r8), intent(out) :: precl(pcols) ! Precipitation rate (m_water / s) |
---|
314 | |
---|
315 | ! |
---|
316 | !---------------------------Local workspace----------------------------- |
---|
317 | ! |
---|
318 | |
---|
319 | ! Integers for loops |
---|
320 | |
---|
321 | integer i,k ! Longitude, level indices |
---|
322 | |
---|
323 | ! Physical Constants - Many of these may be model dependent |
---|
324 | |
---|
325 | real(r8) gravit ! Gravity |
---|
326 | real(r8) rair ! Gas constant for dry air |
---|
327 | real(r8) cpair ! Specific heat of dry air |
---|
328 | real(r8) latvap ! Latent heat of vaporization |
---|
329 | real(r8) rh2o ! Gas constant for water vapor |
---|
330 | real(r8) epsilo ! Ratio of gas constant for dry air to that for vapor |
---|
331 | real(r8) zvir ! Constant for virtual temp. calc. =(rh2o/rair) - 1 |
---|
332 | real(r8) a ! Reference Earth's Radius (m) |
---|
333 | real(r8) omega ! Reference rotation rate of the Earth (s^-1) |
---|
334 | real(r8) pi ! pi |
---|
335 | |
---|
336 | ! Simple Physics Specific Constants |
---|
337 | |
---|
338 | !++++++++ |
---|
339 | real(r8) Tsurf(pcols) ! Sea Surface Temperature (constant for tropical cyclone) |
---|
340 | !++++++++ Tsurf needs to be dependent on latitude for the |
---|
341 | ! moist baroclinic wave test 4-3 with simple-physics, adjust |
---|
342 | |
---|
343 | real(r8) SST_tc ! Sea Surface Temperature for tropical cyclone test |
---|
344 | real(r8) T0 ! Control temp for calculation of qsat |
---|
345 | real(r8) e0 ! Saturation vapor pressure at T0 for calculation of qsat |
---|
346 | real(r8) rhow ! Density of Liquid Water |
---|
347 | real(r8) p0 ! Constant for calculation of potential temperature |
---|
348 | real(r8) Cd0 ! Constant for calculating Cd from Smith and Vogl 2008 |
---|
349 | real(r8) Cd1 ! Constant for calculating Cd from Smith and Vogl 2008 |
---|
350 | real(r8) Cm ! Constant for calculating Cd from Smith and Vogl 2008 |
---|
351 | real(r8) v20 ! Threshold wind speed for calculating Cd from Smith and Vogl 2008 |
---|
352 | real(r8) C ! Drag coefficient for sensible heat and evaporation |
---|
353 | real(r8) T00 ! Horizontal mean T at surface for moist baro test |
---|
354 | real(r8) u0 ! Zonal wind constant for moist baro test |
---|
355 | real(r8) latw ! halfwidth for for baro test |
---|
356 | real(r8) eta0 ! Center of jets (hybrid) for baro test |
---|
357 | real(r8) etav ! Auxiliary variable for baro test |
---|
358 | real(r8) q0 ! Maximum specific humidity for baro test |
---|
359 | |
---|
360 | ! Physics Tendency Arrays |
---|
361 | real(r8) dtdt(pcols,pver) ! Temperature tendency |
---|
362 | real(r8) dqdt(pcols,pver) ! Specific humidity tendency |
---|
363 | real(r8) dudt(pcols,pver) ! Zonal wind tendency |
---|
364 | real(r8) dvdt(pcols,pver) ! Meridional wind tendency |
---|
365 | |
---|
366 | ! Temporary variables for tendency calculations |
---|
367 | |
---|
368 | real(r8) tmp ! Temporary |
---|
369 | real(r8) qsat ! Saturation vapor pressure |
---|
370 | real(r8) qsats ! Saturation vapor pressure of SST |
---|
371 | |
---|
372 | ! Variables for Boundary Layer Calculation |
---|
373 | |
---|
374 | real(r8) wind(pcols) ! Magnitude of Wind |
---|
375 | real(r8) Cd(pcols) ! Drag coefficient for momentum |
---|
376 | real(r8) Km(pcols,pver+1) ! Eddy diffusivity for boundary layer calculations |
---|
377 | real(r8) Ke(pcols,pver+1) ! Eddy diffusivity for boundary layer calculations |
---|
378 | real(r8) rho ! Density at lower/upper interface |
---|
379 | real(r8) za(pcols) ! Heights at midpoints of first model level |
---|
380 | real(r8) dlnpint ! Used for calculation of heights |
---|
381 | real(r8) pbltop ! Top of boundary layer |
---|
382 | real(r8) pblconst ! Constant for the calculation of the decay of diffusivity |
---|
383 | real(r8) CA(pcols,pver) ! Matrix Coefficents for PBL Scheme |
---|
384 | real(r8) CC(pcols,pver) ! Matrix Coefficents for PBL Scheme |
---|
385 | real(r8) CE(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
---|
386 | real(r8) CAm(pcols,pver) ! Matrix Coefficents for PBL Scheme |
---|
387 | real(r8) CCm(pcols,pver) ! Matrix Coefficents for PBL Scheme |
---|
388 | real(r8) CEm(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
---|
389 | real(r8) CFu(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
---|
390 | real(r8) CFv(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
---|
391 | real(r8) CFt(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
---|
392 | real(r8) CFq(pcols,pver+1) ! Matrix Coefficents for PBL Scheme |
---|
393 | |
---|
394 | |
---|
395 | ! Variable for Dry Mass Adjustment, this dry air adjustment is necessary to |
---|
396 | ! conserve the mass of the dry air |
---|
397 | |
---|
398 | real(r8) qini(pcols,pver) ! Initial specific humidity |
---|
399 | |
---|
400 | !=============================================================================== |
---|
401 | ! |
---|
402 | ! Physical Constants - MAY BE MODEL DEPENDENT |
---|
403 | ! |
---|
404 | !=============================================================================== |
---|
405 | gravit = 9.80616_r8 ! Gravity (9.80616 m/s^2) |
---|
406 | rair = 287.0_r8 ! Gas constant for dry air: 287 J/(kg K) |
---|
407 | cpair = 1.0045e3_r8 ! Specific heat of dry air: here we use 1004.5 J/(kg K) |
---|
408 | latvap = 2.5e6_r8 ! Latent heat of vaporization (J/kg) |
---|
409 | rh2o = 461.5_r8 ! Gas constant for water vapor: 461.5 J/(kg K) |
---|
410 | epsilo = rair/rh2o ! Ratio of gas constant for dry air to that for vapor |
---|
411 | zvir = (rh2o/rair) - 1._r8 ! Constant for virtual temp. calc. =(rh2o/rair) - 1 is approx. 0.608 |
---|
412 | a = 6371220.0_r8 ! Reference Earth's Radius (m) |
---|
413 | omega = 7.29212d-5 ! Reference rotation rate of the Earth (s^-1) |
---|
414 | pi = 4._r8*atan(1._r8) ! pi |
---|
415 | |
---|
416 | !=============================================================================== |
---|
417 | ! |
---|
418 | ! Local Constants for Simple Physics |
---|
419 | ! |
---|
420 | !=============================================================================== |
---|
421 | C = 0.0011_r8 ! From Smith and Vogl 2008 |
---|
422 | SST_tc = 302.15_r8 ! Constant Value for SST for tropical cyclone test |
---|
423 | T0 = 273.16_r8 ! control temp for calculation of qsat |
---|
424 | e0 = 610.78_r8 ! saturation vapor pressure at T0 for calculation of qsat |
---|
425 | rhow = 1000.0_r8 ! Density of Liquid Water |
---|
426 | Cd0 = 0.0007_r8 ! Constant for Cd calc. Smith and Vogl 2008 |
---|
427 | Cd1 = 0.000065_r8 ! Constant for Cd calc. Smith and Vogl 2008 |
---|
428 | Cm = 0.002_r8 ! Constant for Cd calc. Smith and Vogl 2008 |
---|
429 | v20 = 20.0_r8 ! Threshold wind speed for calculating Cd from Smith and Vogl 2008 |
---|
430 | p0 = 100000.0_r8 ! Constant for potential temp calculation |
---|
431 | pbltop = 85000._r8 ! Top of boundary layer |
---|
432 | pblconst = 10000._r8 ! Constant for the calculation of the decay of diffusivity |
---|
433 | T00 = 288.0_r8 ! Horizontal mean T at surface for moist baro test |
---|
434 | u0 = 35.0_r8 ! Zonal wind constant for moist baro test |
---|
435 | latw = 2.0_r8*pi/9.0_r8 ! Halfwidth for for baro test |
---|
436 | eta0 = 0.252_r8 ! Center of jets (hybrid) for baro test |
---|
437 | etav = (1._r8-eta0)*0.5_r8*pi ! Auxiliary variable for baro test |
---|
438 | q0 = 0.021 ! Maximum specific humidity for baro test |
---|
439 | |
---|
440 | !=============================================================================== |
---|
441 | ! |
---|
442 | ! Definition of local arrays |
---|
443 | ! |
---|
444 | !=============================================================================== |
---|
445 | ! |
---|
446 | ! Calculate hydrostatic height za of the lowest model level |
---|
447 | ! |
---|
448 | do i=1,pcols |
---|
449 | dlnpint = log(ps(i)) - log(pint(i,pver)) ! ps(i) is identical to pint(i,pver+1), note: this is the correct sign (corrects typo in JAMES paper) |
---|
450 | za(i) = rair/gravit*t(i,pver)*(1._r8+zvir*q(i,pver))*0.5_r8*dlnpint |
---|
451 | end do |
---|
452 | ! |
---|
453 | ! Set Initial Specific Humidity - For dry mass adjustment at the end |
---|
454 | ! |
---|
455 | qini(:pcols,:pver) = q(:pcols,:pver) |
---|
456 | !-------------------------------------------------------------- |
---|
457 | ! Set Sea Surface Temperature (constant for tropical cyclone) |
---|
458 | ! Tsurf needs to be dependent on latitude for the |
---|
459 | ! moist baroclinic wave test 4-3 with simple-physics |
---|
460 | !-------------------------------------------------------------- |
---|
461 | if (test .eq. 1) then ! moist baroclinic wave with simple-physics |
---|
462 | do i=1,pcols |
---|
463 | Tsurf(i) = (T00 + pi*u0/rair * 1.5_r8 * sin(etav) * (cos(etav))**0.5_r8 * & |
---|
464 | ((-2._r8*(sin(lat(i)))**6 * ((cos(lat(i)))**2 + 1._r8/3._r8) + 10._r8/63._r8)* & |
---|
465 | u0 * (cos(etav))**1.5_r8 + & |
---|
466 | (8._r8/5._r8*(cos(lat(i)))**3 * ((sin(lat(i)))**2 + 2._r8/3._r8) - pi/4._r8)*a*omega*0.5_r8 ))/ & |
---|
467 | (1._r8+zvir*q0*exp(-(lat(i)/latw)**4)) |
---|
468 | |
---|
469 | end do |
---|
470 | else |
---|
471 | do i=1,pcols ! constant SST for the tropical cyclone test case |
---|
472 | Tsurf(i) = SST_tc |
---|
473 | end do |
---|
474 | end if |
---|
475 | |
---|
476 | !=============================================================================== |
---|
477 | ! |
---|
478 | ! Set initial physics time tendencies and precipitation field to zero |
---|
479 | ! |
---|
480 | !=============================================================================== |
---|
481 | dtdt(:pcols,:pver) = 0._r8 ! initialize temperature tendency with zero |
---|
482 | dqdt(:pcols,:pver) = 0._r8 ! initialize specific humidity tendency with zero |
---|
483 | dudt(:pcols,:pver) = 0._r8 ! initialize zonal wind tendency with zero |
---|
484 | dvdt(:pcols,:pver) = 0._r8 ! initialize meridional wind tendency with zero |
---|
485 | precl(:pcols) = 0._r8 ! initialize precipitation rate with zero |
---|
486 | |
---|
487 | !=============================================================================== |
---|
488 | ! |
---|
489 | ! Large-Scale Condensation and Precipitation Rate |
---|
490 | ! |
---|
491 | !=============================================================================== |
---|
492 | ! |
---|
493 | ! Calculate Tendencies |
---|
494 | ! |
---|
495 | do k=1,pver |
---|
496 | do i=1,pcols |
---|
497 | qsat = epsilo*e0/pmid(i,k)*exp(-latvap/rh2o*((1._r8/t(i,k))-1._r8/T0)) ! saturation specific humidity |
---|
498 | out_i(i,llm+1-k)=q(i,k)-qsat |
---|
499 | if (q(i,k) > qsat) then ! saturated? |
---|
500 | tmp = 1._r8/dtime*(q(i,k)-qsat)/(1._r8+(latvap/cpair)*(epsilo*latvap*qsat/(rair*t(i,k)**2))) |
---|
501 | dtdt(i,k) = dtdt(i,k)+latvap/cpair*tmp |
---|
502 | dqdt(i,k) = dqdt(i,k)-tmp |
---|
503 | precl(i) = precl(i) + tmp*pdel(i,k)/(gravit*rhow) ! precipitation rate, computed via a vertical integral |
---|
504 | ! corrected in version 1.3 |
---|
505 | end if |
---|
506 | end do |
---|
507 | end do |
---|
508 | ! |
---|
509 | ! Update moisture and temperature fields from Large-Scale Precipitation Scheme |
---|
510 | ! |
---|
511 | do k=1,pver |
---|
512 | do i=1,pcols |
---|
513 | t(i,k) = t(i,k) + dtdt(i,k)*dtime ! update the state variables T and q |
---|
514 | q(i,k) = q(i,k) + dqdt(i,k)*dtime |
---|
515 | end do |
---|
516 | end do |
---|
517 | |
---|
518 | IF (test==0) return |
---|
519 | !=============================================================================== |
---|
520 | ! Send variables to history file - THIS PROCESS WILL BE MODEL SPECIFIC |
---|
521 | ! |
---|
522 | ! note: The variables, as done in many GCMs, are written to the history file |
---|
523 | ! after the moist physics process. This ensures that the moisture fields |
---|
524 | ! are somewhat in equilibrium. |
---|
525 | !=============================================================================== |
---|
526 | ! call diag_phys_writeout(state) ! This is for CESM/CAM |
---|
527 | |
---|
528 | !=============================================================================== |
---|
529 | ! |
---|
530 | ! Turbulent mixing coefficients for the PBL mixing of horizontal momentum, |
---|
531 | ! sensible heat and latent heat |
---|
532 | ! |
---|
533 | ! We are using Simplified Ekman theory to compute the diffusion coefficients |
---|
534 | ! Kx for the boundary-layer mixing. The Kx values are calculated at each time step |
---|
535 | ! and in each column. |
---|
536 | ! |
---|
537 | !=============================================================================== |
---|
538 | ! |
---|
539 | ! Compute magnitude of the wind and drag coeffcients for turbulence scheme: |
---|
540 | ! they depend on the conditions at the lowest model level and stay constant |
---|
541 | ! up to the 850 hPa level. Above this level the coefficients are decreased |
---|
542 | ! and tapered to zero. At the 700 hPa level the strength of the K coefficients |
---|
543 | ! is about 10% of the maximum strength. |
---|
544 | ! |
---|
545 | do i=1,pcols |
---|
546 | wind(i) = sqrt(u(i,pver)**2+v(i,pver)**2) ! wind magnitude at the lowest level |
---|
547 | end do |
---|
548 | do i=1,pcols |
---|
549 | Ke(i,pver+1) = C*wind(i)*za(i) |
---|
550 | if( wind(i) .lt. v20) then |
---|
551 | Cd(i) = Cd0+Cd1*wind(i) |
---|
552 | Km(i,pver+1) = Cd(i)*wind(i)*za(i) |
---|
553 | else |
---|
554 | Cd(i) = Cm |
---|
555 | Km(i,pver+1) = Cm*wind(i)*za(i) |
---|
556 | endif |
---|
557 | end do |
---|
558 | |
---|
559 | do k=1,pver |
---|
560 | do i=1,pcols |
---|
561 | if( pint(i,k) .ge. pbltop) then |
---|
562 | Km(i,k) = Km(i,pver+1) ! constant Km below 850 hPa level |
---|
563 | Ke(i,k) = Ke(i,pver+1) ! constant Ke below 850 hPa level |
---|
564 | else |
---|
565 | Km(i,k) = Km(i,pver+1)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2) ! Km tapered to 0 |
---|
566 | Ke(i,k) = Ke(i,pver+1)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2) ! Ke tapered to 0 |
---|
567 | end if |
---|
568 | end do |
---|
569 | end do |
---|
570 | |
---|
571 | |
---|
572 | !=============================================================================== |
---|
573 | ! Update the state variables u, v, t, q with the surface fluxes at the |
---|
574 | ! lowest model level, this is done with an implicit approach |
---|
575 | ! see Reed and Jablonowski (JAMES, 2012) |
---|
576 | ! |
---|
577 | ! Sea Surface Temperature Tsurf is constant for tropical cyclone test 5-1 |
---|
578 | ! Tsurf needs to be dependent on latitude for the |
---|
579 | ! moist baroclinic wave test 4-3 with simple-physics, adjust |
---|
580 | !=============================================================================== |
---|
581 | |
---|
582 | do i=1,pcols |
---|
583 | qsats = epsilo*e0/ps(i)*exp(-latvap/rh2o*((1._r8/Tsurf(i))-1._r8/T0)) ! saturation specific humidity at the surface |
---|
584 | dudt(i,pver) = dudt(i,pver) + (u(i,pver) & |
---|
585 | /(1._r8+Cd(i)*wind(i)*dtime/za(i))-u(i,pver))/dtime |
---|
586 | dvdt(i,pver) = dvdt(i,pver) + (v(i,pver) & |
---|
587 | /(1._r8+Cd(i)*wind(i)*dtime/za(i))-v(i,pver))/dtime |
---|
588 | u(i,pver) = u(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i)) |
---|
589 | v(i,pver) = v(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i)) |
---|
590 | dtdt(i,pver) = dtdt(i,pver) +((t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i)) & |
---|
591 | /(1._r8+C*wind(i)*dtime/za(i))-t(i,pver))/dtime |
---|
592 | t(i,pver) = (t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i)) & |
---|
593 | /(1._r8+C*wind(i)*dtime/za(i)) |
---|
594 | dqdt(i,pver) = dqdt(i,pver) +((q(i,pver)+C*wind(i)*qsats*dtime/za(i)) & |
---|
595 | /(1._r8+C*wind(i)*dtime/za(i))-q(i,pver))/dtime |
---|
596 | q(i,pver) = (q(i,pver)+C*wind(i)*qsats*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i)) |
---|
597 | end do |
---|
598 | !=============================================================================== |
---|
599 | |
---|
600 | |
---|
601 | !=============================================================================== |
---|
602 | ! Boundary layer mixing, see Reed and Jablonowski (JAMES, 2012) |
---|
603 | !=============================================================================== |
---|
604 | ! Calculate Diagonal Variables for Implicit PBL Scheme |
---|
605 | ! |
---|
606 | do k=1,pver-1 |
---|
607 | do i=1,pcols |
---|
608 | rho = (pint(i,k+1)/(rair*(t(i,k+1)+t(i,k))/2.0_r8)) |
---|
609 | CAm(i,k) = rpdel(i,k)*dtime*gravit*gravit*Km(i,k+1)*rho*rho & |
---|
610 | /(pmid(i,k+1)-pmid(i,k)) |
---|
611 | CCm(i,k+1) = rpdel(i,k+1)*dtime*gravit*gravit*Km(i,k+1)*rho*rho & |
---|
612 | /(pmid(i,k+1)-pmid(i,k)) |
---|
613 | CA(i,k) = rpdel(i,k)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho & |
---|
614 | /(pmid(i,k+1)-pmid(i,k)) |
---|
615 | CC(i,k+1) = rpdel(i,k+1)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho & |
---|
616 | /(pmid(i,k+1)-pmid(i,k)) |
---|
617 | end do |
---|
618 | end do |
---|
619 | do i=1,pcols |
---|
620 | CAm(i,pver) = 0._r8 |
---|
621 | CCm(i,1) = 0._r8 |
---|
622 | CEm(i,pver+1) = 0._r8 |
---|
623 | CA(i,pver) = 0._r8 |
---|
624 | CC(i,1) = 0._r8 |
---|
625 | CE(i,pver+1) = 0._r8 |
---|
626 | CFu(i,pver+1) = 0._r8 |
---|
627 | CFv(i,pver+1) = 0._r8 |
---|
628 | CFt(i,pver+1) = 0._r8 |
---|
629 | CFq(i,pver+1) = 0._r8 |
---|
630 | end do |
---|
631 | do i=1,pcols |
---|
632 | do k=pver,1,-1 |
---|
633 | CE(i,k) = CC(i,k)/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) |
---|
634 | CEm(i,k) = CCm(i,k)/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1)) |
---|
635 | CFu(i,k) = (u(i,k)+CAm(i,k)*CFu(i,k+1)) & |
---|
636 | /(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1)) |
---|
637 | CFv(i,k) = (v(i,k)+CAm(i,k)*CFv(i,k+1)) & |
---|
638 | /(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1)) |
---|
639 | CFt(i,k) = ((p0/pmid(i,k))**(rair/cpair)*t(i,k)+CA(i,k)*CFt(i,k+1)) & |
---|
640 | /(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) |
---|
641 | CFq(i,k) = (q(i,k)+CA(i,k)*CFq(i,k+1)) & |
---|
642 | /(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) |
---|
643 | end do |
---|
644 | end do |
---|
645 | |
---|
646 | ! |
---|
647 | ! Calculate the updated temperature, specific humidity and horizontal wind |
---|
648 | ! |
---|
649 | ! First we need to calculate the updates at the top model level |
---|
650 | ! |
---|
651 | do i=1,pcols |
---|
652 | dudt(i,1) = dudt(i,1)+(CFu(i,1)-u(i,1))/dtime |
---|
653 | dvdt(i,1) = dvdt(i,1)+(CFv(i,1)-v(i,1))/dtime |
---|
654 | u(i,1) = CFu(i,1) |
---|
655 | v(i,1) = CFv(i,1) |
---|
656 | dtdt(i,1) = dtdt(i,1)+(CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair)-t(i,1))/dtime ! corrected in version 1.3 |
---|
657 | t(i,1) = CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair) |
---|
658 | dqdt(i,1) = dqdt(i,1)+(CFq(i,1)-q(i,1))/dtime |
---|
659 | q(i,1) = CFq(i,1) |
---|
660 | end do |
---|
661 | ! |
---|
662 | ! Loop over the remaining level |
---|
663 | ! |
---|
664 | do i=1,pcols |
---|
665 | do k=2,pver |
---|
666 | dudt(i,k) = dudt(i,k)+(CEm(i,k)*u(i,k-1)+CFu(i,k)-u(i,k))/dtime |
---|
667 | dvdt(i,k) = dvdt(i,k)+(CEm(i,k)*v(i,k-1)+CFv(i,k)-v(i,k))/dtime |
---|
668 | u(i,k) = CEm(i,k)*u(i,k-1)+CFu(i,k) |
---|
669 | v(i,k) = CEm(i,k)*v(i,k-1)+CFv(i,k) |
---|
670 | dtdt(i,k) = dtdt(i,k)+((CE(i,k)*t(i,k-1) & |
---|
671 | *(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k)) & |
---|
672 | *(pmid(i,k)/p0)**(rair/cpair)-t(i,k))/dtime |
---|
673 | t(i,k) = (CE(i,k)*t(i,k-1)*(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k)) & |
---|
674 | *(pmid(i,k)/p0)**(rair/cpair) |
---|
675 | dqdt(i,k) = dqdt(i,k)+(CE(i,k)*q(i,k-1)+CFq(i,k)-q(i,k))/dtime |
---|
676 | q(i,k) = CE(i,k)*q(i,k-1)+CFq(i,k) |
---|
677 | end do |
---|
678 | end do |
---|
679 | |
---|
680 | !=============================================================================== |
---|
681 | ! |
---|
682 | ! Dry Mass Adjustment - THIS PROCESS WILL BE MODEL SPECIFIC |
---|
683 | ! |
---|
684 | ! note: Care needs to be taken to ensure that the model conserves the total |
---|
685 | ! dry air mass. Add your own routine here. |
---|
686 | !=============================================================================== |
---|
687 | ! call physics_dme_adjust(state, tend, qini, dtime) ! This is for CESM/CAM |
---|
688 | |
---|
689 | return |
---|
690 | end subroutine simple_physics |
---|
691 | |
---|
692 | |
---|
693 | |
---|
694 | |
---|
695 | END MODULE physics_dcmip_mod |
---|