Line data Source code
1 : module stepon
2 :
3 : !----------------------------------------------------------------------
4 : ! stepon provides the interface layer that allows the different dynamical
5 : ! cores to be called from different locations in the time loop. It also
6 : ! provides a standard interface that is called from the higher level CAM
7 : ! component run methods while leaving non-standardized dycore interface
8 : ! methods to be called from this layer. Ideally only the run methods
9 : ! which allow flexibility in the dynamics/physics calling sequence should
10 : ! remain. The init and finalize methods should be removed and their
11 : ! functionality incorporated in the dycore init and finalize.
12 : !----------------------------------------------------------------------
13 :
14 : use shr_kind_mod, only: r8 => shr_kind_r8
15 :
16 : use spmd_utils, only: mpicom, iam, masterproc
17 : use cam_control_mod, only: initial_run, moist_physics
18 : use ppgrid, only: begchunk, endchunk
19 : use physconst, only: zvir, cappa
20 :
21 : use physics_types, only: physics_state, physics_tend
22 :
23 : use dyn_comp, only: dyn_import_t, dyn_export_t, initial_mr
24 : use dynamics_vars, only: t_fvdycore_state, t_fvdycore_grid
25 : use dyn_internal_state, only: get_dyn_state, get_dyn_state_grid
26 :
27 : use cam_logfile, only: iulog
28 : use cam_abortutils, only: endrun
29 : use perf_mod, only: t_startf, t_stopf, t_barrierf
30 :
31 : use aerosol_properties_mod, only: aerosol_properties
32 : use aerosol_state_mod, only: aerosol_state
33 : use microp_aero, only: aerosol_state_object, aerosol_properties_object
34 :
35 : implicit none
36 : private
37 : save
38 :
39 : public :: &
40 : stepon_init, &! Initialization
41 : stepon_run1, &! run method phase 1
42 : stepon_run2, &! run method phase 2
43 : stepon_run3, &! run method phase 3
44 : stepon_final ! Finalization
45 :
46 : integer :: pdt ! Physics time step
47 : real(r8) :: dtime ! Physics time step
48 : real(r8) :: te0 ! Total energy before dynamics
49 :
50 : ! for fv_out
51 : logical, parameter :: fv_monitor=.true. ! Monitor Mean/Max/Min fields
52 : ! This is CPU-time comsuming;
53 : ! set it to false for production runs
54 : real (r8) :: ptop
55 :
56 : class(aerosol_properties), pointer :: aero_props_obj => null()
57 : logical :: aerosols_transported = .false.
58 :
59 : !=========================================================================================
60 : contains
61 : !=========================================================================================
62 :
63 1536 : subroutine stepon_init(dyn_in, dyn_out)
64 :
65 : use constituents, only: pcnst
66 : use time_manager, only: get_step_size
67 : use physconst, only: rair, cpair
68 : use cam_thermo, only: cam_thermo_calc_kappav
69 : use inic_analytic, only: analytic_ic_active
70 : use cam_initfiles, only: scale_dry_air_mass
71 :
72 : type (dyn_import_t) :: dyn_in ! Dynamics import container
73 : type (dyn_export_t) :: dyn_out ! Dynamics export container
74 :
75 : ! local variables:
76 : type (t_fvdycore_grid), pointer :: grid
77 :
78 : integer :: im, km
79 : integer :: ifirstxy, ilastxy, jfirstxy, jlastxy
80 : integer :: i,k,j,m ! longitude, level, latitude and tracer indices
81 : logical :: nlres = .false. ! true => restart or branch run
82 :
83 : integer :: ks
84 1536 : real (r8), pointer :: ak(:)
85 1536 : real (r8), pointer :: bk(:)
86 :
87 1536 : real(r8), allocatable :: delpdryxy(:,:,:)
88 1536 : real(r8), allocatable :: cap3vi(:,:,:), cappa3v(:,:,:)
89 :
90 : !----------------------------------------------------------------------------
91 :
92 768 : if (.not. initial_run) nlres=.true.
93 :
94 1536 : grid => get_dyn_state_grid()
95 1536 : im = grid%im
96 1536 : km = grid%km
97 :
98 :
99 1536 : ifirstxy = grid%ifirstxy
100 1536 : ilastxy = grid%ilastxy
101 1536 : jfirstxy = grid%jfirstxy
102 1536 : jlastxy = grid%jlastxy
103 :
104 1536 : ks = grid%ks
105 1536 : ptop = grid%ptop
106 1536 : ak => grid%ak
107 1536 : bk => grid%bk
108 :
109 1536 : pdt = get_step_size() ! Physics time step
110 1536 : dtime = pdt
111 :
112 6144 : do j = jfirstxy, jlastxy
113 116736 : do i=ifirstxy, ilastxy
114 115200 : dyn_in%pe(i,1,j) = ptop
115 : enddo
116 : enddo
117 :
118 1536 : if ( nlres) then ! restart or branch run
119 : !
120 : ! read_restart_dynamics delivers phis, ps, u3s, v3s, delp, pt
121 : ! in XY decomposition
122 :
123 : !
124 : ! Do not recalculate delta pressure (delp) if this is a restart run.
125 : ! Re. SJ Lin: The variable "delp" (pressure thikness for a Lagrangian
126 : ! layer) must be in the restart file. This is because delp will be
127 : ! modified "after" the physics update (to account for changes in water
128 : ! vapor), and it can not be reproduced by surface pressure and the
129 : ! ETA coordinate's a's and b's.
130 :
131 : !$omp parallel do private(i,j,k)
132 3072 : do j = jfirstxy, jlastxy
133 76800 : do k=1, km
134 1845504 : do i=ifirstxy, ilastxy
135 1843200 : dyn_in%pe(i,k+1,j) = dyn_in%pe(i,k,j) + dyn_in%delp(i,j,k)
136 : enddo
137 : enddo
138 : enddo
139 : else
140 :
141 : ! Initial run --> generate pe and delp from the surface pressure
142 :
143 : !$omp parallel do private(i,j,k)
144 3072 : do j = jfirstxy, jlastxy
145 79104 : do k=1,km+1
146 1903104 : do i=ifirstxy, ilastxy
147 1900800 : dyn_in%pe(i,k,j) = ak(k) + bk(k) * dyn_in%ps(i,j)
148 : enddo
149 : enddo
150 : enddo
151 :
152 : !$omp parallel do private(i,j,k)
153 25344 : do k = 1, km
154 99072 : do j = jfirstxy, jlastxy
155 1867776 : do i= ifirstxy, ilastxy
156 1843200 : dyn_in%delp(i,j,k) = dyn_in%pe(i,k+1,j) - dyn_in%pe(i,k,j)
157 : enddo
158 : enddo
159 : enddo
160 : endif
161 :
162 : !----------------------------------------------------------
163 : ! Check total dry air mass; set to 982.22 mb if initial run
164 : ! Print out diagnostic message if restart run
165 : !----------------------------------------------------------
166 :
167 1536 : if (scale_dry_air_mass /= 0.0_r8) then
168 : call dryairm( grid, .true., dyn_in%ps, dyn_in%tracer, &
169 1536 : dyn_in%delp, dyn_in%pe, nlres )
170 : endif
171 :
172 1536 : if (grid%iam < grid%npes_xy) then
173 :
174 7680 : allocate( cappa3v(ifirstxy:ilastxy,jfirstxy:jlastxy,km) )
175 7680 : allocate( cap3vi(ifirstxy:ilastxy,jfirstxy:jlastxy,km+1) )
176 1536 : if (grid%high_alt) then
177 0 : call cam_thermo_calc_kappav( dyn_in%tracer, cappa3v )
178 :
179 : !$omp parallel do private(i,j,k)
180 0 : do k=2,km
181 0 : do j=jfirstxy,jlastxy
182 0 : do i=ifirstxy,ilastxy
183 0 : cap3vi(i,j,k) = 0.5_r8*(cappa3v(i,j,k-1)+cappa3v(i,j,k))
184 : enddo
185 : enddo
186 : enddo
187 0 : cap3vi(:,:,1) = 1.5_r8 * cappa3v(:,:,1) - 0.5_r8 * cappa3v(:,:,2)
188 0 : cap3vi(:,:,km+1) = 1.5_r8 * cappa3v(:,:,km) - 0.5_r8 * cappa3v(:,:,km-1)
189 : else
190 3737088 : cappa3v = rair/cpair
191 3853824 : cap3vi = rair/cpair
192 : endif
193 :
194 : ! Initialize pk, edge pressure to the cappa power. Do this with constituent dependent cappa
195 :
196 : !$omp parallel do private(i,j,k)
197 52224 : do k = 1, km+1
198 204288 : do j = jfirstxy, jlastxy
199 3852288 : do i = ifirstxy, ilastxy
200 3801600 : dyn_in%pk(i,j,k) = dyn_in%pe(i,k,j)**cap3vi(i,j,k)
201 : enddo
202 : enddo
203 : enddo
204 :
205 : ! Generate pkz, the conversion factor betw pt and t3
206 :
207 : call pkez(1, im, km, jfirstxy, jlastxy, &
208 : 1, km, ifirstxy, ilastxy, dyn_in%pe, &
209 1536 : dyn_in%pk, cappa3v, ks, dyn_out%peln, dyn_out%pkz, .false., grid%high_alt )
210 :
211 1536 : deallocate( cappa3v, cap3vi )
212 :
213 : endif
214 :
215 1536 : if (initial_run) then
216 :
217 : ! Compute pt for initial run: scaled virtual potential temperature
218 : ! defined as (virtual temp deg K)/pkz. pt will be written to restart (SJL)
219 :
220 : !$omp parallel do private(i,j,k)
221 25344 : do k = 1, km
222 99072 : do j = jfirstxy, jlastxy
223 1867776 : do i = ifirstxy, ilastxy
224 0 : dyn_in%pt(i,j,k) = dyn_in%t3(i,j,k)* &
225 0 : (1._r8 + zvir*dyn_in%tracer(i,j,k,1)) &
226 1843200 : /dyn_in%pkz(i,j,k)
227 : enddo
228 : enddo
229 : enddo
230 :
231 : !----------------------------------------------------------------
232 : ! Convert mixing ratios initialized as dry to moist for dynamics
233 : !----------------------------------------------------------------
234 :
235 : ! on initial time step, dry mixing ratio advected constituents have been
236 : ! initialized to dry mixing ratios. dynpkg expects moist m.r. so convert here.
237 :
238 : ! first calculate delpdry. The set_pdel_state subroutine
239 : ! is called after the dynamics in d_p_coupling to set more variables.
240 : ! This is not in tracers.F90 because it is only used by LR dynamics.
241 3840 : allocate (delpdryxy(ifirstxy:ilastxy,jfirstxy:jlastxy,1:km))
242 25344 : do k = 1, km
243 99072 : do j = jfirstxy, jlastxy
244 1867776 : do i = ifirstxy, ilastxy
245 5308416 : delpdryxy(i,j,k) = dyn_in%delp(i,j,k)* &
246 7151616 : (1._r8 - dyn_in%tracer(i,j,k,1))
247 : enddo
248 : enddo
249 : enddo
250 185088 : do m = 1,pcnst
251 185088 : if (initial_mr(m) == 'dry') then
252 304128 : do k=1, km
253 1188864 : do j = jfirstxy, jlastxy
254 22413312 : do i = ifirstxy, ilastxy
255 0 : dyn_in%tracer(i,j,k,m) = &
256 : dyn_in%tracer(i,j,k,m)* &
257 22118400 : delpdryxy(i,j,k)/dyn_in%delp(i,j,k)
258 : end do
259 : end do
260 : end do
261 : end if
262 : end do
263 768 : deallocate (delpdryxy)
264 :
265 : end if
266 :
267 : ! get aerosol properties
268 1536 : aero_props_obj => aerosol_properties_object()
269 :
270 1536 : if (associated(aero_props_obj)) then
271 : ! determine if there are transported aerosol contistuents
272 1536 : aerosols_transported = aero_props_obj%number_transported()>0
273 : end if
274 :
275 3072 : end subroutine stepon_init
276 :
277 : !=========================================================================================
278 :
279 32256 : subroutine stepon_run1( dtime_out, phys_state, phys_tend, pbuf2d, &
280 : dyn_in, dyn_out )
281 :
282 : ! Phase 1 run of FV dynamics. Run the dynamics, and couple to physics.
283 :
284 1536 : use dp_coupling, only: d_p_coupling
285 : use dyn_comp, only: dyn_run
286 :
287 : use physics_buffer, only: physics_buffer_desc
288 : use advect_tend, only: compute_adv_tends_xyz
289 :
290 : ! arguments
291 : real(r8), intent(out) :: dtime_out ! Time-step
292 : type(physics_state), intent(inout) :: phys_state(begchunk:endchunk)
293 : type(physics_tend), intent(inout) :: phys_tend(begchunk:endchunk)
294 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
295 : type(dyn_import_t) :: dyn_in ! Dynamics import container
296 : type(dyn_export_t) :: dyn_out ! Dynamics export container
297 :
298 : type(T_FVDYCORE_STATE), pointer :: dyn_state
299 :
300 : integer :: rc
301 :
302 : integer :: c
303 : class(aerosol_state), pointer :: aero_state_obj
304 16128 : nullify(aero_state_obj)
305 :
306 16128 : dtime_out = dtime
307 16128 : dyn_state => get_dyn_state()
308 :
309 : ! Dump state variables to IC file
310 16128 : call t_barrierf('sync_diag_dynvar_ic', mpicom)
311 16128 : call t_startf ('diag_dynvar_ic')
312 : call diag_dynvar_ic (dyn_state%grid, dyn_out%phis, dyn_out%ps, &
313 16128 : dyn_out%t3, dyn_out%u3s, dyn_out%v3s, dyn_out%tracer )
314 16128 : call t_stopf ('diag_dynvar_ic')
315 :
316 16128 : call t_startf ('comp_adv_tends1')
317 16128 : call compute_adv_tends_xyz(dyn_state%grid, dyn_in%tracer )
318 16128 : call t_stopf ('comp_adv_tends1')
319 : !
320 : !--------------------------------------------------------------------------
321 : ! Perform finite-volume dynamics -- this dynamical core contains some
322 : ! yet to be published algorithms. Its use in the CAM is
323 : ! for software development purposes only.
324 : ! Please contact S.-J. Lin (Shian-Jiann.Lin@noaa.gov)
325 : ! if you plan to use this mudule for scientific purposes. Contact S.-J. Lin
326 : ! or Will Sawyer (sawyer@gmao.gsfc.nasa.gov) if you plan to modify the
327 : ! software.
328 : !--------------------------------------------------------------------------
329 :
330 : !----------------------------------------------------------
331 : ! For 2-D decomposition, phisxy is input to dynpkg, and the other
332 : ! xy variables are output. Some are computed through direct
333 : ! transposes, and others are derived.
334 : !----------------------------------------------------------
335 16128 : call t_barrierf('sync_dyn_run', mpicom)
336 16128 : call t_startf ('dyn_run')
337 : call dyn_run(ptop, pdt, te0, &
338 16128 : dyn_state, dyn_in, dyn_out, rc )
339 16128 : if ( rc /= 0 ) then
340 0 : write(iulog,*) "STEPON_RUN: dyn_run returned bad error code", rc
341 0 : write(iulog,*) "Quitting."
342 0 : call endrun
343 : endif
344 16128 : call t_stopf ('dyn_run')
345 :
346 16128 : call t_startf ('comp_adv_tends2')
347 16128 : call compute_adv_tends_xyz(dyn_state%grid, dyn_out%tracer )
348 16128 : call t_stopf ('comp_adv_tends2')
349 :
350 : !----------------------------------------------------------
351 : ! Move data into phys_state structure.
352 : !----------------------------------------------------------
353 16128 : call t_barrierf('sync_d_p_coupling', mpicom)
354 16128 : call t_startf('d_p_coupling')
355 16128 : call d_p_coupling(dyn_state%grid, phys_state, phys_tend, pbuf2d, dyn_out)
356 16128 : call t_stopf('d_p_coupling')
357 :
358 : !----------------------------------------------------------
359 : ! update aerosol state object from CAM physics state constituents
360 : !----------------------------------------------------------
361 16128 : if (aerosols_transported) then
362 :
363 0 : do c = begchunk,endchunk
364 0 : aero_state_obj => aerosol_state_object(c)
365 : ! pass number mass or number mixing ratios of aerosol constituents
366 : ! to aerosol state object
367 0 : call aero_state_obj%set_transported(phys_state(c)%q)
368 : end do
369 :
370 : end if
371 :
372 : !EOC
373 16128 : end subroutine stepon_run1
374 :
375 : !-----------------------------------------------------------------------
376 :
377 : !-----------------------------------------------------------------------
378 : !BOP
379 : ! !ROUTINE: stepon_run2 -- second phase run method
380 : !
381 : ! !INTERFACE:
382 29184 : subroutine stepon_run2( phys_state, phys_tend, dyn_in, dyn_out )
383 : ! !USES:
384 16128 : use dp_coupling, only: p_d_coupling
385 : !
386 : ! !INPUT/OUTPUT PARAMETERS:
387 : !
388 : type(physics_state), intent(inout) :: phys_state(begchunk:endchunk)
389 : type(physics_tend), intent(inout) :: phys_tend(begchunk:endchunk)
390 : type (dyn_import_t), intent(inout) :: dyn_in ! Dynamics import container
391 : type (dyn_export_t), intent(inout) :: dyn_out ! Dynamics export container
392 :
393 : type (T_FVDYCORE_GRID), pointer :: grid
394 :
395 : integer :: c
396 : class(aerosol_state), pointer :: aero_state_obj
397 :
398 : !
399 : ! !DESCRIPTION:
400 : !
401 : ! Second phase run method. Couple from physics to dynamics.
402 : !
403 : !EOP
404 : !-----------------------------------------------------------------------
405 : !BOC
406 :
407 : !-----------------------------------------------------------------------
408 :
409 : !----------------------------------------------------------
410 : ! update physics state with aerosol constituents
411 : !----------------------------------------------------------
412 14592 : nullify(aero_state_obj)
413 :
414 14592 : if (aerosols_transported) then
415 0 : do c = begchunk,endchunk
416 0 : aero_state_obj => aerosol_state_object(c)
417 : ! get mass or number mixing ratios of aerosol constituents
418 0 : call aero_state_obj%get_transported(phys_state(c)%q)
419 : end do
420 : end if
421 :
422 : !----------------------------------------------------------
423 : ! Update dynamics variables using phys_state & phys_tend.
424 : ! 2-D decomposition: Compute ptxy and q3xy; for ideal
425 : ! physics, scale ptxy by (old) pkzxy; then transpose to yz variables
426 : ! 1-D decomposition: Compute dudt, dvdt, pt and q3; for ideal physics,
427 : ! scale pt by old pkz.
428 : ! Call uv3s_update to update u3s and v3s from dudt and dvdt.
429 : ! Call p_d_adjust to update pt, q3, pe, delp, ps, piln, pkz and pk.
430 : ! For adiabatic case, transpose to yz variables.
431 : !----------------------------------------------------------
432 14592 : grid => get_dyn_state_grid()
433 :
434 14592 : call t_barrierf('sync_p_d_coupling', mpicom)
435 14592 : call t_startf ('p_d_coupling')
436 : call p_d_coupling(grid, phys_state, phys_tend, &
437 14592 : dyn_in, dtime, zvir, cappa, ptop)
438 14592 : call t_stopf ('p_d_coupling')
439 :
440 : !EOC
441 14592 : end subroutine stepon_run2
442 :
443 : !-----------------------------------------------------------------------
444 :
445 29184 : subroutine stepon_run3(dtime, cam_out, phys_state, &
446 : dyn_in, dyn_out )
447 : ! !USES:
448 14592 : use time_manager, only: get_curr_date
449 : use fv_prints, only: fv_out
450 : use camsrfexch, only: cam_out_t
451 : !
452 : ! !INPUT PARAMETERS:
453 : !
454 : type(physics_state), intent(in):: phys_state(begchunk:endchunk)
455 : real(r8), intent(in) :: dtime ! Time-step
456 : type (dyn_import_t), intent(inout) :: dyn_in ! Dynamics import container
457 : type (dyn_export_t), intent(inout) :: dyn_out ! Dynamics export container
458 : !
459 : ! !INPUT/OUTPUT PARAMETERS:
460 : !
461 : type(cam_out_t), intent(inout) :: cam_out(begchunk:endchunk)
462 : !
463 : ! !DESCRIPTION:
464 : !
465 : ! Final run phase of dynamics. Some printout and time index updates.
466 : !
467 : ! !HISTORY:
468 : ! 2005.09.16 Kluzek Creation
469 : ! 2006.04.13 Sawyer Removed shift_time_indices (not needed in FV)
470 : !
471 : !EOP
472 : !-----------------------------------------------------------------------
473 : !BOC
474 : !
475 : ! !LOCAL VARIABLES:
476 : !
477 :
478 : type(t_fvdycore_state), pointer :: state
479 : type(t_fvdycore_grid), pointer :: grid
480 : integer :: ncdate ! current date in integer format [yyyymmdd]
481 : integer :: ncsec ! time of day relative to current date [seconds]
482 : integer :: yr, mon, day ! year, month, day components of a date
483 : integer :: ncsecp
484 : integer :: freq_diag
485 :
486 : !----------------------------------------------------------
487 : ! Monitor max/min/mean of selected fields
488 : !
489 : ! SEE BELOW **** SEE BELOW **** SEE BELOW
490 :
491 : ! Beware that fv_out uses both dynamics and physics instantiations.
492 : ! However, I think that they are used independently, so that the
493 : ! answers are correct. Still, this violates the notion that the
494 : ! physics state is no longer active after p_d_coupling.
495 : !----------------------------------------------------------
496 14592 : call get_curr_date(yr, mon, day, ncsec)
497 14592 : ncdate = yr*10000 + mon*100 + day
498 14592 : ncsecp = ncsec + pdt ! step complete, but nstep not incremented yet
499 :
500 14592 : state => get_dyn_state()
501 14592 : freq_diag = state%check_dt
502 :
503 14592 : if (fv_monitor .and. mod(ncsecp, freq_diag) == 0) then
504 768 : grid => state%grid
505 :
506 768 : call t_barrierf('sync_fv_out', mpicom)
507 768 : call t_startf('fv_out')
508 : call fv_out(grid, dyn_out%pk, dyn_out%pt, &
509 : ptop, dyn_out%ps, dyn_out%tracer, &
510 : dyn_out%delp, dyn_out%pe, cam_out, &
511 768 : phys_state, ncdate, ncsecp, moist_physics)
512 768 : call t_stopf('fv_out')
513 : endif
514 :
515 : !EOC
516 14592 : end subroutine stepon_run3
517 :
518 : !-----------------------------------------------------------------------
519 :
520 : !-----------------------------------------------------------------------
521 : !BOP
522 : ! !ROUTINE: stepon_final --- Dynamics finalization
523 : !
524 : ! !INTERFACE:
525 1536 : subroutine stepon_final(dyn_in, dyn_out)
526 :
527 : ! !PARAMETERS:
528 : type (dyn_import_t), intent(out) :: dyn_in ! Dynamics import container
529 : type (dyn_export_t), intent(out) :: dyn_out ! Dynamics export container
530 : !
531 : ! !DESCRIPTION:
532 : !
533 : ! Deallocate data needed for dynamics. Finalize any dynamics specific
534 : ! files or subroutines.
535 : !
536 : !EOP
537 : !-----------------------------------------------------------------------
538 : !BOC
539 :
540 : !!! Not yet ready for the call to dyn_final
541 : !!! call dyn_final( RESTART_FILE, dyn_state, dyn_in, dyn_out )
542 : !EOC
543 14592 : end subroutine stepon_final
544 :
545 : !-----------------------------------------------------------------------
546 :
547 : end module stepon
|
