Line data Source code
1 : module zonal_mean_mod
2 : !======================================================================
3 : !
4 : ! Purpose: Compute and make use of Zonal Mean values on physgrid
5 : !
6 : ! This module implements 3 data structures for the spectral analysis
7 : ! and synthesis of zonal mean values based on m=0 spherical harmonics.
8 : !
9 : ! ZonalMean_t: For the analysis/synthesis of zonal mean values
10 : ! on a 2D grid of points distributed over the
11 : ! surface of a sphere.
12 : ! ZonalProfile_t: For the analysis/synthesis of zonal mean values
13 : ! on a meridional grid that spans the latitudes
14 : ! from SP to NP
15 : ! ZonalAverage_t: To calculate zonal mean values via a simple
16 : ! area weighted bin-averaging of 2D grid points
17 : ! assigned to each latitude band.
18 : !
19 : ! NOTE: The weighting of the Zonal Profiles values is scaled such
20 : ! that ZonalMean_t amplitudes can be used to evaluate values
21 : ! on the ZonalProfile_t grid and vice-versa.
22 : !
23 : ! The ZonalMean_t computes global integrals to compute basis
24 : ! amplitudes. For distributed environments the cost of these
25 : ! can be reduced using the The ZonalAverage_t data structures.
26 : !
27 : ! USAGE:
28 : !
29 : ! (1) Compute Zonal mean amplitudes and synthesize values on 2D/3D physgrid
30 : !
31 : ! Usage: type(ZonalMean_t):: ZM
32 : ! =========================================
33 : ! call ZM%init(nbas)
34 : ! ------------------
35 : ! - Initialize the data structure with 'nbas' basis functions
36 : ! for the given physgrid latitudes and areas.
37 : !
38 : ! Arguments:
39 : ! integer ,intent(in):: nbas -Number of m=0 spherical harmonics
40 : !
41 : ! call ZM%calc_amps(Gdata,Bamp)
42 : ! -----------------------------
43 : ! - For the initialized ZonalMean_t; Given Gdata() values on the physgrid,
44 : ! compute the zonal mean basis amplitudes Bamp().
45 : !
46 : ! Interface: 2D data on the physgrid
47 : ! real(r8),intent(in ):: Gdata(pcols,begchunk:endchunk)
48 : ! real(r8),intent(out):: Bamp (nbas)
49 : !
50 : ! Interface: 3D data on the physgrid
51 : ! real(r8),intent(in ):: Gdata(pcols,pver,begchunk:endchunk)
52 : ! real(r8),intent(out):: Bamp (nbas,pver)
53 : !
54 : ! call ZM%eval_grid(Bamp,Gdata)
55 : ! -----------------------------
56 : ! - For the initialized ZonalMean_t; Given Bamp() zonal mean basis
57 : ! amplitudes, compute the Gdata() values on the physgrid.
58 : !
59 : ! Interface: 2D data on the physgrid
60 : ! real(r8),intent(in ):: Bamp (nbas)
61 : ! real(r8),intent(out):: Gdata(pcols,begchunk:endchunk)
62 : !
63 : ! Interface: 3D data on the physgrid
64 : ! real(r8),intent(in ):: Bamp (nbas,pver)
65 : ! real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk)
66 : !
67 : !
68 : ! (2) Compute Zonal mean amplitudes and synthesize values on Zonal profile grid
69 : !
70 : ! Usage: type(ZonalProfile_t):: ZP
71 : ! =========================================
72 : ! call ZP%init(lats,area,nlat,nbas,GEN_GAUSSLATS=.true.)
73 : ! ------------------------------------------------------
74 : ! - Initialize the data structure for the given number of
75 : ! latitudes. Either use the given Latitudes and weights,
76 : ! or OPTIONALLY create profile gridpoints and associated
77 : ! area weights from SP to NP. Then initialize 'nbas' basis
78 : ! functions for the profile gridpoints.
79 : ! If the user supplies the lats/area values, the area values must
80 : ! be correctly scaled such that the global area adds up to 4PI.
81 : ! Otherwise, the ampitudes between ZonalProfile_t and ZonalMean_t
82 : ! are not interchangable.
83 : !
84 : ! Arguments:
85 : ! real(r8),intent(inout):: lats(:) - Latitudes of meridional grid.
86 : ! real(r8),intent(inout):: area(:) - Area of each meridional gridpoint.
87 : ! integer ,intent(in) :: nlat - Number of meridional gridpoints.
88 : ! integer ,intent(in) :: nbas - Number of m=0 spherical harmonics
89 : ! logical ,intent(in),optional:: GEN_GAUSLATS - Flag to generate
90 : ! lats/areas values.
91 : !
92 : ! call ZP%calc_amps(Zdata,Bamp)
93 : ! -----------------------------
94 : ! - Given Zdata() on the Zonal profile grid, compute the
95 : ! zonal basis amplitudes Bamp().
96 : !
97 : ! Interface: 1D data on (nlat) grid
98 : ! real(r8),intent(in ):: Zdata(nlat) - Meridional Profile data
99 : ! real(r8),intent(out):: Bamp (nbas) - Zonal Basis Amplitudes
100 : !
101 : ! Interface: 2D data on (nlat,pver) grid
102 : ! real(r8),intent(in ):: Zdata(nlat,pver) - Meridional Profile data
103 : ! real(r8),intent(out):: Bamp (nbas,pver) - Zonal Basis Amplitudes
104 : !
105 : ! call ZP%eval_grid(Bamp,Zdata)
106 : ! -----------------------------
107 : ! - Given Bamp() zonal basis amplitudes, evaluate the Zdata()
108 : ! values on the Zonal profile grid.
109 : !
110 : ! Interface: 1D data on (nlat) grid
111 : ! real(r8),intent(in ):: Bamp (nbas) - Zonal Basis Amplitudes
112 : ! real(r8),intent(out):: Zdata(nlat) - Meridional Profile data
113 : !
114 : ! Interface: 2D data on (nlat,pver) grid
115 : ! real(r8),intent(in ):: Bamp (nbas,pver) - Zonal Basis Amplitudes
116 : ! real(r8),intent(out):: Zdata(nlat,pver) - Meridional Profile data
117 : !
118 : ! (3) Compute Zonal mean averages (FASTER/LESS-ACCURATE) on Zonal profile grid
119 : ! (For the created zonal profile, just bin average area weighted
120 : ! 2D/3D physgrid grid values)
121 : !
122 : ! Usage: type(ZonalAverage_t):: ZA
123 : ! =========================================
124 : ! call ZA%init(lats,area,nlat,GEN_GAUSSLATS=.true.)
125 : ! --------------------------------------------------
126 : ! - Given the latitude/area for the nlat meridional gridpoints, initialize
127 : ! the ZonalAverage data structure for computing bin-averaging of physgrid
128 : ! values. It is assumed that the domain of these gridpoints of the
129 : ! profile span latitudes from SP to NP.
130 : ! The optional GEN_GAUSSLATS flag allows for the generation of Gaussian
131 : ! latitude gridpoints. The generated grid over-writes the given values
132 : ! lats and area passed by the user.
133 : !
134 : ! Arguments:
135 : ! real(r8),intent(inout):: lats(nlat) - Latitudes of meridional grid.
136 : ! real(r8),intent(inout):: area(nlat) - Area of meridional gridpoints.
137 : ! integer ,intent(in):: nlat - Number of meridional gridpoints
138 : ! logical,intent(in),optional:: GEN_GAUSLATS - Flag to generate
139 : ! lats/areas values.
140 : !
141 : ! call ZA%binAvg(Gdata,Zdata)
142 : ! ---------------------------
143 : ! - For the initialized ZonalAverage_t; Given Gdata() on the physgrid,
144 : ! compute bin averages and return Zdata() on the Zonal profile grid.
145 : !
146 : ! Interface: 2D data on the physgrid
147 : ! real(r8),intent(out):: Gdata(pcols,begchunk:endchunk)
148 : ! real(r8),intent(out):: Zdata(nlat)
149 : !
150 : ! Interface: 3D data on the physgrid
151 : ! real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk)
152 : ! real(r8),intent(out):: Zdata(nlat,pver)
153 : !
154 : !======================================================================
155 :
156 : use shr_kind_mod, only: r8=>SHR_KIND_R8
157 : use phys_grid, only: get_ncols_p, get_rlat_p, get_wght_all_p, get_nlcols_p
158 : use ppgrid, only: begchunk, endchunk, pcols
159 : use shr_reprosum_mod,only: shr_reprosum_calc
160 : use cam_abortutils, only: endrun, handle_allocate_error
161 : use spmd_utils, only: mpicom
162 : use physconst, only: pi
163 : use phys_grid, only: ngcols_p => num_global_phys_cols
164 : use cam_logfile, only: iulog
165 :
166 : implicit none
167 : private
168 :
169 : public :: ZonalMean_t
170 : public :: ZonalProfile_t
171 : public :: ZonalAverage_t
172 :
173 : ! Type definitions
174 : !-------------------
175 : type ZonalMean_t
176 : private
177 : integer :: nbas
178 : real(r8),allocatable:: area (:,:)
179 : real(r8),allocatable:: basis(:,:,:)
180 : real(r8),allocatable:: map (:,:)
181 : contains
182 : procedure,pass:: init => init_ZonalMean
183 : generic,public:: calc_amps => calc_ZonalMean_2Damps, &
184 : calc_ZonalMean_3Damps
185 : generic,public:: eval_grid => eval_ZonalMean_2Dgrid, &
186 : eval_ZonalMean_3Dgrid
187 : procedure,private,pass:: calc_ZonalMean_2Damps
188 : procedure,private,pass:: calc_ZonalMean_3Damps
189 : procedure,private,pass:: eval_ZonalMean_2Dgrid
190 : procedure,private,pass:: eval_ZonalMean_3Dgrid
191 : procedure, pass :: final => final_ZonalMean
192 : end type ZonalMean_t
193 :
194 : type ZonalProfile_t
195 : private
196 : integer :: nlat
197 : integer :: nbas
198 : real(r8),allocatable:: area (:)
199 : real(r8),allocatable:: basis(:,:)
200 : real(r8),allocatable:: map (:,:)
201 : contains
202 : procedure,pass:: init => init_ZonalProfile
203 : generic,public:: calc_amps => calc_ZonalProfile_1Damps, &
204 : calc_ZonalProfile_2Damps
205 : generic,public:: eval_grid => eval_ZonalProfile_1Dgrid, &
206 : eval_ZonalProfile_2Dgrid
207 : procedure,private,pass:: calc_ZonalProfile_1Damps
208 : procedure,private,pass:: calc_ZonalProfile_2Damps
209 : procedure,private,pass:: eval_ZonalProfile_1Dgrid
210 : procedure,private,pass:: eval_ZonalProfile_2Dgrid
211 : procedure, pass :: final => final_ZonalProfile
212 : end type ZonalProfile_t
213 :
214 : type ZonalAverage_t
215 : private
216 : integer :: nlat
217 : real(r8),allocatable:: area (:)
218 : real(r8),allocatable:: a_norm (:)
219 : real(r8),allocatable:: area_g (:,:)
220 : integer ,allocatable:: idx_map(:,:)
221 : contains
222 : procedure,pass:: init => init_ZonalAverage
223 : generic,public:: binAvg => calc_ZonalAverage_2DbinAvg, &
224 : calc_ZonalAverage_3DbinAvg
225 : procedure,private,pass:: calc_ZonalAverage_2DbinAvg
226 : procedure,private,pass:: calc_ZonalAverage_3DbinAvg
227 : procedure, pass :: final => final_ZonalAverage
228 : end type ZonalAverage_t
229 :
230 : real(r8), parameter :: halfPI = 0.5_r8*pi
231 : real(r8), parameter :: twoPI = 2.0_r8*pi
232 : real(r8), parameter :: fourPI = 4.0_r8*pi
233 : real(r8), parameter :: qrtrPI = 0.25_r8*pi
234 : real(r8), parameter :: invSqrt4pi = 1._r8/sqrt(fourPI)
235 :
236 : contains
237 : !=======================================================================
238 0 : subroutine init_ZonalMean(this,I_nbas)
239 : !
240 : ! init_ZonalMean: Initialize the ZonalMean data structures for the
241 : ! physics grid. It is assumed that the domain
242 : ! of these gridpoints spans the surface of the sphere.
243 : ! The representation of basis functions is
244 : ! normalized w.r.t integration over the sphere.
245 : !=====================================================================
246 : !
247 : ! Passed Variables
248 : !------------------
249 : class(ZonalMean_t) :: this
250 : integer ,intent(in):: I_nbas
251 : !
252 : ! Local Values
253 : !--------------
254 0 : real(r8),allocatable:: Clats(:,:)
255 0 : real(r8),allocatable:: Bcoef(:)
256 0 : real(r8),allocatable:: Csum (:,:)
257 0 : real(r8),allocatable:: Cvec (:)
258 0 : real(r8),allocatable:: Bsum (:,:)
259 0 : real(r8),allocatable:: Bnorm(:)
260 0 : real(r8),allocatable:: Bcov (:,:)
261 : real(r8):: area(pcols),rlat
262 :
263 : integer :: nn,n2,nb,lchnk,ncols,cc
264 : integer :: cnum,Cvec_len
265 :
266 : integer :: nlcols, count, astat
267 : character(len=*), parameter :: subname = 'init_ZonalMean'
268 :
269 0 : if (I_nbas<1) then
270 0 : call endrun('ZonalMean%init: ERROR I_nbas must be greater than 0')
271 : end if
272 :
273 : ! Allocate space
274 : !-----------------
275 0 : if(allocated(this%area )) deallocate(this%area)
276 0 : if(allocated(this%basis)) deallocate(this%basis)
277 0 : if(allocated(this%map )) deallocate(this%map)
278 :
279 0 : this%nbas = I_nbas
280 0 : allocate(this%area (pcols,begchunk:endchunk), stat=astat)
281 0 : call handle_allocate_error(astat, subname, 'this%area')
282 0 : allocate(this%basis(pcols,begchunk:endchunk,I_nbas), stat=astat)
283 0 : call handle_allocate_error(astat, subname, 'this%basis')
284 0 : allocate(this%map (I_nbas,I_nbas), stat=astat)
285 0 : call handle_allocate_error(astat, subname, 'this%map')
286 0 : this%area (:,:) = 0._r8
287 0 : this%basis(:,:,:) = 0._r8
288 0 : this%map (:,:) = 0._r8
289 :
290 0 : Cvec_len = 0
291 0 : do nn= 1,this%nbas
292 0 : do n2=nn,this%nbas
293 0 : Cvec_len = Cvec_len + 1
294 : end do
295 : end do
296 :
297 0 : nlcols = get_nlcols_p()
298 :
299 0 : allocate(Clats(pcols,begchunk:endchunk), stat=astat)
300 0 : call handle_allocate_error(astat, subname, 'Clats')
301 0 : allocate(Bcoef(I_nbas/2+1), stat=astat)
302 0 : call handle_allocate_error(astat, subname, 'Bcoef')
303 0 : allocate(Csum (nlcols, Cvec_len), stat=astat)
304 0 : call handle_allocate_error(astat, subname, 'Csum')
305 0 : allocate(Cvec (Cvec_len), stat=astat)
306 0 : call handle_allocate_error(astat, subname, 'Cvec')
307 0 : allocate(Bsum (nlcols, I_nbas), stat=astat)
308 0 : call handle_allocate_error(astat, subname, 'Bsum')
309 0 : allocate(Bnorm(I_nbas), stat=astat)
310 0 : call handle_allocate_error(astat, subname, 'Bnorm')
311 0 : allocate(Bcov (I_nbas,I_nbas), stat=astat)
312 0 : call handle_allocate_error(astat, subname, 'Bcov')
313 :
314 0 : Bsum(:,:) = 0._r8
315 0 : Csum(:,:) = 0._r8
316 :
317 : ! Save a copy of the area weights for each ncol gridpoint
318 : ! and convert Latitudes to SP->NP colatitudes in radians
319 : !-------------------------------------------------------
320 0 : do lchnk=begchunk,endchunk
321 0 : ncols = get_ncols_p(lchnk)
322 0 : call get_wght_all_p(lchnk, ncols, area)
323 0 : do cc = 1,ncols
324 0 : rlat=get_rlat_p(lchnk,cc)
325 0 : this%area(cc,lchnk) = area(cc)
326 0 : Clats (cc,lchnk) = rlat + halfPI
327 : end do
328 : end do
329 :
330 : ! Add first basis for the mean values.
331 : !------------------------------------------
332 0 : this%basis(:,begchunk:endchunk,1) = invSqrt4pi
333 :
334 : ! Loop over the remaining basis functions
335 : !---------------------------------------
336 0 : do nn=2,this%nbas
337 0 : nb = nn-1
338 :
339 : ! Generate coefs for the basis
340 : !------------------------------
341 0 : call sh_gen_basis_coefs(nb,0,Bcoef)
342 :
343 : ! Create basis for the coefs at each ncol gridpoint
344 : !---------------------------------------------------
345 0 : do lchnk=begchunk,endchunk
346 0 : ncols = get_ncols_p(lchnk)
347 0 : do cc = 1,ncols
348 0 : call sh_create_basis(nb,0,Clats(cc,lchnk),Bcoef,this%basis(cc,lchnk,nn))
349 : end do
350 : end do
351 : end do ! nn=2,this%nbas
352 :
353 : ! Numerically normalize the basis funnctions
354 : !--------------------------------------------------------------
355 0 : do nn=1,this%nbas
356 0 : count = 0
357 0 : do lchnk=begchunk,endchunk
358 0 : ncols = get_ncols_p(lchnk)
359 0 : do cc = 1,ncols
360 0 : count=count+1
361 0 : Bsum(count,nn) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk)
362 : end do
363 : end do
364 : end do ! nn=1,this%nbas
365 :
366 0 : call shr_reprosum_calc(Bsum, Bnorm, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom)
367 :
368 0 : do nn=1,this%nbas
369 0 : do lchnk=begchunk,endchunk
370 0 : ncols = get_ncols_p(lchnk)
371 0 : this%basis(:ncols,lchnk,nn) = this%basis(:ncols,lchnk,nn)/sqrt(Bnorm(nn))
372 : end do
373 : end do ! nn=1,this%nbas
374 :
375 : ! Compute covariance matrix for basis functions
376 : ! (Yes, they are theoretically orthonormal, but lets make sure)
377 : !---------------------------------------------------------------
378 0 : cnum = 0
379 0 : do nn= 1,this%nbas
380 0 : do n2=nn,this%nbas
381 0 : cnum = cnum + 1
382 0 : count = 0
383 0 : do lchnk=begchunk,endchunk
384 0 : ncols = get_ncols_p(lchnk)
385 0 : do cc = 1,ncols
386 0 : count=count+1
387 0 : Csum(count,cnum) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,n2)*this%area(cc,lchnk)
388 : end do
389 : end do
390 :
391 : end do
392 : end do
393 :
394 0 : call shr_reprosum_calc(Csum, Cvec, count, nlcols, Cvec_len, gbl_count=ngcols_p, commid=mpicom)
395 :
396 0 : cnum = 0
397 0 : do nn= 1,this%nbas
398 0 : do n2=nn,this%nbas
399 0 : cnum = cnum + 1
400 0 : Bcov(nn,n2) = Cvec(cnum)
401 0 : Bcov(n2,nn) = Cvec(cnum)
402 : end do
403 : end do
404 :
405 : ! Invert to get the basis amplitude map
406 : !--------------------------------------
407 0 : call Invert_Matrix(Bcov,this%nbas,this%map)
408 :
409 : ! End Routine
410 : !------------
411 0 : deallocate(Clats)
412 0 : deallocate(Bcoef)
413 0 : deallocate(Csum )
414 0 : deallocate(Cvec )
415 0 : deallocate(Bsum )
416 0 : deallocate(Bnorm)
417 0 : deallocate(Bcov )
418 :
419 0 : end subroutine init_ZonalMean
420 : !=======================================================================
421 :
422 :
423 : !=======================================================================
424 0 : subroutine calc_ZonalMean_2Damps(this,I_Gdata,O_Bamp)
425 : !
426 : ! calc_ZonalMean_2Damps: Given 2D data values for the ncol gridpoints,
427 : ! compute the zonal mean basis amplitudes.
428 : !=====================================================================
429 : !
430 : ! Passed Variables
431 : !------------------
432 : class(ZonalMean_t) :: this
433 : real(r8),intent(in ) :: I_Gdata(pcols,begchunk:endchunk)
434 : real(r8),intent(out) :: O_Bamp(:)
435 : !
436 : ! Local Values
437 : !--------------
438 0 : real(r8),allocatable :: Csum(:,:)
439 0 : real(r8),allocatable :: Gcov(:)
440 : integer :: nn,n2,ncols,lchnk,cc
441 : integer :: nlcols, count, astat
442 :
443 : character(len=*), parameter :: subname = 'calc_ZonalMean_2Damps'
444 :
445 0 : nlcols = get_nlcols_p()
446 :
447 0 : allocate(Gcov(this%nbas), stat=astat)
448 0 : call handle_allocate_error(astat, subname, 'Gcov')
449 0 : allocate(Csum(nlcols, this%nbas), stat=astat)
450 0 : call handle_allocate_error(astat, subname, 'Csum')
451 0 : Csum(:,:) = 0._r8
452 :
453 : ! Compute Covariance with input data and basis functions
454 : !--------------------------------------------------------
455 0 : do nn= 1,this%nbas
456 0 : count = 0
457 0 : do lchnk=begchunk,endchunk
458 0 : ncols = get_ncols_p(lchnk)
459 0 : do cc = 1,ncols
460 0 : count=count+1
461 0 : Csum(count,nn) = I_Gdata(cc,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk)
462 : end do
463 : end do
464 : end do
465 :
466 0 : call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom)
467 :
468 : ! Multiply by map to get the amplitudes
469 : !-------------------------------------------
470 0 : do nn=1,this%nbas
471 0 : O_Bamp(nn) = 0._r8
472 0 : do n2=1,this%nbas
473 0 : O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2)
474 : end do
475 : end do
476 :
477 : ! End Routine
478 : !------------
479 0 : deallocate(Csum)
480 0 : deallocate(Gcov)
481 :
482 0 : end subroutine calc_ZonalMean_2Damps
483 : !=======================================================================
484 :
485 :
486 : !=======================================================================
487 0 : subroutine calc_ZonalMean_3Damps(this,I_Gdata,O_Bamp)
488 : !
489 : ! calc_ZonalMean_3Damps: Given 3D data values for the ncol,nlev gridpoints,
490 : ! compute the zonal mean basis amplitudes.
491 : !=====================================================================
492 : !
493 : ! Passed Variables
494 : !------------------
495 : class(ZonalMean_t) :: this
496 : real(r8),intent(in ):: I_Gdata(:,:,begchunk:)
497 : real(r8),intent(out):: O_Bamp (:,:)
498 : !
499 : ! Local Values
500 : !--------------
501 0 : real(r8),allocatable:: Csum (:,:)
502 0 : real(r8),allocatable:: Gcov (:)
503 : integer:: nn,n2,ncols,lchnk,cc
504 : integer:: Nsum,ns,ll
505 : integer :: nlcols, count, astat
506 :
507 : integer :: nlev
508 : character(len=*), parameter :: subname = 'calc_ZonalMean_3Damps'
509 :
510 0 : nlev = size(I_Gdata,dim=2)
511 :
512 0 : nlcols = get_nlcols_p()
513 0 : allocate(Gcov(this%nbas), stat=astat)
514 0 : call handle_allocate_error(astat, subname, 'Gcov')
515 0 : allocate(Csum(nlcols, this%nbas), stat=astat)
516 0 : call handle_allocate_error(astat, subname, 'Csum')
517 :
518 0 : Csum(:,:) = 0._r8
519 0 : O_Bamp(:,:) = 0._r8
520 :
521 : ! Compute Covariance with input data and basis functions
522 : !--------------------------------------------------------
523 0 : do ll= 1,nlev
524 :
525 0 : Csum(:,:) = 0._r8
526 0 : Gcov(:) = 0._r8
527 :
528 0 : do nn= 1,this%nbas
529 0 : count = 0
530 0 : do lchnk=begchunk,endchunk
531 0 : ncols = get_ncols_p(lchnk)
532 0 : do cc = 1,ncols
533 0 : count=count+1
534 0 : Csum(count,nn) = I_Gdata(cc,ll,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk)
535 : end do
536 : end do
537 : end do
538 :
539 0 : call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom)
540 :
541 : ! Multiply by map to get the amplitudes
542 : !-------------------------------------------
543 0 : do nn=1,this%nbas
544 0 : O_Bamp(nn,ll) = 0._r8
545 0 : do n2=1,this%nbas
546 0 : O_Bamp(nn,ll) = O_Bamp(nn,ll) + this%map(n2,nn)*Gcov(n2)
547 : end do
548 : end do
549 :
550 : end do
551 :
552 : ! End Routine
553 : !------------
554 0 : deallocate(Csum)
555 0 : deallocate(Gcov)
556 :
557 0 : end subroutine calc_ZonalMean_3Damps
558 : !=======================================================================
559 :
560 :
561 : !=======================================================================
562 0 : subroutine eval_ZonalMean_2Dgrid(this,I_Bamp,O_Gdata)
563 : !
564 : ! eval_ZonalMean_2Dgrid: Given the zonal mean basis amplitudes,
565 : ! compute 2D data values for the ncol gridpoints.
566 : !=====================================================================
567 : !
568 : ! Passed Variables
569 : !------------------
570 : class(ZonalMean_t) :: this
571 : real(r8),intent(in ):: I_Bamp (:)
572 : real(r8),intent(out):: O_Gdata(pcols,begchunk:endchunk)
573 : !
574 : ! Local Values
575 : !--------------
576 : integer:: nn,ncols,lchnk,cc
577 :
578 0 : O_Gdata(:,:) = 0._r8
579 :
580 : ! Construct grid values from basis amplitudes.
581 : !--------------------------------------------------
582 :
583 0 : do nn=1,this%nbas
584 0 : do lchnk=begchunk,endchunk
585 0 : ncols = get_ncols_p(lchnk)
586 0 : do cc = 1,ncols
587 0 : O_Gdata(cc,lchnk) = O_Gdata(cc,lchnk) + (I_Bamp(nn)*this%basis(cc,lchnk,nn))
588 : end do
589 : end do
590 : end do
591 :
592 0 : end subroutine eval_ZonalMean_2Dgrid
593 : !=======================================================================
594 :
595 :
596 : !=======================================================================
597 0 : subroutine eval_ZonalMean_3Dgrid(this,I_Bamp,O_Gdata)
598 : !
599 : ! eval_ZonalMean_3Dgrid: Given the zonal mean basis amplitudes,
600 : ! compute 3D data values for the ncol,nlev gridpoints.
601 : !=====================================================================
602 : !
603 : ! Passed Variables
604 : !------------------
605 : class(ZonalMean_t) :: this
606 : real(r8),intent(in ):: I_Bamp (:,:)
607 : real(r8),intent(out):: O_Gdata(:,:,begchunk:)
608 : !
609 : ! Local Values
610 : !--------------
611 : integer:: nn,ncols,lchnk,cc
612 : integer:: ll
613 :
614 : integer :: nlev
615 0 : nlev = size(O_Gdata,dim=2)
616 :
617 0 : O_Gdata(:,:,:) = 0._r8
618 :
619 : ! Construct grid values from basis amplitudes.
620 : !--------------------------------------------------
621 :
622 0 : do ll = 1,nlev
623 0 : do nn=1,this%nbas
624 0 : do lchnk=begchunk,endchunk
625 0 : ncols = get_ncols_p(lchnk)
626 0 : do cc = 1,ncols
627 0 : O_Gdata(cc,ll,lchnk) = O_Gdata(cc,ll,lchnk) + (I_Bamp(nn,ll)*this%basis(cc,lchnk,nn))
628 : end do
629 : end do
630 : end do
631 : end do
632 :
633 0 : end subroutine eval_ZonalMean_3Dgrid
634 : !=======================================================================
635 :
636 : !=======================================================================
637 1536 : subroutine final_ZonalMean(this)
638 : class(ZonalMean_t) :: this
639 :
640 1536 : if(allocated(this%area )) deallocate(this%area)
641 1536 : if(allocated(this%basis)) deallocate(this%basis)
642 1536 : if(allocated(this%map )) deallocate(this%map)
643 :
644 1536 : end subroutine final_ZonalMean
645 : !=======================================================================
646 :
647 : !=======================================================================
648 0 : subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS)
649 : !
650 : ! init_ZonalProfile: Initialize the ZonalProfile data structure for the
651 : ! given nlat gridpoints. It is assumed that the domain
652 : ! of these gridpoints of the profile span latitudes
653 : ! from SP to NP.
654 : ! The representation of basis functions functions is
655 : ! normalized w.r.t integration over the sphere so that
656 : ! when configured for tha same number of basis functions,
657 : ! the calculated amplitudes are interchangable with
658 : ! those for the ZonalMean_t class.
659 : !
660 : ! The optional GEN_GAUSSLATS flag allows for the
661 : ! generation of Gaussian latitudes. The generated grid
662 : ! over-writes the values of IO_lats/IO_area passed by
663 : ! the user.
664 : !=====================================================================
665 : !
666 : ! Passed Variables
667 : !------------------
668 : class(ZonalProfile_t) :: this
669 : real(r8) ,intent(inout):: IO_lats(:)
670 : real(r8) ,intent(inout):: IO_area(:)
671 : integer ,intent(in):: I_nlat
672 : integer ,intent(in):: I_nbas
673 : logical,optional,intent(in):: GEN_GAUSSLATS
674 : !
675 : ! Local Values
676 : !--------------
677 0 : real(r8),allocatable:: Clats(:)
678 0 : real(r8),allocatable:: Bcoef(:)
679 0 : real(r8),allocatable:: Bcov (:,:)
680 : real(r8):: Bnorm
681 : integer :: ii,nn,n2,nb,ierr, astat
682 : logical :: generate_lats
683 :
684 : character(len=*), parameter :: subname = 'init_ZonalProfile'
685 :
686 0 : generate_lats = .false.
687 :
688 0 : if (present(GEN_GAUSSLATS)) then
689 0 : generate_lats = GEN_GAUSSLATS
690 : end if
691 :
692 : ! Allocate space
693 : !-----------------
694 0 : if(allocated(this%area )) deallocate(this%area)
695 0 : if(allocated(this%basis)) deallocate(this%basis)
696 0 : if(allocated(this%map )) deallocate(this%map)
697 :
698 0 : this%nlat = I_nlat
699 0 : this%nbas = I_nbas
700 0 : allocate(this%area (I_nlat), stat=astat)
701 0 : call handle_allocate_error(astat, subname, 'this%area')
702 0 : allocate(this%basis(I_nlat,I_nbas), stat=astat)
703 0 : call handle_allocate_error(astat, subname, 'this%basis')
704 0 : allocate(this%map (I_nbas,I_nbas), stat=astat)
705 0 : call handle_allocate_error(astat, subname, 'this%map')
706 :
707 0 : allocate(Clats(I_nlat), stat=astat)
708 0 : call handle_allocate_error(astat, subname, 'Clats')
709 0 : allocate(Bcoef(I_nbas/2+1), stat=astat)
710 0 : call handle_allocate_error(astat, subname, 'Bcoef')
711 0 : allocate(Bcov (I_nbas,I_nbas), stat=astat)
712 0 : call handle_allocate_error(astat, subname, 'Bcov')
713 :
714 : ! Optionally create the Latitude Gridpoints
715 : ! and their associated area weights. Otherwise
716 : ! they need to be supplied by the user.
717 : !-----------------------------------------------
718 0 : if(generate_lats) then
719 :
720 : ! Create a Gaussian grid from SP to NP
721 : !--------------------------------------
722 0 : call sh_create_gaus_grid(I_nlat,Clats,IO_area,ierr)
723 0 : if (ierr/=0) then
724 0 : call endrun('init_ZonalProfile: Error creating Gaussian grid')
725 : end if
726 :
727 : ! Convert generated colatitudes SP->NP to Lats and convert
728 : ! to degrees and scale the area for global 2D integrals
729 : !-----------------------------------------------------------
730 0 : do nn=1,I_nlat
731 0 : IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8
732 0 : IO_area(nn) = IO_area(nn)*twoPI
733 : end do
734 : else
735 : ! Convert Latitudes to SP->NP colatitudes in radians
736 : !----------------------------------------------------
737 0 : do nn=1,I_nlat
738 0 : Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8
739 : end do
740 : endif
741 :
742 : ! Copy the area weights for each nlat
743 : ! gridpoint to the data structure
744 : !---------------------------------------
745 0 : this%area(1:I_nlat) = IO_area(1:I_nlat)
746 :
747 : ! Add first basis for the mean values.
748 : !------------------------------------------
749 0 : this%basis(:,1) = invSqrt4pi
750 : Bnorm = 0._r8
751 0 : do ii=1,I_nlat
752 0 : Bnorm = Bnorm + (this%basis(ii,1)*this%basis(ii,1)*this%area(ii))
753 : end do
754 0 : this%basis(:,1) = this%basis(:,1)/sqrt(Bnorm)
755 :
756 : ! Loop over the remaining basis functions
757 : !---------------------------------------
758 0 : do nn=2,I_nbas
759 0 : nb = nn-1
760 :
761 : ! Generate coefs for the basis
762 : !------------------------------
763 0 : call sh_gen_basis_coefs(nb,0,Bcoef)
764 :
765 : ! Create an un-normalized basis for the
766 : ! coefs at each nlat gridpoint
767 : !---------------------------------------
768 0 : do ii=1,I_nlat
769 0 : call sh_create_basis(nb,0,Clats(ii),Bcoef,this%basis(ii,nn))
770 : end do
771 :
772 : ! Numerically normalize the basis funnction
773 : !--------------------------------------------------------------
774 : Bnorm = 0._r8
775 0 : do ii=1,I_nlat
776 0 : Bnorm = Bnorm + (this%basis(ii,nn)*this%basis(ii,nn)*this%area(ii))
777 : end do
778 0 : this%basis(:,nn) = this%basis(:,nn)/sqrt(Bnorm)
779 :
780 : end do ! nn=1,I_nbas
781 :
782 : ! Compute covariance matrix for basis functions
783 : ! (Yes, they are theoretically orthonormal, but lets make sure)
784 : !--------------------------------------------------------------
785 0 : do nn=1,I_nbas
786 0 : do n2=1,I_nbas
787 0 : Bcov(nn,n2) = 0._r8
788 0 : do ii=1,I_nlat
789 0 : Bcov(nn,n2) = Bcov(nn,n2) + (this%basis(ii,nn)*this%basis(ii,n2)*this%area(ii))
790 : end do
791 : end do
792 : end do
793 :
794 : ! Invert to get the basis amplitude map
795 : !--------------------------------------
796 0 : call Invert_Matrix(Bcov,I_nbas,this%map)
797 :
798 : ! End Routine
799 : !------------
800 0 : deallocate(Clats)
801 0 : deallocate(Bcoef)
802 0 : deallocate(Bcov )
803 :
804 0 : end subroutine init_ZonalProfile
805 : !=======================================================================
806 :
807 :
808 : !=======================================================================
809 0 : subroutine calc_ZonalProfile_1Damps(this,I_Zdata,O_Bamp)
810 : !
811 : ! calc_ZonalProfile_1Damps: Given 1D data values for the nlat zonal
812 : ! profiles gridpoints, compute the zonal
813 : ! profile basis amplitudes.
814 : !=====================================================================
815 : !
816 : ! Passed Variables
817 : !------------------
818 : class(ZonalProfile_t):: this
819 : real(r8),intent(in ):: I_Zdata(:)
820 : real(r8),intent(out):: O_Bamp (:)
821 : !
822 : ! Local Values
823 : !--------------
824 0 : real(r8),allocatable:: Gcov(:)
825 : integer:: ii,nn,n2, astat
826 : character(len=*), parameter :: subname = 'calc_ZonalProfile_1Damps'
827 :
828 : ! Compute Covariance with input data and basis functions
829 : !--------------------------------------------------------
830 0 : allocate(Gcov(this%nbas), stat=astat)
831 0 : call handle_allocate_error(astat, subname, 'Gcov')
832 0 : do nn=1,this%nbas
833 0 : Gcov(nn) = 0._r8
834 0 : do ii=1,this%nlat
835 0 : Gcov(nn) = Gcov(nn) + (I_Zdata(ii)*this%basis(ii,nn)*this%area(ii))
836 : end do
837 : end do
838 :
839 : ! Multiply by map to get the amplitudes
840 : !-------------------------------------------
841 0 : do nn=1,this%nbas
842 0 : O_Bamp(nn) = 0._r8
843 0 : do n2=1,this%nbas
844 0 : O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2)
845 : end do
846 : end do
847 :
848 0 : deallocate(Gcov)
849 :
850 0 : end subroutine calc_ZonalProfile_1Damps
851 : !=======================================================================
852 :
853 :
854 : !=======================================================================
855 0 : subroutine calc_ZonalProfile_2Damps(this,I_Zdata,O_Bamp)
856 : !
857 : ! calc_ZonalProfile_2Damps: Given 2D data values for the nlat,nlev zonal
858 : ! profiles gridpoints, compute the zonal
859 : ! profile basis amplitudes.
860 : !=====================================================================
861 : !
862 : ! Passed Variables
863 : !------------------
864 : class(ZonalProfile_t):: this
865 : real(r8),intent(in ):: I_Zdata(:,:)
866 : real(r8),intent(out):: O_Bamp (:,:)
867 : !
868 : ! Local Values
869 : !--------------
870 0 : real(r8),allocatable:: Gcov(:,:)
871 : integer:: ii,nn,n2,ilev
872 :
873 : integer :: nlev, astat
874 : character(len=*), parameter :: subname = 'calc_ZonalProfile_2Damps'
875 :
876 0 : nlev = size(I_Zdata,dim=2)
877 :
878 : ! Compute Covariance with input data and basis functions
879 : !--------------------------------------------------------
880 0 : allocate(Gcov(this%nbas,nlev), stat=astat)
881 0 : call handle_allocate_error(astat, subname, 'Gcov')
882 0 : do ilev=1,nlev
883 0 : do nn=1,this%nbas
884 0 : Gcov(nn,ilev) = 0._r8
885 0 : do ii=1,this%nlat
886 0 : Gcov(nn,ilev) = Gcov(nn,ilev) + (I_Zdata(ii,ilev)*this%basis(ii,nn)*this%area(ii))
887 : end do
888 : end do
889 : end do
890 :
891 : ! Multiply by map to get the amplitudes
892 : !-------------------------------------------
893 0 : do ilev=1,nlev
894 0 : do nn=1,this%nbas
895 0 : O_Bamp(nn,ilev) = 0._r8
896 0 : do n2=1,this%nbas
897 0 : O_Bamp(nn,ilev) = O_Bamp(nn,ilev) + this%map(n2,nn)*Gcov(n2,ilev)
898 : end do
899 : end do
900 : end do
901 0 : deallocate(Gcov)
902 :
903 0 : end subroutine calc_ZonalProfile_2Damps
904 : !=======================================================================
905 :
906 :
907 : !=======================================================================
908 0 : subroutine eval_ZonalProfile_1Dgrid(this,I_Bamp,O_Zdata)
909 : !
910 : ! eval_ZonalProfile_1Dgrid: Given the zonal profile basis amplitudes,
911 : ! compute 1D data values for the nlat gridpoints.
912 : !=====================================================================
913 : !
914 : ! Passed Variables
915 : !------------------
916 : class(ZonalProfile_t):: this
917 : real(r8),intent(in ):: I_Bamp (:)
918 : real(r8),intent(out):: O_Zdata(:)
919 : !
920 : ! Local Values
921 : !--------------
922 : integer:: ii,nn
923 :
924 : ! Construct grid values from basis amplitudes.
925 : !--------------------------------------------------
926 0 : O_Zdata(1:this%nlat) = 0._r8
927 0 : do nn=1,this%nbas
928 0 : do ii=1,this%nlat
929 0 : O_Zdata(ii) = O_Zdata(ii) + (I_Bamp(nn)*this%basis(ii,nn))
930 : end do
931 : end do
932 :
933 0 : end subroutine eval_ZonalProfile_1Dgrid
934 : !=======================================================================
935 :
936 :
937 : !=======================================================================
938 0 : subroutine eval_ZonalProfile_2Dgrid(this,I_Bamp,O_Zdata)
939 : !
940 : ! eval_ZonalProfile_2Dgrid: Given the zonal profile basis amplitudes,
941 : ! compute 2D data values for the nlat,nlev gridpoints.
942 : !=====================================================================
943 : !
944 : ! Passed Variables
945 : !------------------
946 : class(ZonalProfile_t):: this
947 : real(r8),intent(in ):: I_Bamp (:,:)
948 : real(r8),intent(out):: O_Zdata(:,:)
949 : !
950 : ! Local Values
951 : !--------------
952 : integer:: ii,nn,ilev
953 :
954 : integer :: nlev
955 :
956 0 : nlev = size(I_Bamp,dim=2)
957 :
958 : ! Construct grid values from basis amplitudes.
959 : !--------------------------------------------------
960 0 : O_Zdata(1:this%nlat,1:nlev) = 0._r8
961 0 : do nn=1,this%nbas
962 0 : do ilev=1,nlev
963 0 : do ii=1,this%nlat
964 0 : O_Zdata(ii,ilev) = O_Zdata(ii,ilev) + (I_Bamp(nn,ilev)*this%basis(ii,nn))
965 : end do
966 : end do
967 : end do
968 :
969 0 : end subroutine eval_ZonalProfile_2Dgrid
970 : !=======================================================================
971 :
972 : !=======================================================================
973 0 : subroutine final_ZonalProfile(this)
974 : class(ZonalProfile_t) :: this
975 :
976 0 : if(allocated(this%area )) deallocate(this%area)
977 0 : if(allocated(this%basis)) deallocate(this%basis)
978 0 : if(allocated(this%map )) deallocate(this%map)
979 :
980 0 : end subroutine final_ZonalProfile
981 : !=======================================================================
982 :
983 : !=======================================================================
984 0 : subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS)
985 : !
986 : ! init_ZonalAverage: Initialize the ZonalAverage data structure for the
987 : ! given nlat gridpoints. It is assumed that the domain
988 : ! of these gridpoints of the profile span latitudes
989 : ! from SP to NP.
990 : !
991 : ! The optional GEN_GAUSSLATS flag allows for the
992 : ! generation of Gaussian latitudes. The generated grid
993 : ! over-writes the values of IO_lats/IO_area passed by
994 : ! the user.
995 : !=====================================================================
996 : !
997 : ! Passed Variables
998 : !------------------
999 : class(ZonalAverage_t) :: this
1000 : real(r8) ,intent(inout):: IO_lats(:)
1001 : real(r8) ,intent(inout):: IO_area(:)
1002 : integer ,intent(in):: I_nlat
1003 : logical,optional,intent(in):: GEN_GAUSSLATS
1004 : !
1005 : ! Local Values
1006 : !--------------
1007 0 : real(r8),allocatable:: Clats (:)
1008 0 : real(r8),allocatable:: Glats (:,:)
1009 0 : real(r8),allocatable:: BinLat(:)
1010 0 : real(r8),allocatable:: Asum (:,:)
1011 0 : real(r8),allocatable:: Anorm (:)
1012 : real(r8):: area(pcols),rlat
1013 : integer :: nn,jj,ierr, astat
1014 : integer :: ncols,lchnk,cc,jlat
1015 : integer :: nlcols, count
1016 : logical :: generate_lats
1017 : character(len=*), parameter :: subname = 'init_ZonalAverage'
1018 :
1019 0 : generate_lats = .false.
1020 :
1021 0 : if (present(GEN_GAUSSLATS)) then
1022 0 : generate_lats = GEN_GAUSSLATS
1023 : end if
1024 :
1025 0 : nlcols = get_nlcols_p()
1026 :
1027 : ! Allocate space
1028 : !-----------------
1029 0 : if(allocated(this%area )) deallocate(this%area)
1030 0 : if(allocated(this%a_norm )) deallocate(this%a_norm)
1031 0 : if(allocated(this%area_g )) deallocate(this%area_g)
1032 0 : if(allocated(this%idx_map)) deallocate(this%idx_map)
1033 :
1034 0 : this%nlat = I_nlat
1035 0 : allocate(this%area (I_nlat), stat=astat)
1036 0 : call handle_allocate_error(astat, subname, 'this%area')
1037 0 : allocate(this%a_norm (I_nlat), stat=astat)
1038 0 : call handle_allocate_error(astat, subname, 'this%a_norm')
1039 0 : allocate(this%area_g (pcols,begchunk:endchunk), stat=astat)
1040 0 : call handle_allocate_error(astat, subname, 'this%area_g')
1041 0 : allocate(this%idx_map(pcols,begchunk:endchunk), stat=astat)
1042 0 : call handle_allocate_error(astat, subname, 'this%idx_map')
1043 :
1044 0 : allocate(Clats (I_nlat), stat=astat)
1045 0 : call handle_allocate_error(astat, subname, 'Clats')
1046 0 : allocate(BinLat(I_nlat+1), stat=astat)
1047 0 : call handle_allocate_error(astat, subname, 'BinLat')
1048 0 : allocate(Glats (pcols,begchunk:endchunk), stat=astat)
1049 0 : call handle_allocate_error(astat, subname, 'Glats')
1050 0 : allocate(Asum (nlcols,I_nlat), stat=astat)
1051 0 : call handle_allocate_error(astat, subname, 'Asum')
1052 0 : allocate(Anorm (I_nlat), stat=astat)
1053 0 : call handle_allocate_error(astat, subname, 'Anorm')
1054 :
1055 : ! Optionally create the Latitude Gridpoints
1056 : ! and their associated area weights. Otherwise
1057 : ! they need to be supplied by the user.
1058 : !-----------------------------------------------
1059 0 : if(generate_lats) then
1060 :
1061 : ! Create a Gaussin grid from SP to NP
1062 : !--------------------------------------
1063 0 : call sh_create_gaus_grid(this%nlat,Clats,IO_area,ierr)
1064 0 : if (ierr/=0) then
1065 0 : call endrun('init_ZonalAverage: Error creating Gaussian grid')
1066 : end if
1067 :
1068 : ! Convert generated colatitudes SP->NP to Lats and convert
1069 : ! to degrees and scale the area for global 2D integrals
1070 : !-----------------------------------------------------------
1071 0 : do nn=1,this%nlat
1072 0 : IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8
1073 0 : IO_area(nn) = IO_area(nn)*twoPI
1074 : end do
1075 : else
1076 : ! Convert Latitudes to SP->NP colatitudes in radians
1077 : !----------------------------------------------------
1078 0 : do nn=1,this%nlat
1079 0 : Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8
1080 : end do
1081 : endif
1082 :
1083 : ! Copy the Lat grid area weights to the data structure
1084 : !-----------------------------------------------------
1085 0 : this%area(1:this%nlat) = IO_area(1:this%nlat)
1086 :
1087 : ! Save a copy of the area weights for each 2D gridpoint
1088 : ! and convert Latitudes to SP->NP colatitudes in radians
1089 : !-------------------------------------------------------
1090 0 : do lchnk=begchunk,endchunk
1091 0 : ncols = get_ncols_p(lchnk)
1092 0 : call get_wght_all_p(lchnk, ncols, area)
1093 0 : do cc = 1,ncols
1094 0 : rlat=get_rlat_p(lchnk,cc)
1095 0 : this%area_g(cc,lchnk) = area(cc)
1096 0 : Glats (cc,lchnk) = rlat + halfPI
1097 : end do
1098 : end do
1099 :
1100 : ! Set boundaries for Latitude bins
1101 : !-----------------------------------
1102 0 : BinLat(1) = 0._r8
1103 0 : BinLat(this%nlat+1) = pi
1104 0 : do nn=2,this%nlat
1105 0 : BinLat(nn) = (Clats(nn-1)+Clats(nn))/2._r8
1106 : end do
1107 :
1108 : ! Loop over 2D gridpoints and determine its lat bin index
1109 : !---------------------------------------------------------
1110 0 : do lchnk=begchunk,endchunk
1111 0 : ncols = get_ncols_p(lchnk)
1112 0 : do cc = 1,ncols
1113 0 : jlat = -1
1114 0 : if((Glats(cc,lchnk)<=BinLat(2)).and. &
1115 : (Glats(cc,lchnk)>=BinLat(1)) ) then
1116 0 : jlat = 1
1117 0 : elseif((Glats(cc,lchnk)>=BinLat(this%nlat) ).and. &
1118 0 : (Glats(cc,lchnk)<=BinLat(this%nlat+1)) ) then
1119 0 : jlat = this%nlat
1120 : else
1121 0 : do jj=2,(this%nlat-1)
1122 0 : if((Glats(cc,lchnk)>BinLat(jj )).and. &
1123 0 : (Glats(cc,lchnk)<=BinLat(jj+1)) ) then
1124 0 : jlat = jj
1125 0 : exit
1126 : endif
1127 : end do
1128 : endif
1129 0 : if (jlat<1) then
1130 0 : call endrun('ZonalAverage init ERROR: jlat not in range')
1131 : endif
1132 0 : this%idx_map(cc,lchnk) = jlat
1133 : end do
1134 : end do
1135 :
1136 : ! Initialize 2D Area sums for each bin
1137 : !--------------------------------------
1138 0 : Asum(:,:) = 0._r8
1139 0 : Anorm(:) = 0._r8
1140 0 : count = 0
1141 0 : do lchnk=begchunk,endchunk
1142 0 : ncols = get_ncols_p(lchnk)
1143 0 : do cc = 1,ncols
1144 0 : jlat = this%idx_map(cc,lchnk)
1145 0 : count=count+1
1146 0 : Asum(count,jlat) = this%area_g(cc,lchnk)
1147 : end do
1148 : end do
1149 :
1150 0 : call shr_reprosum_calc(Asum, Anorm, count, nlcols, I_nlat, gbl_count=ngcols_p, commid=mpicom)
1151 :
1152 0 : this%a_norm = Anorm
1153 :
1154 0 : if (.not.all(Anorm(:)>0._r8)) then
1155 0 : write(iulog,*) 'init_ZonalAverage -- ERROR in Anorm values: '
1156 0 : do jlat = 1,I_nlat
1157 0 : if (.not.Anorm(jlat)>0._r8) then
1158 0 : write(iulog,*) ' Anorm(',jlat,'): ', Anorm(jlat)
1159 : endif
1160 : end do
1161 0 : call endrun('init_ZonalAverage -- ERROR in Anorm values')
1162 : end if
1163 :
1164 : ! End Routine
1165 : !------------
1166 0 : deallocate(Clats)
1167 0 : deallocate(BinLat)
1168 0 : deallocate(Glats)
1169 0 : deallocate(Asum)
1170 0 : deallocate(Anorm)
1171 :
1172 0 : end subroutine init_ZonalAverage
1173 : !=======================================================================
1174 :
1175 :
1176 : !=======================================================================
1177 0 : subroutine calc_ZonalAverage_2DbinAvg(this,I_Gdata,O_Zdata)
1178 : !
1179 : ! calc_ZonalAverage_2DbinAvg: Given 2D data values for ncol gridpoints,
1180 : ! compute the nlat area weighted binAvg profile
1181 : !=====================================================================
1182 : !
1183 : ! Passed Variables
1184 : !------------------
1185 : class(ZonalAverage_t):: this
1186 : real(r8),intent(in ):: I_Gdata(pcols,begchunk:endchunk)
1187 : real(r8),intent(out):: O_Zdata(:)
1188 : !
1189 : ! Local Values
1190 : !--------------
1191 0 : real(r8),allocatable:: Asum (:,:)
1192 : integer:: nn,ncols,lchnk,cc,jlat
1193 : integer :: nlcols, count, astat
1194 : character(len=*), parameter :: subname = 'calc_ZonalAverage_2DbinAvg'
1195 :
1196 0 : nlcols = get_nlcols_p()
1197 :
1198 :
1199 : ! Initialize Zonal profile
1200 : !---------------------------
1201 0 : allocate(Asum(nlcols,this%nlat), stat=astat)
1202 0 : call handle_allocate_error(astat, subname, 'Asum')
1203 0 : Asum(:,:) = 0._r8
1204 :
1205 0 : O_Zdata(1:this%nlat) = 0._r8
1206 :
1207 : ! Compute area-weighted sums
1208 : !-----------------------------
1209 0 : count = 0
1210 0 : do lchnk=begchunk,endchunk
1211 0 : ncols = get_ncols_p(lchnk)
1212 0 : do cc = 1,ncols
1213 0 : jlat = this%idx_map(cc,lchnk)
1214 0 : count=count+1
1215 0 : Asum(count,jlat) = I_Gdata(cc,lchnk)*this%area_g(cc,lchnk)
1216 : end do
1217 : end do
1218 :
1219 0 : call shr_reprosum_calc(Asum,O_Zdata,count, nlcols, this%nlat,gbl_count=ngcols_p, commid=mpicom)
1220 :
1221 : ! Divide by area norm to get the averages
1222 : !-----------------------------------------
1223 0 : do nn=1,this%nlat
1224 0 : O_Zdata(nn) = O_Zdata(nn)/this%a_norm(nn)
1225 : end do
1226 :
1227 0 : deallocate(Asum)
1228 :
1229 0 : end subroutine calc_ZonalAverage_2DbinAvg
1230 : !=======================================================================
1231 :
1232 :
1233 : !=======================================================================
1234 0 : subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata)
1235 : !
1236 : ! calc_ZonalAverage_3DbinAvg: Given 3D data values for ncol,nlev gridpoints,
1237 : ! compute the nlat,nlev area weighted binAvg profile
1238 : !=====================================================================
1239 : !
1240 : ! Passed Variables
1241 : !------------------
1242 : class(ZonalAverage_t):: this
1243 : real(r8),intent(in ):: I_Gdata(:,:,begchunk:)
1244 : real(r8),intent(out):: O_Zdata(:,:)
1245 : !
1246 : ! Local Values
1247 : !--------------
1248 0 : real(r8),allocatable:: Gsum(:)
1249 0 : real(r8),allocatable:: Asum(:,:)
1250 : integer:: nn,ncols,lchnk,cc,jlat
1251 : integer:: Nsum,ilev,ns
1252 :
1253 : integer :: nlev
1254 : integer :: nlcols, count, astat
1255 : character(len=*), parameter :: subname = 'calc_ZonalAverage_3DbinAvg'
1256 :
1257 0 : nlev = size(I_Gdata,dim=2)
1258 0 : nlcols = get_nlcols_p()
1259 :
1260 : ! Initialize Zonal profile
1261 : !---------------------------
1262 0 : Nsum = this%nlat*nlev
1263 0 : allocate(Gsum(Nsum), stat=astat)
1264 0 : call handle_allocate_error(astat, subname, 'Gsum')
1265 0 : allocate(Asum(nlcols,Nsum), stat=astat)
1266 0 : call handle_allocate_error(astat, subname, 'Asum')
1267 0 : Asum(:,:) = 0._r8
1268 :
1269 0 : O_Zdata(1:this%nlat,1:nlev) = 0._r8
1270 :
1271 : ! Compute area-weighted sums
1272 : !-----------------------------
1273 0 : do ilev = 1,nlev
1274 0 : count = 0
1275 0 : do lchnk=begchunk,endchunk
1276 0 : ncols = get_ncols_p(lchnk)
1277 0 : do cc = 1,ncols
1278 0 : jlat = this%idx_map(cc,lchnk)
1279 0 : ns = jlat + (ilev-1)*this%nlat
1280 0 : count=count+1
1281 0 : Asum(count,ns) = I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk)
1282 : end do
1283 : end do
1284 : end do
1285 :
1286 0 : call shr_reprosum_calc(Asum,Gsum, count, nlcols, Nsum, gbl_count=ngcols_p, commid=mpicom)
1287 :
1288 : ! Divide by area norm to get the averages
1289 : !-----------------------------------------
1290 0 : do ilev = 1,nlev
1291 0 : do nn = 1,this%nlat
1292 0 : ns = nn + (ilev-1)*this%nlat
1293 0 : O_Zdata(nn,ilev) = Gsum(ns)/this%a_norm(nn)
1294 : end do
1295 : end do
1296 :
1297 0 : deallocate(Gsum)
1298 0 : deallocate(Asum)
1299 :
1300 0 : end subroutine calc_ZonalAverage_3DbinAvg
1301 : !=======================================================================
1302 :
1303 : !=======================================================================
1304 1536 : subroutine final_ZonalAverage(this)
1305 : class(ZonalAverage_t) :: this
1306 :
1307 1536 : if(allocated(this%area )) deallocate(this%area)
1308 1536 : if(allocated(this%a_norm )) deallocate(this%a_norm)
1309 1536 : if(allocated(this%area_g )) deallocate(this%area_g)
1310 1536 : if(allocated(this%idx_map)) deallocate(this%idx_map)
1311 :
1312 1536 : end subroutine final_ZonalAverage
1313 : !=======================================================================
1314 :
1315 :
1316 : !=======================================================================
1317 0 : subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat)
1318 : !
1319 : ! Invert_Matrix: Given the NbasxNbas matrix, calculate and return
1320 : ! the inverse of the matrix.
1321 : !
1322 : ! Implemented with the LAPACK DGESV routine.
1323 : !
1324 : !====================================================================
1325 : !
1326 : ! Passed Variables
1327 : !------------------
1328 : real(r8), intent(inout) :: I_Mat(:,:) ! input matrix contains P*L*U
1329 : ! decomposition on output
1330 : integer, intent(in) :: Nbas
1331 : real(r8), intent(out) :: O_InvMat(:,:)
1332 : !
1333 : ! Local Values
1334 : !-------------
1335 0 : integer, allocatable :: Indx(:) ! pivot indices
1336 : integer :: astat, ii
1337 : character(len=*), parameter :: subname = 'Invert_Matrix'
1338 : character(len=80) :: msg
1339 :
1340 : external DGESV
1341 :
1342 : ! Allocate work space
1343 : !---------------------
1344 0 : allocate(Indx(Nbas), stat=astat)
1345 0 : call handle_allocate_error(astat, subname, 'Indx')
1346 :
1347 : ! Initialize the inverse array with the identity matrix
1348 : !-------------------------------------------------------
1349 0 : O_InvMat(:,:) = 0._r8
1350 0 : do ii=1,Nbas
1351 0 : O_InvMat(ii,ii) = 1._r8
1352 : end do
1353 :
1354 0 : call DGESV(Nbas, Nbas, I_Mat, Nbas, Indx, O_InvMat, Nbas, astat)
1355 :
1356 0 : if (astat < 0) then
1357 0 : write(msg, '(a, i1, a)') 'argument # ', abs(astat), ' has an illegal value'
1358 0 : call endrun(subname//': DGESV error return: '//msg)
1359 0 : else if (astat > 0) then
1360 0 : call endrun(subname//': DGESV error return: matrix is singular')
1361 : end if
1362 :
1363 0 : deallocate(Indx)
1364 :
1365 0 : end subroutine Invert_Matrix
1366 : !=======================================================================
1367 :
1368 : !=======================================================================
1369 : ! legacy spherepack routines
1370 : !=======================================================================
1371 0 : subroutine sh_gen_basis_coefs(nn,mm,cp)
1372 : !
1373 : ! spherepack alfk
1374 : !
1375 : ! dimension of real cp(nn/2 + 1)
1376 : ! arguments
1377 : !
1378 : ! purpose computes fourier coefficients in the trigonometric series
1379 : ! representation of the normalized associated
1380 : ! legendre function pbar(nn,mm,theta) for use by
1381 : ! sh_gen_basis_coefs in calculating pbar(nn,mm,theta).
1382 : !
1383 : ! first define the normalized associated
1384 : ! legendre functions
1385 : !
1386 : ! pbar(mm,nn,theta) = sqrt((2*nn+1)*factorial(nn-mm)
1387 : ! /(2*factorial(nn+mm)))*sin(theta)**mm/(2**nn*
1388 : ! factorial(nn)) times the (nn+mm)th derivative of
1389 : ! (x**2-1)**nn with respect to x=cos(theta)
1390 : !
1391 : ! where theta is colatitude.
1392 : !
1393 : ! then subroutine sh_gen_basis_coefs computes the coefficients
1394 : ! cp(k) in the following trigonometric
1395 : ! expansion of pbar(m,n,theta).
1396 : !
1397 : ! 1) for n even and m even, pbar(mm,nn,theta) =
1398 : ! .5*cp(1) plus the sum from k=1 to k=nn/2
1399 : ! of cp(k+1)*cos(2*k*th)
1400 : !
1401 : ! 2) for nn even and mm odd, pbar(mm,nn,theta) =
1402 : ! the sum from k=1 to k=nn/2 of
1403 : ! cp(k)*sin(2*k*th)
1404 : !
1405 : ! 3) for n odd and m even, pbar(mm,nn,theta) =
1406 : ! the sum from k=1 to k=(nn+1)/2 of
1407 : ! cp(k)*cos((2*k-1)*th)
1408 : !
1409 : ! 4) for nn odd and mm odd, pbar(mm,nn,theta) =
1410 : ! the sum from k=1 to k=(nn+1)/2 of
1411 : ! cp(k)*sin((2*k-1)*th)
1412 : !
1413 : ! arguments
1414 : !
1415 : ! on input nn
1416 : ! nonnegative integer specifying the degree of
1417 : ! pbar(nn,mm,theta)
1418 : !
1419 : ! mm
1420 : ! is the order of pbar(nn,mm,theta). mm can be
1421 : ! any integer however cp is computed such that
1422 : ! pbar(nn,mm,theta) = 0 if abs(m) is greater
1423 : ! than nn and pbar(nn,mm,theta) = (-1)**mm*
1424 : ! pbar(nn,-mm,theta) for negative mm.
1425 : !
1426 : ! on output cp
1427 : ! array of length (nn/2)+1
1428 : ! which contains the fourier coefficients in
1429 : ! the trigonometric series representation of
1430 : ! pbar(nn,mm,theta)
1431 : !
1432 : ! special conditions none
1433 : !
1434 : ! algorithm the highest order coefficient is determined in
1435 : ! closed form and the remainig coefficients are
1436 : ! determined as the solution of a backward
1437 : ! recurrence relation.
1438 : !
1439 : !=====================================================================
1440 : !
1441 : ! Passed Variables
1442 : !------------------
1443 : integer ,intent(in ):: nn
1444 : integer ,intent(in ):: mm
1445 : real(r8),intent(out):: cp(nn/2+1)
1446 : !
1447 : ! Local Values
1448 : !----------------
1449 : real(r8):: fnum,fnmh
1450 : real(r8):: pm1
1451 : real(r8):: t1,t2
1452 : real(r8):: fden
1453 : real(r8):: cp2
1454 : real(r8):: fnnp1
1455 : real(r8):: fnmsq
1456 : real(r8):: fk
1457 : real(r8):: a1,b1,C1
1458 : integer :: ma,nmms2,nex
1459 : integer :: ii,jj
1460 :
1461 : real(r8),parameter:: SC10=1024._r8
1462 : real(r8),parameter:: SC20=SC10*SC10
1463 : real(r8),parameter:: SC40=SC20*SC20
1464 :
1465 0 : cp(1) = 0._r8
1466 0 : ma = abs(mm)
1467 0 : if(ma>nn) return
1468 :
1469 0 : if((nn-1)<0) then
1470 0 : cp(1) = sqrt(2._r8)
1471 0 : return
1472 0 : elseif((nn-1)==0) then
1473 0 : if(ma/=0) then
1474 0 : cp(1) = sqrt(.75_r8)
1475 0 : if(mm==-1) cp(1) = -cp(1)
1476 : else
1477 0 : cp(1) = sqrt(1.5_r8)
1478 : endif
1479 0 : return
1480 : else
1481 0 : if(mod(nn+ma,2)/=0) then
1482 0 : nmms2 = (nn-ma-1)/2
1483 0 : fnum = nn + ma + 2
1484 0 : fnmh = nn - ma + 2
1485 0 : pm1 = -1._r8
1486 : else
1487 0 : nmms2 = (nn-ma)/2
1488 0 : fnum = nn + ma + 1
1489 0 : fnmh = nn - ma + 1
1490 0 : pm1 = 1._r8
1491 : endif
1492 : endif
1493 :
1494 0 : t1 = 1._r8/SC20
1495 0 : nex = 20
1496 0 : fden = 2._r8
1497 0 : if(nmms2>=1) then
1498 0 : do ii = 1,nmms2
1499 0 : t1 = fnum*t1/fden
1500 0 : if (t1>SC20) then
1501 0 : t1 = t1/SC40
1502 0 : nex = nex + 40
1503 : endif
1504 0 : fnum = fnum + 2._r8
1505 0 : fden = fden + 2._r8
1506 : end do
1507 : endif
1508 :
1509 0 : if(mod(ma/2,2)/=0) then
1510 0 : t1 = -t1/2._r8**(nn-1-nex)
1511 : else
1512 0 : t1 = t1/2._r8**(nn-1-nex)
1513 : endif
1514 0 : t2 = 1._r8
1515 0 : if(ma/=0) then
1516 0 : do ii = 1,ma
1517 0 : t2 = fnmh*t2/ (fnmh+pm1)
1518 0 : fnmh = fnmh + 2._r8
1519 : end do
1520 : endif
1521 :
1522 0 : cp2 = t1*sqrt((nn+.5_r8)*t2)
1523 0 : fnnp1 = nn*(nn+1)
1524 0 : fnmsq = fnnp1 - 2._r8*ma*ma
1525 :
1526 0 : if((mod(nn,2)==0).and.(mod(ma,2)==0)) then
1527 0 : jj = 1+(nn+1)/2
1528 : else
1529 0 : jj = (nn+1)/2
1530 : endif
1531 :
1532 0 : cp(jj) = cp2
1533 0 : if(mm<0) then
1534 0 : if(mod(ma,2)/=0) cp(jj) = -cp(jj)
1535 : endif
1536 0 : if(jj<=1) return
1537 :
1538 0 : fk = nn
1539 0 : a1 = (fk-2._r8)*(fk-1._r8) - fnnp1
1540 0 : b1 = 2._r8* (fk*fk-fnmsq)
1541 0 : cp(jj-1) = b1*cp(jj)/a1
1542 :
1543 0 : jj = jj - 1
1544 0 : do while(jj>1)
1545 0 : fk = fk - 2._r8
1546 0 : a1 = (fk-2._r8)*(fk-1._r8) - fnnp1
1547 0 : b1 = -2._r8*(fk*fk-fnmsq)
1548 0 : c1 = (fk+1._r8)*(fk+2._r8) - fnnp1
1549 0 : cp(jj-1) = -(b1*cp(jj)+c1*cp(jj+1))/a1
1550 0 : jj = jj - 1
1551 : end do
1552 :
1553 : end subroutine sh_gen_basis_coefs
1554 : !=======================================================================
1555 :
1556 : !=======================================================================
1557 0 : subroutine sh_create_basis(nn,mm,theta,cp,pb)
1558 : !
1559 : ! spherepack lfpt
1560 : !
1561 : ! dimension of
1562 : ! arguments
1563 : ! cp((nn/2)+1)
1564 : !
1565 : ! purpose routine sh_create_basis uses coefficients computed by
1566 : ! routine sh_gen_basis_coefs to compute the
1567 : ! normalized associated legendre function pbar(nn,mm,theta)
1568 : ! at colatitude theta.
1569 : !
1570 : ! arguments
1571 : !
1572 : ! on input nn
1573 : ! nonnegative integer specifying the degree of
1574 : ! pbar(nn,mm,theta)
1575 : ! mm
1576 : ! is the order of pbar(nn,mm,theta). mm can be
1577 : ! any integer however pbar(nn,mm,theta) = 0
1578 : ! if abs(mm) is greater than nn and
1579 : ! pbar(nn,mm,theta) = (-1)**mm*pbar(nn,-mm,theta)
1580 : ! for negative mm.
1581 : !
1582 : ! theta
1583 : ! colatitude in radians
1584 : !
1585 : ! cp
1586 : ! array of length (nn/2)+1
1587 : ! containing coefficients computed by routine
1588 : ! sh_gen_basis_coefs
1589 : !
1590 : ! on output pb
1591 : ! variable containing pbar(n,m,theta)
1592 : !
1593 : ! special conditions calls to routine sh_create_basis must be preceded by an
1594 : ! appropriate call to routine sh_gen_basis_coefs.
1595 : !
1596 : ! algorithm the trigonometric series formula used by
1597 : ! routine sh_create_basis to calculate pbar(nn,mm,theta) at
1598 : ! colatitude theta depends on mm and nn as follows:
1599 : !
1600 : ! 1) for nn even and mm even, the formula is
1601 : ! .5*cp(1) plus the sum from k=1 to k=n/2
1602 : ! of cp(k)*cos(2*k*theta)
1603 : ! 2) for nn even and mm odd. the formula is
1604 : ! the sum from k=1 to k=nn/2 of
1605 : ! cp(k)*sin(2*k*theta)
1606 : ! 3) for nn odd and mm even, the formula is
1607 : ! the sum from k=1 to k=(nn+1)/2 of
1608 : ! cp(k)*cos((2*k-1)*theta)
1609 : ! 4) for nn odd and mm odd, the formula is
1610 : ! the sum from k=1 to k=(nn+1)/2 of
1611 : ! cp(k)*sin((2*k-1)*theta)
1612 : !
1613 : !=====================================================================
1614 : integer, intent(in) :: nn,mm
1615 : real(r8), intent(in) :: theta
1616 : real(r8), intent(in) :: cp(:)
1617 : real(r8), intent(out) :: pb
1618 :
1619 : real(r8) :: cdt
1620 : real(r8) :: sdt
1621 : real(r8) :: ct
1622 : real(r8) :: st
1623 : real(r8) :: summ
1624 : real(r8) :: cth
1625 :
1626 : integer:: ma,nmod,mmod,kdo
1627 : integer:: kp1,kk
1628 :
1629 0 : pb = 0._r8
1630 0 : ma = abs(mm)
1631 0 : if(ma>nn) return
1632 :
1633 0 : if(nn<=0) then
1634 0 : if(ma<=0) then
1635 0 : pb = sqrt(.5_r8)
1636 0 : return
1637 : endif
1638 : endif
1639 :
1640 0 : nmod = mod(nn,2)
1641 0 : mmod = mod(ma,2)
1642 :
1643 0 : if(nmod<=0) then
1644 0 : if(mmod<=0) then
1645 0 : kdo = nn/2 + 1
1646 0 : cdt = cos(theta+theta)
1647 0 : sdt = sin(theta+theta)
1648 0 : ct = 1._r8
1649 0 : st = 0._r8
1650 0 : summ = .5_r8*cp(1)
1651 0 : do kp1 = 2,kdo
1652 0 : cth = cdt*ct - sdt*st
1653 0 : st = sdt*ct + cdt*st
1654 0 : ct = cth
1655 0 : summ = summ + cp(kp1)*ct
1656 : end do
1657 0 : pb = summ
1658 0 : return
1659 : endif
1660 0 : kdo = nn/2
1661 0 : cdt = cos(theta+theta)
1662 0 : sdt = sin(theta+theta)
1663 0 : ct = 1._r8
1664 0 : st = 0._r8
1665 0 : summ = 0._r8
1666 0 : do kk = 1,kdo
1667 0 : cth = cdt*ct - sdt*st
1668 0 : st = sdt*ct + cdt*st
1669 0 : ct = cth
1670 0 : summ = summ + cp(kk)*st
1671 : end do
1672 0 : pb = summ
1673 0 : return
1674 : endif
1675 :
1676 0 : kdo = (nn+1)/2
1677 0 : if(mmod<=0) then
1678 0 : cdt = cos(theta+theta)
1679 0 : sdt = sin(theta+theta)
1680 0 : ct = cos(theta)
1681 0 : st = -sin(theta)
1682 0 : summ = 0._r8
1683 0 : do kk = 1,kdo
1684 0 : cth = cdt*ct - sdt*st
1685 0 : st = sdt*ct + cdt*st
1686 0 : ct = cth
1687 0 : summ = summ + cp(kk)*ct
1688 : end do
1689 0 : pb = summ
1690 0 : return
1691 : endif
1692 :
1693 0 : cdt = cos(theta+theta)
1694 0 : sdt = sin(theta+theta)
1695 0 : ct = cos(theta)
1696 0 : st = -sin(theta)
1697 0 : summ = 0._r8
1698 0 : do kk = 1,kdo
1699 0 : cth = cdt*ct - sdt*st
1700 0 : st = sdt*ct + cdt*st
1701 0 : ct = cth
1702 0 : summ = summ + cp(kk)*st
1703 : end do
1704 0 : pb = summ
1705 :
1706 : end subroutine sh_create_basis
1707 : !=======================================================================
1708 :
1709 : !=======================================================================
1710 0 : subroutine sh_create_gaus_grid(nlat,theta,wts,ierr)
1711 : !
1712 : ! spherepack gaqd
1713 : ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1714 : ! . .
1715 : ! . copyright (c) 2001 by ucar .
1716 : ! . .
1717 : ! . university corporation for atmospheric research .
1718 : ! . .
1719 : ! . all rights reserved .
1720 : ! . .
1721 : ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1722 : !
1723 : ! February 2002
1724 : !
1725 : ! gauss points and weights are computed using the fourier-newton
1726 : ! described in "on computing the points and weights for
1727 : ! gauss-legendre quadrature", paul n. swarztrauber, siam journal
1728 : ! on scientific computing (DOI 10.1137/S1064827500379690).
1729 : ! This routine is faster and more accurate than older program
1730 : ! with the same name.
1731 : !
1732 : ! computes the nlat gaussian colatitudes and weights
1733 : ! in double precision. the colatitudes are in radians and lie in the
1734 : ! in the interval (0,pi).
1735 : !
1736 : ! input parameters
1737 : !
1738 : ! nlat the number of gaussian colatitudes in the interval (0,pi)
1739 : ! (between the two poles). nlat must be greater than zero.
1740 : !
1741 : ! output parameters
1742 : !
1743 : ! theta a double precision array with length nlat
1744 : ! containing the gaussian colatitudes in
1745 : ! increasing radians on the interval (0,pi).
1746 : !
1747 : ! wts a double precision array with lenght nlat
1748 : ! containing the gaussian weights.
1749 : !
1750 : ! ierror = 0 no errors
1751 : ! = 1 if nlat<=0
1752 : !
1753 : !===================================================================
1754 : !
1755 : ! Passed variables
1756 : !-----------------
1757 : integer ,intent(in ) :: nlat
1758 : real(r8),intent(out) :: theta(nlat)
1759 : real(r8),intent(out) :: wts(nlat)
1760 : integer ,intent(out) :: ierr
1761 : !
1762 : ! Local Values
1763 : !-------------
1764 : real(r8):: sgnd
1765 : real(r8):: xx,dtheta,dthalf
1766 : real(r8):: cmax,zprev,zlast,zero,zhold,pb,dpb,dcor,summ,cz
1767 : integer :: mnlat,ns2,nhalf,nix,it,ii
1768 :
1769 : real(r8), parameter :: eps = epsilon(1._r8)
1770 :
1771 : ! check work space length
1772 : !------------------------
1773 0 : if(nlat<=0) then
1774 0 : ierr = 1
1775 0 : return
1776 : endif
1777 0 : ierr = 0
1778 :
1779 : ! compute weights and points analytically when nlat=1,2
1780 : !-------------------------------------------------------
1781 0 : if(nlat==1) then
1782 0 : theta(1) = acos(0._r8)
1783 0 : wts (1) = 2._r8
1784 0 : return
1785 0 : elseif(nlat==2) then
1786 0 : xx = sqrt(1._r8/3._r8)
1787 0 : theta(1) = acos( xx)
1788 0 : theta(2) = acos(-xx)
1789 0 : wts (1) = 1._r8
1790 0 : wts (2) = 1._r8
1791 0 : return
1792 : endif
1793 :
1794 : ! Proceed for nlat > 2
1795 : !----------------------
1796 0 : mnlat = mod(nlat,2)
1797 0 : ns2 = nlat/2
1798 0 : nhalf = (nlat+1)/2
1799 :
1800 0 : call sh_fourier_coefs_dp(nlat,cz,theta(ns2+1),wts(ns2+1))
1801 :
1802 0 : dtheta = halfPI/nhalf
1803 0 : dthalf = dtheta/2._r8
1804 0 : cmax = .2_r8*dtheta
1805 :
1806 : ! estimate first point next to theta = pi/2
1807 : !-------------------------------------------
1808 0 : if(mnlat/=0) then
1809 0 : zero = halfPI - dtheta
1810 0 : zprev = halfPI
1811 0 : nix = nhalf - 1
1812 : else
1813 0 : zero = halfPI - dthalf
1814 0 : nix = nhalf
1815 : endif
1816 :
1817 0 : do while(nix/=0)
1818 : dcor = huge(1._r8)
1819 : it = 0
1820 0 : do while (abs(dcor) > eps*abs(zero))
1821 : it = it + 1
1822 : ! newton iterations
1823 : !-----------------------
1824 0 : call sh_legp_dlegp_theta(nlat,zero,cz,theta(ns2+1),wts(ns2+1),pb,dpb)
1825 0 : dcor = pb/dpb
1826 0 : if(dcor.ne.0._r8) then
1827 0 : sgnd = dcor/abs(dcor)
1828 : else
1829 : sgnd = 1._r8
1830 : endif
1831 0 : dcor = sgnd*min(abs(dcor),cmax)
1832 0 : zero = zero - dcor
1833 : end do
1834 :
1835 0 : theta(nix) = zero
1836 0 : zhold = zero
1837 :
1838 : ! wts(nix) = (nlat+nlat+1)/(dpb*dpb)
1839 : ! yakimiw's formula permits using old pb and dpb
1840 : !--------------------------------------------------
1841 0 : wts(nix) = (nlat+nlat+1)/ (dpb+pb*dcos(zlast)/dsin(zlast))**2
1842 0 : nix = nix - 1
1843 0 : if(nix==nhalf-1) zero = 3._r8*zero - pi
1844 0 : if(nix<nhalf-1) zero = zero + zero - zprev
1845 0 : zprev = zhold
1846 : end do
1847 :
1848 : ! extend points and weights via symmetries
1849 : !-------------------------------------------
1850 0 : if(mnlat/=0) then
1851 0 : theta(nhalf) = halfPI
1852 0 : call sh_legp_dlegp_theta(nlat,halfPI,cz,theta(ns2+1),wts(ns2+1),pb,dpb)
1853 0 : wts(nhalf) = (nlat+nlat+1)/ (dpb*dpb)
1854 : endif
1855 :
1856 0 : do ii = 1,ns2
1857 0 : wts (nlat-ii+1) = wts(ii)
1858 0 : theta(nlat-ii+1) = pi - theta(ii)
1859 : end do
1860 :
1861 : summ = 0._r8
1862 0 : do ii = 1,nlat
1863 0 : summ = summ + wts(ii)
1864 : end do
1865 0 : do ii = 1,nlat
1866 0 : wts(ii) = 2._r8*wts(ii)/summ
1867 : end do
1868 :
1869 : end subroutine sh_create_gaus_grid
1870 : !=======================================================================
1871 :
1872 : !=======================================================================
1873 0 : subroutine sh_fourier_coefs_dp(nn,cz,cp,dcp)
1874 : !
1875 : ! spherepack cpdp
1876 : !
1877 : ! computes the fourier coefficients of the legendre
1878 : ! polynomial p_n^0 and its derivative.
1879 : ! nn is the degree and nn/2 or (nn+1)/2
1880 : ! coefficients are returned in cp depending on whether
1881 : ! nn is even or odd. The same number of coefficients
1882 : ! are returned in dcp. For nn even the constant
1883 : ! coefficient is returned in cz.
1884 : !=====================================================================
1885 : !
1886 : ! Passed variables
1887 : !-----------------
1888 : integer, intent(in) :: nn
1889 : real(r8), intent(out) :: cz
1890 : real(r8), intent(out) :: cp(nn/2+1)
1891 : real(r8), intent(out) :: dcp(nn/2+1)
1892 : !
1893 : ! Local Values
1894 : !--------------
1895 : real(r8):: t1,t2,t3,t4
1896 : integer :: ncp,jj
1897 :
1898 0 : ncp = (nn+1)/2
1899 0 : t1 = -1._r8
1900 0 : t2 = nn + 1._r8
1901 0 : t3 = 0._r8
1902 0 : t4 = nn + nn + 1._r8
1903 0 : if(mod(nn,2)==0) then
1904 0 : cp(ncp) = 1._r8
1905 0 : do jj = ncp,2,-1
1906 0 : t1 = t1 + 2._r8
1907 0 : t2 = t2 - 1._r8
1908 0 : t3 = t3 + 1._r8
1909 0 : t4 = t4 - 2._r8
1910 0 : cp(jj-1) = (t1*t2)/ (t3*t4)*cp(jj)
1911 : end do
1912 0 : t1 = t1 + 2._r8
1913 0 : t2 = t2 - 1._r8
1914 0 : t3 = t3 + 1._r8
1915 0 : t4 = t4 - 2._r8
1916 0 : cz = (t1*t2)/ (t3*t4)*cp(1)
1917 0 : do jj = 1,ncp
1918 0 : dcp(jj) = (jj+jj)*cp(jj)
1919 : end do
1920 : else
1921 0 : cp(ncp) = 1._r8
1922 0 : do jj = ncp - 1,1,-1
1923 0 : t1 = t1 + 2._r8
1924 0 : t2 = t2 - 1._r8
1925 0 : t3 = t3 + 1._r8
1926 0 : t4 = t4 - 2._r8
1927 0 : cp(jj) = (t1*t2)/ (t3*t4)*cp(jj+1)
1928 : end do
1929 0 : do jj = 1,ncp
1930 0 : dcp(jj) = (jj+jj-1)*cp(jj)
1931 : end do
1932 : endif
1933 :
1934 0 : end subroutine sh_fourier_coefs_dp
1935 : !=======================================================================
1936 :
1937 : !=======================================================================
1938 0 : subroutine sh_legp_dlegp_theta(nn,theta,cz,cp,dcp,pb,dpb)
1939 : !
1940 : ! spherepack tpdp
1941 : !
1942 : ! computes pn(theta) and its derivative dpb(theta) with
1943 : ! respect to theta
1944 : !=====================================================================
1945 : !
1946 : ! Passed variables
1947 : !------------------
1948 : integer, intent(in) :: nn
1949 : real(r8), intent(in) :: theta
1950 : real(r8), intent(in) :: cz
1951 : real(r8), intent(in) :: cp (nn/2+1)
1952 : real(r8), intent(in) :: dcp(nn/2+1)
1953 : real(r8), intent(out) :: pb
1954 : real(r8), intent(out) :: dpb
1955 : !
1956 : ! Local Values
1957 : !--------------
1958 : real(r8):: cdt,sdt,cth,sth,chh
1959 : integer :: kdo,kk
1960 :
1961 0 : cdt = dcos(theta+theta)
1962 0 : sdt = dsin(theta+theta)
1963 0 : if(mod(nn,2)==0) then
1964 : ! n even
1965 : !----------
1966 0 : kdo = nn/2
1967 0 : pb = .5_r8*cz
1968 0 : dpb = 0._r8
1969 0 : if(nn>0) then
1970 : cth = cdt
1971 : sth = sdt
1972 0 : do kk = 1,kdo
1973 0 : pb = pb + cp(kk)*cth
1974 0 : dpb = dpb - dcp(kk)*sth
1975 0 : chh = cdt*cth - sdt*sth
1976 0 : sth = sdt*cth + cdt*sth
1977 0 : cth = chh
1978 : end do
1979 : endif
1980 : else
1981 : ! n odd
1982 : !-----------
1983 0 : kdo = (nn+1)/2
1984 0 : pb = 0._r8
1985 0 : dpb = 0._r8
1986 0 : cth = dcos(theta)
1987 0 : sth = dsin(theta)
1988 0 : do kk = 1,kdo
1989 0 : pb = pb + cp(kk)*cth
1990 0 : dpb = dpb - dcp(kk)*sth
1991 0 : chh = cdt*cth - sdt*sth
1992 0 : sth = sdt*cth + cdt*sth
1993 0 : cth = chh
1994 : end do
1995 : endif
1996 :
1997 0 : end subroutine sh_legp_dlegp_theta
1998 : !=======================================================================
1999 :
2000 0 : end module zonal_mean_mod
|
