1module verticalCoord_mod
2 ! MODULE verticalCoord_mod (prefix='vco' category='7. Low-level data objects')
3 !
4 !:Purpose: Derived type and procedures related to the vertical levels.
5 ! The derived type includes a pointer to the associated VGRID
6 ! descriptor.
7 !
8 use midasMpi_mod
9 use Vgrid_Descriptors
10 use varNameList_mod
11 use utilities_mod
12
13 implicit none
14
15 private
16
17 ! public derived type
18 public :: struct_vco
19 ! public variables
20 public :: vco_ip1_other, vco_maxNumLevels
21 ! public procedures
22 public :: vco_setupFromFile, vco_getNumLev, vco_equal, vco_deallocate, vco_mpiBcast
23 public :: vco_subsetOrNot, vco_levelMatchingList
24
25 integer, parameter :: maxNumOtherLevels = 20
26 integer :: vco_ip1_other(maxNumOtherLevels)
27
28 integer, parameter :: vco_maxNumLevels = 200
29
30 type struct_vco
31 logical :: initialized=.false.
32 integer :: Vcode = -1
33 integer :: nlev_T = 0
34 integer :: nlev_M = 0
35 integer :: nlev_Other(vnl_numvarmaxOther) = 0
36 integer :: nlev_Depth = 0
37 integer :: ip1_seaLevel
38 integer :: ip1_sfc ! ip1 value for the surface (hybrid = 1)
39 integer :: ip1_T_2m ! ip1 value for the 2m thermodynamic level
40 integer :: ip1_M_10m ! ip1 value for the 10m momentum level
41 integer, pointer :: ip1_T(:) => null()
42 integer, pointer :: ip1_M(:) => null() ! encoded IP1 levels (Thermo/Moment)
43 integer, pointer :: ip1_depth(:) => null() ! encoded IP1 levels (Ocean depth levels)
44 type(vgrid_descriptor) :: vgrid
45 logical :: vgridPresent
46 real(8), pointer :: depths(:) => null()
47 end type struct_vco
48
49contains
50
51 !--------------------------------------------------------------------------
52 ! vco_allocateIp1
53 !--------------------------------------------------------------------------
54 subroutine vco_allocateIp1(vco)
55 !
56 !:Purpose: Allocate the ip1 arrays of a vertical coordinate object.
57 !
58 implicit none
59
60 ! Arguments:
61 type(struct_vco), pointer, intent(inout) :: vco ! Vertical coordinate object
62
63 ! Locals:
64 integer :: numLev
65
66 numLev = vco_getNumLev(vco,'MM')
67 if (numLev > 0) allocate (vco%ip1_M(numLev))
68
69 numLev = vco_getNumLev(vco,'TH')
70 if (numLev > 0) allocate (vco%ip1_T(numLev))
71
72 end subroutine vco_allocateIp1
73
74 !--------------------------------------------------------------------------
75 ! vco_SetupFromFile
76 !--------------------------------------------------------------------------
77 subroutine vco_SetupFromFile(vco, templatefile, etiket_opt, beSilent_opt)
78 !
79 !:Purpose: Initialize vertical coordinate object with information from
80 ! a standard file.
81 !
82 implicit none
83
84 ! Arguments:
85 type(struct_vco),pointer, intent(inout) :: vco ! Vertical coordinate object
86 character(len=*), intent(in) :: templatefile ! Template file
87 character(len=*), optional, intent(in) :: etiket_opt ! Optional argument etiket
88 logical, optional, intent(in) :: beSilent_opt ! Optional argument beSilent
89
90 ! Locals:
91 logical :: beSilent
92 character(len=12) :: etiket
93 integer :: nultemplate,ierr
94 integer, parameter :: maxNumRecords = 1000
95 integer :: recordIndex, numRecords, ikeys(maxNumRecords)
96 integer :: fnom,fstouv,fstfrm,fclos,fstprm,fstinl
97 integer :: ip1_sfc
98 character(len=10) :: blk_S
99 logical :: fileExists, atmFieldFound, sfcFieldFound, oceanFieldFound
100 integer :: IP1kind
101 real :: vertCoordValue
102 integer :: ideet, inpas, dateStamp_origin, ini, inj, ink, inbits, idatyp
103 integer :: ip1, ip2, ip3, ig1, ig2, ig3, ig4, iswa, ilng, idltf, iubc
104 integer :: iextra1, iextra2, iextra3
105 character(len=2) :: typvar
106 character(len=4) :: nomvar
107 character(len=1) :: grtyp
108
109 if ( associated(vco) ) then
110 call utl_abort('vco_setupFromFile: the supplied vco pointer is not null!')
111 end if
112
113 allocate(vco)
114 call convip(vco%ip1_seaLevel, 0.0, 0, 2, blk_s, .false.)
115
116 if (present(etiket_opt)) then
117 etiket = etiket_opt
118 else
119 etiket = ' '
120 end if
121
122 if (present(beSilent_opt)) then
123 beSilent = beSilent_opt
124 else
125 beSilent = .false.
126 end if
127
128 ! Check if template file exists
129 if (mmpi_myid == 0 .and. .not. beSilent) then
130 write(*,*) 'vco_setupFromFile: Template File = ', trim(templatefile)
131 end if
132 inquire(file=templatefile,exist=fileExists)
133 if ( .not. fileExists )then
134 write(*,*) 'vco_setupFromFile: Template File = ', trim(templatefile)
135 call utl_abort('vco_setupFromFile: CANNOT FIND TEMPLATE FILE!')
136 end if
137
138 ! First priority, check if vgrid descriptor record present - if so, atmospheric fields
139 atmFieldFound = utl_varNamePresentInFile('!!',fileName_opt=trim(templatefile))
140
141 ! If not atmospheric field, we need to examine data records in the template file
142 if (.not. atmFieldFound) then
143
144 if (.not. beSilent) then
145 write(*,*) 'vco_setupFromFile: No vgrid descriptor found, examine data records in template file'
146 end if
147
148 ! open the template file
149 nultemplate=0
150 ierr=fnom(nultemplate,templatefile,'RND+OLD+R/O',0)
151 if ( ierr == 0 ) then
152 ierr = fstouv(nultemplate,'RND+OLD')
153 else
154 write(*,*) 'vco_setupFromFile: Template File = ', trim(templatefile)
155 call utl_abort('vco_setupFromFile: CANNOT OPEN TEMPLATE FILE!')
156 end if
157
158 ierr = fstinl(nultemplate,ini,inj,ink,-1,etiket,-1,-1,-1,' ', &
159 ' ',ikeys,numRecords,maxNumRecords)
160 if (ikeys(1) <= 0) then
161 call utl_abort('vco_setupFromFile: Could not find any records in the supplied file')
162 end if
163 if (.not. beSilent) then
164 write(*,*) 'vco_setupFromFile: number of records found = ', numRecords
165 end if
166
167 ! check records for ocean or surface fields
168 sfcFieldFound = .false.
169 oceanFieldFound = .false.
170 record_loop: do recordIndex = 1, numRecords
171 ierr = fstprm(ikeys(recordIndex), dateStamp_origin, ideet, inpas, ini, inj, &
172 ink, inbits, idatyp, ip1, ip2, ip3, &
173 typvar, nomvar, etiket, grtyp, ig1, ig2, ig3, ig4, &
174 iswa, ilng, idltf, iubc, iextra1, iextra2, iextra3)
175
176 ! ignore any variables not present in varnamelist_mod
177 if (.not. vnl_varnameIsValid(trim(nomvar))) cycle record_loop
178
179 ! check for record with ocean data
180 call convip(ip1, vertCoordValue, Ip1Kind, -1, blk_s, .false.)
181 if (Ip1Kind == 0 .and. &
182 vnl_varKindFromVarname(trim(nomvar)) == 'OC') then
183 oceanFieldFound = .true.
184 exit record_loop
185 end if
186
187 ! check for record with surface data (and that are not ocean variables)
188 call convip(ip1_sfc, 1.0, 5, 2, blk_s, .false.)
189 if (ip1 == 0 .or. ip1 == ip1_sfc .or. ip1 == vco%ip1_seaLevel .or. Ip1Kind == 3) then
190 sfcFieldFound = .true.
191 exit record_loop
192 end if
193
194 ! something else was found, abort (not sure abort is necessary)
195 write(*,*) 'vco_setupFromFile: found record with unknown vertical coordinate'
196 write(*,*) 'varName = ', trim(nomvar), ' typvar = ', typvar, ' ip1 = ', ip1
197 call utl_abort('vco_setupFromFile: found a non-surface field')
198
199 end do record_loop
200
201 ierr = fstfrm(nultemplate)
202 ierr = fclos (nultemplate)
203
204 end if
205
206 ! call appropriate setup routine based on what was found in template file
207 if (atmFieldFound) then
208 call vco_setupAtmFromFile(vco, templatefile, etiket, beSilent)
209 else if (oceanFieldFound) then
210 call vco_setupOceanFromFile(vco, templatefile, etiket, beSilent)
211 else if (sfcFieldFound) then
212 call vco_setupSfcFromFile(vco, beSilent)
213 else
214 call utl_abort('vco_setupFromFile: could not setup vco from template file')
215 end if
216
217 end subroutine vco_setupFromFile
218
219 !--------------------------------------------------------------------------
220 ! vco_setupAtmFromFile
221 !--------------------------------------------------------------------------
222 subroutine vco_setupAtmFromFile(vco, templatefile, etiket, beSilent)
223 !
224 !:Purpose: Initialize vertical coordinate object with information from
225 ! a standard file. Use vgrid descriptor for atmospheric fields.
226 !
227 implicit none
228
229 ! Arguments:
230 type(struct_vco), pointer, intent(inout) :: vco ! Vertical coordinate object
231 character(len=*), intent(in) :: templatefile ! Template file
232 character(len=*), intent(in) :: etiket
233 logical, intent(in) :: beSilent
234
235 ! Locals:
236 integer :: Vcode, jlev, nlevMatched, stat, nultemplate, ierr, ikey
237 integer :: fnom, fstouv, fstfrm, fclos, fstinf
238 integer :: vgd_nlev_M, vgd_nlev_T, ip1_sfc
239 integer :: ni, nj, nk, varListIndex, IP1kind
240 integer, pointer :: vgd_ip1_M(:), vgd_ip1_T(:)
241 logical :: ip1_found
242 character(len=10) :: blk_S
243 character(len=4) :: nomvar_T, nomvar_M, nomvar_Other
244 character(len=10) :: IP1string
245 real :: otherVertCoordValue
246
247 ! Open the template file
248 nultemplate = 0
249 ierr = fnom(nultemplate,templatefile,'RND+OLD+R/O',0)
250 if (ierr == 0) then
251 ierr = fstouv(nultemplate,'RND+OLD')
252 else
253 write(*,*) 'vco_setupFromFile: Template File = ', trim(templatefile)
254 call utl_abort('vco_setupAtmFromFile: CANNOT OPEN TEMPLATE FILE!')
255 end if
256
257 ! Try creating vgrid descriptor
258 stat = vgd_new(vco%vgrid, unit = nultemplate, format = "fst", ip1 = -1, ip2 = -1)
259 if (stat == VGD_OK) then
260 vco%vgridPresent = .true.
261 else
262 call utl_abort('vco_setupAtmFromFile: !! record exists, but not able to create descriptor object')
263 end if
264
265 ! Print out vertical structure
266 if (mmpi_myid == 0 .and. .not. beSilent) then
267 call flush(6) ! possibly needed so vgd_print output appears correctly in listing
268 stat = vgd_print(vco%vgrid)
269 if (stat /= VGD_OK)then
270 call utl_abort('vco_setupAtmFromFile: ERROR with vgd_print')
271 end if
272 end if
273
274 ! Get version of the vertical coordinate
275 stat = vgd_get(vco%vgrid, key = 'ig_1 - vertical coord code', value = Vcode)
276 if ( stat /= VGD_OK ) then
277 call utl_abort('vco_setupAtmFromFile: problem with vgd_get: key= ig_1 - vertical coord code')
278 end if
279 if (Vcode /= 5002 .and. Vcode /= 5005 .and. Vcode /= 5100 .and. Vcode /= 21001) then
280 call utl_abort('vco_setupAtmFromFile: Invalid Vcode. Currently only 5002, 5005, 5100 and 21001 supported.')
281 end if
282 vco%Vcode = Vcode
283
284 ! Get vgrid values for ip1
285 stat = vgd_get(vco%vgrid, key='vipm - vertical levels (m)', value = vgd_ip1_m)
286 stat = vgd_get(vco%vgrid, key='vipt - vertical ip1 levels (t)', value = vgd_ip1_t)
287
288 vgd_nlev_M = size(vgd_ip1_M)
289 vgd_nlev_T = size(vgd_ip1_T)
290
291 ! Set the number of vertical levels and allocate vco arrays
292 vco%nlev_T = 0
293 nomvar_T = 'TH '
294 do jlev = 1, vgd_nlev_T
295 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vgd_ip1_T(jlev), -1, -1, ' ', nomvar_T)
296 if (ikey > 0) vco%nlev_T = vco%nlev_T + 1
297 end do
298 if (vco%nlev_T == 0) then
299 if (mmpi_myid == 0 .and. .not. beSilent) then
300 write(*,*) 'vco_setupAtmFromFile: TH not found looking for TT to get nlev_T'
301 end if
302 nomvar_T = 'TT '
303 do jlev = 1, vgd_nlev_T
304 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vgd_ip1_T(jlev), -1, -1, ' ', nomvar_T)
305 if (ikey > 0) vco%nlev_T = vco%nlev_T + 1
306 end do
307 end if
308 if (vco%nlev_T == 0 .and. .not. beSilent) then
309 write(*,*)
310 write(*,*) 'vco_setupAtmFromFile: Could not find a valid thermodynamic variable in the template file!'
311 else if (vco%nlev_T > vco_maxNumLevels) then
312 write(*,*)
313 write(*,*) 'nlev_T = ',vco%nlev_T
314 write(*,*) 'vco_maxNumLevels = ',vco_maxNumLevels
315 call utl_abort('vco_setupAtmFromFile: nlev_T is greater than vco_maxNumLevels!')
316 end if
317
318 vco%nlev_M = 0
319 nomvar_M = 'MM '
320 do jlev = 1, vgd_nlev_M
321 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vgd_ip1_M(jlev), -1, -1, ' ', nomvar_M)
322 if (ikey > 0) vco%nlev_M = vco%nlev_M + 1
323 end do
324 if (vco%nlev_M == 0) then
325 if (mmpi_myid == 0 .and. .not. beSilent) then
326 write(*,*) 'vco_setupAtmFromFile: MM not found looking for UU to get nlev_M'
327 end if
328 nomvar_M = 'UU '
329 do jlev = 1, vgd_nlev_M
330 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vgd_ip1_M(jlev), -1, -1, ' ', nomvar_M)
331 if (ikey > 0) vco%nlev_M = vco%nlev_M + 1
332 end do
333 end if
334 if (vco%nlev_M == 0) then
335 if (mmpi_myid == 0 .and. .not. beSilent) then
336 write(*,*) 'vco_setupAtmFromFile: UU not found looking for PP to get nlev_M'
337 end if
338 nomvar_M = 'PP '
339 do jlev = 1, vgd_nlev_M
340 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vgd_ip1_M(jlev), -1, -1, ' ', nomvar_M)
341 if (ikey > 0) vco%nlev_M = vco%nlev_M + 1
342 end do
343 end if
344 if (vco%nlev_M == 0 .and. .not. beSilent) then
345 write(*,*)
346 write(*,*) 'vco_setupAtmFromFile: Could not find a valid momentum variable in the template file!'
347 else if (vco%nlev_M > vco_maxNumLevels) then
348 write(*,*)
349 write(*,*) 'nlev_M = ',vco%nlev_M
350 write(*,*) 'vco_maxNumLevels = ',vco_maxNumLevels
351 call utl_abort('vco_setupAtmFromFile: nlev_M is greater than vco_maxNumLevels!')
352 end if
353
354 do jlev = 1, maxNumOtherLevels
355 otherVertCoordValue = real(jlev)
356 IP1kind = 3
357 call convip_plus(vco_ip1_other(jlev), otherVertCoordValue, IP1kind, +2, IP1string, .false.)
358 end do
359 vco%nlev_Other(:) = 0
360 do varListIndex = 1, vnl_numvarmaxOther
361 nomvar_Other = vnl_varNameListOther(varListIndex)
362 do jlev = 1, maxNumOtherLevels
363 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vco_ip1_other(jlev), -1, -1, ' ', nomvar_Other)
364 if (ikey > 0) vco%nlev_Other(varListIndex) = vco%nlev_Other(varListIndex) + 1
365 end do
366 if (mmpi_myid == 0 .and. .not. beSilent) then
367 if (vco%nlev_Other(varListIndex) == 0) then
368 write(*,*)
369 write(*,*) 'vco_setupAtmFromFile: Found no levels in template file for OTHER type variable ', nomvar_Other
370 else
371 write(*,*) 'vco_setupAtmFromFile: Found ', vco%nlev_Other(varListIndex), &
372 ' levels in template file for OTHER type variable ', nomvar_Other
373 end if
374 end if
375 end do
376
377 if (vco%nlev_M == 0 .and. vco%nlev_T == 0) then
378 call utl_abort('vco_setupAtmFromFile: they were no valid momentum and thermodynamic variables in the template file!')
379 end if
380
381 if (mmpi_myid == 0 .and. .not.beSilent) then
382 write(*,*) 'vco_setupAtmFromFile: nlev_M, nlev_T=',vco%nlev_M,vco%nlev_T
383 end if
384
385 call vco_allocateIp1(vco)
386
387 ! Match up ip1 values from file and vgrid for momentum levels
388 nlevMatched = 0
389 do jlev = 1, vgd_nlev_M
390 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vgd_ip1_M(jlev), -1, -1, ' ', nomvar_M)
391 if (ikey > 0) then
392 nlevMatched = nlevMatched + 1
393 if (nlevMatched > vco%nlev_M) then
394 call utl_abort('vco_setupAtmFromFile: Problem with consistency between vgrid descriptor and template file (momentum)')
395 end if
396 vco%ip1_M(nlevMatched) = vgd_ip1_M(jlev)
397 end if
398 end do
399 if (nlevMatched /= vco%nlev_M) then
400 write(*,*) 'vco_setupAtmFromFile: nlevMatched = ', nlevMatched, ', nlev_M = ', vco%nlev_M
401 call utl_abort('vco_setupAtmFromFile: Problem with consistency between vgrid descriptor and template file (momentum)')
402 end if
403
404 ! Match up ip1 values from file and vgrid for thermo levels
405 nlevMatched = 0
406 do jlev = 1, vgd_nlev_T
407 ikey = fstinf(nultemplate, ni, nj, nk, -1 ,etiket, vgd_ip1_T(jlev), -1, -1, ' ', nomvar_T)
408 if (ikey > 0) then
409 nlevMatched = nlevMatched + 1
410 if (nlevMatched > vco%nlev_T) then
411 call utl_abort('vco_setupAtmFromFile: Problem with consistency between vgrid descriptor and template file (thermo)')
412 end if
413 vco%ip1_T(nlevMatched) = vgd_ip1_T(jlev)
414 end if
415 end do
416 if (nlevMatched /= vco%nlev_T) then
417 write(*,*) 'vco_setupAtmFromFile: nlevMatched = ', nlevMatched, ', nlev_T = ', vco%nlev_T
418 call utl_abort('vco_setupAtmFromFile: Problem with consistency between vgrid descriptor and template file (thermo)')
419 end if
420
421 ! determine IP1 of sfc (hyb=1.0)
422 if (Vcode == 5002 .or. Vcode == 5005 .or. Vcode == 5100) then
423 call convip(ip1_sfc, 1.0, 5, 2, blk_s, .false.)
424 else if (Vcode == 21001) then
425 call convip(ip1_sfc, 0.0, 21, 2, blk_s, .false.)
426 end if
427 ip1_found = .false.
428 do jlev = 1, vgd_nlev_T
429 if (ip1_sfc == vgd_ip1_T(jlev)) then
430 ip1_found = .true.
431 vco%ip1_sfc = vgd_ip1_T(jlev)
432 end if
433 end do
434 if (.not. ip1_found) then
435 write(*,*) 'vco_setupAtmFromFile: Could not find IP1=',ip1_sfc
436 call utl_abort('vco_setupAtmFromFile: No surface level found in Vgrid!!!')
437 else
438 if (mmpi_myid == 0 .and. .not. beSilent) write(*,*) 'vco_setupAtmFromFile: Set surface level IP1=', vco%ip1_sfc
439 end if
440
441 ! determine IP1s of 2m and 10m levels
442 call set_2m_10m_levels(vco)
443
444 vco%initialized = .true.
445
446 ierr = fstfrm(nultemplate)
447 ierr = fclos (nultemplate)
448
449 end subroutine vco_setupAtmFromFile
450
451 !--------------------------------------------------------------------------
452 ! vco_setupOceanFromFile
453 !--------------------------------------------------------------------------
454 subroutine vco_setupOceanFromFile(vco,templatefile,etiket,beSilent)
455 !
456 !:Purpose: Initialize vertical coordinate object with information from
457 ! a standard file. For Ocean fields on depth levels.
458 !
459 implicit none
460
461 ! Arguments:
462 type(struct_vco), pointer, intent(inout) :: vco ! Vertical coordinate object
463 character(len=*), intent(in) :: templatefile ! Template file
464 character(len=*), intent(in) :: etiket
465 logical, intent(in) :: beSilent
466
467 ! Locals:
468 integer :: nultemplate, ierr
469 integer, parameter :: maxNumRecords = 500
470 integer :: recordIndex, numRecords, ikeys(maxNumRecords)
471 integer :: fnom,fstouv,fstfrm,fclos,fstprm,fstinl
472 character(len=10) :: blk_S
473 integer :: IP1kind
474 real :: vertCoordValue
475 integer, parameter :: maxNumDepthLevels = 200
476 real(8) :: depths(maxNumDepthLevels)
477 integer :: ip1_depth(maxNumDepthLevels)
478 integer :: ideet, inpas, dateStamp_origin, ini, inj, ink, inbits, idatyp
479 integer :: ip1, ip2, ip3, ig1, ig2, ig3, ig4, iswa, ilng, idltf, iubc
480 integer :: iextra1, iextra2, iextra3
481 character(len=2) :: typvar
482 character(len=4) :: nomvar
483 character(len=1) :: grtyp
484
485 if (.not. beSilent) write(*,*) 'vco_setupOceanFromFile: found ocean fields'
486
487 ! initialize some components of vco
488 vco%vgridPresent = .false.
489 vco%nlev_T = 0
490 vco%nlev_M = 0
491 vco%Vcode = 0
492 vco%initialized = .true.
493
494 ! open the template file
495 nultemplate = 0
496 ierr = fnom(nultemplate,templatefile,'RND+OLD+R/O',0)
497 if ( ierr == 0 ) then
498 ierr = fstouv(nultemplate,'RND+OLD')
499 else
500 write(*,*) 'vco_setupFromFile: Template File = ', trim(templatefile)
501 call utl_abort('vco_setupOceanFromFile: CANNOT OPEN TEMPLATE FILE!')
502 end if
503
504 ierr = fstinl(nultemplate,ini,inj,ink,-1,etiket,-1,-1,-1,' ', &
505 ' ',ikeys,numRecords,maxNumRecords)
506 if ( ikeys(1) <= 0 ) then
507 call utl_abort('vco_setupOceanFromFile: Could not find any records ' // &
508 'in the supplied file')
509 end if
510 record_loop: do recordIndex = 1, numRecords
511 ierr = fstprm(ikeys(recordIndex), dateStamp_origin, ideet, inpas, ini, inj, &
512 ink, inbits, idatyp, ip1, ip2, ip3, &
513 typvar, nomvar, etiket, grtyp, ig1, ig2, ig3, ig4, &
514 iswa, ilng, idltf, iubc, iextra1, iextra2, iextra3)
515
516 ! ignore any variables not present in varnamelist_mod
517 if (.not. vnl_varnameIsValid(trim(nomvar))) cycle record_loop
518
519 ! check for record with ocean data on depth levels
520 call convip(ip1, vertCoordValue, Ip1Kind, -1, blk_s, .false.)
521 if ( Ip1Kind == 0 .and. vertCoordValue >= 0.0 .and. &
522 vnl_varKindFromVarname(trim(nomvar)) == 'OC' .and. &
523 vnl_varLevelFromVarname(trim(nomvar)) == 'DP' ) then
524 ! check if we've NOT already recorded this depth level
525 if ( .not. any(vertCoordValue == depths(1:vco%nLev_depth)) ) then
526 vco%nLev_depth = vco%nLev_depth + 1
527 depths(vco%nLev_depth) = vertCoordValue
528 ip1_depth(vco%nLev_depth) = ip1
529 if (mmpi_myid == 0 .and. .not.beSilent) then
530 write(*,*) 'vco_setupOceanFromFile: found ocean record: nLev_depth = ', &
531 vco%nLev_depth, 'varName = ', trim(nomvar), &
532 ', value = ', depths(vco%nLev_depth)
533 end if
534 end if
535 cycle record_loop
536 end if
537
538 end do record_loop
539
540 ierr = fstfrm(nultemplate)
541 ierr = fclos (nultemplate)
542
543 ! Allocate object arrays and copy in depth information
544 allocate(vco%depths(vco%nLev_depth))
545 allocate(vco%ip1_depth(vco%nLev_depth))
546 vco%depths(:) = depths(1:vco%nLev_depth)
547 vco%ip1_depth(:) = ip1_depth(1:vco%nLev_depth)
548
549 ! Check if ocean depth levels are in correct order (ascending in value)
550 if ( vco%nLev_depth > 1 ) then
551 if ( any(vco%depths(2:vco%nLev_depth)-vco%depths(1:(vco%nLev_depth-1)) < 0.0) ) then
552 call utl_abort('vco_setupOceanFromFile: some depth levels not in ascending order')
553 end if
554 end if
555
556 end subroutine vco_setupOceanFromFile
557
558 !--------------------------------------------------------------------------
559 ! vco_setupSfcFromFile
560 !--------------------------------------------------------------------------
561 subroutine vco_setupSfcFromFile(vco,beSilent)
562 !
563 !:Purpose: Initialize vertical coordinate object with information from
564 ! a standard file. For surface only fields.
565 !
566 implicit none
567
568 ! Arguments:
569 type(struct_vco), pointer, intent(inout) :: vco ! Vertical coordinate object
570 logical, intent(in) :: beSilent
571
572 if (.not. beSilent) write(*,*) 'vco_setupSfcFromFile: found surface fields'
573
574 ! initialize some components of vco
575 vco%vgridPresent = .false.
576 vco%nlev_T = 0
577 vco%nlev_M = 0
578 vco%Vcode = 0
579 vco%initialized = .true.
580
581 end subroutine vco_setupSfcFromFile
582
583 !--------------------------------------------------------------------------
584 ! vco_deallocate
585 !--------------------------------------------------------------------------
586 subroutine vco_deallocate(vco)
587 !
588 !:Purpose: Deallocate vertical coordinate object
589 !
590 implicit none
591
592 ! Arguments:
593 type(struct_vco), pointer, intent(inout) :: vco ! Vertical coordinate object
594
595 ! Locals:
596 integer :: stat
597
598 if ( vco%vgridPresent ) then
599 deallocate(vco%ip1_M)
600 deallocate(vco%ip1_T)
601 stat = vgd_free(vco%vgrid)
602 end if
603
604 if ( vco%nLev_depth > 0 ) then
605 deallocate(vco%depths)
606 deallocate(vco%ip1_depth)
607 end if
608 nullify(vco)
609
610 end subroutine vco_deallocate
611
612 !--------------------------------------------------------------------------
613 ! vco_getNumLev
614 !--------------------------------------------------------------------------
615 function vco_getNumLev(vco,varLevel,varName_opt) result(nlev)
616 !
617 !:Purpose: get number of vertical levels
618 !
619 implicit none
620
621 ! Arguments:
622 type(struct_vco), pointer, intent(in) :: vco ! Vertical coordinate object
623 character(len=*), intent(in) :: varLevel ! 'TH', 'MM', 'SF', 'SFMM', 'SFTH', 'DP', 'SS' or 'OT'
624 character(len=*), optional, intent(in) :: varName_opt ! only needed for varLevel='OT'
625 ! Result:
626 integer :: nlev
627
628 ! Locals:
629 integer :: varListIndex
630
631 if (varLevel == 'MM') then
632 nlev = vco%nlev_M
633 else if (varLevel == 'TH') then
634 nlev = vco%nlev_T
635 else if (varLevel == 'SF' .or. varLevel == 'SFTH' .or. &
636 varLevel == 'SFMM' .or. varLevel == 'SS') then
637 nlev = 1
638 else if (varLevel == 'OT') then
639 if (.not. present(varName_opt)) then
640 call utl_abort('vco_getNumLev: varName must be specified for varLevel=OT')
641 end if
642 varListIndex = vnl_varListIndexOther(varName_opt)
643 nlev = vco%nlev_Other(varListIndex)
644 else if (varLevel == 'DP') then
645 nlev = vco%nlev_depth
646 else
647 call utl_abort('vco_getNumLev: Unknown variable type! ' // varLevel)
648 end if
649
650 end function vco_getNumLev
651
652 !--------------------------------------------------------------------------
653 ! vco_mpiBcast
654 !--------------------------------------------------------------------------
655 subroutine vco_mpiBcast(vco)
656 !
657 !:Purpose: MPI broadcast of vertical coordinate object
658 !
659 implicit none
660
661 ! Arguments:
662 type(struct_vco), pointer, intent(inout) :: vco ! vertical coordinate object
663
664 ! Locals:
665 integer :: ierr, vgd_nlev_M, vgd_nlev_T
666 integer :: vgdig1, vgdig2, vgdig3, vgdig4, vgdip1, vgdip2, vgdip3, vgddate
667 integer :: vgdtable_dim1, vgdtable_dim2, vgdtable_dim3
668 character(len=12) :: vgdetik
669 real(8), pointer :: vgdtable(:,:,:)
670
671 nullify(vgdtable)
672
673 write(*,*) 'vco_mpiBcast: starting'
674
675 if (mmpi_myid > 0) then
676 if (.not.associated(vco)) then
677 allocate(vco)
678 else
679 call utl_abort('vco_mpiBcast: vco must be nullified for mpi task id > 0')
680 end if
681 end if
682
683 call rpn_comm_bcast(vco%initialized , 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
684 call rpn_comm_bcast(vco%vgridPresent, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
685 call rpn_comm_bcast(vco%nlev_T , 1, 'MPI_INTEGER', 0, 'GRID', ierr)
686 call rpn_comm_bcast(vco%nlev_M , 1, 'MPI_INTEGER', 0, 'GRID', ierr)
687 call rpn_comm_bcast(vco%ip1_sfc , 1, 'MPI_INTEGER', 0, 'GRID', ierr)
688 call rpn_comm_bcast(vco%ip1_T_2m , 1, 'MPI_INTEGER', 0, 'GRID', ierr)
689 call rpn_comm_bcast(vco%ip1_M_10m , 1, 'MPI_INTEGER', 0, 'GRID', ierr)
690 call rpn_comm_bcast(vco%nlev_depth , 1, 'MPI_INTEGER', 0, 'GRID', ierr)
691 call rpn_comm_bcast(vco%Vcode , 1, 'MPI_INTEGER', 0, 'GRID', ierr)
692 call rpn_comm_bcast(vco%nlev_other, vnl_numvarmaxOther, 'MPI_INTEGER', 0, 'GRID', ierr)
693 if (vco%nLev_depth > 0) then
694 if (mmpi_myid > 0) then
695 allocate(vco%ip1_depth(vco%nlev_depth))
696 allocate(vco%depths(vco%nlev_depth))
697 end if
698 call rpn_comm_bcast(vco%ip1_depth , vco%nlev_depth, 'MPI_INTEGER', 0, 'GRID', ierr)
699 call rpn_comm_bcast(vco%depths , vco%nlev_depth, 'MPI_REAL8' , 0, 'GRID', ierr)
700 end if
701 if (vco%vgridPresent) then
702 if (mmpi_myid == 0) then
703 vgd_nlev_M = size(vco%ip1_M)
704 vgd_nlev_T = size(vco%ip1_T)
705 end if
706 call rpn_comm_bcast(vgd_nlev_M, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
707 call rpn_comm_bcast(vgd_nlev_T, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
708 if (mmpi_myid > 0) then
709 allocate(vco%ip1_M(vgd_nlev_M))
710 allocate(vco%ip1_T(vgd_nlev_T))
711 end if
712 if (vgd_nlev_M > 0) call rpn_comm_bcast(vco%ip1_M, vgd_nlev_M, 'MPI_INTEGER', 0, 'GRID', ierr)
713 if (vgd_nlev_T > 0) call rpn_comm_bcast(vco%ip1_T, vgd_nlev_T, 'MPI_INTEGER', 0, 'GRID', ierr)
714 end if
715
716 ! now do bcast for vgrid object
717 if (vco%vgridPresent) then
718
719 if (mmpi_myid == 0) then
720 ierr = vgd_get(vco%vgrid,'VTBL',vgdtable)
721 vgdtable_dim1 = size(vgdtable,1)
722 vgdtable_dim2 = size(vgdtable,2)
723 vgdtable_dim3 = size(vgdtable,3)
724 ierr = vgd_get(vco%vgrid,'DATE',vgddate)
725 ierr = vgd_get(vco%vgrid,'ETIK',vgdetik)
726 ierr = vgd_get(vco%vgrid,'IG_1',vgdig1)
727 ierr = vgd_get(vco%vgrid,'IG_2',vgdig2)
728 ierr = vgd_get(vco%vgrid,'IG_3',vgdig3)
729 ierr = vgd_get(vco%vgrid,'IG_4',vgdig4)
730 ierr = vgd_get(vco%vgrid,'IP_1',vgdip1)
731 ierr = vgd_get(vco%vgrid,'IP_2',vgdip2)
732 ierr = vgd_get(vco%vgrid,'IP_3',vgdip3)
733 end if
734
735 ! 3D table of real*8
736 call rpn_comm_bcast(vgdtable_dim1, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
737 call rpn_comm_bcast(vgdtable_dim2, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
738 call rpn_comm_bcast(vgdtable_dim3, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
739 if (mmpi_myid > 0) allocate(vgdtable(vgdtable_dim1, vgdtable_dim2, vgdtable_dim3))
740 call rpn_comm_bcast(vgdtable, size(vgdtable), 'MPI_REAL8', 0, 'GRID', ierr)
741 ! others
742 call rpn_comm_bcast(vgddate, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
743 call rpn_comm_bcastc(vgdetik, len(vgdetik), 'MPI_CHARACTER', 0, 'GRID', ierr)
744 call rpn_comm_bcast(vgdig1, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
745 call rpn_comm_bcast(vgdig2, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
746 call rpn_comm_bcast(vgdig3, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
747 call rpn_comm_bcast(vgdig4, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
748 call rpn_comm_bcast(vgdip1, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
749 call rpn_comm_bcast(vgdip2, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
750 call rpn_comm_bcast(vgdip3, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
751
752 if (mmpi_myid > 0) then
753 ierr = vgd_new(vco%vgrid,vgdtable)
754 ierr = vgd_put(vco%vgrid,'DATE',vgddate)
755 ierr = vgd_put(vco%vgrid,'ETIK',vgdetik)
756 ierr = vgd_put(vco%vgrid,'IG_1',vgdig1)
757 ierr = vgd_put(vco%vgrid,'IG_2',vgdig2)
758 ierr = vgd_put(vco%vgrid,'IG_3',vgdig3)
759 ierr = vgd_put(vco%vgrid,'IG_4',vgdig4)
760 ierr = vgd_put(vco%vgrid,'IP_1',vgdip1)
761 ierr = vgd_put(vco%vgrid,'IP_2',vgdip2)
762 ierr = vgd_put(vco%vgrid,'IP_3',vgdip3)
763 end if
764
765 deallocate(vgdtable)
766
767 end if
768
769 write(*,*) 'vco_mpiBcast: done'
770
771 end subroutine vco_mpiBcast
772
773 !--------------------------------------------------------------------------
774 ! vco_equal
775 !--------------------------------------------------------------------------
776 function vco_equal(vco1,vco2) result(equal)
777 !
778 !:Purpose: Compare two vertical grid object and provide a logical result if they are equal or not
779 !
780 implicit none
781
782 ! Arguments:
783 type(struct_vco), pointer, intent(in) :: vco1 ! vertical coordinate object one
784 type(struct_vco), pointer, intent(in) :: vco2 ! vertical coordinate object two
785 ! Result:
786 logical :: equal
787
788 equal = .true.
789
790 equal = equal .and. (vco1%Vcode == vco2%Vcode)
791 if (.not. equal) then
792 write(*,*) 'vco_equal: Vcode not equal'
793 return
794 end if
795 if ( vco1%vgridPresent .and. vco2%vgridPresent ) then
796 equal = equal .and. (vco1%vgrid == vco2%vgrid)
797 if (.not. equal) then
798 write(*,*) 'vco_equal: vgrid not equal'
799 return
800 end if
801 end if
802
803 ! Even if vgrid defined, not enough just to compare vgrid, must compare everything
804 equal = equal .and. (vco1%nlev_T == vco2%nlev_T)
805 if (.not. equal) then
806 write(*,*) 'vco_equal: nlev_T not equal', vco1%nlev_T, vco2%nlev_T
807 return
808 end if
809 equal = equal .and. (vco1%nlev_M == vco2%nlev_M)
810 if (.not. equal) then
811 write(*,*) 'vco_equal: nlev_M not equal', vco1%nlev_M, vco2%nlev_M
812 return
813 end if
814 if (vco1%vgridPresent .and. vco2%vgridPresent .and. &
815 vco1%nlev_T > 0 .and. vco2%nlev_T > 0) then
816 equal = equal .and. all(vco1%ip1_T(:) == vco2%ip1_T(:))
817 if (.not. equal) then
818 write(*,*) 'vco_equal: ip1_T not equal'
819 return
820 end if
821 equal = equal .and. all(vco1%ip1_M(:) == vco2%ip1_M(:))
822 if (.not. equal) then
823 write(*,*) 'vco_equal: ip1_M not equal'
824 return
825 end if
826 equal = equal .and. (vco1%ip1_sfc == vco2%ip1_sfc)
827 if (.not. equal) then
828 write(*,*) 'vco_equal: ip1_sfc not equal'
829 return
830 end if
831 if (vco1%Vcode == 5002 .or. vco1%Vcode == 5005 .or. vco1%Vcode == 5100) then
832 equal = equal .and. hybridCoefEqualOrNot(vco1, vco2)
833 if (.not. equal) then
834 write(*,*) 'vco_equal: hybrid parameters are not equal'
835 return
836 end if
837 end if
838 end if
839
840 ! For ocean fields, check depth levels
841 if (vco1%nLev_depth > 0) then
842 equal = equal .and. all(vco1%depths(:) == vco2%depths(:))
843 if (.not. equal) then
844 write(*,*) 'vco_equal: ocean depth levels are not equal'
845 return
846 end if
847 end if
848
849 end function vco_equal
850
851 !--------------------------------------------------------------------------
852 ! vco_subsetOrNot
853 !--------------------------------------------------------------------------
854 function vco_subsetOrNot(vco_template, vco_full) result(subset)
855 !
856 !:Purpose: This function determines if vco_template is a subset of vco_full.
857 !
858 implicit none
859
860 ! Arguments:
861 type(struct_vco), pointer, intent(in) :: vco_full ! vertical coordinate object full
862 type(struct_vco), pointer, intent(in) :: vco_template ! vertical coordinate object template
863 ! Result:
864 logical :: subset
865
866 ! Locals:
867 integer, allocatable :: THlevelWanted(:), MMlevelWanted(:)
868
869 !
870 !- Compare the vCode
871 !
872 if (vco_template%Vcode /= vco_full%Vcode) then
873 subset = .false.
874 return
875 end if
876
877 !
878 !- Check if there are any thermo or momentum levels
879 !
880 if ( (vco_template%nlev_T == 0) .and. (vco_template%nlev_M == 0) ) then
881 subset = .false.
882 return
883 end if
884
885 !
886 !- Compare the IP1s
887 !
888 allocate(THlevelWanted(vco_template%nlev_T))
889 allocate(MMlevelWanted(vco_template%nlev_M))
890
891 call vco_levelMatchingList( THlevelWanted, MMlevelWanted, & ! OUT
892 vco_template, vco_full ) ! IN
893
894 if (any(THlevelWanted == -1) .or. any(MMlevelWanted == -1)) then
895 subset = .false.
896 return
897 end if
898
899 deallocate(MMlevelWanted)
900 deallocate(THlevelWanted)
901
902 !
903 !- For hybrid coordinates, compare additional grid parameters
904 !
905 if ( .not. hybridCoefEqualOrNot(vco_template, vco_full) ) then
906 subset = .false.
907 return
908 end if
909
910 !
911 !- When reaching this point, we assume that we have a subset
912 !
913 subset = .true.
914
915 end function vco_subsetOrNot
916
917 !--------------------------------------------------------------------------
918 ! hybridCoefEqualOrNot
919 !--------------------------------------------------------------------------
920 function hybridCoefEqualOrNot(vco1, vco2) result(equal)
921 !
922 !:Purpose: To compare two vertical coordinate hybrid coefficient object
923 !
924 implicit none
925
926 ! Arguments:
927 type(struct_vco), pointer, intent(in) :: vco1 ! vertical coordinate object one
928 type(struct_vco), pointer, intent(in) :: vco2 ! vertical coordinate object two
929 ! Result:
930 logical :: equal
931
932 ! Locals:
933 real(8) :: ptop1, ptop2
934 real(8), pointer :: coefA1(:), coefA2(:)
935 real(8), pointer :: coefB1(:), coefB2(:)
936 real :: coefR11, coefR12
937 real :: coefR21, coefR22
938 integer :: stat, levIndex
939
940 if (vco1%Vcode == 5002) then
941 !- Ptop
942 stat = vgd_get(vco1%vgrid,key='PTOP - top level pressure',value=ptop1)
943 stat = vgd_get(vco2%vgrid,key='PTOP - top level pressure',value=ptop2)
944 if (ptop1 /= ptop2) then
945 equal = .false.
946 return
947 end if
948 end if
949
950 if (vco1%Vcode == 5002 .or. vco1%Vcode == 5005 .or. vco1%Vcode == 5100) then
951
952 !- Pref
953 stat = vgd_get(vco1%vgrid,key='PREF - reference pressure',value=ptop1)
954 stat = vgd_get(vco2%vgrid,key='PREF - reference pressure',value=ptop2)
955 if (ptop1 /= ptop2) then
956 equal = .false.
957 return
958 end if
959 !- R-coef 1
960 stat = vgd_get(vco1%vgrid,key='RC_1 - first R-coef value',value=coefR11)
961 stat = vgd_get(vco2%vgrid,key='RC_1 - first R-coef value',value=coefR12)
962 if (coefR11 /= coefR12) then
963 equal = .false.
964 return
965 end if
966 !- R-coef 2
967 stat = vgd_get(vco1%vgrid,key='RC_2 - second R-coef value',value=coefR21)
968 stat = vgd_get(vco2%vgrid,key='RC_2 - second R-coef value',value=coefR22)
969 if (coefR21 /= coefR22) then
970 equal = .false.
971 return
972 end if
973 !- A
974 nullify(coefA1)
975 nullify(coefA2)
976 stat = vgd_get(vco1%vgrid,key='CA_M - vertical A coefficient (m)',value=coefA1)
977 stat = vgd_get(vco2%vgrid,key='CA_M - vertical A coefficient (m)',value=coefA2)
978 if ( size(coefA1) /= size(coefA2) ) then
979 equal = .false.
980 return
981 end if
982 do levIndex = 1, size(coefA1)
983 if (coefA1(levIndex) /= coefA2(levIndex)) then
984 equal = .false.
985 return
986 end if
987 end do
988 !- B
989 nullify(coefB1)
990 nullify(coefB2)
991 stat = vgd_get(vco1%vgrid,key='CB_M - vertical B coefficient (m)',value=coefB1)
992 stat = vgd_get(vco2%vgrid,key='CB_M - vertical B coefficient (m)',value=coefB2)
993 if ( size(coefB1) /= size(coefB2) ) then
994 equal = .false.
995 return
996 end if
997 do levIndex = 1, size(coefB1)
998 if (coefB1(levIndex) /= coefB2(levIndex)) then
999 equal = .false.
1000 return
1001 end if
1002 end do
1003
1004 end if
1005
1006 !
1007 !- When reaching this point, we assume that they are equal
1008 !
1009 equal = .true.
1010
1011 end function hybridCoefEqualOrNot
1012
1013 !--------------------------------------------------------------------------
1014 ! vco_levelMatchingList
1015 !--------------------------------------------------------------------------
1016 subroutine vco_levelMatchingList(THmatchingList, MMmatchingList, vco1, vco2)
1017 !
1018 !:Purpose: This subroutine returns arrays of array indices of the levels (ip1s) in vco2
1019 ! corresponding with the levels (ip1s) in vco1
1020 !
1021 implicit none
1022
1023 ! Arguments:
1024 type(struct_vco), pointer, intent(in) :: vco1 ! vertical coordinate object one
1025 type(struct_vco), pointer, intent(in) :: vco2 ! vertical coordinate object two
1026 integer, intent(out) :: THmatchingList(vco1%nlev_T) ! TH matching list
1027 integer, intent(out) :: MMmatchingList(vco1%nlev_M) ! MM matching list
1028
1029 ! Locals:
1030 integer :: levIndex1, levIndex2
1031
1032 !
1033 !- Do momentum levels...
1034 !
1035 MMmatchingList(:) = -1
1036 do levIndex1 = 1, vco1%nlev_M
1037 do levIndex2 = 1, vco2%nlev_M
1038 if ( (vco2%ip1_M(levIndex2) == vco1%ip1_M(levIndex1)) .or. &
1039 (vco2%ip1_M(levIndex2) == vco2%ip1_M_10m .and. &
1040 vco1%ip1_M(levIndex1) == vco1%ip1_M_10m) ) then
1041 MMmatchingList(levIndex1) = levIndex2
1042 exit
1043 end if
1044 end do
1045 end do
1046
1047 !
1048 !- Do thermo levels...
1049 !
1050 THmatchingList(:) = -1
1051 do levIndex1 = 1, vco1%nlev_T
1052 do levIndex2 = 1, vco2%nlev_T
1053 if ( (vco2%ip1_T(levIndex2) == vco1%ip1_T(levIndex1)) .or. &
1054 (vco2%ip1_T(levIndex2) == vco2%ip1_T_2m .and. &
1055 vco1%ip1_T(levIndex1) == vco1%ip1_T_2m) ) then
1056 THmatchingList(levIndex1) = levIndex2
1057 exit
1058 end if
1059 end do
1060 end do
1061
1062 end subroutine vco_levelMatchingList
1063
1064 !--------------------------------------------------------------------------
1065 ! set_2m_10m_levels
1066 !--------------------------------------------------------------------------
1067 subroutine set_2m_10m_levels(vco)
1068 !
1069 !:Purpose: To set 2-m and 10-m levels
1070 !
1071 implicit none
1072
1073 ! Arguments:
1074 type(struct_vco), pointer, intent(in) :: vco ! vertical coordinate object
1075
1076 ! Locals:
1077 character(len=10) :: blk_s
1078
1079 if (vco%Vcode == 5002) then
1080 vco%ip1_T_2m = vco%ip1_sfc
1081 vco%ip1_M_10m = vco%ip1_sfc
1082 else if (vco%Vcode == 5005 .or. vco%Vcode == 5100) then
1083 call convip(vco%ip1_T_2m , 1.5, 4, 2, blk_s, .false.)
1084 call convip(vco%ip1_M_10m, 10.0, 4, 2, blk_s, .false.)
1085 else
1086 vco%ip1_T_2m = -1
1087 vco%ip1_M_10m = -1
1088 end if
1089
1090 end subroutine set_2m_10m_levels
1091
1092end module verticalCoord_mod
