Line data Source code
1 :
2 : module qbo
3 : !--------------------------------------------------------------------
4 : ! This module performes a relaxation towards either a cyclic idealized
5 : ! QBO sequence derived from observations (here:28months) or towards
6 : ! observed equatorial wind data
7 : !
8 : ! Author: Katja Matthes
9 : ! Date: February 2005
10 : !
11 : ! implementation into WACCM standard version, K. Matthes, October 2005
12 : ! 3 possibilities:
13 : ! default: no QBO relaxation
14 : ! options: 1. cyclic QBO time series, currently a 28month sequence available: qbocyclic28months.nc
15 : ! 2. observed QBO time series, currently CCMVal QBO series, extended: qboseries_ext.nc
16 : ! 3. fixed QBOe or QBOw phase, currently : qboeast.nc, qbowest.nc
17 : !
18 : !--------------------------------------------------------------------
19 : ! modified by Anne Smith, December 2009
20 : !
21 : ! Now accepts alternative input in the form of Fourier coefficients.
22 : ! In this case, the QBO wind is calculated from the expression
23 : ! u_qbo(k,n) = ubar_qbo(k) + sum_n[fcos_qbo(k,n) * cos(ffreq_qbo(n)*(cday-cday_ref))
24 : ! + fsin_qbo(k,n) * sin(ffreq_qbo(n)*(cday-cday_ref))]
25 : ! where cday is day of model run
26 : ! cday_ref is a reference date so that historical data line up correctly
27 : ! k is level index
28 : ! n is coefficient index
29 : ! sum_n is the sum over n
30 : !---------------------------------------------------------------------
31 :
32 : use shr_kind_mod, only: r8 => shr_kind_r8
33 : use spmd_utils, only: masterproc
34 : use ppgrid, only: pcols, pver
35 : use time_manager, only: get_curr_date, get_curr_calday
36 : use phys_grid, only: get_rlat_p
37 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
38 : use cam_abortutils, only: endrun
39 : use cam_logfile, only: iulog
40 : use bnddyi_mod, only: bnddyi
41 :
42 : implicit none
43 :
44 : private
45 : save
46 :
47 : !---------------------------------------------------------------------
48 : ! Public methods
49 : !---------------------------------------------------------------------
50 : public :: qbo_readnl ! read namelist
51 : public :: qbo_init ! initialize qbo package
52 : public :: qbo_timestep_init ! interpolate to current time
53 : public :: qbo_relax ! relax zonal mean wind
54 :
55 : !---------------------------------------------------------------------
56 : ! Private module data
57 : !---------------------------------------------------------------------
58 : integer,parameter :: qbo_dypm = 30 ! days in a cyclic qbo "month"
59 :
60 : integer :: nm, np ! Indices for previous, next month
61 : integer :: np1 = 0 ! tempory for next time index of dataset
62 : integer :: ktop ! Model layers within qbo region
63 : integer :: kbot ! Model layers within qbo region
64 : integer :: timesiz ! size of time dimension on dataset
65 : integer :: levsiz ! size of lev dimension on dataset
66 : integer :: qbo_mons ! length of cyclic qbo in months
67 : integer :: delt ! time step in seconds
68 : real(r8) :: qbo_days ! length of cyclic qbo in days
69 : integer, allocatable :: date_qbo(:) ! Date on qbo dataset (YYYYMMDD) (timesiz)
70 : integer, allocatable :: secd_qbo(:) ! seconds of date (0-86399) (timesiz)
71 :
72 : real(r8) :: cdaym ! dataset calendar day previous month
73 : real(r8) :: cdayp ! dataset calendar day next month
74 : real(r8),allocatable :: u_qbo(:,:) ! qbo winds(pver,timesiz)
75 : real(r8),allocatable :: tauz(:)
76 : real(r8) :: u_tstep(pver) ! qbo winds for this time step
77 :
78 : integer :: coefsiz ! size of coefficient dimension on fft dataset
79 : real(r8) :: cday_ref ! reference day for fft input
80 : real(r8),allocatable :: ubar_qbo(:) ! qbo mean winds(pver)
81 : real(r8),allocatable :: fcos_qbo(:,:) ! qbo cosine coefficients(pver,coefsiz)
82 : real(r8),allocatable :: fsin_qbo(:,:) ! qbo sine coefficients(pver,coefsiz)
83 : real(r8),allocatable :: ffreq_qbo(:) ! frequencies for expanding qbo coefficients (coefsiz)
84 :
85 : ! Index into physics buffer for zonal mean zonal wind
86 : integer :: uzm_idx = -1
87 :
88 : !
89 : ! Options for controlling QBO relaxation
90 : !
91 : character(len=256) :: qbo_forcing_file = ""
92 : logical,public :: qbo_use_forcing = .FALSE. ! .TRUE. => this package is active
93 : logical :: qbo_cyclic = .FALSE. ! .TRUE. => assume cyclic qbo data
94 :
95 : logical :: has_monthly_data = .TRUE. ! .TRUE. => data file has monthly winds
96 : ! .FALSE.=> data file has fft coefficients
97 :
98 :
99 : contains
100 :
101 : !================================================================================================
102 :
103 1536 : subroutine qbo_readnl(nlfile)
104 :
105 : use namelist_utils, only: find_group_name
106 : use units, only: getunit, freeunit
107 : use mpishorthand
108 :
109 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
110 :
111 : ! Local variables
112 : integer :: unitn, ierr
113 : character(len=*), parameter :: subname = 'qbo_readnl'
114 :
115 : namelist /qbo_nl/ qbo_use_forcing, qbo_forcing_file, qbo_cyclic
116 :
117 1536 : if (masterproc) then
118 2 : unitn = getunit()
119 2 : open( unitn, file=trim(nlfile), status='old' )
120 2 : call find_group_name(unitn, 'qbo_nl', status=ierr)
121 2 : if (ierr == 0) then
122 2 : read(unitn, qbo_nl, iostat=ierr)
123 2 : if (ierr /= 0) then
124 0 : call endrun(subname // ':: ERROR reading namelist')
125 : end if
126 : end if
127 2 : close(unitn)
128 2 : call freeunit(unitn)
129 :
130 : ! Check to make sure forcing file was set if qbo forcing is on.
131 2 : if (qbo_use_forcing .and. trim(qbo_forcing_file) == "") then
132 : call endrun(subname // ':: qbo forcing is on but no forcing &
133 0 : &file was set.')
134 : end if
135 : end if
136 :
137 : #ifdef SPMD
138 1536 : call mpibcast (qbo_forcing_file, len(qbo_forcing_file ), mpichar, 0, mpicom)
139 1536 : call mpibcast (qbo_use_forcing, 1, mpilog, 0, mpicom)
140 1536 : call mpibcast (qbo_cyclic, 1, mpilog, 0, mpicom)
141 : #endif
142 :
143 1536 : end subroutine qbo_readnl
144 :
145 : !================================================================================================
146 :
147 :
148 1536 : subroutine qbo_init
149 : !---------------------------------------------------------------------
150 : ! initialize qbo module
151 : !---------------------------------------------------------------------
152 : use physics_buffer, only: pbuf_get_index
153 : use time_manager, only : get_step_size
154 : use ref_pres, only : pref_mid
155 : use phys_control, only : phys_getopts
156 : use cam_history, only : addfld, add_default
157 : use wrap_nf
158 : #if (defined SPMD )
159 : use mpishorthand
160 : #endif
161 :
162 : !---------------------------------------------------------------------
163 : ! Local workspace
164 : !---------------------------------------------------------------------
165 : real(r8), parameter :: hPa2Pa = 100._r8
166 :
167 : integer :: yr, mon, day ! components of a date
168 : integer :: ncsec ! current time of day [seconds]
169 : integer :: ncid ! netcdf ID for input dataset
170 : integer :: astat ! allocate status
171 :
172 : integer :: k,kk ! vertical interpolation indexes
173 : integer :: levdimid ! netcdf id for level dimension
174 : integer :: timdimid ! netcdf id for time dimension
175 : integer :: levid ! netcdf id for level variable
176 : integer :: dateid ! netcdf id for date variable
177 : integer :: secdid ! netcdf id for seconds variable
178 : integer :: uqboid ! netcdf id for qbo wind variable
179 :
180 : real(r8) :: fac1, fac2 ! interpolation factors
181 : real(r8) :: calday ! current calendar day
182 : real(r8) :: cday ! current day (in cycle, if cycling)
183 1536 : real(r8), allocatable :: p_inp(:) ! qbo pressure levels(levsiz)
184 1536 : real(r8), allocatable :: u_inp(:,:) ! input qbo winds(levsiz,timesiz) have to be sorted
185 : ! from top to bottom like winds in model
186 :
187 : integer :: fftdimid ! netcdf id for fft coefficient dimension
188 : integer :: ubarqboid ! netcdf id for mean wind variable
189 : integer :: cosqboid ! netcdf id for cosine coefficient variable
190 : integer :: sinqboid ! netcdf id for sine coefficient variable
191 : integer :: freqqboid ! netcdf id for fft frequency variable
192 : integer :: refdatid ! netcdf id for reference date
193 :
194 : integer :: ref_date ! date corresponding to beginning of QBO data
195 : integer :: yr_ref ! year for reference day
196 1536 : real(r8), allocatable :: ubar_inp(:) ! input time-mean winds for fft input (levsiz)
197 1536 : real(r8), allocatable :: cosq_inp(:,:) ! input cosine coefficients for fft input (levsiz,coefsiz)
198 1536 : real(r8), allocatable :: sinq_inp(:,:) ! input sine coefficients for fft input (levsiz,coefsiz)
199 :
200 : logical :: found
201 :
202 : logical :: history_waccm
203 :
204 1536 : if( .not. qbo_use_forcing ) then
205 : return
206 : end if
207 :
208 0 : call phys_getopts(history_waccm_out=history_waccm)
209 :
210 : !---------------------------------------------------------------------
211 : ! SPMD: Master does all the work of reading files. Sends needed info to slaves
212 : !---------------------------------------------------------------------
213 0 : if( masterproc ) then
214 : !---------------------------------------------------------------------
215 : ! Get qbo file
216 : !---------------------------------------------------------------------
217 0 : call wrap_open( qbo_forcing_file, NF90_NOERR, ncid )
218 0 : write(iulog,*) 'qbo_init: successfully opened ',trim(qbo_forcing_file)
219 : !---------------------------------------------------------------------
220 : ! Figure out if the file contains qbo winds by month or fft coefficients
221 : ! by looking for the variable DATE
222 : !---------------------------------------------------------------------
223 0 : call wrap_inq_varid( ncid, 'date' , dateid, abort=.false. )
224 :
225 0 : if (dateid > 0) then
226 0 : has_monthly_data=.true.
227 : else
228 0 : has_monthly_data=.false.
229 : end if
230 0 : write(iulog,*) 'has_monthly_data=', has_monthly_data
231 :
232 : !---------------------------------------------------------------------
233 : ! Get and check dimension info
234 : !---------------------------------------------------------------------
235 0 : if (has_monthly_data) then
236 0 : call wrap_inq_dimid( ncid, 'level' , levdimid )
237 0 : call wrap_inq_dimid( ncid, 'time', timdimid )
238 :
239 0 : call wrap_inq_dimlen( ncid, levdimid, levsiz )
240 0 : call wrap_inq_dimlen( ncid, timdimid, timesiz )
241 :
242 0 : call wrap_inq_varid( ncid, 'date' , dateid )
243 0 : call wrap_inq_varid( ncid, 'secs' , secdid )
244 0 : call wrap_inq_varid( ncid, 'qbo' , uqboid )
245 0 : call wrap_inq_varid( ncid, 'level', levid )
246 : else
247 0 : call wrap_inq_dimid( ncid, 'level' , levdimid )
248 0 : call wrap_inq_dimid( ncid, 'ncoef', fftdimid )
249 :
250 0 : call wrap_inq_dimlen( ncid, levdimid, levsiz )
251 0 : call wrap_inq_dimlen( ncid, fftdimid, coefsiz )
252 :
253 0 : call wrap_inq_varid( ncid, 'ref_date', refdatid )
254 :
255 0 : call wrap_inq_varid( ncid, 'ubar' , ubarqboid )
256 0 : call wrap_inq_varid( ncid, 'cosqbo' , cosqboid )
257 0 : call wrap_inq_varid( ncid, 'sinqbo' , sinqboid )
258 0 : call wrap_inq_varid( ncid, 'freqqbo', freqqboid )
259 0 : call wrap_inq_varid( ncid, 'level' , levid )
260 : end if
261 : end if
262 :
263 :
264 : !---------------------------------------------------------------------
265 : ! Broadcast the logical flag has_monthly_data
266 : !---------------------------------------------------------------------
267 : #if (defined SPMD )
268 0 : call mpibcast( has_monthly_data, 1, mpilog, 0, mpicom )
269 : #endif
270 :
271 0 : if (.not.has_monthly_data) then
272 : !---------------------------------------------------------------------
273 : ! Broadcast coefsiz
274 : !---------------------------------------------------------------------
275 : #if (defined SPMD )
276 0 : call mpibcast( coefsiz, 1, mpiint, 0, mpicom )
277 : #endif
278 : end if
279 :
280 0 : if (has_monthly_data) then
281 : #if (defined SPMD )
282 : !---------------------------------------------------------------------
283 : ! Broadcast the time size to all tasks for case with monthly input
284 : !---------------------------------------------------------------------
285 0 : call mpibcast( timesiz, 1, mpiint, 0, mpicom )
286 : #endif
287 :
288 0 : delt = get_step_size()
289 0 : if( qbo_cyclic ) then
290 0 : qbo_mons = timesiz
291 0 : qbo_days = qbo_mons*qbo_dypm
292 : end if
293 :
294 0 : allocate( date_qbo(timesiz), secd_qbo(timesiz), stat=astat )
295 0 : if( astat /= 0 ) then
296 0 : write(iulog,*) 'qbo_init: failed to allocate date_qbo ... secd_qbo; error = ',astat
297 0 : call endrun
298 : end if
299 : end if
300 : !---------------------------------------------------------------------
301 : ! Get input variables
302 : !---------------------------------------------------------------------
303 : master_proc : &
304 0 : if( masterproc ) then
305 : !---------------------------------------------------------------------
306 : ! Allocate arrays, depending on type of input data: monthly or fft
307 : !---------------------------------------------------------------------
308 0 : if (has_monthly_data) then
309 0 : allocate( u_inp(levsiz,timesiz), p_inp(levsiz), stat=astat )
310 0 : if( astat /= 0 ) then
311 0 : write(iulog,*) 'qbo_init: failed to allocate u_inp ... p_inp; error = ',astat
312 0 : call endrun
313 : end if
314 0 : astat = nf90_get_var( ncid, uqboid, u_inp )
315 0 : if (astat/=NF90_NOERR) then
316 0 : write(iulog,*) "QBO: NF90_GET_VAR: error reading varid =", uqboid
317 0 : call endrun
318 : end if
319 0 : call wrap_get_var_realx( ncid, levid , p_inp )
320 : else
321 0 : allocate( p_inp(levsiz), stat=astat )
322 0 : if( astat /= 0 ) then
323 0 : write(iulog,*) 'qbo_init: failed to allocate p_inp; error = ',astat
324 0 : call endrun
325 : end if
326 0 : allocate( ubar_inp(levsiz), stat=astat )
327 0 : if( astat /= 0 ) then
328 0 : write(iulog,*) 'qbo_init: failed to allocate ubar_inp; error = ',astat
329 0 : call endrun
330 : end if
331 0 : allocate( cosq_inp(levsiz,coefsiz), stat=astat )
332 0 : if( astat /= 0 ) then
333 0 : write(iulog,*) 'qbo_init: failed to allocate cosq_inp; error = ',astat
334 0 : call endrun
335 : end if
336 0 : allocate( sinq_inp(levsiz,coefsiz), stat=astat )
337 0 : if( astat /= 0 ) then
338 0 : write(iulog,*) 'qbo_init: failed to allocate sinq_inp; error = ',astat
339 0 : call endrun
340 : end if
341 0 : allocate( ffreq_qbo(coefsiz), stat=astat )
342 0 : if( astat /= 0 ) then
343 0 : write(iulog,*) 'qbo_init: failed to allocate ffreq_qbo; error = ',astat
344 0 : call endrun
345 : end if
346 0 : call wrap_get_var_realx( ncid, levid , p_inp )
347 0 : call wrap_get_var_realx( ncid, ubarqboid, ubar_inp )
348 0 : astat = nf90_get_var( ncid, cosqboid, cosq_inp )
349 0 : if (astat/=NF90_NOERR) then
350 0 : write(iulog,*) "QBO: NF90_GET_VAR: error reading varid =", cosqboid
351 0 : call endrun
352 : end if
353 0 : astat = nf90_get_var( ncid, sinqboid, sinq_inp )
354 0 : if (astat/=NF90_NOERR) then
355 0 : write(iulog,*) "QBO: NF90_GET_VAR: error reading varid =", sinqboid
356 0 : call endrun
357 : end if
358 0 : call wrap_get_var_realx( ncid, freqqboid, ffreq_qbo )
359 0 : call wrap_get_scalar_int(ncid, refdatid, ref_date )
360 :
361 0 : call bnddyi( ref_date, 0, cday_ref )
362 0 : yr_ref = ref_date/10000
363 0 : cday_ref = cday_ref + yr_ref*365._r8
364 :
365 : end if
366 :
367 : !---------------------------------------------------------------------
368 : ! Convert from millibars to pascals
369 : !---------------------------------------------------------------------
370 0 : p_inp(:) = p_inp(:)*hPa2Pa
371 0 : write(iulog,*) 'qbo_init: p_inp', p_inp/hPa2Pa
372 :
373 0 : if (has_monthly_data) then
374 0 : call wrap_get_var_int( ncid, dateid, date_qbo )
375 0 : call wrap_get_var_int( ncid, secdid, secd_qbo )
376 :
377 0 : write(iulog,*) 'qbo_init: u_inp', u_inp(:,1)
378 : end if
379 :
380 : !---------------------------------------------------------------------
381 : ! Find first model layer within qbo range
382 : !---------------------------------------------------------------------
383 0 : do ktop = 1,pver
384 0 : if( pref_mid(ktop) >= p_inp(1) ) then
385 : exit
386 : end if
387 : end do
388 0 : write(iulog,*) 'qbo_init: ktop = ', ktop, pref_mid(ktop)/hPa2Pa
389 :
390 : !---------------------------------------------------------------------
391 : ! Find last model layer within qbo range
392 : !---------------------------------------------------------------------
393 0 : do kbot = pver,ktop,-1
394 0 : if( pref_mid(kbot) <= p_inp(levsiz) ) then
395 : exit
396 : end if
397 : end do
398 0 : write(iulog,*) 'qbo_init: kbot = ', kbot, pref_mid(kbot)/hPa2Pa
399 :
400 0 : if (has_monthly_data) then
401 0 : allocate( u_qbo(ktop:kbot,timesiz), stat=astat )
402 0 : if( astat /= 0 ) then
403 0 : write(iulog,*) 'qbo_init: failed to allocate u_qbo; error = ',astat
404 0 : call endrun
405 : end if
406 : else
407 0 : allocate( ubar_qbo(ktop:kbot), stat=astat )
408 0 : if( astat /= 0 ) then
409 0 : write(iulog,*) 'qbo_init: failed to allocate ubar_qbo; error = ',astat
410 0 : call endrun
411 : end if
412 0 : allocate( fcos_qbo(ktop:kbot,coefsiz), stat=astat )
413 0 : if( astat /= 0 ) then
414 0 : write(iulog,*) 'qbo_init: failed to allocate fcos_qbo; error = ',astat
415 0 : call endrun
416 : end if
417 0 : allocate( fsin_qbo(ktop:kbot,coefsiz), stat=astat )
418 0 : if( astat /= 0 ) then
419 0 : write(iulog,*) 'qbo_init: failed to allocate fsin_qbo; error = ',astat
420 0 : call endrun
421 : end if
422 : end if
423 : !---------------------------------------------------------------------
424 : ! Vertically interpolate input winds to model reference pressures
425 : !---------------------------------------------------------------------
426 0 : do k = ktop,kbot
427 0 : do kk = 1,levsiz-1
428 0 : if( p_inp(kk+1) > pref_mid(k) ) then
429 : exit
430 : end if
431 : end do
432 0 : fac1 = (pref_mid(k) - p_inp(kk+1)) / (p_inp(kk) - p_inp(kk+1))
433 0 : fac2 = (p_inp(kk) - pref_mid(k) ) / (p_inp(kk) - p_inp(kk+1))
434 0 : if (has_monthly_data) then
435 0 : u_qbo(k,:) = u_inp(kk,:)*fac1 + u_inp(kk+1,:)*fac2
436 : else
437 0 : ubar_qbo(k) = ubar_inp(kk)*fac1 + ubar_inp(kk+1)*fac2
438 0 : fcos_qbo(k,:) = cosq_inp(kk,:)*fac1 + cosq_inp(kk+1,:)*fac2
439 0 : fsin_qbo(k,:) = sinq_inp(kk,:)*fac1 + sinq_inp(kk+1,:)*fac2
440 : endif
441 : end do
442 :
443 0 : if (has_monthly_data) then
444 0 : deallocate( u_inp, p_inp )
445 0 : write(iulog,*) 'qbo_init: u', u_qbo(ktop:kbot,1)
446 : else
447 0 : deallocate( p_inp )
448 0 : deallocate( ubar_inp )
449 0 : deallocate( cosq_inp )
450 0 : deallocate( sinq_inp )
451 0 : write(iulog,*) 'qbo_init: fcos_qbo', fcos_qbo(ktop:kbot,1)
452 : end if
453 :
454 :
455 : !---------------------------------------------------------------------
456 : ! Dates are not used for cyclic QBO
457 : !---------------------------------------------------------------------
458 0 : if( has_monthly_data .and. qbo_cyclic ) then
459 0 : secd_qbo(:) = 0
460 0 : date_qbo(:) = 0
461 : end if
462 : end if master_proc
463 :
464 : #if (defined SPMD )
465 : !---------------------------------------------------------------------
466 : ! Broadcast the vertical limits
467 : !---------------------------------------------------------------------
468 0 : call mpibcast( ktop, 1, mpiint, 0, mpicom )
469 0 : call mpibcast( kbot, 1, mpiint, 0, mpicom )
470 : #endif
471 :
472 0 : if( .not. masterproc ) then
473 0 : if (has_monthly_data) then
474 0 : allocate( u_qbo(ktop:kbot,timesiz), stat=astat )
475 0 : if( astat /= 0 ) then
476 0 : write(iulog,*) 'qbo_init: failed to allocate u_qbo; error = ',astat
477 0 : call endrun
478 : end if
479 : else
480 0 : allocate( ubar_qbo(ktop:kbot), stat=astat )
481 0 : if( astat /= 0 ) then
482 0 : write(iulog,*) 'qbo_init: failed to allocate ubar_qbo; error = ',astat
483 0 : call endrun
484 : end if
485 0 : allocate( fcos_qbo(ktop:kbot,coefsiz), stat=astat )
486 0 : if( astat /= 0 ) then
487 0 : write(iulog,*) 'qbo_init: failed to allocate fcos_qbo; error = ',astat
488 0 : call endrun
489 : end if
490 0 : allocate( fsin_qbo(ktop:kbot,coefsiz), stat=astat )
491 0 : if( astat /= 0 ) then
492 0 : write(iulog,*) 'qbo_init: failed to allocate fsin_qbo; error = ',astat
493 0 : call endrun
494 : end if
495 0 : allocate( ffreq_qbo(coefsiz), stat=astat )
496 0 : if( astat /= 0 ) then
497 0 : write(iulog,*) 'qbo_init: failed to allocate ffreq_qbo; error = ',astat
498 0 : call endrun
499 : end if
500 : end if
501 : end if
502 :
503 : !---------------------------------------------------------------------
504 : ! Broadcast input data to all tasks
505 : !---------------------------------------------------------------------
506 0 : kk = kbot-ktop+1
507 0 : if (has_monthly_data) then
508 : #if (defined SPMD )
509 0 : call mpibcast( u_qbo , (kbot-ktop+1)*timesiz, mpir8, 0, mpicom )
510 0 : call mpibcast( date_qbo, timesiz, mpiint, 0, mpicom )
511 0 : call mpibcast( secd_qbo, timesiz, mpiint, 0, mpicom )
512 : #endif
513 : else
514 : #if (defined SPMD )
515 0 : call mpibcast( cday_ref , 1, mpir8, 0, mpicom )
516 : #endif
517 : #if (defined SPMD )
518 0 : call mpibcast( ubar_qbo , (kbot-ktop+1), mpir8, 0, mpicom )
519 0 : call mpibcast( fcos_qbo , (kbot-ktop+1)*coefsiz, mpir8, 0, mpicom )
520 0 : call mpibcast( fsin_qbo , (kbot-ktop+1)*coefsiz, mpir8, 0, mpicom )
521 : #endif
522 : #if (defined SPMD )
523 0 : call mpibcast( ffreq_qbo , coefsiz, mpir8, 0, mpicom )
524 : #endif
525 : endif
526 :
527 : !---------------------------------------------------------------------
528 : ! setup vertical factor
529 : !---------------------------------------------------------------------
530 0 : allocate( tauz(ktop-1:kbot+1), stat=astat )
531 0 : if( astat /= 0 ) then
532 0 : write(iulog,*) 'qbo_init: failed to allocate tauz; error = ',astat
533 0 : call endrun
534 : end if
535 :
536 0 : tauz(ktop-1) = 2._r8
537 0 : tauz(kbot+1) = 2._r8
538 0 : tauz(ktop:kbot) = 1._r8
539 :
540 : !---------------------------------------------------------------------
541 : ! Get current day of year and date
542 : !---------------------------------------------------------------------
543 0 : calday = get_curr_calday()
544 0 : call get_curr_date( yr, mon, day, ncsec )
545 0 : if (masterproc) write(iulog,*) 'qbo_init: get_curr_date = ', yr,mon,day,ncsec
546 :
547 : !---------------------------------------------------------------------
548 : ! Set current day in run or in qbo cycle
549 : !---------------------------------------------------------------------
550 0 : cday = calday + yr*365._r8
551 0 : if (masterproc) write(iulog,*) 'qbo_init: cday = ', cday
552 0 : if( has_monthly_data .and. qbo_cyclic ) then
553 0 : cday = mod( cday,qbo_days )
554 : end if
555 0 : if (masterproc) write(iulog,*) 'qbo_init: cday = ', cday
556 :
557 0 : if( has_monthly_data ) then
558 0 : if( qbo_cyclic ) then
559 : !---------------------------------------------------------------------
560 : ! Set up cyclic qbo
561 : !---------------------------------------------------------------------
562 : ! Set past and future month indexes
563 : !---------------------------------------------------------------------
564 0 : if( cday == 0._r8 ) then
565 0 : nm = qbo_mons - 1
566 0 : np = qbo_mons
567 : else
568 0 : nm = cday / qbo_dypm
569 0 : np = mod( nm,qbo_mons ) + 1
570 0 : if( nm == 0 ) then
571 0 : nm = qbo_mons
572 : end if
573 : end if
574 0 : if (masterproc) write(iulog,*) 'qbo_init: nm,np = ', nm,np
575 : !---------------------------------------------------------------------
576 : ! Set past and future days for data, generate day for cyclic data
577 : !---------------------------------------------------------------------
578 0 : cdayp = np * qbo_dypm
579 0 : cdaym = nm * qbo_dypm
580 : else
581 : !---------------------------------------------------------------------
582 : ! Set up noncyclic qbo
583 : !---------------------------------------------------------------------
584 0 : found = .false.
585 0 : do nm = 1, timesiz-1
586 0 : np = nm+1
587 0 : call bnddyi( date_qbo(nm), secd_qbo(nm), cdaym )
588 0 : call bnddyi( date_qbo(np), secd_qbo(np), cdayp )
589 0 : yr = date_qbo(nm)/10000
590 0 : cdaym = cdaym + yr*365._r8
591 0 : yr = date_qbo(np)/10000
592 0 : cdayp = cdayp + yr*365._r8
593 0 : if (masterproc) write(iulog,*) 'qbo_init: nm, date_qbo(nm), date_qbo(np), cdaym, cdayp = ', &
594 0 : nm, date_qbo(nm), date_qbo(np), cdaym, cdayp
595 0 : if( cday >= cdaym .and. cday < cdayp ) then
596 : found = .true.
597 : exit
598 : end if
599 : end do
600 0 : if( .not. found ) then
601 0 : call endrun( 'QBO_INIT: failed to find bracketing dates' )
602 : end if
603 : end if
604 0 : if (masterproc) write(iulog,*) 'qbo_init: cdaym,cdayp = ', cdaym,cdayp
605 : endif
606 :
607 : !---------------------------------------------------------------------
608 : ! Initialize output buffer for two fields: QBO forcing wind and wind tendency of qbo relaxation
609 : !----------------------------------------------------------------------
610 : !
611 : !----------------------------------------------------------------------
612 0 : call addfld ('QBOTEND',(/ 'lev' /), 'A','M/S/S','Wind tendency from QBO relaxation')
613 0 : call addfld ('QBO_U0', (/ 'lev' /), 'A','M/S','Specified wind used for QBO')
614 :
615 0 : if (history_waccm) then
616 0 : call add_default ('QBOTEND', 1, ' ')
617 0 : call add_default ('QBO_U0', 1, ' ')
618 : end if
619 :
620 : !----------------------------------------------------------------------
621 : ! Get zonal mean zonal wind index in pbuf.
622 : !----------------------------------------------------------------------
623 0 : uzm_idx = pbuf_get_index("UZM")
624 :
625 0 : if (masterproc) write(iulog,*) 'end of qbo_init'
626 :
627 1536 : end subroutine qbo_init
628 :
629 : !================================================================================================
630 :
631 16128 : subroutine qbo_timestep_init
632 : !-----------------------------------------------------------------------
633 : !
634 : ! Purpose: Interpolate QBO zonal wind to current time
635 : !
636 : ! Method: Linear interpolation between dates on QBO file,
637 : ! vertically and horizontally
638 : !
639 : ! Author:
640 : !
641 : !-----------------------------------------------------------------------
642 :
643 : !-----------------------------------------------------------------------
644 : ! Local variables
645 : !-----------------------------------------------------------------------
646 : integer :: k ! level index
647 : integer :: yr, mon, day ! components of a date
648 : integer :: yrl, monl, dayl ! components of a date
649 : integer :: ncdate ! current date in integer format [yyyymmdd]
650 : integer :: ncsec ! current time of day [seconds]
651 : integer :: ncsecl ! current time of day [seconds]
652 :
653 : real(r8) :: fact1, fact2 ! time interpolation factors
654 : real(r8) :: calday ! day of year at end of present time step
655 : real(r8) :: caldayl ! day of year at begining of present time step
656 : real(r8) :: cday ! day within qbo period at end of present time step
657 : real(r8) :: cdayl ! day within qbo period at begining of present time step
658 : real(r8) :: deltat ! time difference (days) between cdaym and cdayp
659 :
660 : integer :: n ! coefficient index
661 : real(r8) :: ccc ! cosine term for expanding coefficients
662 : real(r8) :: sss ! sine term for expanding coefficients
663 :
664 : logical :: new_interval ! flag for new time interval
665 :
666 : has_qbo_forcing : &
667 16128 : if( qbo_use_forcing ) then
668 : !-----------------------------------------------------------------------
669 : ! Get current day of year and date
670 : !-----------------------------------------------------------------------
671 0 : caldayl = get_curr_calday( -delt )
672 0 : call get_curr_date( yrl, monl, dayl, ncsecl, -delt )
673 0 : calday = get_curr_calday()
674 0 : call get_curr_date( yr, mon, day, ncsec )
675 0 : ncdate = yr*10000 + mon*100 + day
676 : #ifdef QBO_DIAGS
677 : write(iulog,*) 'qbo_timestep_init: calday = ', calday
678 : write(iulog,*) 'qbo_timestep_init: ncdate = ', ncdate
679 : #endif
680 :
681 : !-----------------------------------------------------------------------
682 : ! Set current day in run or in qbo cycle
683 : !-----------------------------------------------------------------------
684 0 : cday = calday + yr*365._r8
685 0 : cdayl = caldayl + yrl*365._r8
686 : #ifdef QBO_DIAGS
687 : write(iulog,*) 'qbo_timestep_init: cday = ', cday
688 : #endif
689 :
690 0 : if( has_monthly_data ) then
691 : !-----------------------------------------------------------------------
692 : ! Time interpolation for cases with monthly input data
693 : !-----------------------------------------------------------------------
694 0 : if( .not. qbo_cyclic ) then
695 0 : new_interval = cday > cdayp
696 : else
697 0 : cday = mod( cday,qbo_days )
698 0 : cdayl = mod( cdayl,qbo_days )
699 0 : if( cday > cdayl ) then
700 0 : new_interval = cday > cdayp
701 : else
702 : new_interval = .true.
703 : end if
704 : end if
705 : #ifdef QBO_DIAGS
706 : write(iulog,*) 'qbo_timestep_init: cday = ', cday
707 : #endif
708 : !-----------------------------------------------------------------------
709 : ! If model time is past current forward timeslice, then switch to next one.
710 : ! If qbo_cyclic = .true. interpolation between end and beginning of data (np == 1).
711 : ! Note that np is never 1 when qbo_cyclic is .false.
712 : !-----------------------------------------------------------------------
713 : next_interval : &
714 0 : if( new_interval ) then
715 :
716 : !-----------------------------------------------------------------------
717 : ! Increment index of future time sample
718 : !-----------------------------------------------------------------------
719 0 : if( qbo_cyclic ) then
720 0 : np1 = mod( np,qbo_mons ) + 1
721 : else
722 0 : np1 = np + 1
723 : end if
724 0 : if( np1 > timesiz ) then
725 0 : call endrun ('QBO_TIMESTEP_INIT: Attempt to go past end of QBO data')
726 : end if
727 : !-----------------------------------------------------------------------
728 : ! Set indexes into u table
729 : !-----------------------------------------------------------------------
730 0 : nm = np
731 0 : np = np1
732 : #ifdef QBO_DIAGS
733 : write(iulog,*) 'qbo_timestep_init: nm,np = ', nm,np
734 : write(iulog,*) 'qbo_timestep_init: date_qbo(np), secd_qbo(np) = ', date_qbo(np), secd_qbo(np)
735 : #endif
736 :
737 : !-----------------------------------------------------------------------
738 : ! Set past and future days for data, generate day for cyclic data
739 : !-----------------------------------------------------------------------
740 0 : cdaym = cdayp
741 0 : if( qbo_cyclic ) then
742 0 : cdayp = np * qbo_dypm
743 : else
744 0 : call bnddyi( date_qbo(np), secd_qbo(np), cdayp )
745 0 : yr = date_qbo(np)/10000
746 0 : cdayp = cdayp + yr*365._r8
747 : end if
748 : #ifdef QBO_DIAGS
749 : write(iulog,*) 'qbo_timestep_init: cdaym,cdayp = ', cdaym,cdayp
750 : #endif
751 0 : if( np /= 1 .and. cday > cdayp ) then
752 0 : write(iulog,*) 'qbo_timestep_init: Input qbo for date',date_qbo(np),' sec ',secd_qbo(np), &
753 0 : 'does not exceed model date',ncdate,' sec ',ncsec,' Stopping.'
754 0 : call endrun
755 : end if
756 : end if next_interval
757 :
758 : !-----------------------------------------------------------------------
759 : ! Determine time interpolation factors. Account for December-January
760 : ! interpolation if dataset is being cycled.
761 : !-----------------------------------------------------------------------
762 0 : if( qbo_cyclic .and. np == 1 ) then ! Dec-Jan interpolation
763 0 : deltat = cdayp + qbo_days - cdaym
764 0 : if (cday > cdayp) then ! We are in December
765 0 : fact1 = (cdayp + qbo_days - cday)/deltat
766 0 : fact2 = (cday - cdaym)/deltat
767 : else ! We are in January
768 0 : fact1 = (cdayp - cday)/deltat
769 0 : fact2 = (cday + qbo_days - cdaym)/deltat
770 : end if
771 : else
772 0 : deltat = cdayp - cdaym
773 0 : fact1 = (cdayp - cday )/deltat
774 0 : fact2 = (cday - cdaym)/deltat
775 : end if
776 : #ifdef QBO_DIAGS
777 : write(iulog,*) 'qbo_timestep_init: fact1,fact2 = ', fact1, fact2
778 : #endif
779 :
780 : !-----------------------------------------------------------------------
781 : ! Time interpolation
782 : !-----------------------------------------------------------------------
783 0 : do k = ktop, kbot
784 0 : u_tstep(k) = u_qbo(k,nm)*fact1 + u_qbo(k,np)*fact2
785 : end do
786 0 : if( ktop > 1 ) then
787 0 : u_tstep(ktop-1) = u_tstep(ktop)
788 : end if
789 0 : if( kbot < pver ) then
790 0 : u_tstep(kbot+1) = u_tstep(kbot)
791 : end if
792 :
793 : else
794 :
795 : !-----------------------------------------------------------------------
796 : ! Wind at this timestep for fft input data
797 : !-----------------------------------------------------------------------
798 : !-----------------------------------------------------------------------
799 : ! Set past and future days for data, generate winds for current day from fft data
800 : !-----------------------------------------------------------------------
801 :
802 0 : do k = ktop, kbot
803 0 : u_tstep(k) = ubar_qbo(k)
804 : end do
805 0 : do n=1,coefsiz
806 0 : ccc = cos(ffreq_qbo(n)*(cday-cday_ref))
807 0 : sss = sin(ffreq_qbo(n)*(cday-cday_ref))
808 0 : do k = ktop, kbot
809 0 : u_tstep(k) = u_tstep(k) + fcos_qbo(k,n)*ccc + fsin_qbo(k,n)*sss
810 : end do
811 : end do
812 0 : if( ktop > 1 ) then
813 0 : u_tstep(ktop-1) = u_tstep(ktop)
814 : end if
815 0 : if( kbot < pver ) then
816 0 : u_tstep(kbot+1) = u_tstep(kbot)
817 : end if
818 : end if
819 :
820 : #ifdef QBO_DIAGS
821 : write(iulog,*) 'qbo_timestep_init: u_tstep ', u_tstep(ktop:kbot)
822 : #endif
823 :
824 : end if has_qbo_forcing
825 :
826 1536 : end subroutine qbo_timestep_init
827 :
828 : !================================================================================================
829 :
830 14810880 : subroutine qbo_relax( state, pbuf, ptend )
831 : use physics_buffer, only: physics_buffer_desc, pbuf_get_field
832 : use cam_history, only: outfld
833 : !------------------------------------------------------------------------
834 : ! relax zonal mean wind towards qbo sequence
835 : !------------------------------------------------------------------------
836 :
837 : !--------------------------------------------------------------------------------
838 : ! ... dummy arguments
839 : !--------------------------------------------------------------------------------
840 : type(physics_state), intent(in) :: state ! Physics state variables
841 : type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer
842 : type(physics_ptend), intent(out) :: ptend ! individual parameterization tendencies
843 :
844 : !--------------------------------------------------------------------------------
845 : ! Local variables
846 : !--------------------------------------------------------------------------------
847 : integer :: lchnk ! chunk identifier
848 : integer :: ncol ! number of atmospheric columns
849 : integer :: i, k ! loop indices
850 : integer :: kl, ku ! loop indices
851 :
852 : real(r8) :: tauxi(pcols) ! latitudes in radians for present chunk
853 : real(r8) :: tauzz
854 : real(r8) :: u
855 : real(r8) :: rlat ! latitudes in radians for present chunk
856 : real(r8) :: crelax ! relaxation constant
857 :
858 : real(r8), parameter :: tconst = 10._r8 ! relaxation time constant in days
859 : real(r8), parameter :: tconst1 = tconst * 86400._r8
860 : real(r8) :: qbo_u0(pcols,pver) ! QBO wind used for driving parameterization
861 :
862 : ! Zonal mean zonal wind from pbuf.
863 72960 : real(r8), pointer :: uzm(:,:)
864 :
865 72960 : lchnk = state%lchnk
866 72960 : ncol = state%ncol
867 :
868 72960 : call physics_ptend_init(ptend, state%psetcols, 'qbo', lu=.true.)
869 :
870 : has_qbo_forcing : &
871 72960 : if( qbo_use_forcing ) then
872 :
873 0 : call pbuf_get_field(pbuf, uzm_idx, uzm)
874 :
875 0 : kl = max( 1,ktop-1 )
876 0 : ku = min( pver,kbot+1 )
877 : !--------------------------------------------------------------------------------
878 : ! get latitude in radians for present chunk
879 : !--------------------------------------------------------------------------------
880 0 : do i = 1,ncol
881 0 : rlat = get_rlat_p( lchnk, i )
882 0 : tauxi(i) = tconst1*taux( rlat )
883 : end do
884 :
885 0 : qbo_u0(:,:) = 0._r8
886 :
887 0 : do k = kl,ku
888 0 : tauzz = tauz(k)
889 0 : u = u_tstep(k)
890 0 : do i = 1,ncol
891 : !--------------------------------------------------------------------------------
892 : ! determine relaxation constant
893 : !--------------------------------------------------------------------------------
894 0 : crelax = tauxi(i)*tauzz
895 0 : if( crelax /= 0._r8 ) then
896 0 : crelax = 1._r8 / crelax
897 : !--------------------------------------------------------------------------------
898 : ! do relaxation of zonal mean wind
899 : !--------------------------------------------------------------------------------
900 :
901 0 : if(u < 50.0_r8) then
902 0 : ptend%u(i,k) = crelax * (u - uzm(i,k))
903 : end if
904 : end if
905 : !--------------------------------------------------------------------------------
906 : ! variable representing QBO wind
907 : !--------------------------------------------------------------------------------
908 0 : if((u < 50.0_r8) .and. (crelax /= 0._r8)) then
909 0 : qbo_u0(i,k) = u/tauzz/tauxi(i)*tconst1
910 : end if
911 : end do
912 : end do
913 :
914 : !--------------------------------------------------------------------------------
915 : !output tendency of relaxation to monthly ('h1') output file
916 : !--------------------------------------------------------------------------------
917 0 : call outfld( 'QBOTEND', ptend%u(:,:), pcols, lchnk )
918 :
919 : !--------------------------------------------------------------------------------
920 : !output specified QBO wind to h0 output file
921 : !--------------------------------------------------------------------------------
922 0 : call outfld( 'QBO_U0', qbo_u0, pcols, lchnk )
923 : end if has_qbo_forcing
924 :
925 145920 : end subroutine qbo_relax
926 :
927 : !================================================================================================
928 :
929 0 : function taux( rlat )
930 : !------------------------------------------------------------------------
931 : ! calculates relaxation constant in latitude
932 : !------------------------------------------------------------------------
933 :
934 : !------------------------------------------------------------------------
935 : ! ... dummy arguments
936 : !------------------------------------------------------------------------
937 : real(r8), intent(in) :: rlat ! latitude in radians for present chunk
938 :
939 : !------------------------------------------------------------------------
940 : ! ... local variables
941 : !------------------------------------------------------------------------
942 : real(r8), parameter :: factor = 1._r8/(2._r8*0.174532925_r8*0.174532925_r8)
943 : real(r8) :: alat ! abs rlat
944 :
945 : !------------------------------------------------------------------------
946 : ! ... function declaration
947 : !------------------------------------------------------------------------
948 : real(r8) :: taux ! relaxation constant in latitude
949 :
950 :
951 0 : alat = abs( rlat )
952 0 : if( alat <= .035_r8 ) then
953 : !------------------------------------------------------------------------
954 : ! rlat=0.035 (latitude in radians): rlat*180/pi=2 degrees, around equator full relaxation
955 : !------------------------------------------------------------------------
956 : taux = 1._r8
957 0 : else if( alat <= .384_r8) then
958 : !------------------------------------------------------------------------
959 : ! from 6 to 22 degrees latitude weakening of relaxation with Gaussian distribution
960 : ! half width=10° => in radians: 0.174532925
961 : !------------------------------------------------------------------------
962 0 : taux = exp( rlat*rlat*factor )
963 : else
964 : !------------------------------------------------------------------------
965 : ! other latitudes no relaxation
966 : !------------------------------------------------------------------------
967 : taux = 0._r8
968 : end if
969 :
970 72960 : end function taux
971 :
972 : end module qbo
|
