1module lamAnalysisGridTransforms_mod
2 ! MODULE lamAnalysisGridTransforms_mod (prefix='lgt' category='7. Low-level data objects')
3 !
4 !:Purpose: Performs some horizontal grid-point variable transforms
5 ! for the limited-area computational analysis grids (extended and
6 ! non-extended).
7 !
8 use earthConstants_mod
9 use mathPhysConstants_mod
10 use horizontalCoord_mod
11 use verticalCoord_mod
12 use midasMpi_mod
13 use utilities_mod
14 use Vgrid_Descriptors
15 implicit none
16 save
17 private
18
19 ! public procedures
20 public :: lgt_setupFromHCO, lgt_mach, lgt_mach_r4
21 public :: lgt_PsiChiToUV, lgt_PsiChiToUVAdj, lgt_UVToVortDiv
22 public :: lgt_createLamTemplateGrids
23
24 ! Definition of some parameters characterizing the geometry of
25 ! the Limited-Area (LA) analysis grid and associated metric factors
26 type :: struct_glmf
27 real(8), allocatable :: rlat (:) ! Latitudes of Scalar gridpoints
28 real(8), allocatable :: rlon (:) ! Longitudes of Scalar gridpoints
29 real(8) :: rdlon ! Latitude difference of gridpoints
30 real(8) :: rdlat ! Longitude differences of gridpoints
31 real(8), allocatable :: cos2 (:) ! Grid metric for Psi-Chi to U-V
32 real(8), allocatable :: cos2h(:) ! Grid metric for Psi-Chi to U-V
33 real(8), allocatable :: cos2vd (:) ! Grid metric for U-V to Vort-Div
34 real(8), allocatable :: cos2hvd(:) ! Grid metric for U-V to Vort-Div
35 real(8), allocatable :: idmuh(:) ! Grid metric for U-V to Vort-Div
36 real(8), allocatable :: idmu (:) ! Grid metric for U-V to Vort-Div
37 real(8) :: dx ! Grid-point spacing (uniform)
38 real(8), allocatable :: conphy (:) ! to go from Wind Images to true winds
39 real(8), allocatable :: conima (:) ! to go from true winds to Wind Images
40 end type struct_glmf
41
42 type(struct_hco), pointer :: hco_core => null()
43 type(struct_hco), pointer :: hco_ext => null()
44 type(struct_glmf):: glmf
45
46 integer :: ni_ext, nj_ext ! With Extension for Bi-Fourrier
47 integer :: ni_core, nj_core ! Without Extension for Bi-Fourrier
48 integer :: ext_i , ext_j ! Extension gridpoints
49
50 integer :: ni_ext_per, nj_ext_per
51
52 integer :: istart, iend, jstart, jend
53
54 integer :: LonPerPE, LonPerPEmax, myLonBeg, myLonEnd
55 integer :: LatPerPE, LatPerPEmax, myLatBeg, myLatEnd
56
57 logical :: initialized = .false.
58
59contains
60
61 !--------------------------------------------------------------------------
62 ! lgt_SetupFromHCO
63 !--------------------------------------------------------------------------
64 subroutine lgt_SetupFromHCO( hco_ext_in, hco_core_opt )
65 implicit none
66
67 ! Arguments:
68 type(struct_hco), pointer, intent(in) :: hco_ext_in
69 type(struct_hco), pointer, optional, intent(in) :: hco_core_opt
70
71 ! Locals:
72 real(8), allocatable :: rlath (:) ! Latitudes of half grid-points of gridpoints in lat-direction
73 real(8), allocatable :: rrcos (:) ! 1.0/Cos(Latitudes of gridpoints)
74 real(8), allocatable :: rdmu (:) ! Differences of mu=sin(lat)
75 real(8), allocatable :: rdmuh (:) ! Differences of muh=sin(lath)
76 real(8), allocatable :: r1mmu2 (:) ! (1.-mu**2)
77 real(8), allocatable :: r1mmu2h(:) ! (1.-muh**2)
78
79 real(8) :: dlon_test, dlat_test, dlon_ref, dlat_ref
80
81 integer :: i, j
82
83 logical, save :: firstCall = .true.
84
85 ! Ensure subroutine only runs one time during program execution
86 if (firstCall) then
87 firstCall = .false.
88 else
89 return
90 end if
91
92 write(*,*)
93 write(*,*) 'lgt_SetupFromHCO: Starting...'
94
95 initialized = .true.
96
97 !
98 !- 1. Check the hco structures
99 !
100 hco_ext => hco_ext_in
101 if( present(hco_core_opt) ) then
102 hco_core => hco_core_opt
103 else
104 hco_core => hco_ext_in
105 end if
106
107 if ( (.not. hco_core%initialized) .or. &
108 (.not. hco_ext%initialized) ) then
109 write(*,*)
110 write(*,*) 'lgt_SetupFromHCO: At least one hco structure was not initilzed'
111 write(*,*) 'hco_core = ', hco_core%initialized
112 write(*,*) 'hco_ext = ', hco_ext%initialized
113 call utl_abort('lgt_SetupFromHCO: abort')
114 end if
115
116 if ( (hco_core%global) .or. &
117 (hco_ext%global) ) then
118 write(*,*)
119 write(*,*) 'lgt_SetupFromHCO: At least one hco structure is from a global grid, skipping rest of setup'
120 write(*,*) 'hco_core = ', hco_core%global
121 write(*,*) 'hco_ext = ', hco_ext%global
122 return
123 end if
124
125 !
126 !- 2. Dimension settings and Memory allocation
127 !
128
129 !- 2.1 Dimensions
130 ni_core = hco_core%ni
131 nj_core = hco_core%nj
132
133 ni_ext = hco_ext%ni
134 nj_ext = hco_ext%nj
135
136 ext_i = ni_ext - ni_core
137 ext_j = nj_ext - nj_core
138
139 if ( ext_i == 0 .and. ext_j == 0 ) then
140 write(*,*)
141 write(*,*) 'lgt_SetupFromHCO: LAM core and extended grids are identical'
142 else if ( ext_i < 10 .or. ext_j < 10 ) then
143 write(*,*)
144 write(*,*) 'lgt_SetupFromHCO: LAM domain extension is less than 10 gridpoints'
145 write(*,*) ' ext_i = ', ext_i,' ext_j = ', ext_j
146 end if
147
148 ni_ext_per = ni_ext + 1 ! Fields will be periodic (i = 1 repeated) at ni_ext+1
149 nj_ext_per = nj_ext + 1 ! Fields will be periodic (j = 1 repeated) at nj_ext+1
150
151 !- 2.2 First Make sure we have a uniform grid in x and y direction
152
153 !- 2.2.1 Analysis core grid
154 dlon_ref = hco_core%lon(2) - hco_core%lon(1)
155 do i = 2, ni_core
156 dlon_test = hco_core%lon(i) - hco_core%lon(i-1)
157 if ( (dlon_test - dlon_ref) > dlon_ref/100.0d0 ) then
158 write(*,*)
159 write(*,*) 'lgt_SetupFromHCO: Core grid spacing is not uniform in x-direction'
160 write(*,*) ' i = ', i
161 write(*,*) ' dlon = ', dlon_test
162 write(*,*) ' dlon ref = ', dlon_ref
163 call utl_abort('lgt_SetupFromHCO')
164 end if
165 end do
166
167 dlat_ref = hco_core%lat(2) - hco_core%lat(1)
168 do j = 2, nj_core
169 dlat_test = hco_core%lat(j) - hco_core%lat(j-1)
170 if ( (dlat_test - dlat_ref) > dlat_ref/100.0d0 ) then
171 write(*,*)
172 write(*,*) 'lgt_SetupFromHCO: Core grid spacing is not uniform in x-direction'
173 write(*,*) ' j = ', j
174 write(*,*) ' dlat = ', dlon_test
175 write(*,*) ' dlat ref = ', dlon_ref
176 call utl_abort('lgt_SetupFromHCO')
177 end if
178 end do
179
180 !- 2.2.2 Extended Analysis grid
181 dlon_ref = hco_ext%lon(2) - hco_ext%lon(1)
182 do i = 2, ni_ext
183 dlon_test = hco_ext%lon(i) - hco_ext%lon(i-1)
184 if ( (dlon_test - dlon_ref) > dlon_ref/100.0d0 ) then
185 write(*,*)
186 write(*,*) 'lgt_SetupFromHCO: Extended grid spacing is not uniform in x-direction'
187 write(*,*) ' i = ', i
188 write(*,*) ' dlon = ', dlon_test
189 write(*,*) ' dlon ref = ', dlon_ref
190 call utl_abort('lgt_SetupFromHCO')
191 end if
192 end do
193
194 dlat_ref = hco_ext%lat(2) - hco_ext%lat(1)
195 do j = 2, nj_ext
196 dlat_test = hco_ext%lat(j) - hco_ext%lat(j-1)
197 if ( (dlat_test - dlat_ref) > dlat_ref/100.0d0 ) then
198 write(*,*)
199 write(*,*) 'lgt_SetupFromHCO: Extended grid spacing is not uniform in x-direction'
200 write(*,*) ' j = ', j
201 write(*,*) ' dlat = ', dlon_test
202 write(*,*) ' dlat ref = ', dlon_ref
203 call utl_abort('lgt_SetupFromHCO')
204 end if
205 end do
206
207 !- 2.3 Dimensions for variables needed to be symmetric
208 istart = -4 ! 5 gridpoints padding in West direction
209 iend = ni_ext + 4 ! 4 gridpoints padding in East direction
210 jstart = -4 ! 5 gridpoints padding in South direction
211 jend = nj_ext + 4 ! 4 gridpoints padding in North direction
212
213 !- 2.4 Allocations
214 allocate(glmf%rlon(istart:iend))
215 allocate(glmf%rlat(jstart:jend))
216 allocate(rlath (jstart:jend))
217 allocate(rrcos (jstart:jend))
218 allocate(rdmu (jstart:jend))
219 allocate(rdmuh (jstart:jend))
220 allocate(r1mmu2 (jstart:jend))
221 allocate(r1mmu2h (jstart:jend))
222
223 !
224 !- 3. Set Lat-Lon of the computational grid
225 !
226
227 ! 3.1 Set (constant) Grid spacing
228 glmf%rdlon = (hco_ext%lon(2) - hco_ext%lon(1))
229 glmf%rdlat = (hco_ext%lat(2) - hco_ext%lat(1))
230
231 !- 3.2 Lat-Lon
232
233 ! 3.2.1 Extract the lat-lon from the extended grid
234 glmf%rlon(1:ni_ext) = hco_ext%lon(1:ni_ext)
235 glmf%rlat(1:nj_ext) = hco_ext%lat(1:nj_ext)
236
237 ! 3.2.2 Extend to the full computational grid (larger than the extended grid!)
238
239 ! West
240 do i = istart, 0
241 glmf%rlon(i) = glmf%rlon(1) + (i-1) * glmf%rdlon
242 end do
243 ! East
244 do i = ni_ext+1, iend
245 glmf%rlon(i) = glmf%rlon(ni_ext) + (i-ni_ext) * glmf%rdlon
246 end do
247 ! North
248 do j = nj_ext+1, jend
249 glmf%rlat(j) = glmf%rlat(nj_ext) + (j-nj_ext) * glmf%rdlat
250 end do
251 ! South
252 do j = jstart, 0
253 glmf%rlat(j) = glmf%rlat(1) + (j-1) * glmf%rdlat
254 end do
255
256 !- 3.2 Half Lat-Lon
257 do j = jstart+1, jend-1
258 rlath(j) = ( glmf%rlat(j) + glmf%rlat(j+1) ) / 2.0d0
259 end do
260
261 !
262 !- 4. Set Metric Factors
263 !
264
265 !- 4.1 Compute local factors
266 do j = jstart+1, jend-2
267 rrcos (j) = 1.0d0 / cos(glmf%rlat (j))
268 rdmu (j) = sin(glmf%rlat (j+1)) - sin(glmf%rlat (j))
269 rdmuh (j) = sin(rlath(j+1)) - sin(rlath(j))
270 r1mmu2 (j) = (cos(glmf%rlat (j)))**2
271 r1mmu2h(j) = (cos(rlath(j)))**2
272 end do
273
274 !- 4.2 Bi-periodize and symmetrize Metric coefficients
275 call lgt_mach(rdmu (1:nj_ext), & ! INOUT
276 1, nj_ext,1) ! IN
277 call lgt_mach(rdmuh (1:nj_ext), & ! INOUT
278 1, nj_ext,1) ! IN
279 call lgt_mach(r1mmu2 (1:nj_ext), & ! INOUT
280 1, nj_ext,1) ! IN
281 call lgt_mach(r1mmu2h(1:nj_ext), & ! INOUT
282 1, nj_ext,1) ! IN
283
284 call symmetrize_coef(rdmu ) ! INOUT
285 call symmetrize_coef(rdmuh ) ! INOUT
286 call symmetrize_coef(r1mmu2 ) ! INOUT
287 call symmetrize_coef(r1mmu2h) ! INOUT
288
289 !- 4.3 Compute global factors
290 glmf%dx = 1.d0 / (ec_ra * glmf%rdlon)
291
292 allocate(glmf%cos2 ( 0:nj_ext+1))
293 allocate(glmf%cos2h ( 0:nj_ext+1))
294 do j = 0, nj_ext+1
295 glmf%cos2 (j) = r1mmu2 (j) / (ec_ra * rdmuh(j-1))
296 glmf%cos2h(j) = r1mmu2h(j) / (ec_ra * rdmu (j ))
297 end do
298
299 allocate(glmf%idmuh (0:nj_ext-1))
300 allocate(glmf%idmu (1:nj_ext ))
301 allocate(glmf%cos2vd (1:nj_ext ))
302 allocate(glmf%cos2hvd(1:nj_ext ))
303
304 do j = 1, nj_ext
305 glmf%idmuh (j-1) = 1.d0 / rdmuh(j-1)
306 glmf%idmu (j) = 1.d0 / rdmu (j)
307 glmf%cos2vd (j) = 1.d0 / r1mmu2 (j)
308 glmf%cos2hvd(j) = 1.d0 / r1mmu2h(j)
309 end do
310
311 !
312 ! Conversion of wind images to physical winds and vice-versa
313 ! N.B.: Those are geometrical factors of the COMPUTATIONAL grid
314 ! ==> only computational latitude variation...
315 !
316 allocate(glmf%conphy(nj_ext))
317 allocate(glmf%conima(nj_ext))
318
319 call lgt_mach(rrcos(1:nj_ext), & ! INOUT
320 1, nj_ext,1) ! IN
321
322 do j = 1, nj_ext
323 glmf%conphy(j) = rrcos(j) ! to go from Wind Images to true winds
324 glmf%conima(j) = 1.0d0 / glmf%conphy(j) ! to go from true winds to Wind Images
325 end do
326
327 !
328 !- 5. MPI partitionning
329 !
330 if ( mmpi_nprocs /= 0 ) then
331 call mmpi_setup_lonbands(ni_ext, & ! IN
332 lonPerPE, lonPerPEmax, myLonBeg, myLonEnd ) ! OUT
333
334 call mmpi_setup_latbands(nj_ext, & ! IN
335 latPerPE, latPerPEmax, myLatBeg, myLatEnd ) ! OUT
336 else
337 ! This option is needed for the biper program
338 lonPerPE = ni_ext
339 myLonBeg = 1
340 myLonEnd = ni_ext
341 latPerPE = nj_ext
342 myLatBeg = 1
343 myLatEnd = nj_ext
344 write(*,*)
345 write(*,*) 'WARNING: the module will be executed in NO-MPI MODE!'
346 end if
347
348 !
349 !- 6. Ending
350 !
351 deallocate(rlath )
352 deallocate(rrcos )
353 deallocate(rdmu )
354 deallocate(rdmuh )
355 deallocate(r1mmu2 )
356 deallocate(r1mmu2h)
357
358 write(*,*)
359 write(*,*) 'lgt_SetupFromHCO: Done!'
360
361 end subroutine lgt_SetupFromHCO
362
363 !--------------------------------------------------------------------------
364 ! symmetrize_coef
365 !--------------------------------------------------------------------------
366 subroutine symmetrize_coef(coef_inout)
367 !
368 !:Purpose: Extend symmetrically Metric coefficients.
369 !
370 implicit none
371
372 ! Arguments:
373 real(8), intent(inout) :: coef_inout(jstart:jend)
374
375 ! Locals:
376 integer :: j
377
378 do j = jstart, 0
379 coef_inout(j) = coef_inout(nj_ext+j)
380 end do
381
382 do j = nj_ext+1, jend
383 coef_inout(j) = coef_inout(j-nj_ext)
384 end do
385
386 end subroutine symmetrize_coef
387
388 !--------------------------------------------------------------------------
389 ! lgt_PsiChiToUV
390 !--------------------------------------------------------------------------
391 subroutine lgt_PsiChiToUV(psi, chi, uphy, vphy, nk)
392 implicit none
393
394 ! Arguments:
395 integer, intent(in) :: nk
396 real(8), intent(in) :: psi(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
397 real(8), intent(in) :: chi(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
398 real(8), intent(out) :: uphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
399 real(8), intent(out) :: vphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
400
401 ! Locals:
402 real(8), allocatable :: psi_ext(:,:,:)
403 real(8), allocatable :: chi_ext(:,:,:)
404 real(8), allocatable :: uimg(:,:,:)
405 real(8), allocatable :: vimg(:,:,:)
406 real(8), allocatable :: uimgs(:,:,:)
407 real(8), allocatable :: vimgs(:,:,:)
408 integer :: i,j,k
409
410 if ( hco_ext%global ) then
411 call utl_abort('lgt_PsiChiToUV: Not compatible with global grid')
412 endif
413
414 if ( .not. initialized ) then
415 call utl_abort('lgt_PsiChiToUV: AnalysisGrid not initialized')
416 endif
417
418 allocate(psi_ext( (myLonBeg-1):(myLonEnd+1), (myLatBeg-1):(myLatEnd+1), nk))
419 allocate(chi_ext( (myLonBeg-1):(myLonEnd+1), (myLatBeg-1):(myLatEnd+1), nk))
420 allocate(uimg ( myLonBeg:myLonEnd , myLatBeg:myLatEnd , nk))
421 allocate(vimg ( myLonBeg:myLonEnd , myLatBeg:myLatEnd , nk))
422 allocate(uimgs ( (myLonBeg-1):myLonEnd , myLatBeg:(myLatEnd+1) , nk))
423 allocate(vimgs ( myLonBeg:(myLonEnd+1) , (myLatBeg-1):myLatEnd , nk))
424
425 !
426 !- 1. Symmetrize
427 !
428 call symmetrize( psi_ext, & ! OUT
429 psi, myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
430 call symmetrize( chi_ext, & ! OUT
431 chi, myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
432
433 !
434 !- 2. Compute Wind on staggered grid
435 !
436
437 !$OMP PARALLEL DO PRIVATE (k,j,i)
438 do k = 1, nk
439 !- 2.1 u-wind component
440 do j = myLatBeg, myLatEnd+1
441 do i = myLonBeg-1, myLonEnd
442 uimgs(i,j,k) = glmf%dx * ( chi_ext(i+1,j,k) - chi_ext(i,j,k) ) &
443 - glmf%cos2(j) * ( psi_ext(i,j,k) - psi_ext(i,j-1,k) )
444 end do
445 end do
446 !- 2.2 v-wind component
447 do j = myLatBeg-1, myLatEnd
448 do i = myLonBeg, myLonEnd+1
449 vimgs(i,j,k) = glmf%cos2h(j) * ( chi_ext(i,j+1,k) - chi_ext(i,j,k) ) &
450 + glmf%dx * ( psi_ext(i,j,k) - psi_ext(i-1,j,k) )
451 end do
452 end do
453 end do
454 !$OMP END PARALLEL DO
455
456 !
457 !- 3. Move to collocated (scalar) grid
458 !
459 call uvStagToColloc( uimgs, vimgs, & ! IN
460 uimg , vimg , & ! OUT
461 myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
462
463 !
464 !- 4. Convert Wind images to Physical (true) winds
465 !
466 !$OMP PARALLEL DO PRIVATE (j)
467 do j = myLatBeg, myLatEnd
468 uphy(:,j,:) = glmf%conphy(j) * uimg(:,j,:)
469 vphy(:,j,:) = glmf%conphy(j) * vimg(:,j,:)
470 end do
471 !$OMP END PARALLEL DO
472
473 deallocate(psi_ext)
474 deallocate(chi_ext)
475 deallocate(uimg)
476 deallocate(vimg)
477 deallocate(uimgs)
478 deallocate(vimgs)
479
480 end subroutine lgt_PsiChiToUV
481
482 !--------------------------------------------------------------------------
483 ! lgt_PsiChiToUVAdj
484 !--------------------------------------------------------------------------
485 subroutine lgt_PsiChiToUVAdj(psi, chi, uphy, vphy, nk)
486 implicit none
487
488 ! Arguments:
489 integer, intent(in) :: nk
490 real(8), intent(out) :: psi(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
491 real(8), intent(out) :: chi(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
492 real(8), intent(in) :: uphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
493 real(8), intent(in) :: vphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
494
495 ! Locals:
496 real(8), allocatable :: psi_ext(:,:,:)
497 real(8), allocatable :: chi_ext(:,:,:)
498 real(8), allocatable :: uimg(:,:,:)
499 real(8), allocatable :: vimg(:,:,:)
500 real(8), allocatable :: uimgs(:,:,:)
501 real(8), allocatable :: vimgs(:,:,:)
502 integer :: i,j,k
503
504 if ( hco_ext%global ) then
505 call utl_abort('lgt_PsiChiToUVAdj: Not compatible with global grid')
506 endif
507
508 if ( .not. initialized ) then
509 call utl_abort('lgt_PsiChiToUV: AnalysisGrid not initialized')
510 endif
511
512 allocate(psi_ext( (myLonBeg-1):(myLonEnd+1), (myLatBeg-1):(myLatEnd+1), nk))
513 allocate(chi_ext( (myLonBeg-1):(myLonEnd+1), (myLatBeg-1):(myLatEnd+1), nk))
514 allocate(uimg ( myLonBeg:myLonEnd , myLatBeg:myLatEnd , nk))
515 allocate(vimg ( myLonBeg:myLonEnd , myLatBeg:myLatEnd , nk))
516 allocate(uimgs ( (myLonBeg-1):myLonEnd , myLatBeg:(myLatEnd+1) , nk))
517 allocate(vimgs ( myLonBeg:(myLonEnd+1) , (myLatBeg-1):myLatEnd , nk))
518
519 !
520 !- 4. Convert Physical (true) winds to Wind images
521 !
522 !$OMP PARALLEL DO PRIVATE (j)
523 do j = myLatBeg, myLatEnd
524 uimg(:,j,:) = glmf%conphy(j) * uphy(:,j,:)
525 vimg(:,j,:) = glmf%conphy(j) * vphy(:,j,:)
526 end do
527 !$OMP END PARALLEL DO
528
529 !
530 !- 3. Move to stagerred grid
531 !
532 uimgs(:,:,:) = 0.d0
533 vimgs(:,:,:) = 0.d0
534
535 call uvStagToCollocAdj( uimgs, vimgs, & ! OUT
536 uimg , vimg , & ! IN
537 myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
538
539 !
540 !- 2. Compute Psi-Chi
541 !
542 chi_ext(:,:,:) = 0.d0
543 psi_ext(:,:,:) = 0.d0
544
545 !$OMP PARALLEL DO PRIVATE (k,j,i)
546 do k = 1, nk
547 !- 2.2 from v-wind component
548 do j = myLatEnd, myLatBeg-1, -1
549 do i = myLonEnd+1, myLonBeg, -1
550 chi_ext(i ,j+1,k) = chi_ext(i ,j+1,k) + vimgs(i,j,k) * glmf%cos2h(j)
551 chi_ext(i ,j ,k) = chi_ext(i ,j ,k) - vimgs(i,j,k) * glmf%cos2h(j)
552 psi_ext(i ,j ,k) = psi_ext(i ,j ,k) + vimgs(i,j,k) * glmf%dx
553 psi_ext(i-1,j ,k) = psi_ext(i-1,j ,k) - vimgs(i,j,k) * glmf%dx
554 end do
555 end do
556 !- 2.1 from u-wind component
557 do j = myLatEnd+1, myLatBeg, -1
558 do i = myLonEnd, myLonBeg-1, -1
559 chi_ext(i+1,j ,k) = chi_ext(i+1,j ,k) + uimgs(i,j,k) * glmf%dx
560 chi_ext(i ,j ,k) = chi_ext(i ,j ,k) - uimgs(i,j,k) * glmf%dx
561 psi_ext(i ,j ,k) = psi_ext(i ,j ,k) - uimgs(i,j,k) * glmf%cos2(j)
562 psi_ext(i ,j-1,k) = psi_ext(i ,j-1,k) + uimgs(i,j,k) * glmf%cos2(j)
563 end do
564 end do
565 end do
566 !$OMP END PARALLEL DO
567
568 !
569 !- 1. De-Symmetrize
570 !
571 chi(:,:,:) = 0.d0
572 psi(:,:,:) = 0.d0
573
574 call symmetrizeAdj( psi_ext, & ! IN
575 psi, & ! OUT
576 myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
577 call symmetrizeAdj( chi_ext, & ! IN
578 chi, & ! OUT
579 myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
580
581 deallocate(psi_ext)
582 deallocate(chi_ext)
583 deallocate(uimg)
584 deallocate(vimg)
585 deallocate(uimgs)
586 deallocate(vimgs)
587
588 end subroutine lgt_PsiChiToUVAdj
589
590 !--------------------------------------------------------------------------
591 ! uvStagToColloc
592 !--------------------------------------------------------------------------
593 subroutine uvStagToColloc(uStag, vStag, uColloc, vColloc, iBeg, iEnd, jBeg, jEnd , nk)
594 implicit none
595
596 ! Arguments:
597 integer, intent(in) :: iBeg
598 integer, intent(in) :: iEnd
599 integer, intent(in) :: jBeg
600 integer, intent(in) :: JEnd
601 integer, intent(in) :: nk
602 real(8), intent(out) :: uColloc(iBeg:iEnd ,jBeg :jEnd ,nk)
603 real(8), intent(out) :: vColloc(iBeg:iEnd ,jBeg :jEnd ,nk)
604 real(8), intent(in) :: uStag (iBeg-1:iEnd,jBeg :jEnd+1,nk)
605 real(8), intent(in) :: vStag (iBeg:iEnd+1,jBeg-1:jEnd ,nk)
606
607 ! Locals:
608 integer :: i,j,k
609
610 !$OMP PARALLEL DO PRIVATE (k,j,i)
611 do k = 1, nk
612 do j = jBeg, jEnd
613 do i = iBeg, iEnd
614 uColloc(i,j,k) = ( uStag(i-1,j ,k) + uStag(i,j,k) ) / 2.0d0
615 vColloc(i,j,k) = ( vStag(i ,j-1,k) + vStag(i,j,k) ) / 2.0d0
616 end do
617 end do
618 end do
619 !$OMP END PARALLEL DO
620
621 end subroutine uvStagToColloc
622
623 !--------------------------------------------------------------------------
624 ! uvStagToCollocAdj
625 !--------------------------------------------------------------------------
626 subroutine uvStagToCollocAdj(uStag, vStag, uColloc, vColloc, iBeg, iEnd, jBeg, jEnd , nk)
627 implicit none
628
629 ! Arguments:
630 integer, intent(in) :: iBeg
631 integer, intent(in) :: iEnd
632 integer, intent(in) :: jBeg
633 integer, intent(in) :: jEnd
634 integer, intent(in) :: nk
635 real(8), intent(in) :: uColloc(iBeg:iEnd ,jBeg :jEnd ,nk)
636 real(8), intent(in) :: vColloc(iBeg:iEnd ,jBeg :jEnd ,nk)
637 real(8), intent(out) :: uStag (iBeg-1:iEnd,jBeg :jEnd+1,nk)
638 real(8), intent(out) :: vStag (iBeg:iEnd+1,jBeg-1:jEnd ,nk)
639
640 ! Locals:
641 integer :: i, j, k
642
643 !$OMP PARALLEL DO PRIVATE (k,j,i)
644 do k = 1, nk
645 do j = jEnd, jBeg, -1
646 do i = iEnd, iBeg, -1
647 uStag(i-1,j ,k) = uStag(i-1,j ,k) + uColloc(i,j,k) / 2.0d0
648 uStag(i ,j ,k) = uStag(i ,j ,k) + uColloc(i,j,k) / 2.0d0
649 vStag(i ,j-1,k) = vStag(i ,j-1,k) + vColloc(i,j,k) / 2.0d0
650 vStag(i ,j ,k) = vStag(i ,j ,k) + vColloc(i,j,k) / 2.0d0
651 end do
652 end do
653 end do
654 !$OMP END PARALLEL DO
655
656 end subroutine uvStagToCollocAdj
657
658 !--------------------------------------------------------------------------
659 ! Symmetrize
660 !--------------------------------------------------------------------------
661 subroutine symmetrize(field_out, field_in, iBeg, iEnd, jBeg, jEnd, nk)
662 !
663 !:Purpose: Extend symmetrically outside 1 grid point all around LAM-boundary
664 ! ready for finite differences
665 !
666 implicit none
667
668 ! Arguments:
669 integer, intent(in) :: iBeg
670 integer, intent(in) :: iEnd
671 integer, intent(in) :: jBeg
672 integer, intent(in) :: jEnd
673 integer, intent(in) :: nk
674 real(8), intent(out):: field_out(iBeg-1:iEnd+1,jBeg-1:jEnd+1, nk)
675 real(8), intent(in) :: field_in(iBeg:iEnd,jBeg:jEnd,nk)
676
677 ! Locals:
678 real(8), allocatable :: field_8(:,:,:)
679
680 integer :: ni,nj
681
682 ni = iEnd-iBeg+1
683 nj = jEnd-jBeg+1
684
685 allocate(field_8(0:ni+1,0:nj+1, nk))
686 field_8(:,:,:) = 0.d0
687
688 field_8(1:ni,1:nj,:) = field_in(iBeg:iEnd,jBeg:jEnd,:)
689
690 call RPN_COMM_xch_halo_8(field_8, & ! INOUT
691 0,ni+1,0,nj+1,ni,nj,nk, & ! IN
692 1,1,.true.,.true.,ni,0) ! IN
693
694 field_out(iBeg-1:iEnd+1,jBeg-1:jEnd+1,:) = field_8(0:ni+1,0:nj+1,:)
695 deallocate(field_8)
696
697 end subroutine symmetrize
698
699 !--------------------------------------------------------------------------
700 ! SymmetrizeAdj
701 !--------------------------------------------------------------------------
702 subroutine symmetrizeAdj(field_in, field_out, iBeg, iEnd, jBeg, jEnd, nk)
703 !
704 !:Purpose: Adjoint of sub. symmetrize.
705 !
706 implicit none
707
708 ! Arguments:
709 integer, intent(in) :: iBeg
710 integer, intent(in) :: iEnd
711 integer, intent(in) :: jBeg
712 integer, intent(in) :: jEnd
713 integer, intent(in) :: nk
714 real(8), intent(in) :: field_in(iBeg-1:iEnd+1,jBeg-1:jEnd+1, nk)
715 real(8), intent(out) :: field_out(iBeg:iEnd,jBeg:jEnd,nk)
716
717 ! Locals:
718 real(8), allocatable :: field_8(:,:,:)
719 integer :: i,j,k,ni,nj
720
721 ni = iEnd-iBeg+1
722 nj = jEnd-jBeg+1
723
724 allocate(field_8(0:ni+1,0:nj+1, nk))
725 field_8(0:ni+1,0:nj+1,:) = field_in(iBeg-1:iEnd+1,jBeg-1:jEnd+1,:)
726
727 call RPN_COMM_adj_halo8(field_8, & ! INOUT
728 0,ni+1,0,nj+1,ni,nj,nk, & ! IN
729 1,1,.true.,.true.,ni,0) ! IN
730
731 !$OMP PARALLEL DO PRIVATE (k,j,i)
732 do k = 1, nk
733 do j = jBeg, jEnd
734 do i = iBeg, iEnd
735 field_out(i,j,k) = field_out(i,j,k) + field_8(i-iBeg+1,j-jBeg+1,k)
736 end do
737 end do
738 end do
739 !$OMP END PARALLEL DO
740
741 deallocate(field_8)
742
743 end subroutine symmetrizeAdj
744
745 !--------------------------------------------------------------------------
746 ! lgt_Mach
747 !--------------------------------------------------------------------------
748 subroutine lgt_mach(gd,ni,nj,nk)
749 !
750 !:Purpose: [to be completed]
751 !
752 !:Arguments:
753 ! :ni: Maximum I-dimension where the input array is assumed to carry
754 ! information. Will be used as I-limit where backward
755 ! derivatives will be evaluated
756 !
757 ! :nj: Maximum J-dimension where the input array is assumed to carry
758 ! information. Will be used as J-limit where backward derivatives
759 ! will be evaluated
760 !
761 implicit none
762
763 ! Arguments:
764 integer, intent(in) :: ni
765 integer, intent(in) :: nj
766 integer, intent(in) :: nk
767 real(8), intent(inout) :: gd(ni,nj,nk)
768
769 ! Locals:
770 integer :: istart, jstart
771 integer :: i,j,k
772 real(8) :: con, xp, yp, a0, a1, b1, b2
773 real(8) :: deriv_istart, deriv_jstart, deriv_i0, deriv_j0, del
774
775 if ( hco_ext%global ) then
776 call utl_abort('lgt_Mach: Not compatible with global grid')
777 endif
778
779 if ( .not. initialized ) then
780 call utl_abort('lgt_Mach: AnalysisGrid not initialized')
781 endif
782
783 if ( (ni /= 1 .and. ni /= ni_ext) .or. &
784 (nj /= 1 .and. nj /= nj_ext) ) then
785 call utl_abort('lgt_Mach: Invalid Dimensions')
786 end if
787
788 !$OMP PARALLEL
789 !$OMP DO PRIVATE (k,j,i,istart,jstart,con,a0,a1,b1,b2,del,deriv_istart,deriv_i0,deriv_jstart,deriv_j0,xp,yp)
790 do k = 1, nk
791 !
792 !- 1. Periodicized in x-direction from ni_core to ni
793 !
794 if ( ni > 1 ) then
795 istart = ni_core - 2 ! I-limit where backward derivatives will be evaluated
796 con = 1.d0 / ( glmf%rdlon * MPC_DEGREES_PER_RADIAN_R8 * 111.d+3)
797 do j = 1, nj
798 a0 = 0.5d0 * ( gd(istart,j,k) + gd(1,j,k) )
799 a1 = 0.5d0 * ( gd(istart,j,k) - gd(1,j,k) )
800 deriv_istart = con * ( gd(istart,j,k) - gd(istart-1,j,k))
801 deriv_i0 = con * ( gd(2,j,k) - gd(1,j,k) )
802 b1 = 0.5d0 * ( deriv_istart - deriv_i0 )
803 b2 = 0.25d0* ( deriv_istart + deriv_i0 )
804 del = real(ni_ext_per-istart,8)
805 do i = istart, ni
806 xp = MPC_PI_R8 * real(i-istart,8) / del
807 gd(i,j,k) = a0 + a1*cos(xp) + b1*sin(xp) + b2*sin(2.d0*xp)
808 end do
809 end do
810 end if
811 !
812 !- 2. Periodicized in y-direction from nj_core to nj
813 !
814 if ( nj > 1 ) then
815 jstart = nj_core - 2
816 con = 1.d0 / (glmf%rdlat * MPC_DEGREES_PER_RADIAN_R8 * 111.d+3)
817 do i = 1, ni
818 a0 = 0.5d0 * ( gd(i,jstart,k) + gd(i,1,k) )
819 a1 = 0.5d0 * ( gd(i,jstart,k) - gd(i,1,k) )
820 deriv_jstart = con * ( gd(i,jstart,k) - gd(i,jstart-1,k))
821 deriv_j0 = con * ( gd(i,2,k) - gd(i,1,k) )
822 b1 = 0.5d0 * ( deriv_jstart - deriv_j0 )
823 b2 = 0.25d0* ( deriv_jstart + deriv_j0 )
824 del = real(nj_ext_per-jstart,8)
825 do j = jstart , nj
826 yp = MPC_PI_R8 * real(j-jstart,8) / del
827 gd(i,j,k) = a0 + a1*cos(yp) + b1*sin(yp) + b2*sin(2.d0*yp)
828 end do
829 end do
830 end if
831
832 end do
833 !$OMP END DO
834 !$OMP END PARALLEL
835
836 end subroutine lgt_mach
837
838 !--------------------------------------------------------------------------
839 ! lgt_Mach_r4
840 !--------------------------------------------------------------------------
841 subroutine lgt_mach_r4(gd,ni,nj,nk)
842 !:Purpose: [to be completed]
843 !
844 !:Arguments:
845 ! :ni: Maximum I-dimension where the input array is assumed to carry
846 ! information. Will be used as I-limit where backward
847 ! derivatives will be evaluated
848 ! :nj: Maximum J-dimension where the input array is assumed to carry
849 ! information. Will be used as J-limit where backward
850 ! derivatives will be evaluated
851 !
852 implicit none
853
854 ! Arguments:
855 integer, intent(in) :: ni
856 integer, intent(in) :: nj
857 integer, intent(in) :: nk
858 real(4), intent(inout) :: gd(ni,nj,nk)
859
860 ! Locals:
861 integer :: istart, jstart
862 integer :: i,j,k
863 real(4) :: con, xp, yp, a0, a1, b1, b2
864 real(4) :: deriv_istart, deriv_jstart, deriv_i0, deriv_j0, del
865
866 if ( hco_ext%global ) then
867 call utl_abort('lgt_Mach_r4: Not compatible with global grid')
868 endif
869
870 if ( .not. initialized ) then
871 call utl_abort('lgt_Mach_r4: AnalysisGrid not initialized')
872 endif
873
874 if ( (ni /= 1 .and. ni /= ni_ext) .or. &
875 (nj /= 1 .and. nj /= nj_ext) ) then
876 call utl_abort('lgt_Mach_r4 : Invalid Dimensions')
877 end if
878
879 !$OMP PARALLEL
880 !$OMP DO PRIVATE (k,j,i,istart,jstart,con,a0,a1,b1,b2,del,deriv_istart,deriv_i0,deriv_jstart,deriv_j0,xp,yp)
881 do k = 1, nk
882 !
883 !- 1. Periodicized in x-direction from ni_core to ni
884 !
885 if ( ni > 1 ) then
886 istart = ni_core - 2 ! I-limit where backward derivatives will be evaluated
887 con = 1.0 / ( glmf%rdlon * MPC_DEGREES_PER_RADIAN_R4 * 111.e+3)
888 do j = 1, nj
889 a0 = 0.50 * ( gd(istart,j,k) + gd(1,j,k) )
890 a1 = 0.50 * ( gd(istart,j,k) - gd(1,j,k) )
891 deriv_istart = con * ( gd(istart,j,k) - gd(istart-1,j,k))
892 deriv_i0 = con * ( gd(2,j,k) - gd(1,j,k) )
893 b1 = 0.50 * ( deriv_istart - deriv_i0 )
894 b2 = 0.250* ( deriv_istart + deriv_i0 )
895 del = real(ni_ext_per-istart,8)
896 do i = istart, ni
897 xp = MPC_PI_R4 * real(i-istart,8) / del
898 gd(i,j,k) = a0 + a1*cos(xp) + b1*sin(xp) + b2*sin(2.0*xp)
899 end do
900 end do
901 end if
902 !
903 !- 2. Periodicized in y-direction from nj_core to nj
904 !
905 if ( nj > 1 ) then
906 jstart = nj_core - 2
907 con = 1.0 / (glmf%rdlat * MPC_DEGREES_PER_RADIAN_R4 * 111.e+3)
908 do i = 1, ni
909 a0 = 0.50 * ( gd(i,jstart,k) + gd(i,1,k) )
910 a1 = 0.50 * ( gd(i,jstart,k) - gd(i,1,k) )
911 deriv_jstart = con * ( gd(i,jstart,k) - gd(i,jstart-1,k))
912 deriv_j0 = con * ( gd(i,2,k) - gd(i,1,k) )
913 b1 = 0.50 * ( deriv_jstart - deriv_j0 )
914 b2 = 0.250* ( deriv_jstart + deriv_j0 )
915 del = real(nj_ext_per-jstart,8)
916 do j = jstart , nj
917 yp = MPC_PI_R4 * real(j-jstart,8) / del
918 gd(i,j,k) = a0 + a1*cos(yp) + b1*sin(yp) + b2*sin(2.0*yp)
919 end do
920 end do
921 end if
922
923 end do
924 !$OMP END DO
925 !$OMP END PARALLEL
926
927 end subroutine lgt_mach_r4
928
929 !--------------------------------------------------------------------------
930 ! lgt_UVToVortDiv
931 !--------------------------------------------------------------------------
932 subroutine lgt_UVToVortDiv(Vorticity, Divergence, uphy, vphy, nk)
933 implicit none
934
935 ! Arguments:
936 integer, intent(in) :: nk
937 real(8), intent(out) :: Vorticity (myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
938 real(8), intent(out) :: Divergence(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
939 real(8), intent(in) :: uphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
940 real(8), intent(in) :: vphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nk)
941
942 ! Locals:
943 real(8), allocatable :: uimg(:,:,:)
944 real(8), allocatable :: vimg(:,:,:)
945 real(8), allocatable :: uimg_sym(:,:,:)
946 real(8), allocatable :: vimg_sym(:,:,:)
947 integer :: i,j,k
948
949 if ( hco_ext%global ) then
950 call utl_abort('lgt_UVToVortDiv: Not compatible with global grid')
951 endif
952
953 if ( .not. initialized ) then
954 call utl_abort('lgt_UVToVortDiv: AnalysisGrid not initialized')
955 endif
956
957 allocate(uimg_sym( (myLonBeg-1):(myLonEnd+1), (myLatBeg-1):(myLatEnd+1), nk))
958 allocate(vimg_sym( (myLonBeg-1):(myLonEnd+1), (myLatBeg-1):(myLatEnd+1), nk))
959 allocate(uimg ( myLonBeg:myLonEnd , myLatBeg:myLatEnd , nk))
960 allocate(vimg ( myLonBeg:myLonEnd , myLatBeg:myLatEnd , nk))
961
962 !
963 !- 2. Convert Physical (true) winds to Wind images
964 !
965 !$OMP PARALLEL DO PRIVATE (j)
966 do j = myLatBeg, myLatEnd
967 uimg(:,j,:) = glmf%conima(j) * uphy(:,j,:)
968 vimg(:,j,:) = glmf%conima(j) * vphy(:,j,:)
969 end do
970 !$OMP END PARALLEL DO
971
972 !
973 !- 3. Symmetrize
974 !
975 call symmetrize( uimg_sym, & ! OUT
976 uimg, myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
977 call symmetrize( vimg_sym, & ! OUT
978 vimg, myLonBeg, myLonEnd, myLatBeg, myLatEnd, nk ) ! IN
979 deallocate(uimg,vimg)
980
981 !
982 !- 4. Compute Vorticity and Divergence
983 !
984
985 !$OMP PARALLEL DO PRIVATE (k,j,i)
986 do k = 1, nk
987 do j = myLatBeg, myLatEnd
988 do i = myLonBeg, myLonEnd
989 vorticity (i,j,k) = &
990 glmf%cos2hvd(j) * glmf%dx * (vimg_sym(i+1,j,k)-vimg_sym(i,j,k)) - &
991 ec_r1sa * glmf%idmu(j) * (uimg_sym(i,j+1,k)-uimg_sym(i,j,k))
992
993 divergence(i,j,k) = &
994 glmf%cos2vd(j) * glmf%dx * (uimg_sym(i,j,k)-uimg_sym(i-1,j,k)) + &
995 ec_r1sa * glmf%idmuh(j-1)* (vimg_sym(i,j,k)-vimg_sym(i,j-1,k))
996 end do
997 end do
998 end do
999 !$OMP END PARALLEL DO
1000
1001 deallocate(uimg_sym)
1002 deallocate(vimg_sym)
1003
1004 end subroutine lgt_UVToVortDiv
1005
1006 !--------------------------------------------------------------------------
1007 ! lgt_createLamTemplateGrids
1008 !--------------------------------------------------------------------------
1009 subroutine lgt_createLamTemplateGrids(templateFileName, hco_core, vco, &
1010 grd_ext_x, grd_ext_y)
1011 implicit none
1012
1013 ! Arguments:
1014 character(len=*), intent(in) :: templateFileName
1015 type(struct_hco), intent(in) :: hco_core
1016 type(struct_vco), intent(in) :: vco
1017 integer , intent(in) :: grd_ext_x
1018 integer , intent(in) :: grd_ext_y
1019
1020 ! Locals:
1021 integer :: ni_ext, nj_ext, i, j, lev, ni, nj, nk
1022 integer :: iun = 0
1023 integer :: ier, fnom, fstouv, fstfrm, fclos, fstecr
1024 real(8), allocatable :: Field2d(:,:)
1025 real(8), allocatable :: lat_ext(:)
1026 real(8), allocatable :: lon_ext(:)
1027 real(4), allocatable :: dummy2D(:,:)
1028 real(8) :: dlat, dlon
1029 real(4) :: work
1030 integer :: dateo,npak,status
1031 integer :: ip1,ip2,ip3,deet,npas,datyp,ig1,ig2,ig3,ig4
1032 integer :: ig1_tictac,ig2_tictac,ig3_tictac,ig4_tictac
1033 character(len=1) :: grtyp
1034 character(len=2) :: typvar
1035 character(len=12) :: etiket
1036
1037 !
1038 !- 1. Opening the output template file
1039 !
1040 ier = fnom(iun, trim(templateFileName), 'RND', 0)
1041 ier = fstouv(iun, 'RND')
1042
1043 npak = -32
1044
1045 !
1046 !- 2. Writing the core grid (Ensemble) template
1047 !
1048
1049 !- 2.1 Tic-Tac
1050 deet = 0
1051 ip1 = hco_core%ig1
1052 ip2 = hco_core%ig2
1053 ip3 = hco_core%ig3
1054 npas = 0
1055 datyp = 1
1056 grtyp = hco_core%grtypTicTac
1057 typvar = 'X'
1058 etiket = 'COREGRID'
1059 dateo = 0
1060
1061 call cxgaig ( grtyp, & ! IN
1062 ig1_tictac, ig2_tictac, ig3_tictac, ig4_tictac, & ! OUT
1063 real(hco_core%xlat1), real(hco_core%xlon1), & ! IN
1064 real(hco_core%xlat2), real(hco_core%xlon2) ) ! IN
1065
1066 ig1 = ig1_tictac
1067 ig2 = ig2_tictac
1068 ig3 = ig3_tictac
1069 ig4 = ig4_tictac
1070
1071 ier = utl_fstecr(hco_core%lon*MPC_DEGREES_PER_RADIAN_R8, npak, &
1072 iun, dateo, deet, npas, hco_core%ni, 1, 1, ip1, &
1073 ip2, ip3, typvar, '>>', etiket, grtyp, ig1, &
1074 ig2, ig3, ig4, datyp, .true.)
1075
1076 ier = utl_fstecr(hco_core%lat*MPC_DEGREES_PER_RADIAN_R8, npak, &
1077 iun, dateo, deet, npas, 1, hco_core%nj, 1, ip1, &
1078 ip2, ip3, typvar, '^^', etiket, grtyp, ig1, &
1079 ig2, ig3, ig4, datyp, .true.)
1080
1081 !- 2.2 2D Field
1082 allocate(Field2d(hco_core%ni,hco_core%nj))
1083 Field2d(:,:) = 10.d0
1084
1085 deet = 0
1086 ip1 = 0
1087 ip2 = 0
1088 ip3 = 0
1089 npas = 0
1090 datyp = 1
1091 grtyp = hco_core%grtyp
1092 typvar = 'A'
1093 etiket = 'COREGRID'
1094 dateo = 0
1095 ig1 = hco_core%ig1
1096 ig2 = hco_core%ig2
1097 ig3 = hco_core%ig3
1098 ig4 = hco_core%ig4
1099
1100 ier = utl_fstecr(Field2d, npak, &
1101 iun, dateo, deet, npas, hco_core%ni, hco_core%nj, 1, ip1, &
1102 ip2, ip3, typvar, 'P0', etiket, grtyp, ig1, &
1103 ig2, ig3, ig4, datyp, .true.)
1104
1105 deallocate(Field2d)
1106
1107 !
1108 !- 3. Create and Write the extended grid (Analysis) template
1109 !
1110 ni_ext = hco_core%ni + grd_ext_x
1111 nj_ext = hco_core%nj + grd_ext_y
1112
1113 !- 3.1 Tic-Tac
1114 allocate(lon_ext(ni_ext))
1115 allocate(lat_ext(nj_ext))
1116
1117 !- Copy core grid info
1118 lon_ext(1:hco_core%ni) = hco_core%lon(:)
1119 lat_ext(1:hco_core%nj) = hco_core%lat(:)
1120
1121 !- Extend the lat lon
1122 dlon = hco_core%lon(2) - hco_core%lon(1)
1123 do i = hco_core%ni + 1, ni_ext
1124 lon_ext(i) = lon_ext(hco_core%ni) + (i - hco_core%ni) * dlon
1125 end do
1126
1127 dlat = hco_core%lat(2) - hco_core%lat(1)
1128 do j = hco_core%nj + 1, nj_ext
1129 lat_ext(j) = lat_ext(hco_core%nj) + (j - hco_core%nj) * dlat
1130 end do
1131
1132 !- Write
1133 deet = 0
1134 ip1 = hco_core%ig1 + 100 ! Must be different from the core grid
1135 ip2 = hco_core%ig2 + 100 ! Must be different from the core grid
1136 if (hco_core%ig3 > 0) then
1137 ip3 = hco_core%ig3 + 100 ! Must be different from the core grid
1138 else
1139 ip3 = 0
1140 end if
1141 npas = 0
1142 datyp = 1
1143 grtyp = hco_core%grtypTicTac
1144 typvar = 'X'
1145 etiket = 'ANALYSIS'
1146 dateo = 0
1147 ig1 = ig1_tictac
1148 ig2 = ig2_tictac
1149 ig3 = ig3_tictac
1150 ig4 = ig4_tictac
1151
1152 ier = utl_fstecr(lon_ext*MPC_DEGREES_PER_RADIAN_R8, npak, &
1153 iun, dateo, deet, npas, ni_ext, 1, 1, ip1, &
1154 ip2, ip3, typvar, '>>', etiket, grtyp, ig1, &
1155 ig2, ig3, ig4, datyp, .true.)
1156
1157 ier = utl_fstecr(lat_ext*MPC_DEGREES_PER_RADIAN_R8, npak, &
1158 iun, dateo, deet, npas, 1, nj_ext, 1, ip1, &
1159 ip2, ip3, typvar, '^^', etiket, grtyp, ig1, &
1160 ig2, ig3, ig4, datyp, .true.)
1161
1162 deallocate(lon_ext)
1163 deallocate(lat_ext)
1164
1165 !- 3.2 2D Field
1166 allocate(Field2d(ni_ext,nj_ext))
1167 Field2d(:,:) = 10.d0
1168
1169 deet = 0
1170 ip1 = 0
1171 ip2 = 0
1172 ip3 = 0
1173 npas = 0
1174 datyp = 1
1175 grtyp = hco_core%grtyp
1176 typvar = 'A'
1177 etiket = 'ANALYSIS'
1178 dateo = 0
1179 ig1 = hco_core%ig1 + 100 ! Must be different from the core grid
1180 ig2 = hco_core%ig2 + 100 ! Must be different from the core grid
1181 if (hco_core%ig3 > 0) then
1182 ig3 = hco_core%ig3 + 100 ! Must be different from the core grid
1183 else
1184 ig3 = 0
1185 end if
1186 ig4 = 0
1187
1188 ier = utl_fstecr(Field2d, npak, &
1189 iun, dateo, deet, npas, ni_ext, nj_ext, 1, ip1, &
1190 ip2, ip3, typvar, 'P0', etiket, grtyp, ig1, &
1191 ig2, ig3, ig4, datyp, .true.)
1192
1193 deallocate(Field2d)
1194
1195 !
1196 !- 4. Write the vertical grid description
1197 !
1198
1199 if (vco%vgridPresent) then
1200 !- 4.1 Write the toc-toc
1201 status = vgd_write(vco%vgrid,iun,'fst')
1202
1203 if ( status /= VGD_OK ) then
1204 call utl_abort('createLamTemplateGrids: ERROR with vgd_write')
1205 end if
1206
1207 !- 4.2 Write a dummy 2D field for each MM and TH levels
1208 npak = -12
1209 dateo = 0
1210 deet = 0
1211 npas = 0
1212 ni = 4
1213 nj = 2
1214 nk = 1
1215 ip2 = 0
1216 ip3 = 0
1217 typvar = 'A'
1218 etiket = 'VERTICALGRID'
1219 grtyp = 'G'
1220 ig1 = 0
1221 ig2 = 0
1222 ig3 = 0
1223 ig4 = 0
1224 datyp = 1
1225
1226 allocate(dummy2D(ni,nj))
1227 dummy2D(:,:) = 0.0
1228
1229 do lev = 1, vco%nlev_M
1230 ip1 = vco%ip1_M(lev)
1231 ier = fstecr(dummy2D, work, npak, iun, dateo, deet, npas, ni, nj, &
1232 nk, ip1, ip2, ip3, typvar, 'MM', etiket, grtyp, &
1233 ig1, ig2, ig3, ig4, datyp, .true.)
1234 end do
1235 do lev = 1, vco%nlev_T
1236 ip1 = vco%ip1_T(lev)
1237 ier = fstecr(dummy2D, work, npak, iun, dateo, deet, npas, ni, nj, &
1238 nk, ip1, ip2, ip3, typvar, 'TH', etiket, grtyp, &
1239 ig1, ig2, ig3, ig4, datyp, .true.)
1240 end do
1241
1242 deallocate(dummy2D)
1243 end if ! vco%vgridPresent
1244
1245 !
1246 !- 5. Closing the output template file
1247 !
1248 ier = fstfrm(iun)
1249 ier = fclos (iun)
1250
1251 end subroutine lgt_createLamTemplateGrids
1252
1253end module lamAnalysisGridTransforms_mod
