1module horizontalCoord_mod
2 ! MODULE horizontalCoord_mod (prefix='hco' category='7. Low-level data objects')
3 !
4 !:Purpose: Derived type and procedures related to the horizontal grid
5 ! coordinate for various grids (global and limited area).
6 !
7 use midasMpi_mod
8 use earthConstants_mod
9 use mathPhysConstants_mod
10 use utilities_mod
11 use varNameList_mod
12 use physicsFunctions_mod
13
14 implicit none
15 save
16 private
17
18 ! Public derived type
19 public :: struct_hco
20
21 ! Public Subroutines
22 public :: hco_SetupFromFile, hco_equal, hco_deallocate, hco_mpiBcast, hco_weight, hco_setupYgrid
23
24 integer, parameter :: maxNumSubGrid = 2
25
26 type :: struct_hco
27 character(len=32) :: gridname
28 logical :: initialized = .false.
29 integer :: ni
30 integer :: nj
31 character(len=1) :: grtyp
32 character(len=1) :: grtypTicTac
33 integer :: ig1
34 integer :: ig2
35 integer :: ig3
36 integer :: ig4
37 integer :: EZscintID = -1
38 integer :: numSubGrid
39 integer :: EZscintIDsubGrids(maxNumSubGrid)
40 real(8), allocatable :: lat(:) ! in radians
41 real(8), allocatable :: lon(:) ! in radians
42 real(4), allocatable :: lat2d_4(:,:) ! in radians
43 real(4), allocatable :: lon2d_4(:,:) ! in radians
44 real(8) :: dlat ! in radians
45 real(8) :: dlon ! in radians
46 real(8) :: maxGridSpacing ! in meter
47 real(8) :: minGridSpacing ! in meter
48 logical :: global
49 logical :: rotated
50 real(8) :: xlat1, xlat1_yan
51 real(8) :: xlon1, xlon1_yan
52 real(8) :: xlat2, xlat2_yan
53 real(8) :: xlon2, xlon2_yan
54 real(4), allocatable :: tictacU(:)
55 end type struct_hco
56
57contains
58
59 !--------------------------------------------------------------------------
60 ! hco_SetupFromFile
61 !--------------------------------------------------------------------------
62 subroutine hco_SetupFromFile(hco, TemplateFile, EtiketName, GridName_opt, &
63 varName_opt)
64 !
65 !:Purpose: to initialize hco structure from a template file
66 !
67 implicit none
68
69 ! Arguments:
70 type(struct_hco), pointer, intent(out) :: hco
71 character(len=*), intent(in) :: TemplateFile
72 character(len=*), intent(in) :: EtiketName
73 character(len=*), optional, intent(in) :: GridName_opt
74 character(len=*), optional, intent(in) :: varName_opt
75
76 ! Locals:
77 real(8), allocatable :: lat_8(:)
78 real(8), allocatable :: lon_8(:)
79 real(8) :: maxDeltaLat, maxDeltaLon, maxGridSpacing, deltaLon, deltaLat
80 real(8) :: minDeltaLat, minDeltaLon, minGridSpacing
81 real(8) :: deltaLon1, deltaLon2, deltaLon3
82 real(8) :: deltaLat1, deltaLat2, deltaLat3
83 real(8), save :: maxGridSpacingPrevious = -1.0d0
84 real(8), save :: minGridSpacingPrevious = -1.0d0
85 real(4) :: xlat1_4, xlon1_4, xlat2_4, xlon2_4
86 real(4) :: xlat1_yan_4, xlon1_yan_4, xlat2_yan_4, xlon2_yan_4
87 integer :: iu_template, numSubGrid, varIndex
88 integer :: fnom, fstlir, fstouv, fstfrm, fclos
89 integer :: ezqkdef, ezget_nsubgrids, ezget_subgridids, ezgprm
90 integer :: key, fstinf, fstprm, ier, EZscintID, EZscintIDsubGrids(maxNumSubGrid)
91 integer :: ni, nj, ni_tictacU, ni_t, nj_t, nlev_t, gdll
92 integer :: dateo, deet, npas, nk, nbits, datyp
93 integer :: ip1, ip2, ip3, swa, lng, dltf, ubc
94 integer :: extra1, extra2, extra3
95 integer :: ig1, ig2, ig3, ig4
96 integer :: ig1_tictac, ig2_tictac, ig3_tictac, ig4_tictac
97 integer :: ni_yy, nj_yy, ig1_yy, ig2_yy, ig3_yy, ig4_yy
98 integer :: latIndex, lonIndex, latIndexBeg, latIndexEnd
99 real(4), parameter :: absMaxLat = 85. ! abs of latitude threshold where to compute minGridSpacing in 3.2
100 logical :: FileExist, global, rotated, foundVarNameInFile
101 character(len=4 ) :: nomvar
102 character(len=2 ) :: typvar
103 character(len=1 ) :: grtyp, grtypTicTac
104 character(len=12) :: etiket
105
106 if(.not.associated(hco)) then
107 allocate(hco)
108 else
109 call utl_abort('hco_setupFromFile: supplied hco must be null')
110 end if
111
112 !
113 !- 1.1 Determine which variable to use for defining the grid
114 !
115 if (present(varName_opt)) then
116
117 ! User specified variable name
118 nomvar = varName_opt
119
120 else
121
122 ! First try to use P0
123 nomvar = 'P0'
124
125 if (.not. utl_varNamePresentInFile(nomvar,fileName_opt=trim(TemplateFile))) then
126 ! P0 not present, look for another suitable variable in the file
127
128 foundVarNameInFile = .false.
129 do varIndex = 1, vnl_numvarmax
130 nomvar = vnl_varNameList(varIndex)
131
132 ! check if variable is in the file
133 if (.not. utl_varNamePresentInFile(nomvar,fileName_opt=trim(TemplateFile))) cycle
134
135 foundVarNameInFile = .true.
136 exit
137
138 end do
139
140 if (.not. foundVarNameInFile) call utl_abort('hco_SetupFromFile: NO variables found in the file!!!')
141
142 end if
143
144 end if
145
146 write(*,*) 'hco_SetupFromFile: defining hco by varname= ', nomvar
147
148
149 !
150 !- 1.2 Open/Check template file
151 !
152 inquire(file=trim(TemplateFile), exist=FileExist)
153
154 if (FileExist) then
155 iu_template = 0
156 ier = fnom(iu_template,trim(TemplateFile),'RND+OLD+R/O',0)
157 if (ier == 0) then
158 write(*,*)
159 write(*,*) 'hco_setupFromFile: Template File =', trim(TemplateFile)
160 ier = fstouv(iu_template,'RND+OLD')
161 else
162 write(*,*)
163 write(*,*) 'hco_SetupFromFile: Error in opening the template grid file'
164 write(*,*) trim(TemplateFile)
165 call utl_abort('hco_SetupFromFile')
166 end if
167 else
168 write(*,*)
169 write(*,*) 'hco_SetupFromFile: template grid file DOES NOT EXIST'
170 write(*,*) trim(TemplateFile)
171 call utl_abort('hco_SetupFromFile')
172 end if
173
174 !
175 !- 2. Get Horizontal grid info
176 !
177
178 !- 2.1 Grid size and grid projection info
179 dateo = -1
180 etiket = EtiketName
181 ip1 = -1
182 ip2 = -1
183 ip3 = -1
184 typvar = ' '
185
186 key = fstinf(iu_template, & ! IN
187 ni, nj, nk, & ! OUT
188 dateo, etiket, ip1, ip2, ip3, typvar, nomvar)! IN
189
190 if (key < 0) then
191 write(*,*)
192 write(*,*) 'hco_SetupFromFile: Unable to find output horiz grid info using = ',nomvar
193 write(*,*) ' with etiket = ',trim(EtiketName)
194 call utl_abort('hco_setupFromFile: unable to setup the structure')
195 end if
196
197 ier = fstprm(key, & ! IN
198 dateo, deet, npas, ni, nj, nk, nbits, & ! OUT
199 datyp, ip1, ip2, ip3, typvar, nomvar, etiket, & ! OUT
200 grtyp, ig1, ig2, ig3, & ! OUT
201 ig4, swa, lng, dltf, ubc, extra1, extra2, extra3) ! OUT
202
203 if (trim(grtyp) == 'G' .and. ig2 == 1) then
204 call utl_abort('hco_setupFromFile: ERROR: due to bug in ezsint, Gaussian grid with ig2=1 no longer supported')
205 end if
206
207 EZscintID = ezqkdef(ni, nj, grtyp, ig1, ig2, ig3, ig4, iu_template) ! IN
208 numSubGrid = 1
209 EZscintIDsubGrids(:) = MPC_missingValue_INT
210
211 allocate(lat_8(1:nj))
212 allocate(lon_8(1:ni))
213
214 allocate(hco%lat2d_4(1:ni,1:nj))
215 allocate(hco%lon2d_4(1:ni,1:nj))
216
217 ier = gdll(EZscintID, & ! IN
218 hco%lat2d_4, hco%lon2d_4) ! OUT
219
220 xlat1_yan_4 = MPC_missingValue_R4
221 xlon1_yan_4 = MPC_missingValue_R4
222 xlat2_yan_4 = MPC_missingValue_R4
223 xlon2_yan_4 = MPC_missingValue_R4
224
225 grtypTicTac = 'X'
226
227 if (mmpi_myid == 0) write(*,*) 'hco_setupFromFile: grtyp, ni, nj, EZscintID = ', grtyp, ni, nj, EZscintID
228
229 !- 2.2 Rotated lat-lon grid
230 if (trim(grtyp) == 'Z') then
231
232 !- 2.2.1 Read the Longitudes
233 dateo = -1
234 etiket = ''
235 ip1 = ig1
236 ip2 = ig2
237 ip3 = ig3
238 typvar = ''
239 nomvar = '>>'
240
241 ier = utl_fstlir(lon_8, & ! OUT
242 iu_template, & ! IN
243 ni_t, nj_t, nlev_t, & ! OUT
244 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
245
246 if (ier < 0) then
247 write(*,*)
248 write(*,*) 'hco_SetupFromFile: Unable to find >> grid descriptors'
249 call utl_abort('hco_setupFromFile')
250 end if
251
252 ! Test if the dimensions are compatible with the grid
253 if (ni_t /= ni .or. nj_t /= 1) then
254 write(*,*)
255 write(*,*) 'hco_SetupFromFile: Incompatible >> grid descriptors !'
256 write(*,*) 'Found :', ni_t, nj_t
257 write(*,*) 'Should be :', ni, 1
258 call utl_abort('hco_setupFromFile')
259 end if
260
261 !- 2.2.2 Read the latitudes
262 dateo = -1
263 etiket = ''
264 ip1 = ig1
265 ip2 = ig2
266 ip3 = ig3
267 typvar = ''
268 nomvar = '^^'
269
270 ier = utl_fstlir(lat_8, & ! OUT
271 iu_template, & ! IN
272 ni_t, nj_t, nlev_t, & ! OUT
273 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
274
275 if (ier < 0) then
276 write(*,*)
277 write(*,*) 'hco_SetupFromFile: Unable to find ^^ grid descriptors'
278 call utl_abort('hco_setupFromFile')
279 end if
280
281 ! Test if the dimensions are compatible with the grid
282 if (ni_t /= 1 .or. nj_t /= nj) then
283 write(*,*)
284 write(*,*) 'hco_SetupFromFile: Incompatible ^^ grid descriptors !'
285 write(*,*) 'Found :', ni_t, nj_t
286 write(*,*) 'Should be :', 1, nj
287 call utl_abort('hco_setupFromFile')
288 end if
289
290 !- 2.2.3 Do we have a rotated grid ?
291 dateo = -1
292 etiket = ''
293 ip1 = ig1
294 ip2 = ig2
295 ip3 = ig3
296 typvar = ''
297 nomvar = '^^'
298
299 key = fstinf(iu_template, & ! IN
300 ni_t, nj_t, nk, & ! OUT
301 dateo, etiket, ip1, ip2, ip3, typvar, nomvar) ! IN
302
303 ier = fstprm(key, & ! IN
304 dateo, deet, npas, ni_t, nj_t, nk, nbits, & ! OUT
305 datyp, ip1, ip2, ip3, typvar, nomvar, etiket, & ! OUT
306 grtypTicTac, ig1_tictac, ig2_tictac, & ! OUT
307 ig3_tictac, ig4_tictac, swa, lng, dltf, & ! OUT
308 ubc, extra1, extra2, extra3) ! OUT
309
310 call cigaxg (grtypTicTac, & ! IN
311 xlat1_4, xlon1_4, xlat2_4, xlon2_4, & ! OUT
312 ig1_tictac, ig2_tictac, ig3_tictac, ig4_tictac) ! IN
313
314 if (xlat1_4 == 0.0 .and. xlat2_4 == 0.0) then
315 rotated = .false.
316 else
317 rotated = .true.
318 end if
319 if (mmpi_myid == 0) then
320 write(*,*) 'hco_setupFromFile: xlat1/2, xlon1/2, rotated = ', &
321 xlat1_4, xlat2_4, xlon1_4, xlon2_4, rotated
322 end if
323
324 !- 2.2.4 Is this a global or a LAM domain ?
325 call global_or_lam(global, & ! OUT
326 lon_8, ni) ! IN
327
328 !- 2.3 B Grid
329 else if (trim(grtyp) == 'B') then
330
331 !- 2.3.1 Find the latitudes and longitudes
332 lon_8(:) = real(hco%lon2d_4(:,nj/2),8)
333 lat_8(:) = real(hco%lat2d_4(1,:),8)
334
335 !- 2.3.2 This grid type is not rotated
336 rotated = .false.
337 xlat1_4 = 0.0
338 xlon1_4 = 180.0
339 xlat2_4 = 0.0
340 xlon2_4 = 180.0
341
342 !- 2.3.3 We know this is a global grid
343 global = .true.
344
345 !- 2.4 Gaussian Grid
346 else if (trim(grtyp) == 'G') then
347
348 !- 2.4.1 Find the latitudes and longitudes
349 lon_8(:) = real(hco%lon2d_4(:,nj/2),8)
350 lat_8(:) = real(hco%lat2d_4(1,:),8)
351
352 !- 2.4.2 This grid type is not rotated
353 rotated = .false.
354 xlat1_4 = 0.0
355 xlon1_4 = 180.0
356 xlat2_4 = 0.0
357 xlon2_4 = 180.0
358
359 !- 2.4.3 We know this is a global grid
360 global = .true.
361
362 !- 2.5 Universal Grid (Yin-Yang) - not fully supported: use at own risk!
363 else if (trim(grtyp) == 'U') then
364
365 !- 2.5.1 Read the tic-tac vector
366 dateo = -1
367 etiket = ' '
368 ip1 = ig1
369 ip2 = ig2
370 ip3 = ig3
371 typvar = 'X'
372 nomvar = '^>'
373
374 ni_tictacU = 5 + 2 * (10 + ni + nj/2)
375 allocate(hco%tictacU(ni_tictacU))
376 ier = fstlir(hco%tictacU, & ! OUT
377 iu_template, & ! IN
378 ni_t, nj_t, nlev_t, & ! OUT
379 dateo, etiket, ip1, ip2, ip3, typvar, nomvar) ! IN
380
381 if (ier < 0) then
382 write(*,*)
383 write(*,*) 'hco_SetupFromFile: Unable to find ^> grid descriptors'
384 call utl_abort('hco_setupFromFile')
385 end if
386
387 ! Test if the dimensions are compatible with the grid
388 if (ni_t /= ni_tictacU .or. nj_t /= 1) then
389 write(*,*)
390 write(*,*) 'hco_SetupFromFile: Incompatible ^> grid descriptors !'
391 write(*,*) 'Found :', ni_t, nj_t
392 write(*,*) 'Should be :', ni_tictacU, 1
393 call utl_abort('hco_setupFromFile')
394 end if
395
396 !- 2.5.1 Initialize latitudes and longitudes to dummy values - should not be used!
397 lon_8(:) = MPC_missingValue_R8
398 lat_8(:) = MPC_missingValue_R8
399
400 !- 2.5.2 Yin-Yan subgrid IDs
401 numSubGrid = ezget_nsubgrids(EZscintID)
402 if (numSubGrid /= 2) then
403 call utl_abort('hco_setupFromFile: CONFUSED! number of sub grids must be 2')
404 end if
405 ier = ezget_subgridids(EZscintID, EZscintIDsubGrids)
406
407 !- 2.5.3 Determine parameters related to Yin and Yan grid rotations
408 rotated = .true. ! since Yin-Yan is made up of 2 grids with different rotations
409
410 ier = ezgprm(EZscintIDsubGrids(1), grtypTicTac, ni_yy, nj_yy, ig1_yy, ig2_yy, ig3_yy, ig4_yy)
411 grtypTicTac = 'E' ! needed since ezgprm returns 'Z', but grtyp for tictac should be 'E'
412 call cigaxg (grtypTicTac, & ! IN
413 xlat1_4, xlon1_4, xlat2_4, xlon2_4, & ! OUT
414 ig1_yy, ig2_yy, ig3_yy, ig4_yy) ! IN
415
416 ier = ezgprm(EZscintIDsubGrids(2), grtypTicTac, ni_yy, nj_yy, ig1_yy, ig2_yy, ig3_yy, ig4_yy)
417 grtypTicTac = 'E' ! needed since ezgprm returns 'Z', but grtyp for tictac should be 'E'
418 call cigaxg (grtypTicTac, & ! IN
419 xlat1_yan_4, xlon1_yan_4, xlat2_yan_4, xlon2_yan_4, & ! OUT
420 ig1_yy, ig2_yy, ig3_yy, ig4_yy) ! IN
421
422 rotated = .true. ! since Yin-Yan is made up of 2 grids with different rotations
423
424 !- 2.5.3 We know this is a global grid
425 global = .true.
426
427 !- 2.5 Irregular structure
428 else if (trim(grtyp) == 'Y') then
429
430 !- 2.6.1 This grid type is not rotated
431 rotated = .false.
432 xlat1_4 = 0.0
433 xlon1_4 = 0.0
434 xlat2_4 = 1.0
435 xlon2_4 = 1.0
436
437 grtypTicTac = 'L'
438
439 !- 2.6.2 Test using first row of longitudes (should work for ORCA grids)
440 lon_8(:) = hco%lon2d_4(:,1)
441 call global_or_lam(global, & ! OUT
442 lon_8, ni) ! IN
443
444 !- 2.6.3 Initialize latitudes and longitudes to dummy values - should not be used!
445 lon_8(:) = MPC_missingValue_R8
446 lat_8(:) = MPC_missingValue_R8
447
448 else
449 write(*,*)
450 write(*,*) 'hco_SetupFromFile: Only grtyp = Z or G or B or U or Y are supported !, grtyp = ', trim(grtyp)
451 call utl_abort('hco_setupFromFile')
452 end if
453
454 !
455 !- 3. Initialized Horizontal Grid Structure
456 !
457 allocate(hco%lat(1:nj))
458 allocate(hco%lon(1:ni))
459
460 if (present(gridName_opt)) then
461 hco%gridname = trim(gridName_opt)
462 else
463 hco%gridname = 'UNDEFINED'
464 end if
465 hco%ni = ni
466 hco%nj = nj
467 hco%grtyp = trim(grtyp)
468 hco%grtypTicTac = trim(grtypTicTac)
469 hco%ig1 = ig1
470 hco%ig2 = ig2
471 hco%ig3 = ig3
472 hco%ig4 = ig4
473 hco%EZscintID = EZscintID
474 hco%numSubGrid = numSubGrid
475 hco%EZscintIDsubGrids(:) = EZscintIDsubGrids(:)
476 hco%lon(:) = lon_8(:) * MPC_RADIANS_PER_DEGREE_R8
477 hco%lat(:) = lat_8(:) * MPC_RADIANS_PER_DEGREE_R8
478 hco%dlon = (lon_8(2) - lon_8(1)) * MPC_RADIANS_PER_DEGREE_R8
479 hco%dlat = (lat_8(2) - lat_8(1)) * MPC_RADIANS_PER_DEGREE_R8
480 hco%global = global
481 hco%rotated = rotated
482 hco%xlat1 = real(xlat1_4,8)
483 hco%xlon1 = real(xlon1_4,8)
484 hco%xlat2 = real(xlat2_4,8)
485 hco%xlon2 = real(xlon2_4,8)
486 hco%xlat1_yan = real(xlat1_yan_4,8)
487 hco%xlon1_yan = real(xlon1_yan_4,8)
488 hco%xlat2_yan = real(xlat2_yan_4,8)
489 hco%xlon2_yan = real(xlon2_yan_4,8)
490 hco%initialized = .true.
491
492 hco%lat2d_4(:,:) = hco%lat2d_4(:,:) * MPC_RADIANS_PER_DEGREE_R8
493 hco%lon2d_4(:,:) = hco%lon2d_4(:,:) * MPC_RADIANS_PER_DEGREE_R8
494
495 deallocate(lat_8)
496 deallocate(lon_8)
497
498 !- 3.1 Compute maxGridSpacing
499
500 latIndexBeg = 1
501 if (trim(grtyp) == 'U') then
502 latIndexEnd = nj / 2
503 else
504 latIndexEnd = nj
505 end if
506
507 maxDeltaLat = 0.0d0
508 do lonIndex = 1, ni - 1
509 do latIndex = latIndexBeg, latIndexEnd - 1
510
511 deltaLat1 = abs(hco%lat2d_4(lonIndex, latIndex) - hco%lat2d_4(lonIndex , latIndex + 1))
512 deltaLat2 = abs(hco%lat2d_4(lonIndex, latIndex) - hco%lat2d_4(lonIndex + 1, latIndex ))
513 deltaLat3 = abs(hco%lat2d_4(lonIndex, latIndex) - hco%lat2d_4(lonIndex + 1, latIndex + 1))
514
515 deltaLat = max(deltaLat1, deltaLat2, deltaLat3)
516 if (deltaLat > maxDeltaLat) maxDeltaLat = deltaLat
517
518 end do
519 end do
520
521 maxDeltaLon = 0.0d0
522 do lonIndex = 1, ni - 1
523 do latIndex = latIndexBeg, latIndexEnd - 1
524
525 deltaLon1 = abs(hco%lon2d_4(lonIndex, latIndex) - hco%lon2d_4(lonIndex , latIndex + 1))
526 deltaLon2 = abs(hco%lon2d_4(lonIndex, latIndex) - hco%lon2d_4(lonIndex + 1, latIndex ))
527 deltaLon3 = abs(hco%lon2d_4(lonIndex, latIndex) - hco%lon2d_4(lonIndex + 1, latIndex + 1))
528
529 if (deltaLon1 > MPC_PI_R8) deltaLon1 = deltaLon1 - 2.0d0 * MPC_PI_R8
530 deltaLon1 = abs(deltaLon1 * cos(hco%lat2d_4(lonIndex,latIndex)))
531 if (deltaLon2 > MPC_PI_R8) deltaLon2 = deltaLon2 - 2.0d0 * MPC_PI_R8
532 deltaLon2 = abs(deltaLon2 * cos(hco%lat2d_4(lonIndex,latIndex)))
533 if (deltaLon3 > MPC_PI_R8) deltaLon3 = deltaLon3 - 2.0d0 * MPC_PI_R8
534 deltaLon3 = abs(deltaLon3 * cos(hco%lat2d_4(lonIndex,latIndex)))
535
536 deltaLon = max(deltaLon1, deltaLon2, deltaLon3)
537 if (deltaLon > maxDeltaLon) maxDeltaLon = deltaLon
538
539 end do
540 end do
541
542 maxGridSpacing = ec_ra * sqrt(2.0d0) * max(maxDeltaLon, maxDeltaLat)
543
544 if (mmpi_myid == 0 .and. maxGridSpacing /= maxGridSpacingPrevious) then
545 maxGridSpacingPrevious = maxGridSpacing
546 write(*,*) 'hco_setupFromFile: maxDeltaLat=', maxDeltaLat * MPC_DEGREES_PER_RADIAN_R8, ' deg'
547 write(*,*) 'hco_setupFromFile: maxDeltaLon=', maxDeltaLon * MPC_DEGREES_PER_RADIAN_R8, ' deg'
548 write(*,*) 'hco_setupFromFile: maxGridSpacing=', maxGridSpacing, ' m'
549 end if
550
551 if (maxGridSpacing > 1.0d6) then
552 call utl_abort('hco_setupFromFile: maxGridSpacing is greater than 1000 km.')
553 end if
554
555 hco%maxGridSpacing = maxGridSpacing
556
557 !- 3.2 Compute minGridSpacing
558
559 minDeltaLat = 1.0d6
560 do lonIndex = 1, ni - 1
561 do latIndex = latIndexBeg, latIndexEnd - 1
562
563 deltaLat1 = abs(hco%lat2d_4(lonIndex, latIndex) - hco%lat2d_4(lonIndex , latIndex + 1))
564 deltaLat2 = abs(hco%lat2d_4(lonIndex, latIndex) - hco%lat2d_4(lonIndex + 1, latIndex ))
565 deltaLat3 = abs(hco%lat2d_4(lonIndex, latIndex) - hco%lat2d_4(lonIndex + 1, latIndex + 1))
566
567 deltaLat = max(deltaLat1, deltaLat2, deltaLat3)
568 if (deltaLat < minDeltaLat) minDeltaLat = deltaLat
569
570 end do
571 end do
572
573 minDeltaLon = 1.0d6
574 do lonIndex = 1, ni - 1
575 do latIndex = latIndexBeg, latIndexEnd - 1
576
577 if(abs(hco%lat2d_4(lonIndex, latIndex)) * MPC_DEGREES_PER_RADIAN_R8 < absMaxLat) then
578
579 deltaLon1 = abs(hco%lon2d_4(lonIndex, latIndex) - hco%lon2d_4(lonIndex , latIndex + 1))
580 deltaLon2 = abs(hco%lon2d_4(lonIndex, latIndex) - hco%lon2d_4(lonIndex + 1, latIndex ))
581 deltaLon3 = abs(hco%lon2d_4(lonIndex, latIndex) - hco%lon2d_4(lonIndex + 1, latIndex + 1))
582
583 if (deltaLon1 > MPC_PI_R8) deltaLon1 = deltaLon1 - 2.0d0 * MPC_PI_R8
584 deltaLon1 = abs(deltaLon1 * cos(hco%lat2d_4(lonIndex,latIndex)))
585 if (deltaLon2 > MPC_PI_R8) deltaLon2 = deltaLon2 - 2.0d0 * MPC_PI_R8
586 deltaLon2 = abs(deltaLon2 * cos(hco%lat2d_4(lonIndex,latIndex)))
587 if (deltaLon3 > MPC_PI_R8) deltaLon3 = deltaLon3 - 2.0d0 * MPC_PI_R8
588 deltaLon3 = abs(deltaLon3 * cos(hco%lat2d_4(lonIndex,latIndex)))
589
590 deltaLon = max(deltaLon1, deltaLon2, deltaLon3)
591 if (deltaLon < minDeltaLon) minDeltaLon = deltaLon
592
593 end if
594
595 end do
596 end do
597
598 minGridSpacing = ec_ra * sqrt(2.0d0) * min(minDeltaLon, minDeltaLat)
599
600 if (mmpi_myid == 0 .and. minGridSpacing /= minGridSpacingPrevious) then
601 minGridSpacingPrevious = minGridSpacing
602 write(*,*) 'hco_setupFromFile: minDeltaLat=', minDeltaLat * MPC_DEGREES_PER_RADIAN_R8, ' deg'
603 write(*,*) 'hco_setupFromFile: minDeltaLon=', minDeltaLon * MPC_DEGREES_PER_RADIAN_R8, ' deg'
604 write(*,*) 'hco_setupFromFile: minGridSpacing=', minGridSpacing, ' m'
605 end if
606
607 if (minGridSpacing > 1.0d6) then
608 call utl_abort('hco_setupFromFile: minGridSpacing is greater than 1000 km.')
609 end if
610
611 hco%minGridSpacing = minGridSpacing
612
613 !
614 !- 4. Close the input file
615 !
616 ier = fstfrm(iu_template)
617 ier = fclos (iu_template)
618
619 end subroutine hco_SetupFromFile
620
621 !--------------------------------------------------------------------------
622 ! Global_or_lam
623 !--------------------------------------------------------------------------
624 subroutine global_or_lam(global, lon, ni)
625 !
626 !:Purpose: to decide if a given grid is global or lam from input longitude array
627 !
628 implicit none
629
630 ! Arguments:
631 integer, intent(in) :: ni
632 real(8), intent(in) :: lon(ni)
633 logical, intent(out) :: global
634
635 ! Locals:
636 real(8) :: dx, next_lon
637
638 dx = lon(2) - lon(1)
639 next_lon = lon(ni) + 1.5d0 * dx
640
641 write(*,*)
642 write(*,*) 'dx = ',dx
643 write(*,*) 'lon(ni) = ',lon(ni)
644 write(*,*) 'next_lon = ',next_lon
645 write(*,*) 'lon(1) = ',lon(1)
646
647 if (next_lon - lon(1) > 360.0d0 .or. &
648 next_lon - lon(1) < 3.0*dx) then
649
650 global = .true.
651 if (lon(1) == lon(ni)) then
652 write(*,*)
653 write(*,*) ' *** Global Grid where i = ni (repetition) '
654 else
655 write(*,*)
656 write(*,*) ' *** Global Grid where i /= ni '
657 end if
658
659 else
660
661 global = .false.
662 write(*,*)
663 write(*,*) ' *** Limited-Area Grid '
664
665 end if
666
667 end subroutine global_or_lam
668
669 !--------------------------------------------------------------------------
670 ! mpiBcast
671 !--------------------------------------------------------------------------
672 subroutine hco_mpiBcast(hco)
673 !
674 !:Purpose: to broadcast hco strucure from MPI task 0 to other tasks
675 !
676 implicit none
677
678 ! Arguments:
679 type(struct_hco), pointer, intent(inout) :: hco
680
681 ! Locals:
682 integer :: ierr
683 integer, external :: ezqkdef
684
685 write(*,*) 'hco_mpiBcast: starting'
686
687 if (mmpi_myid > 0) then
688 if(.not.associated(hco)) then
689 allocate(hco)
690 else
691 call utl_abort('hco_mpiBcast: hco must be nullified for mpi task id > 0')
692 end if
693 end if
694
695 call rpn_comm_bcastc(hco%gridname, len(hco%gridname), 'MPI_CHARACTER', 0, 'GRID', ierr)
696 call rpn_comm_bcast(hco%initialized, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
697 call rpn_comm_bcast(hco%ni, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
698 call rpn_comm_bcast(hco%nj, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
699 call rpn_comm_bcastc(hco%grtyp, len(hco%grtyp), 'MPI_CHARACTER', 0, 'GRID', ierr)
700 call rpn_comm_bcastc(hco%grtypTicTac, len(hco%grtypTicTac), 'MPI_CHARACTER', 0, 'GRID', ierr)
701 call rpn_comm_bcast(hco%ig1, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
702 call rpn_comm_bcast(hco%ig2, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
703 call rpn_comm_bcast(hco%ig3, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
704 call rpn_comm_bcast(hco%ig4, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
705 if (mmpi_myid > 0) then
706 allocate(hco%lat(hco%nj))
707 allocate(hco%lon(hco%ni))
708 allocate(hco%lat2d_4(hco%ni,hco%nj))
709 allocate(hco%lon2d_4(hco%ni,hco%nj))
710 end if
711 call rpn_comm_bcast(hco%lat, size(hco%lat), 'MPI_REAL8', 0, 'GRID', ierr)
712 call rpn_comm_bcast(hco%lon, size(hco%lon), 'MPI_REAL8', 0, 'GRID', ierr)
713 call rpn_comm_bcast(hco%lat2d_4, size(hco%lat2d_4), 'MPI_REAL4', 0, 'GRID', ierr)
714 call rpn_comm_bcast(hco%lon2d_4, size(hco%lon2d_4), 'MPI_REAL4', 0, 'GRID', ierr)
715 call rpn_comm_bcast(hco%dlat, 1, 'MPI_REAL8', 0, 'GRID', ierr)
716 call rpn_comm_bcast(hco%dlon, 1, 'MPI_REAL8', 0, 'GRID', ierr)
717 call rpn_comm_bcast(hco%global, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
718 call rpn_comm_bcast(hco%rotated, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
719 call rpn_comm_bcast(hco%xlat1, 1, 'MPI_REAL8', 0, 'GRID', ierr)
720 call rpn_comm_bcast(hco%xlon1, 1, 'MPI_REAL8', 0, 'GRID', ierr)
721 call rpn_comm_bcast(hco%xlat2, 1, 'MPI_REAL8', 0, 'GRID', ierr)
722 call rpn_comm_bcast(hco%xlon2, 1, 'MPI_REAL8', 0, 'GRID', ierr)
723 call rpn_comm_bcast(hco%xlat1_yan, 1, 'MPI_REAL8', 0, 'GRID', ierr)
724 call rpn_comm_bcast(hco%xlon1_yan, 1, 'MPI_REAL8', 0, 'GRID', ierr)
725 call rpn_comm_bcast(hco%xlat2_yan, 1, 'MPI_REAL8', 0, 'GRID', ierr)
726 call rpn_comm_bcast(hco%xlon2_yan, 1, 'MPI_REAL8', 0, 'GRID', ierr)
727 if (hco%grtyp == 'U') then
728 if (mmpi_myid > 0) then
729 allocate(hco%tictacU(5 + 2 * (10 + hco%ni + hco%nj/2)))
730 end if
731 call rpn_comm_bcast(hco%tictacU, size(hco%tictacU), 'MPI_REAL4', 0, 'GRID', ierr)
732 end if
733
734 if (mmpi_myid > 0) then
735 if (hco%grtyp == 'G' .or. hco%grtyp == 'B') then
736 hco%EZscintID = ezqkdef(hco%ni, hco%nj, hco%grtyp, hco%ig1, hco%ig2, hco%ig3, hco%ig4, 0)
737 else
738 ! special treatment since EZscintID not properly communicated: keep as is
739 write(*,*) 'hco_mpiBcast: Warning! Grid ID for EZSCINT not communicated for grtyp = ', hco%grtyp
740 write(*,*) 'hco_mpiBcast: Warning! Grid ID for EZSCINT equals = ', hco%EZscintID
741 end if
742 end if
743
744 write(*,*) 'hco_mpiBcast: done'
745
746 end subroutine hco_mpiBcast
747
748 !--------------------------------------------------------------------------
749 ! Equal ?
750 !--------------------------------------------------------------------------
751 function hco_equal(hco1, hco2) result(equal)
752 !
753 !:Purpose: to check if two given hco strucures are equal or not
754 !
755 implicit none
756
757 ! Arguments:
758 type(struct_hco), pointer, intent(in) :: hco1
759 type(struct_hco), pointer, intent(in) :: hco2
760 ! Result:
761 logical :: equal
762
763 equal = .true.
764 equal = equal .and. (hco1%ni == hco2%ni)
765 equal = equal .and. (hco1%nj == hco2%nj)
766 if (.not. equal) then
767 write(*,*) 'hco_equal: dimensions not equal ', hco1%ni, hco1%nj, hco2%ni, hco2%nj
768 return
769 end if
770
771 equal = equal .and. (hco1%grtyp == hco2%grtyp)
772 if (.not. equal) then
773 write(*,*) 'hco_equal: grid type not equal'
774 return
775 end if
776
777 equal = equal .and. (hco1%dlat == hco2%dlat)
778 equal = equal .and. (hco1%dlon == hco2%dlon)
779 if (.not. equal) then
780 write(*,*) 'hco_equal: grid spacing not equal'
781 return
782 end if
783
784 if(hco1%grtyp == 'G' .or. hco1%grtyp == 'B') then
785 equal = equal .and. (hco1%ig2 == hco2%ig2)
786 if (.not. equal) then
787 write(*,*) 'hco_equal: Gaussian grid ig2 not equal'
788 return
789 end if
790 end if
791
792 equal = equal .and. (hco1%rotated .eqv. hco2%rotated)
793 equal = equal .and. (hco1%xlat1 == hco2%xlat1)
794 equal = equal .and. (hco1%xlon1 == hco2%xlon1)
795 equal = equal .and. (hco1%xlat2 == hco2%xlat2)
796 equal = equal .and. (hco1%xlon2 == hco2%xlon2)
797 equal = equal .and. (hco1%xlat1_yan == hco2%xlat1_yan)
798 equal = equal .and. (hco1%xlon1_yan == hco2%xlon1_yan)
799 equal = equal .and. (hco1%xlat2_yan == hco2%xlat2_yan)
800 equal = equal .and. (hco1%xlon2_yan == hco2%xlon2_yan)
801 if (.not. equal) then
802 write(*,*) 'hco_equal: rotation not equal: ', hco1%rotated, hco2%rotated
803 write(*,*) 'hco_equal: xlat1: ', hco1%xlat1, hco2%xlat1
804 write(*,*) 'hco_equal: xlat2: ', hco1%xlat2, hco2%xlat2
805 write(*,*) 'hco_equal: xlon1: ', hco1%xlon1, hco2%xlon1
806 write(*,*) 'hco_equal: xlon2: ', hco1%xlon2, hco2%xlon2
807 write(*,*) 'hco_equal: xlat1_yan: ', hco1%xlat1_yan, hco2%xlat1_yan
808 write(*,*) 'hco_equal: xlat2_yan: ', hco1%xlat2_yan, hco2%xlat2_yan
809 write(*,*) 'hco_equal: xlon1_yan: ', hco1%xlon1_yan, hco2%xlon1_yan
810 write(*,*) 'hco_equal: xlon2_yan: ', hco1%xlon2_yan, hco2%xlon2_yan
811 return
812 end if
813
814 equal = equal .and. all(hco1%lat(:) == hco2%lat(:))
815 equal = equal .and. all(hco1%lon(:) == hco2%lon(:))
816 if (.not. equal) then
817 write(*,*) 'hco_equal: lat/lon not equal'
818 return
819 end if
820
821 end function hco_equal
822
823 !--------------------------------------------------------------------------
824 ! hco_deallocate
825 !--------------------------------------------------------------------------
826 subroutine hco_deallocate(hco)
827 implicit none
828
829 ! Arguments:
830 type(struct_hco), pointer, intent(inout) :: hco
831
832 if (allocated(hco % lat)) deallocate(hco % lat)
833 if (allocated(hco % lon)) deallocate(hco % lon)
834 if (allocated(hco % lat2d_4)) deallocate(hco % lat2d_4)
835 if (allocated(hco % lon2d_4)) deallocate(hco % lon2d_4)
836 if (allocated(hco % tictacU)) deallocate(hco % tictacU)
837
838 if (associated(hco)) deallocate(hco)
839 nullify(hco)
840
841 end subroutine hco_deallocate
842
843 !--------------------------------------------------------------------------
844 ! grid_mask
845 !--------------------------------------------------------------------------
846 subroutine grid_mask(F_mask_8,dx,dy,xg,yg,ni,nj)
847 !
848 !:Purpose: 1) Find out where YIN lat lon points are in (YAN) grid with call to smat.
849 ! 2) If they are not outside of Yin grid, put area to zero for those points.
850 !
851 ! Author Qaddouri
852 !
853 implicit none
854
855 ! Arguments:
856 integer, intent(in) :: Ni
857 integer, intent(in) :: Nj
858 real(8) , intent(out) :: F_mask_8(Ni,Nj)
859 real(8), intent(in) :: dx
860 real(8), intent(in) :: dy
861 real(4), intent(in) :: xg(ni)
862 real(4), intent(in) :: yg(nj)
863
864 ! Locals:
865 integer :: lonIndex,latIndex,np_subd
866 real(8) :: poids(ni,nj),x_a_4,y_a_4,sp,sf,sp1,sf1
867 real(4) :: area_4(ni,nj)
868
869 np_subd = 4*ni
870
871 sp = 0.d0
872 sf = 0.d0
873
874 do latIndex = 1, nj
875 y_a_4 = yg(latIndex)
876 do lonIndex = 1, ni
877
878 x_a_4 = xg(lonIndex)-acos(-1.d0)
879
880 area_4(lonIndex,latIndex) = dx*dy*cos(yg(latIndex))
881 poids (lonIndex,latIndex) = yyg_weight (x_a_4,y_a_4,dx,dy,np_subd)
882
883 !Check if poids <0
884 if (poids(lonIndex,latIndex)*(1.d0-poids(lonIndex,latIndex)) > 0.d0) then
885 sp = sp + poids(lonIndex,latIndex)*area_4(lonIndex,latIndex)
886 else if (abs(poids(lonIndex,latIndex)-1.d0) < 1.d-14) then
887 sf = sf + poids(lonIndex,latIndex)*area_4(lonIndex,latIndex)
888 end if
889
890 end do
891 end do
892
893 !Correct and scale poids
894 !-----------------------
895 sp1 = 0.d0
896 sf1 = 0.d0
897
898 do latIndex = 1, nj
899 do lonIndex = 1, ni
900
901 x_a_4 = poids(lonIndex,latIndex)*(2.d0*acos(-1.d0) - sf)/sp
902
903 if (poids(lonIndex,latIndex)*(1.d0-poids(lonIndex,latIndex)) > 0.d0) then
904 poids(lonIndex,latIndex) = min(1.0d0, x_a_4)
905 end if
906 if (poids(lonIndex,latIndex)*(1.0-poids(lonIndex,latIndex)) > 0.d0) then
907 sp1 = sp1 + poids(lonIndex,latIndex)*area_4(lonIndex,latIndex)
908 else if (abs(poids(lonIndex,latIndex)-1.d0) < 1.d-14) then
909 sf1 = sf1 + poids(lonIndex,latIndex)*area_4(lonIndex,latIndex)
910 end if
911
912 end do
913 end do
914
915 !Correct
916 !-------
917 do latIndex = 1, nj
918 do lonIndex = 1, ni
919 x_a_4 = poids(lonIndex,latIndex)*(2.d0*acos(-1.d0) - sf1)/sp1
920
921 if (poids(lonIndex,latIndex)*(1.d0-poids(lonIndex,latIndex)) > 0.d0) then
922 poids(lonIndex,latIndex) = min(1.d0, x_a_4)
923 end if
924
925 end do
926 end do
927
928 F_mask_8 = 0.d0
929 do latIndex=1,nj
930 do lonIndex = 1,ni
931 F_mask_8(lonIndex,latIndex) = poids(lonIndex,latIndex)
932 end do
933 end do
934
935 end subroutine grid_mask
936
937 !--------------------------------------------------------------------------
938 ! inter_curve_boundary_yy
939 !--------------------------------------------------------------------------
940 subroutine inter_curve_boundary_yy (x,y,xi,yi,np)
941 !
942 !:Purpose: compute the intersections between a line and the panel
943 ! (yin or yang) boundary. The line passes through the panel
944 ! center point (0, 0) and the cell center point (x, y).
945 !
946 ! Author A. Qaddouri. October 2016
947 !
948 ! Note: this routine has been taken and adjusted from a routine
949 ! with the same name in the GEM model.
950 !
951 !:Arguments: input: (x, y): longitude, latitude of the cell
952 ! center point,
953 ! np: workspace,
954 ! output: (xi, yi): longitude, latitude of the
955 ! intersection point.
956 implicit none
957
958 ! Arguments:
959 real(8), intent(in) :: x
960 real(8), intent(in) :: y
961 real(8), intent(out) :: xi
962 real(8), intent(out) :: yi
963 integer, intent(in) :: np
964
965 ! Locals:
966 real(8) :: tol, pi, xmin, ymin, xb
967 real(8) :: xc, yc, s1, s2, x1, test, dxs
968 real(8) :: xp1, xp2, xr1, yr1, xr2, yr2
969 integer :: i
970
971 tol = 1.0d-16
972 pi = MPC_PI_R8
973 xmin = -3.d0*pi/4.d0
974 ymin = -pi/4.d0
975 xb = -0.5d0*pi
976
977 xc = x
978 yc = y
979
980 if (x > 0.d0) xc = - x
981 if (y > 0.d0) yc = - y
982
983 if (abs(xc) < tol) then
984 xi = xc
985 yi = ymin
986 else
987
988 s1 = yc / xc
989 s2 = ymin/(xb-xmin)
990
991 x1 = s2*xmin/(s2-s1)
992 if (x1 > xb) then
993 xi = ymin/s1
994 yi = ymin
995 else
996 test = -1.d0
997 dxs = -xb/(np-1)
998 i = 1
999 do while (test < 0.d0)
1000 xp1 = (i-1)*dxs + xb
1001 xp2 = (i)*dxs + xb
1002 xr1 = atan2(sin(ymin),-cos(ymin)*cos(xp1))
1003 yr1 = asin(cos(ymin)*sin(xp1))
1004 xr2 = atan2(sin(ymin),-cos(ymin)*cos(xp2))
1005 yr2 = asin(cos(ymin)*sin(xp2))
1006 s2 = (yr1-yr2)/(xr1-xr2)
1007 xi = (s2*xr2-yr2)/(s2-s1)
1008 yi = s1*xi
1009 test=(xi-xr1)*(xr2-xi)
1010 i = i+1
1011 end do
1012 end if
1013 end if
1014
1015 if (x > 0.d0) xi = - xi
1016 if (y > 0.d0) yi = - yi
1017
1018 end subroutine inter_curve_boundary_yy
1019
1020 !--------------------------------------------------------------------------
1021 ! hco_weight
1022 !--------------------------------------------------------------------------
1023 subroutine hco_weight(hco, weight)
1024 !
1025 !:Purpose: given the horizontal grid definition of the grid,
1026 ! return appropriate weights for individual points (avoiding
1027 ! double counting in the overlap regions in the case of a Yin-Yang grid).
1028 !
1029 !:Author: Abdessamad Qaddouri and Peter Houtekamer
1030 ! October 2016
1031 !
1032 !:Revision: imported code for the Yin-Yang grid on May 2021 from the EnKF library and
1033 ! combined with code for other grid types from the MIDAS library.
1034 !
1035 implicit none
1036
1037 ! Arguments:
1038 type(struct_hco), intent(in) :: hco ! structure with the specification of the horizontal grid
1039 real(8), intent(out) :: weight(:,:) ! weight to be given when computing a horizontal average
1040
1041 ! Locals:
1042 integer :: sindx
1043 integer :: ni,nj
1044 integer :: lonIndex,latIndex,lonIndexP1,latIndexP1
1045 real(8), allocatable :: F_mask_8(:,:), F_mask(:,:)
1046 real(8) :: deg2rad,dx,dy,sum_weight
1047 real(8) :: lon1,lon2,lon3,lat1,lat2,lat3
1048 real(4), allocatable :: xg(:),yg(:)
1049
1050 deg2rad= MPC_RADIANS_PER_DEGREE_R8
1051 sindx = 6
1052
1053 if (trim(hco%grtyp) == 'U') then ! case of a Yin-Yang grid
1054 write(*,*) 'compute weights for Yin_Yang grid'
1055 ni = nint(hco%tictacU(sindx))
1056 nj = nint(hco%tictacU(sindx+1))
1057 allocate (F_mask_8 (ni,nj))
1058 allocate (F_mask(ni,2*nj))
1059 allocate (xg(ni))
1060 allocate (yg(nj))
1061
1062 dx = hco%tictacU(sindx+10+1) -hco%tictacU(sindx+10)
1063 dy= hco%tictacU(sindx+10+ni+1)-hco%tictacU(sindx+10+ni)
1064 dx= deg2rad* dx
1065 dy= deg2rad* dy
1066 do lonIndex=1,ni
1067 xg(lonIndex)=deg2rad* hco%tictacU(sindx+10+lonIndex-1)
1068 end do
1069 do latIndex=1,nj
1070 yg(latIndex)= deg2rad*hco%tictacU(sindx+10+ni+latIndex-1)
1071 weight (:,latIndex) = cos(deg2rad* hco%tictacU(sindx+10+ni+latIndex-1))
1072 weight (:,nj+latIndex)= weight (:,latIndex)
1073 end do
1074 call grid_mask (F_mask_8,dx,dy,xg,yg,ni,nj)
1075 do latIndex=1,nj
1076 F_mask (:,latIndex) = F_mask_8(:,latIndex)
1077 F_mask (:,nj+latIndex)= F_mask (:,latIndex)
1078 end do
1079 do latIndex=1,nj*2
1080 weight(:,latIndex) = weight(:, latIndex) * F_mask(:, latIndex)
1081 end do
1082 sum_weight=sum(weight)
1083 weight = weight/sum_weight
1084 deallocate(F_mask_8)
1085 deallocate(F_mask)
1086 deallocate(xg)
1087 deallocate(yg)
1088 else
1089 write(*,*) 'compute weights for grid type: ',hco%grtyp
1090 do latIndex=1,hco%nj
1091 latIndexP1 = min(hco%nj,latIndex+1)
1092 do lonIndex=1,hco%ni
1093 lonIndexP1 = min(hco%ni,lonIndex+1)
1094 lon1 = hco%lon2d_4(lonIndex,latIndex)
1095 lon2 = hco%lon2d_4(lonIndexP1,latIndex)
1096 lon3 = hco%lon2d_4(lonIndex,latIndexP1)
1097 lat1 = hco%lat2d_4(lonIndex,latIndex)
1098 lat2 = hco%lat2d_4(lonIndexP1,latIndex)
1099 lat3 = hco%lat2d_4(lonIndex,latIndexP1)
1100 dx = phf_calcDistance(lat1, lon1, lat2, lon2)/1000.0D0
1101 dy = phf_calcDistance(lat1, lon1, lat3, lon3)/1000.0D0
1102 weight(lonIndex,latIndex) = dx * dy
1103 end do
1104 end do
1105 sum_weight=sum(weight)
1106 weight = weight/sum_weight
1107 end if
1108
1109 end subroutine hco_weight
1110
1111 !--------------------------------------------------------------------------
1112 ! yyg_weight
1113 !--------------------------------------------------------------------------
1114 real(8) function yyg_weight (x,y,dx,dy,np)
1115 !
1116 !:Purpose: Based on a Draft from Zerroukat (2013) - Evaluates weight
1117 ! for each cell, so that the overlap is computed once.
1118 !
1119 !:Author: Author Abdessamad Qaddouri -- Summer 2014
1120 !
1121 !:Note: this routine has been taken and adjusted from a routine
1122 ! with the same name in the GEM model.
1123 !
1124 implicit none
1125
1126 ! Arguments:
1127 real(8), intent(in) :: x ! x,y: (lon, lat) in radians of the cell center point.
1128 real(8), intent(in) :: y
1129 real(8), intent(in) :: dx ! horizontal grid resolution in radians
1130 real(8), intent(in) :: dy
1131 integer, intent(in) :: np ! working dimension
1132
1133 ! Locals:
1134 real(8) :: pi, xmin, xmax, ymin, ymax, xb1, xb2
1135 real(8) :: t1x, t2x, t1y, dcell, xi, yi, di, dp, df, d
1136
1137 pi = MPC_PI_R8
1138 xmin = -3.d0*pi/4.d0 ; xmax = 3.d0*pi/4.d0
1139 ymin = -pi/4.d0 ; ymax = pi/4.d0
1140 xb1 = -0.5d0*pi ; xb2 = 0.5d0*pi
1141
1142 t1x = (x-xmin)*(xmax-x)
1143 t2x = (x-xb1)*(xb2-x)
1144 t1y = (y-ymin)*(ymax-y)
1145
1146 if (t1x < 0.d0 .or. t1y < 0.d0) then
1147 yyg_weight = 0.0
1148 else if (t2x > 0.d0 .and. t1y > 0.d0) then
1149 yyg_weight = 1.d0
1150 else
1151 dcell = 0.5d0*dsqrt(dx**2 + dy**2)
1152
1153 call inter_curve_boundary_yy (x, y, xi, yi, np)
1154
1155 di = sqrt(xi**2 + yi**2)
1156 dp = sqrt(x**2 + y**2)
1157 df = dp - di
1158 d = min(max(-dcell,df),dcell)
1159 yyg_weight = 0.5d0*(1.d0 - (d/dcell))
1160 end if
1161
1162 end function yyg_weight
1163
1164 !--------------------------------------------------------------------------
1165 ! hco_setupYgrid
1166 !--------------------------------------------------------------------------
1167 subroutine hco_setupYgrid(hco, ni, nj)
1168 !
1169 !:Purpose: to initialize hco structure for a Y grid
1170 !
1171 implicit none
1172
1173 ! Arguments:
1174 type(struct_hco), pointer, intent(inout) :: hco
1175 integer, intent(in) :: ni
1176 integer, intent(in) :: nj
1177
1178 allocate(hco)
1179 if (mmpi_myid == 0) then
1180 hco%initialized = .true.
1181 hco%ni = ni
1182 hco%nj = nj
1183 hco%grtyp = 'Y'
1184 hco%grtypTicTac = 'L'
1185 if (allocated(hco%lat2d_4)) then
1186 deallocate(hco%lat2d_4)
1187 deallocate(hco%lon2d_4)
1188 end if
1189 allocate(hco%lat2d_4(ni, nj))
1190 allocate(hco%lon2d_4(ni, nj))
1191 hco%xlat1 = 0.d0
1192 hco%xlon1 = 0.d0
1193 hco%xlat2 = 1.d0
1194 hco%xlon2 = 1.d0
1195 end if
1196
1197 end subroutine hco_setupYgrid
1198
1199end module horizontalCoord_mod
