Line data Source code
1 : module inic_analytic
2 :
3 : !-----------------------------------------------------------------------
4 : !
5 : ! Purpose: Set analytic initial conditions based on input coordinates
6 : !
7 : !
8 : !-----------------------------------------------------------------------
9 : use cam_logfile, only: iulog
10 : use shr_kind_mod, only: r8 => shr_kind_r8
11 : use cam_abortutils, only: endrun
12 : use shr_sys_mod, only: shr_sys_flush
13 : use inic_analytic_utils, only: analytic_ic_active, analytic_ic_type
14 :
15 : implicit none
16 : private
17 :
18 : public :: analytic_ic_active ! forwarded from init_analytic_utils
19 : public :: analytic_ic_set_ic ! Set analytic initial conditions
20 : public :: dyn_set_inic_col
21 :
22 : interface analytic_ic_set_ic
23 : module procedure dyn_set_inic_cblock
24 : end interface analytic_ic_set_ic
25 :
26 : ! Private module variables
27 : integer :: call_num = 0
28 :
29 : ! Private interface
30 : #ifdef ANALYTIC_IC
31 : interface get_input_shape
32 : module procedure get_input_shape_2d
33 : module procedure get_input_shape_3d
34 : end interface get_input_shape
35 : #endif
36 :
37 : !==============================================================================
38 : CONTAINS
39 : !==============================================================================
40 :
41 0 : subroutine dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint, U, V, T, &
42 : PS, PHIS_IN, PHIS_OUT, Q, m_cnst, mask, verbose)
43 : use cam_initfiles, only: pertlim
44 : #ifdef ANALYTIC_IC
45 : use ic_held_suarez, only: hs94_set_ic
46 : use ic_baroclinic, only: bc_wav_set_ic
47 : use ic_baro_dry_jw06, only: bc_dry_jw06_set_ic
48 : use ic_us_standard_atmosphere, only: us_std_atm_set_ic
49 : #endif
50 : use spmd_utils, only: masterproc
51 : !-----------------------------------------------------------------------
52 : !
53 : ! Purpose: Set analytic initial values for dynamics state variables
54 : !
55 : !-----------------------------------------------------------------------
56 :
57 : ! Dummy arguments
58 : integer , intent(in) :: vcoord ! See dyn_tests_utils
59 : real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol)
60 : real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol)
61 : integer, intent(in) :: glob_ind(:) ! global column index
62 : real(r8), optional, intent(in) :: zint(:,:) ! height at layer interfaces
63 : real(r8), optional, intent(inout) :: U(:,:) ! zonal velocity
64 : real(r8), optional, intent(inout) :: V(:,:) ! meridional velocity
65 : real(r8), optional, intent(inout) :: T(:,:) ! temperature
66 : real(r8), optional, intent(inout) :: PS(:) ! surface pressure
67 : real(r8), optional, intent(in) :: PHIS_IN(:) ! surface geopotential
68 : real(r8), optional, intent(out) :: PHIS_OUT(:) ! surface geopotential
69 : real(r8), optional, intent(inout) :: Q(:,:,:) ! tracer (ncol, lev, m)
70 : integer, optional, intent(in) :: m_cnst(:) ! tracer indices (reqd. if Q)
71 : logical, optional, intent(in) :: mask(:) ! Only init where .true.
72 : logical, optional, intent(in) :: verbose ! For internal use
73 :
74 : ! Local variables
75 : logical :: verbose_use
76 : logical, allocatable :: mask_use(:)
77 : real(r8) :: pertval
78 : integer, allocatable :: rndm_seed(:)
79 : integer :: rndm_seed_sz
80 : integer :: i, k
81 : integer :: ncol, nlev
82 : character(len=*), parameter :: subname = 'DYN_SET_INIC_COL'
83 :
84 : #ifdef ANALYTIC_IC
85 : allocate(mask_use(size(latvals)))
86 : if (present(mask)) then
87 : if (size(mask_use) /= size(mask)) then
88 : call endrun('cnst_init_default: input, mask, is wrong size')
89 : end if
90 : mask_use = mask
91 : else
92 : mask_use = .true.
93 : end if
94 :
95 : if (present(verbose)) then
96 : verbose_use = verbose
97 : else
98 : verbose_use = .true.
99 : end if
100 :
101 : ! Basic size sanity checks
102 : if (size(latvals) /= size(lonvals)) then
103 : call endrun(subname//'latvals and lonvals must have same size')
104 : end if
105 : if (present(U)) then
106 : if (size(U) > 0) then
107 : call check_array_size(U(:,1), 'U', latvals, subname)
108 : else
109 : return
110 : end if
111 : end if
112 : if (present(V)) then
113 : if (size(V) > 0) then
114 : call check_array_size(V(:,1), 'V', latvals, subname)
115 : else
116 : return
117 : end if
118 : end if
119 : if (present(T)) then
120 : if (size(T) > 0) then
121 : call check_array_size(T(:,1), 'T', latvals, subname)
122 : else
123 : return
124 : end if
125 : end if
126 : if (present(PS)) then
127 : if (size(PS) > 0) then
128 : call check_array_size(PS, 'PS', latvals, subname)
129 : else
130 : return
131 : end if
132 : end if
133 : if (present(PHIS_IN)) then
134 : if (size(PHIS_IN) > 0) then
135 : call check_array_size(PHIS_IN, 'PHIS_IN', latvals, subname)
136 : else
137 : return
138 : end if
139 : end if
140 : if (present(PHIS_OUT)) then
141 : if (size(PHIS_OUT) > 0) then
142 : call check_array_size(PHIS_OUT, 'PHIS_OUT', latvals, subname)
143 : else
144 : return
145 : end if
146 : end if
147 : ! Some special checks on the tracer argument
148 : if (present(Q)) then
149 : if (.not. present(m_cnst)) then
150 : call endrun(subname//'m_cnst is required if Q is present')
151 : end if
152 : if (size(Q, 3) /= size(m_cnst, 1)) then
153 : call endrun(subname//': size of m_cnst must match last dimension of Q')
154 : end if
155 : if (size(Q) > 0) then
156 : call check_array_size(Q(:,1,1), 'Q', latvals, subname)
157 : else
158 : return
159 : end if
160 : end if
161 : select case(trim(analytic_ic_type))
162 : case('held_suarez_1994')
163 : call hs94_set_ic(latvals, lonvals, U=U, V=V, T=T, PS=PS, PHIS=PHIS_OUT, &
164 : Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)
165 :
166 : case('moist_baroclinic_wave_dcmip2016', 'dry_baroclinic_wave_dcmip2016')
167 : call bc_wav_set_ic(vcoord, latvals, lonvals, zint=zint, U=U, V=V, T=T, PS=PS, &
168 : PHIS=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)
169 :
170 : case('dry_baroclinic_wave_jw2006')
171 : call bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U=U, V=V, T=T, PS=PS, &
172 : PHIS=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)
173 :
174 : case('us_standard_atmosphere')
175 : call us_std_atm_set_ic(latvals, lonvals, zint=zint, U=U, V=V, T=T, PS=PS, PHIS_IN=PHIS_IN, &
176 : PHIS_OUT=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)
177 :
178 : case default
179 : call endrun(subname//': Unknown analytic_ic_type, "'//trim(analytic_ic_type)//'"')
180 : end select
181 :
182 : ! Maybe peturb T initial conditions
183 : if (present(T) .and. (pertlim /= 0.0_r8)) then
184 :
185 : ! Add random perturbation to temperature if required
186 : if(masterproc .and. verbose_use) then
187 : write(iulog,*) trim(subname), ': Adding random perturbation bounded by +/-', &
188 : pertlim,' to initial temperature field'
189 : end if
190 : call random_seed(size=rndm_seed_sz)
191 : allocate(rndm_seed(rndm_seed_sz))
192 :
193 : ncol = size(T, 1)
194 : nlev = size(T, 2)
195 : do i = 1, ncol
196 : if (mask_use(i)) then
197 : ! seed random_number generator based on global column index
198 : rndm_seed(:) = glob_ind(i)
199 : call random_seed(put=rndm_seed)
200 : do k = 1, nlev
201 : call random_number(pertval)
202 : pertval = 2.0_r8 * pertlim * (0.5_r8 - pertval)
203 : T(i,k) = T(i,k) * (1.0_r8 + pertval)
204 : end do
205 : end if
206 : end do
207 :
208 : deallocate(rndm_seed)
209 : end if
210 :
211 : ! To get different random seeds each time
212 : call_num = call_num + 1
213 : #else
214 0 : call endrun(subname//': analytic initial conditions are not enabled')
215 : #endif
216 :
217 0 : end subroutine dyn_set_inic_col
218 :
219 0 : subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
220 : PS, PHIS_IN, PHIS_OUT, Q, m_cnst, mask)
221 : !-----------------------------------------------------------------------
222 : !
223 : ! Purpose: Set analytic initial values for dynamics state variables
224 : !
225 : !-----------------------------------------------------------------------
226 :
227 : ! Dummy arguments
228 : integer, intent(in) :: vcoord ! See dyn_tests_utils
229 : real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol)
230 : real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol)
231 : integer, intent(in) :: glob_ind(:) ! global column index
232 : real(r8), optional, intent(inout) :: U(:,:,:) ! zonal velocity
233 : real(r8), optional, intent(inout) :: V(:,:,:) ! meridional velocity
234 : real(r8), optional, intent(inout) :: T(:,:,:) ! temperature
235 : real(r8), optional, intent(inout) :: PS(:,:) ! surface pressure
236 : real(r8), optional, intent(in) :: PHIS_IN(:,:)! surface geopotential
237 : real(r8), optional, intent(out) :: PHIS_OUT(:,:)! surface geopotential
238 : real(r8), optional, intent(inout) :: Q(:,:,:,:) ! tracer (ncol,lev,blk,m)
239 : integer, optional, intent(in) :: m_cnst(:) ! tracer indices (reqd. if Q)
240 : logical, optional, intent(in) :: mask(:) ! Only init where .true.
241 :
242 : ! Local variables
243 : real(r8), allocatable :: lat_use(:)
244 : integer :: i, bbeg, bend
245 : integer :: size1, size2, size3
246 : integer :: nblks, blksize
247 : logical :: verbose
248 : character(len=4) :: mname
249 : character(len=*), parameter :: subname = 'DYN_SET_INIC_CBLOCK'
250 :
251 : #ifdef ANALYTIC_IC
252 : verbose = .true. ! So subroutines can report setting variables
253 : ! Figure out what sort of blocks we have, all variables should be the same
254 : size1 = -1
255 : mname = ''
256 : if (present(U)) then
257 : call get_input_shape(U, 'U', mname, size1, size2, size3, subname)
258 : end if
259 : if(present(V)) then
260 : call get_input_shape(V, 'V', mname, size1, size2, size3, subname)
261 : end if
262 : if(present(T)) then
263 : call get_input_shape(T, 'T', mname, size1, size2, size3, subname)
264 : end if
265 : if(present(Q)) then
266 : call get_input_shape(Q(:,:,:,1), 'Q', mname, size1, size2, size3, subname)
267 : end if
268 : ! Need to do all 3-D variables before any 2-D variables
269 : if(present(PS)) then
270 : call get_input_shape(PS, 'PS', mname, size1, size2, size3, subname)
271 : end if
272 : if(present(PHIS_IN)) then
273 : call get_input_shape(PHIS_IN, 'PHIS_IN', mname, size1, size2, size3, subname)
274 : end if
275 : if(present(PHIS_OUT)) then
276 : call get_input_shape(PHIS_OUT, 'PHIS_OUT', mname, size1, size2, size3, subname)
277 : end if
278 : if (size1 < 0) then
279 : call endrun(subname//': No state variables to initialize')
280 : end if
281 : if ((size(latvals) == size1*size3) .and. (size(lonvals) == size1*size3)) then
282 : ! Case: unstructured with blocks in 3rd dim
283 : if (size(glob_ind) /= size(latvals)) then
284 : call endrun(subname//': there must be a global index for every column')
285 : end if
286 : nblks = size3
287 : blksize = size1
288 : bend = 0
289 : do i = 1, nblks
290 : bbeg = bend + 1
291 : bend = bbeg + blksize - 1
292 : if (present(mask)) then
293 : if (size(mask) /= size(latvals)) then
294 : call endrun(subname//': incorrect mask size')
295 : end if
296 : if (present(U)) then
297 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
298 : glob_ind(bbeg:bend), U=U(:,:,i), mask=mask(bbeg:bend), verbose=verbose)
299 : end if
300 : if (present(V)) then
301 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
302 : glob_ind(bbeg:bend), V=V(:,:,i), mask=mask(bbeg:bend), verbose=verbose)
303 : end if
304 : if (present(PS).and.present(PHIS_IN).and.present(T)) then
305 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
306 : glob_ind(bbeg:bend), PS=PS(:,i), PHIS_IN=PHIS_IN(:,i), T=T(:,:,i), &
307 : mask=mask(bbeg:bend), verbose=verbose)
308 : else
309 : if (present(T)) then
310 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
311 : glob_ind(bbeg:bend), T=T(:,:,i), mask=mask(bbeg:bend), verbose=verbose)
312 : end if
313 : if (present(PHIS_OUT)) then
314 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
315 : glob_ind(bbeg:bend), PHIS_OUT=PHIS_OUT(:,i), mask=mask(bbeg:bend), verbose=verbose)
316 : end if
317 : if (present(PS)) then
318 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
319 : glob_ind(bbeg:bend), PS=PS(:,i), mask=mask(bbeg:bend), verbose=verbose)
320 : end if
321 : end if
322 : if (present(Q)) then
323 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
324 : glob_ind(bbeg:bend), Q=Q(:,:,i,:), m_cnst=m_cnst, &
325 : mask=mask(bbeg:bend), verbose=verbose)
326 : end if
327 : else
328 : if (present(U)) then
329 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
330 : glob_ind(bbeg:bend), U=U(:,:,i), verbose=verbose)
331 : end if
332 : if (present(V)) then
333 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
334 : glob_ind(bbeg:bend), V=V(:,:,i), verbose=verbose)
335 : end if
336 : if (present(PS).and.present(PHIS_IN).and.present(T)) then
337 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
338 : glob_ind(bbeg:bend), PHIS_IN=PHIS_IN(:,i),PS=PS(:,i),T=T(:,:,i), &
339 : verbose=verbose)
340 : else
341 : if (present(T)) then
342 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
343 : glob_ind(bbeg:bend), T=T(:,:,i), verbose=verbose)
344 : end if
345 : if (present(PS)) then
346 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
347 : glob_ind(bbeg:bend), PS=PS(:,i), verbose=verbose)
348 : end if
349 : if (present(PHIS_OUT)) then
350 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
351 : glob_ind(bbeg:bend), PHIS_OUT=PHIS_OUT(:,i), verbose=verbose)
352 : end if
353 : end if
354 : if (present(Q)) then
355 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
356 : glob_ind(bbeg:bend), Q=Q(:,:,i,:), m_cnst=m_cnst, &
357 : verbose=verbose)
358 : end if
359 : end if
360 : verbose = .false.
361 : end do
362 : else if ((size(latvals) == size1*size2) .and. (size(lonvals) == size1*size2)) then
363 : ! Case: unstructured with blocks in 2nd dim
364 : if (size(glob_ind) /= size(latvals)) then
365 : call endrun(subname//': there must be a global index for every column')
366 : end if
367 : nblks = size2
368 : blksize = size1
369 : bend = 0
370 : do i = 1, nblks
371 : bbeg = bend + 1
372 : bend = bbeg + blksize - 1
373 : if (present(mask)) then
374 : if (size(mask) /= size(latvals)) then
375 : call endrun(subname//': incorrect mask size')
376 : end if
377 : if (present(U)) then
378 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
379 : glob_ind(bbeg:bend), U=U(:,i,:), mask=mask(bbeg:bend), verbose=verbose)
380 : end if
381 : if (present(V)) then
382 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
383 : glob_ind(bbeg:bend), V=V(:,i,:), mask=mask(bbeg:bend), verbose=verbose)
384 : end if
385 : if (present(PS).and.present(PHIS_IN).and.present(T)) then
386 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
387 : glob_ind(bbeg:bend), PHIS_IN=PHIS_IN(:,i),PS=PS(:,i),T=T(:,i,:), &
388 : mask=mask(bbeg:bend), verbose=verbose)
389 : else
390 : if (present(T)) then
391 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
392 : glob_ind(bbeg:bend), T=T(:,i,:), mask=mask(bbeg:bend), verbose=verbose)
393 : end if
394 : if (present(PS)) then
395 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
396 : glob_ind(bbeg:bend), PS=PS(:,i), mask=mask(bbeg:bend), verbose=verbose)
397 : end if
398 : if (present(PHIS_OUT)) then
399 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
400 : glob_ind(bbeg:bend), PHIS_OUT=PHIS_OUT(:,i), mask=mask(bbeg:bend), verbose=verbose)
401 : end if
402 : end if
403 : if (present(Q)) then
404 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
405 : glob_ind(bbeg:bend), Q=Q(:,i,:,:), m_cnst=m_cnst, &
406 : mask=mask(bbeg:bend), verbose=verbose)
407 : end if
408 : else
409 : if (present(U)) then
410 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
411 : glob_ind(bbeg:bend), U=U(:,i,:), verbose=verbose)
412 : end if
413 : if (present(V)) then
414 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
415 : glob_ind(bbeg:bend), V=V(:,i,:), verbose=verbose)
416 : end if
417 : if (present(PS).and.present(PHIS_IN).and.present(T)) then
418 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
419 : glob_ind(bbeg:bend), PHIS_IN=PHIS_IN(:,i),PS=PS(:,i),T=T(:,i,:), &
420 : verbose=verbose)
421 : else
422 : if (present(T)) then
423 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
424 : glob_ind(bbeg:bend), T=T(:,i,:), verbose=verbose)
425 : end if
426 : if (present(PS)) then
427 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
428 : glob_ind(bbeg:bend), PS=PS(:,i), verbose=verbose)
429 : end if
430 : if (present(PHIS_OUT)) then
431 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
432 : glob_ind(bbeg:bend), PHIS_OUT=PHIS_OUT(:,i), verbose=verbose)
433 : end if
434 : end if
435 : if (present(Q)) then
436 : call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
437 : glob_ind(bbeg:bend), Q=Q(:,i,:,:), m_cnst=m_cnst, &
438 : verbose=verbose)
439 : end if
440 : end if
441 : verbose = .false.
442 : end do
443 : else if ((size(latvals) == size2) .and. (size(lonvals) == size1)) then
444 : ! Case: lon,lat,lev
445 : if (size(glob_ind) /= (size2 * size1)) then
446 : call endrun(subname//': there must be a global index for every column')
447 : end if
448 : nblks = size2
449 : allocate(lat_use(size(lonvals)))
450 : if (present(mask)) then
451 : call endrun(subname//': mask not supported for lon/lat')
452 : else
453 : bend = 0
454 : do i = 1, nblks
455 : bbeg = bend + 1
456 : bend = bbeg + size1 - 1
457 : lat_use = latvals(i)
458 : if (present(U)) then
459 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
460 : U=U(:,i,:), verbose=verbose)
461 : end if
462 : if (present(V)) then
463 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
464 : V=V(:,i,:), verbose=verbose)
465 : end if
466 : if (present(PS).and.present(PHIS_IN).and.present(T)) then
467 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
468 : PS=PS(:,i),T=T(:,i,:),PHIS_IN=PHIS_IN(:,i), verbose=verbose)
469 : else
470 : if (present(T)) then
471 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
472 : T=T(:,i,:), verbose=verbose)
473 : end if
474 : if (present(PS)) then
475 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
476 : PS=PS(:,i), verbose=verbose)
477 : end if
478 : if (present(PHIS_OUT)) then
479 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
480 : PHIS_OUT=PHIS_OUT(:,i), verbose=verbose)
481 : end if
482 : end if
483 : if (present(Q)) then
484 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
485 : Q=Q(:,i,:,:), m_cnst=m_cnst, verbose=verbose)
486 : end if
487 : verbose = .false.
488 : end do
489 : end if
490 : deallocate(lat_use)
491 : else if ((size(latvals) == size3) .and. (size(lonvals) == size1)) then
492 : if (size(glob_ind) /= (size3 * size1)) then
493 : call endrun(subname//': there must be a global index for every column')
494 : end if
495 : ! Case: lon,lev,lat
496 : nblks = size3
497 : allocate(lat_use(size(lonvals)))
498 : if (present(mask)) then
499 : call endrun(subname//': mask not supported for lon/lat')
500 : else
501 : bend = 0
502 : do i = 1, nblks
503 : bbeg = bend + 1
504 : bend = bbeg + size1 - 1
505 : lat_use = latvals(i)
506 : if (present(U)) then
507 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
508 : U=U(:,:,i), verbose=verbose)
509 : end if
510 : if (present(V)) then
511 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
512 : V=V(:,:,i), verbose=verbose)
513 : end if
514 : if (present(PS).and.present(PHIS_IN).and.present(T)) then
515 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
516 : T=T(:,:,i),PS=PS(:,i), PHIS_IN=PHIS_IN(:,i), verbose=verbose)
517 : else
518 : if (present(T)) then
519 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
520 : T=T(:,:,i), verbose=verbose)
521 : end if
522 : if (present(PS)) then
523 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
524 : PS=PS(:,i), verbose=verbose)
525 : end if
526 : if (present(PHIS_OUT)) then
527 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
528 : PHIS_OUT=PHIS_OUT(:,i), verbose=verbose)
529 : end if
530 : end if
531 : if (present(Q)) then
532 : call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
533 : Q=Q(:,:,i,:), m_cnst=m_cnst, verbose=verbose)
534 : end if
535 : verbose = .false.
536 : end do
537 : end if
538 : deallocate(lat_use)
539 : else
540 : call endrun(subname//': Unknown state variable layout')
541 : end if
542 : #else
543 0 : call endrun(subname//': analytic initial conditions are not enabled')
544 : #endif
545 0 : end subroutine dyn_set_inic_cblock
546 :
547 : #ifdef ANALYTIC_IC
548 : subroutine get_input_shape_2d(array, aname, sname, size1, size2, size3, es)
549 : real(r8), intent(in) :: array(:,:)
550 : character(len=*), intent(in) :: aname
551 : character(len=*), intent(inout) :: sname
552 : integer, intent(inout) :: size1
553 : integer, intent(inout) :: size2
554 : integer, intent(inout) :: size3
555 : character(len=*), intent(in) :: es
556 :
557 : if ((size1 < 0) .and. (size(array) == 0)) then
558 : ! The shape has not yet been set, set it to zero
559 : size1 = 0
560 : size2 = 0
561 : size3 = 0
562 : sname = trim(aname)
563 : else if (size1 < 0) then
564 : ! The shape has not yet been set, set it
565 : size1 = size(array, 1)
566 : size2 = size(array, 2)
567 : size3 = 1
568 : sname = trim(aname)
569 : else if ((size1 * size2 * size3) > 0) then
570 : ! For 2-D variables, the second dimension is always the block size
571 : ! However, since the shape may have been set by a 3-D variable, we
572 : ! need to pass either possibility
573 : if ( (size1 /= size(array, 1)) .or. &
574 : ((size2 /= size(array, 2)) .and. (size3 /= size(array, 2)))) then
575 : call endrun(trim(es)//': shape of '//trim(aname)//' does not match shape of '//trim(sname))
576 : end if
577 : ! No else, we cannot compare to zero size master array
578 : end if
579 :
580 : end subroutine get_input_shape_2d
581 :
582 : subroutine get_input_shape_3d(array, aname, sname, size1, size2, size3, es)
583 : real(r8), intent(in) :: array(:,:,:)
584 : character(len=*), intent(in) :: aname
585 : character(len=*), intent(inout) :: sname
586 : integer, intent(inout) :: size1
587 : integer, intent(inout) :: size2
588 : integer, intent(inout) :: size3
589 : character(len=*), intent(in) :: es
590 :
591 : if ((size1 < 0) .and. (size(array) == 0)) then
592 : ! The shape has not yet been set, set it to zero
593 : size1 = 0
594 : size2 = 0
595 : size3 = 0
596 : sname = trim(aname)
597 : else if (size1 < 0) then
598 : ! The shape has not yet been set, set it
599 : size1 = size(array, 1)
600 : size2 = size(array, 2)
601 : size3 = size(array, 3)
602 : sname = trim(aname)
603 : else if ((size1 * size2 * size3) > 0) then
604 : ! We have a shape, make sure array matches it
605 : if ((size1 /= size(array, 1)) .or. (size2 /= size(array, 2)) .or. (size3 /= size(array, 3))) then
606 : call endrun(trim(es)//': shape of '//trim(aname)//' does not match shape of '//trim(sname))
607 : end if
608 : end if
609 : ! No else, we cannot compare to zero size master array
610 : end subroutine get_input_shape_3d
611 :
612 : subroutine check_array_size(array, aname, check, subname)
613 : real(r8), intent(in) :: array(:)
614 : character(len=*), intent(in) :: aname
615 : real(r8), intent(in) :: check(:)
616 : character(len=*), intent(in) :: subname
617 :
618 : if (size(array, 1) /= size(check, 1)) then
619 : call endrun(trim(subname)//': '//trim(aname)//' has the wrong first dimension')
620 : end if
621 :
622 : end subroutine check_array_size
623 : #endif
624 :
625 : end module inic_analytic
|
