forked from NCAR/ccpp-physics
-
Notifications
You must be signed in to change notification settings - Fork 37
/
Copy pathcires_ugwp.F90
448 lines (377 loc) · 19.6 KB
/
cires_ugwp.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
!> \file cires_ugwp.F90
!! This file contains the Unified Gravity Wave Physics (UGWP) scheme by Valery Yudin (University of Colorado, CIRES)
!! See Valery Yudin's presentation at 2017 NGGPS PI meeting:
!! Gravity waves (GWs): Mesoscale GWs transport momentum, energy (heat) , and create eddy mixing in the whole atmosphere domain; Breaking and dissipating GWs deposit: (a) momentum; (b) heat (energy); and create (c) turbulent mixing of momentum, heat, and tracers
!! To properly incorporate GW effects (a-c) unresolved by DYCOREs we need GW physics
!! "Unified": a) all GW effects due to both dissipation/breaking; b) identical GW solvers for all GW sources; c) ability to replace solvers.
!! Unified Formalism:
!! 1. GW Sources: Stochastic and physics based mechanisms for GW-excitations in the lower atmosphere, calibrated by the high-res analyses/forecasts, and observations (3 types of GW sources: orography, convection, fronts/jets).
!! 2. GW Propagation: Unified solver for "propagation, dissipation and breaking" excited from all type of GW sources.
!! 3. GW Effects: Unified representation of GW impacts on the "resolved" flow for all sources (energy-balanced schemes for momentum, heat and mixing).
!! https://www.weather.gov/media/sti/nggps/Presentations%202017/02%20NGGPS_VYUDIN_2017_.pdf
module cires_ugwp
use machine, only: kind_phys
use cires_ugwpv0_module, only: knob_ugwp_version, cires_ugwpv0_mod_init, cires_ugwpv0_mod_finalize
use ugwp_driver_v0
use gwdps, only: gwdps_run
use cires_ugwp_triggers
implicit none
private
public cires_ugwp_init, cires_ugwp_run, cires_ugwp_finalize
logical :: is_initialized = .False.
contains
! ------------------------------------------------------------------------
! CCPP entry points for CIRES Unified Gravity Wave Physics (UGWP) scheme v0
! ------------------------------------------------------------------------
!>\defgroup cires_ugwp_run_mod CIRES Unified Gravity Wave Physics v0 Module
!> @{
!>@ The subroutine initializes the CIRES UGWP V0.
!> \section arg_table_cires_ugwp_init Argument Table
!! \htmlinclude cires_ugwp_init.html
!!
subroutine cires_ugwp_init (me, master, nlunit, input_nml_file, logunit, &
fn_nml2, lonr, latr, levs, ak, bk, dtp, cdmbgwd, cgwf, &
pa_rf_in, tau_rf_in, con_p0, gwd_opt,do_ugwp, errmsg, errflg)
!---- initialization of cires_ugwp
implicit none
integer, intent (in) :: me
integer, intent (in) :: master
integer, intent (in) :: nlunit
character(len=*), intent (in) :: input_nml_file(:)
integer, intent (in) :: logunit
integer, intent (in) :: lonr
integer, intent (in) :: levs
integer, intent (in) :: latr
real(kind=kind_phys), intent (in) :: ak(:), bk(:)
real(kind=kind_phys), intent (in) :: dtp
real(kind=kind_phys), intent (in) :: cdmbgwd(:), cgwf(:) ! "scaling" controls for "old" GFS-GW schemes
real(kind=kind_phys), intent (in) :: pa_rf_in, tau_rf_in
real(kind=kind_phys), intent (in) :: con_p0
integer, intent(in) :: gwd_opt
logical, intent (in) :: do_ugwp
character(len=*), intent (in) :: fn_nml2
!character(len=*), parameter :: fn_nml='input.nml'
integer :: ios
logical :: exists
real :: dxsg
integer :: k
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
if (is_initialized) return
! Consistency checks
if (gwd_opt/=1) then
write(errmsg,'(*(a))') "Logic error: namelist choice of gravity wave &
& drag is different from cires_ugwp scheme"
errflg = 1
return
end if
if (do_ugwp .or. cdmbgwd(3) > 0.0) then
call cires_ugwpv0_mod_init (me, master, nlunit, input_nml_file, logunit, &
fn_nml2, lonr, latr, levs, ak, bk, con_p0, dtp, &
cdmbgwd(1:2), cgwf, pa_rf_in, tau_rf_in)
else
write(errmsg,'(*(a))') "Logic error: cires_ugwp_init called but do_ugwp is false and cdmbgwd(3) <= 0"
errflg = 1
return
end if
if (.not.knob_ugwp_version==0) then
write(errmsg,'(*(a))') 'Logic error: CCPP only supports version zero of UGWP'
errflg = 1
return
end if
is_initialized = .true.
end subroutine cires_ugwp_init
! -----------------------------------------------------------------------
! finalize of cires_ugwp (_finalize)
! -----------------------------------------------------------------------
!>@brief The subroutine finalizes the CIRES UGWP
#if 0
!> \section arg_table_cires_ugwp_finalize Argument Table
!! \htmlinclude cires_ugwp_finalize.html
!!
#endif
subroutine cires_ugwp_finalize(errmsg, errflg)
implicit none
!
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
if (.not.is_initialized) return
call cires_ugwpv0_mod_finalize()
is_initialized = .false.
end subroutine cires_ugwp_finalize
! -----------------------------------------------------------------------
! originally from ugwp_driver_v0.f
! driver of cires_ugwp (_driver)
! -----------------------------------------------------------------------
! driver is called after pbl & before chem-parameterizations
! -----------------------------------------------------------------------
! order = dry-adj=>conv=mp-aero=>radiation -sfc/land- chem -> vertdiff-> [rf-gws]=> ion-re
! -----------------------------------------------------------------------
!>@brief These subroutines and modules execute the CIRES UGWP Version 0.
!> \section gen_cires_ugwp CIRES UGWP V0 Scheme General Algorithm
!! The physics of Non-Orographic Gravity Waves (NGWs) in the UGWP framework
!!(Yudin et al. 2018 \cite yudin_et_al_2018) is represented by four GW-solvers, introduced in
!!Lindzen (1981) \cite lindzen_1981, Hines (1997) \cite hines_1997, Alexander
!!and Dunkerton (1999) \cite alexander_and_dunkerton_1999, and Scinocca (2003) \cite scinocca_2003.
!!A major modification of these GW solvers was introduced with the addition of the
!!background dissipation of temperature and winds to the saturation criteria for wave breaking.
!!This feature is important in the mesosphere and thermosphere for WAM applications and it
!!considers appropriate scale-dependent dissipation of waves near the model top lid providing
!!the momentum and energy conservation in the vertical column physics (Shaw and
!!Shepherd (2009) \cite shaw_and_shepherd_2009). In the UGWP-v0 scheme, a modification of the
!!Scinocca (2003) \cite scinocca_2003 algorithm for NGWs with non-hydrostatic and rotational
!!effects for GW propagations and background dissipation is contained in the subroutine
!!fv3_ugwp_solv2_v0. Future development plans for the UGWP scheme include additional GW-solvers
!!to be implemented along with physics-based triggering of waves and stochastic approaches
!!for selection of GW modes characterized by horizontal phase velocities, azimuthal
!!directions and magnitude of the vertical momentum flux (VMF).
!!
!! In UGWP-v0, the specification for the VMF function is adopted from the
!! GEOS-5 global atmosphere model of GMAO NASA/GSFC, as described in
!! Molod et al. (2015) \cite molod_et_al_2015 and employed in the MERRRA-2
!! reanalysis (Gelaro et al., 2017 \cite gelaro_et_al_2017). The Fortran
!! subroutine slat_geos5_tamp_v0() describes the latitudinal shape of
!! VMF-function as displayed in Figure 3 of Molod et al. (2015)
!! \cite molod_et_al_2015. It shows that the enhanced values of
!! VMF in the equatorial region gives opportunity to simulate the
!! QBO-like oscillations in the equatorial zonal winds and lead to more
!! realistic simulations of the equatorial dynamics in GEOS-5 operational
!! and MERRA-2 reanalysis products. For the first vertically extended
!! version of FV3GFS in the stratosphere and mesosphere, this simplified
!! function of VMF allows us to tune the model climate and to evaluate
!! multi-year simulations of FV3GFS with the MERRA-2 and ERA-5 reanalysis
!! products, along with temperature, ozone, and water vapor observations
!! of current satellite missions. After delivery of the UGWP-code, the
!! EMC group developed and tested approach to modulate the zonal mean
!! NGW forcing by 3D-distributions of the total precipitation as a proxy
!! for the excitation of NGWs by convection and the vertically-integrated
!! (surface - tropopause) Turbulent Kinetic Energy (TKE). The verification
!! scores with updated NGW forcing, as reported elsewhere by EMC researchers,
!! display noticeable improvements in the forecast scores produced by
!! FV3GFS configuration extended into the mesosphere.
!!
!> \section arg_table_cires_ugwp_run Argument Table
!! \htmlinclude cires_ugwp_run.html
!!
! \section det_cires_ugwp CIRES UGWP V0 Scheme Detailed Algorithm
subroutine cires_ugwp_run(do_ugwp, me, master, im, levs, ntrac, dtp, kdt, lonr, &
oro, oro_uf, hprime, nmtvr, oc, theta, sigma, gamma, elvmax, clx, oa4, &
do_tofd, ldiag_ugwp, cdmbgwd, xlat, xlat_d, sinlat, coslat, &
area, ugrs, vgrs, tgrs, qgrs, prsi, prsl, prslk, phii, phil, &
del, kpbl, dusfcg, dvsfcg, gw_dudt, gw_dvdt, gw_dtdt, gw_kdis, &
tau_tofd, tau_mtb, tau_ogw, tau_ngw, zmtb, zlwb, zogw, &
dusfc_ms, dvsfc_ms, dusfc_bl, dvsfc_bl, &
dudt_ogw, dtauy2d_ms, dtaux2d_bl, dtauy2d_bl, &
dudt_mtb, dudt_tms, du3dt_mtb, du3dt_ogw, du3dt_tms, &
dudt, dvdt, dtdt, rdxzb, con_g, con_pi, con_cp, con_rd, con_rv, con_fvirt, &
con_omega, rain, ntke, q_tke, dqdt_tke, lprnt, ipr, &
dtend, dtidx, index_of_x_wind, index_of_y_wind, index_of_temperature, &
index_of_process_orographic_gwd, index_of_process_nonorographic_gwd, &
ldiag3d, lssav, flag_for_gwd_generic_tend, errmsg, errflg)
implicit none
! interface variables
integer, intent(in) :: me, master, im, levs, ntrac, kdt, lonr, nmtvr
integer, intent(in), dimension(:) :: kpbl
real(kind=kind_phys), intent(in), dimension(:) :: oro, oro_uf, hprime, oc, theta, sigma, gamma
logical, intent(in) :: flag_for_gwd_generic_tend
! elvmax is intent(in) for CIRES UGWP, but intent(inout) for GFS GWDPS
real(kind=kind_phys), intent(inout), dimension(:) :: elvmax
real(kind=kind_phys), intent(in), dimension(:, :) :: clx, oa4
real(kind=kind_phys), intent(in), dimension(:) :: xlat, xlat_d, sinlat, coslat, area
real(kind=kind_phys), intent(in), dimension(:, :) :: del, ugrs, vgrs, tgrs, prsl, prslk, phil
real(kind=kind_phys), intent(in), dimension(:, :) :: prsi, phii
real(kind=kind_phys), intent(in), dimension(:,:,:):: qgrs
real(kind=kind_phys), intent(in) :: dtp, cdmbgwd(:)
logical, intent(in) :: do_ugwp, do_tofd, ldiag_ugwp
real(kind=kind_phys), intent(out), dimension(:) :: dusfcg, dvsfcg
real(kind=kind_phys), intent(out), dimension(:) :: zmtb, zlwb, zogw, rdxzb
real(kind=kind_phys), intent(out), dimension(:) :: tau_mtb, tau_ogw, tau_tofd, tau_ngw
real(kind=kind_phys), intent(out), dimension(:, :):: gw_dudt, gw_dvdt, gw_dtdt, gw_kdis
real(kind=kind_phys), intent(out), dimension(:, :):: dudt_mtb, dudt_ogw, dudt_tms
real(kind=kind_phys), intent(out), dimension(:) :: dusfc_ms, dvsfc_ms, dusfc_bl, dvsfc_bl
real(kind=kind_phys), intent(out), dimension(:, :) :: dtauy2d_ms
real(kind=kind_phys), intent(out), dimension(:, :) :: dtaux2d_bl, dtauy2d_bl
! dtend is only allocated if ldiag=.true.
real(kind=kind_phys), optional, intent(inout) :: dtend(:,:,:)
integer, intent(in) :: dtidx(:,:), &
index_of_x_wind, index_of_y_wind, index_of_temperature, &
index_of_process_orographic_gwd, index_of_process_nonorographic_gwd
logical, intent(in) :: ldiag3d, lssav
! These arrays only allocated if ldiag_ugwp = .true.
real(kind=kind_phys), intent(inout), dimension(:,:) :: du3dt_mtb, du3dt_ogw, du3dt_tms
real(kind=kind_phys), intent(inout), dimension(:, :):: dudt, dvdt, dtdt
real(kind=kind_phys), intent(in) :: con_g, con_pi, con_cp, con_rd, con_rv, con_fvirt, con_omega
real(kind=kind_phys), intent(in), dimension(:) :: rain
integer, intent(in) :: ntke
real(kind=kind_phys), intent(in), dimension(:,:) :: q_tke, dqdt_tke
logical, intent(in) :: lprnt
integer, intent(in) :: ipr
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
! local variables
integer :: i, k, idtend
real(kind=kind_phys), dimension(im) :: sgh30
real(kind=kind_phys), dimension(im, levs) :: Pdvdt, Pdudt
real(kind=kind_phys), dimension(im, levs) :: Pdtdt, Pkdis
real(kind=kind_phys), dimension(im, levs) :: ed_dudt, ed_dvdt, ed_dtdt
! from ugwp_driver_v0.f -> cires_ugwp_initialize.F90 -> module ugwp_wmsdis_init
real(kind=kind_phys), parameter :: tamp_mpa=30.e-3
! switches that activate impact of OGWs and NGWs (WL* how to deal with them? *WL)
real(kind=kind_phys), parameter :: pogw=1., pngw=1., pked=1.
real(kind=kind_phys), dimension(:,:), allocatable :: tke
real(kind=kind_phys), dimension(:), allocatable :: turb_fac, tem
real(kind=kind_phys) :: rfac, tx1
! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
! 1) ORO stationary GWs
! ------------------
! wrap everything in a do_ugwp 'if test' in order not to break the namelist functionality
if (do_ugwp) then ! calling revised old GFS gravity wave drag
! topo paras
! w/ orographic effects
if(nmtvr == 14)then
! calculate sgh30 for TOFD
sgh30 = abs(oro - oro_uf)
! w/o orographic effects
else
sgh30 = 0.
endif
zlwb(:) = 0.
call GWDPS_V0(im, levs, lonr, do_tofd, Pdvdt, Pdudt, Pdtdt, Pkdis, &
ugrs, vgrs, tgrs, qgrs(:,:,1), kpbl, prsi,del,prsl, prslk, phii, phil, &
dtp, kdt, sgh30, hprime, oc, oa4, clx, theta, sigma, gamma, elvmax, &
dusfcg, dvsfcg, xlat_d, sinlat, coslat, area, cdmbgwd(1:2), &
me, master, rdxzb, con_g, con_omega, zmtb, zogw, tau_mtb, tau_ogw, &
tau_tofd, dudt_mtb, dudt_ogw, dudt_tms)
else ! calling old GFS gravity wave drag as is
do k=1,levs
do i=1,im
Pdvdt(i,k) = 0.0
Pdudt(i,k) = 0.0
Pdtdt(i,k) = 0.0
Pkdis(i,k) = 0.0
enddo
enddo
if (cdmbgwd(1) > 0.0 .or. cdmbgwd(2) > 0.0) then
call gwdps_run(im, levs, Pdvdt, Pdudt, Pdtdt, &
ugrs, vgrs, tgrs, qgrs(:,:,1), &
kpbl, prsi, del, prsl, prslk, phii, phil, dtp, kdt, &
hprime, oc, oa4, clx, theta, sigma, gamma, &
elvmax, dusfcg, dvsfcg, dudt_ogw, dtauy2d_ms, &
dtaux2d_bl, dtauy2d_bl, dusfc_ms, dvsfc_ms, &
dusfc_bl, dvsfc_bl, &
con_g, con_cp, con_rd, con_rv, lonr, &
nmtvr, cdmbgwd, me, lprnt, ipr, rdxzb, ldiag_ugwp, &
errmsg, errflg)
if (errflg/=0) return
endif
tau_mtb = 0.0 ; tau_ogw = 0.0 ; tau_tofd = 0.0
if (ldiag_ugwp) then
du3dt_mtb = 0.0 ; du3dt_ogw = 0.0 ; du3dt_tms= 0.0
endif
endif ! do_ugwp
if(ldiag3d .and. lssav .and. .not. flag_for_gwd_generic_tend) then
idtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd)
if(idtend>=1) then
dtend(:,:,idtend) = dtend(:,:,idtend) + Pdudt*dtp
endif
idtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd)
if(idtend>=1) then
dtend(:,:,idtend) = dtend(:,:,idtend) + Pdvdt*dtp
endif
idtend = dtidx(index_of_temperature,index_of_process_orographic_gwd)
if(idtend>=1) then
dtend(:,:,idtend) = dtend(:,:,idtend) + Pdtdt*dtp
endif
endif
if (cdmbgwd(3) > 0.0) then
! 2) non-stationary GW-scheme with GMAO/MERRA GW-forcing
call slat_geos5_tamp_v0(im, tamp_mpa, xlat_d, tau_ngw)
if (abs(1.0-cdmbgwd(3)) > 1.0e-6) then
if (cdmbgwd(4) > 0.0) then
allocate(turb_fac(im))
do i=1,im
turb_fac(i) = 0.0
enddo
if (ntke > 0) then
allocate(tke(im,levs))
allocate(tem(im))
tke(:,:) = q_tke(:,:) + dqdt_tke(:,:) * dtp
tem(:) = 0.0
do k=1,(levs+levs)/3
do i=1,im
turb_fac(i) = turb_fac(i) + del(i,k) * tke(i,k)
tem(i) = tem(i) + del(i,k)
enddo
enddo
do i=1,im
turb_fac(i) = turb_fac(i) / tem(i)
enddo
deallocate(tke)
deallocate(tem)
endif
rfac = 86400000 / dtp
do i=1,im
tx1 = cdmbgwd(4)*min(10.0, max(turb_fac(i),rain(i)*rfac))
tau_ngw(i) = tau_ngw(i) * max(0.1, min(5.0, tx1))
enddo
deallocate(turb_fac)
endif
do i=1,im
tau_ngw(i) = tau_ngw(i) * cdmbgwd(3)
enddo
endif
call fv3_ugwp_solv2_v0(im, levs, dtp, tgrs, ugrs, vgrs,qgrs(:,:,1), &
prsl, prsi, phil, xlat_d, sinlat, coslat, &
gw_dudt, gw_dvdt, gw_dtdt, gw_kdis, tau_ngw, &
me, master, kdt)
do k=1,levs
do i=1,im
gw_dtdt(i,k) = pngw*gw_dtdt(i,k) + pogw*Pdtdt(i,k)
gw_dudt(i,k) = pngw*gw_dudt(i,k) + pogw*Pdudt(i,k)
gw_dvdt(i,k) = pngw*gw_dvdt(i,k) + pogw*Pdvdt(i,k)
gw_kdis(i,k) = pngw*gw_kdis(i,k) + pogw*Pkdis(i,k)
! accumulation of tendencies for CCPP to replicate EMC-physics updates (!! removed in latest code commit to VLAB)
!dudt(i,k) = dudt(i,k) + gw_dudt(i,k)
!dvdt(i,k) = dvdt(i,k) + gw_dvdt(i,k)
!dtdt(i,k) = dtdt(i,k) + gw_dtdt(i,k)
enddo
enddo
else
do k=1,levs
do i=1,im
gw_dtdt(i,k) = Pdtdt(i,k)
gw_dudt(i,k) = Pdudt(i,k)
gw_dvdt(i,k) = Pdvdt(i,k)
gw_kdis(i,k) = Pkdis(i,k)
enddo
enddo
endif
if (pogw == 0.0) then
tau_mtb = 0. ; tau_ogw = 0. ; tau_tofd = 0.
dudt_mtb = 0. ; dudt_ogw = 0. ; dudt_tms = 0.
endif
if(ldiag3d .and. lssav .and. .not. flag_for_gwd_generic_tend) then
idtend = dtidx(index_of_x_wind,index_of_process_nonorographic_gwd)
if(idtend>=1) then
dtend(:,:,idtend) = dtend(:,:,idtend) + (gw_dudt - Pdudt)*dtp
endif
idtend = dtidx(index_of_y_wind,index_of_process_nonorographic_gwd)
if(idtend>=1) then
dtend(:,:,idtend) = dtend(:,:,idtend) + (gw_dvdt - Pdvdt)*dtp
endif
idtend = dtidx(index_of_temperature,index_of_process_nonorographic_gwd)
if(idtend>=1) then
dtend(:,:,idtend) = dtend(:,:,idtend) + (gw_dtdt - Pdtdt)*dtp
endif
endif
end subroutine cires_ugwp_run
!> @}
end module cires_ugwp