| blob | 8a0a5573d6d2e0d6451dc3afae110d42e2d22418 |
1 ! FCOPY, ICOPY, etc.: C pre-processor macros for C->N and N->C
2 ! horizontal interpolation. These exist solely to reduce code
3 ! complexity.
5 ! Copy from C array to N array:
6 #define FCOPY(N,C,i,j,k) \
7 N(k,i-nits+1)=W1(i,j)*C(II(i,j) ,JJ(i,j) ,k) \
8 +W2(i,j)*C(II(i,j)+1,JJ(i,j) ,k) \
9 +W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k) \
10 +W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k)
12 ! Copy from C array to (k,used) indexed N array:
13 #define uFCOPY(N,C,i,j,k) \
14 N(k,used)=W1(i,j)*C(II(i,j) ,JJ(i,j) ,k) \
15 +W2(i,j)*C(II(i,j)+1,JJ(i,j) ,k) \
16 +W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k) \
17 +W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k)
19 ! ICOPY: Grab from C array, but do not assign:
20 #define ICOPY(C,i,j,k) \
21 (W1(i,j)*C(II(i,j) ,JJ(i,j) ,k) \
22 + W2(i,j)*C(II(i,j)+1,JJ(i,j) ,k) \
23 + W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k) \
24 + W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k))
26 ! IKJCOPY: grab from IKJ C array, do not assign:
27 #define IKJCOPY(C,i,k,j) \
28 (W1(i,j)*C(II(i,j) ,k,JJ(i,j) ) \
29 + W2(i,j)*C(II(i,j)+1,k,JJ(i,j) ) \
30 + W3(i,j)*C(II(i,j)+a,k,JJ(i,j)-1) \
31 + W4(i,j)*C(II(i,j)+a,k,JJ(i,j)+1))
33 ! ICOPY2D: grab from IJ C array, do not assign:
34 #define ICOPY2D(C,i,j) \
35 (W1(i,j)*C(II(i,j) ,JJ(i,j) ) \
36 + W2(i,j)*C(II(i,j)+1,JJ(i,j) ) \
37 + W3(i,j)*C(II(i,j)+a,JJ(i,j)-1) \
38 + W4(i,j)*C(II(i,j)+a,JJ(i,j)+1))
40 ! Copying from N array to C array:
41 #define UPCOPY(C,N,i,j,k,ni,nj) \
42 C(k,i-istart+1)=(N(ni,nj+2,k) \
43 + N(ni-a ,nj+1,k)+ N(ni+1-a,nj+1,k) \
44 + N(ni-1,nj ,k)+ N(ni,nj ,k) + N(ni+1,nj ,k) \
45 + N(ni-a ,nj-1,k)+ N(ni+1-a,nj-1,k) \
46 + N(ni,nj-2,k) \
47 ) / 9
49 ! Average to C points from N points without assignment:
50 #define NGRAB(N,ni,nj,nk) \
51 (N(ni,nj+2,nk) \
52 + N(ni-a ,nj+1,nk)+ N(ni+1-a,nj+1,nk)\
53 + N(ni-1,nj ,nk) + N(ni,nj ,nk) + N(ni+1,nj ,nk)\
54 + N(ni-a ,nj-1,nk)+ N(ni+1-a,nj-1,nk)\
55 + N(ni,nj-2,nk) \
56 ) / 9
58 ! Average to C points from N points without assignment on an IKJ grid:
59 #define NGRABIKJ(N,ni,nk,nj) \
60 (N(ni,nk,nj+2) \
61 + N(ni-a ,nk,nj+1)+ N(ni+1-a,nk,nj+1) \
62 + N(ni-1,nk,nj ) + N(ni,nk,nj ) + N(ni+1,nk,nj )\
63 + N(ni-a ,nk,nj-1)+ N(ni+1-a,nk,nj-1)\
64 + N(ni,nk,nj-2) \
65 ) / 9
67 ! Average to C points from N points without assignment, no vertical levels:
68 #define NGRAB2D(N,ni,nj) \
69 (N(ni,nj+2) \
70 + N(ni-a ,nj+1)+ N(ni+1-a,nj+1)\
71 + N(ni-1,nj ) + N(ni,nj ) + N(ni+1,nj )\
72 + N(ni-a ,nj-1)+ N(ni+1-a,nj-1)\
73 + N(ni,nj-2) \
74 ) / 9
76 ! Copying from N array to I array:
77 #define N2ICOPY(C,N,i,j,k,ni,nj) \
78 C(i,j,k)=( N(ni,nj+2,k) \
79 + N(ni-a ,nj+1,k)+ N(ni+1-a,nj+1,k) \
80 + N(ni-1,nj ,k)+ N(ni,nj ,k) + N(ni+1,nj ,k) \
81 + N(ni-a ,nj-1,k)+ N(ni+1-a,nj-1,k) \
82 + N(ni,nj-2,k) \
83 ) / 9
85 ! Maximum value from N array to C array:
86 #define N2I_SET_MAX(C,N,i,j,k,ni,nj) \
87 C(i,j,k)=max(C(i,j,k),max(N(ni,nj+2,k), \
88 max(N(ni-1,nj ,k),max(N(ni,nj ,k), max(N(ni+1,nj ,k),\
89 max(N(ni-a ,nj-1,k),max(N(ni+1-a,nj-1,k),\
90 N(ni,nj-2,k) )))))))
92 ! Maximum value from N array to C array:
93 #define N2I_SET_MAX_IJ(C,N,i,j,ni,nj) \
94 C(i,j)=max(C(i,j), max(N(ni,nj+2),\
95 max(N(ni-1,nj ),max(N(ni,nj ), max(N(ni+1,nj ), \
96 max(N(ni-a ,nj-1),max(N(ni+1-a,nj-1), \
97 N(ni,nj-2) )))))))
100 module module_interp_nmm
101 use module_model_constants, only: g, R_D, p608
102 implicit none
104 private
105 public :: find_kpres
107 public :: nmm_interp_pd, nmm_keep_pd, nmm_method_linear
109 public :: c2b_fulldom, c2n_fulldom, n2c_fulldom
110 public :: n2c_max2d, n2c_max3d
111 public :: c2b_fulldom_new, n2c_fulldom_new
112 public :: c2b_mass, c2n_mass, n2c_mass
113 public :: c2n_massikj, n2c_massikj
114 public :: c2b_copy3d, c2n_copy3d, n2c_copy3d
116 public :: c2b_copy2d, c2n_copy2d, n2c_copy2d, c2n_copy2d_nomask
117 public :: c2b_near2d, c2n_near2d, n2c_near2d
118 public :: c2b_inear2d, c2n_inear2d, n2c_inear2d
120 public :: c2n_near3dikj, c2n_sst, c2n_near3d
122 ! unimplemented: public :: nmm_method_quadratic, nmm_method_spline_log, nmm_method_copy, nmm_method_spline
124 ! integer, parameter :: nmm_method_copy = 5
125 ! integer, parameter :: nmm_method_spline_log = 4
126 ! integer, parameter :: nmm_method_spline = 3
127 ! integer, parameter :: nmm_method_quadratic = 2
129 integer, parameter :: nmm_method_linear = 1
131 integer, parameter :: nmm_interp_pd = 1
132 integer, parameter :: nmm_keep_pd = 0
134 ! Constants from the original base_state_parent:
135 REAL, PARAMETER :: LAPSR=6.5E-3, GI=1./G,D608=0.608
136 REAL, PARAMETER :: COEF3=287.05*GI*LAPSR, COEF2=-1./COEF3
137 REAL, PARAMETER :: TRG=2.0*R_D*GI,LAPSI=1.0/LAPSR
138 REAL, PARAMETER :: P_REF=103000.
140 integer, parameter :: EConst=0, ECopy=1, EExtrap=2 ! MUST match module_dm
142 contains
143 ! ********************************************************************
144 ! subs *_NEAR2D -- horizontally interpolate via nearest neighbor
145 ! method, but do not vertically interpolate. Only handles 2D H grid
146 ! arrays.
147 ! ********************************************************************
149 subroutine c2b_near2d (inear,jnear,cfield, &
150 fbxs,fbxe,fbys,fbye, &
151 cims, cime, cjms, cjme, &
152 nids, nide, njds, njde, &
153 nims, nime, njms, njme, &
154 nits, nite, njts, njte)
155 implicit none
156 integer, parameter :: bdyw = 1
158 integer, intent(in):: &
159 cims, cime, cjms, cjme, &
160 nids, nide, njds, njde, &
161 nims, nime, njms, njme, &
162 nits, nite, njts, njte
164 integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
165 real, intent(in) :: cfield(cims:cime,cjms:cjme)
167 ! Output field boundary info:
168 real,dimension(nims:nime,bdyw) :: fbys,fbye
169 real,dimension(njms:njme,bdyw) :: fbxs,fbxe
171 integer :: j,i,a,nx,nz,ck,k,j1,j2,add
172 real :: weight
174 nx=min(nide-1,nite)-max(nids,nits)+1
176 j1=max(njts-1,njds)
177 if(mod(j1,2)/=1) j1=j1+1
179 j2=min(njte+1,njde-1)
180 if(mod(j2,2)/=1) j2=j2-1
182 if(nits==1) then
183 i=1
184 do j=j1,j2,2
185 fbxs(j,1)=cfield(inear(i,j),jnear(i,j))
186 enddo
187 endif
188 if(nite>=nide-1) then
189 i=nide-1
190 do j=j1,j2,2
191 fbxe(j,1)=cfield(inear(i,j),jnear(i,j))
192 enddo
193 endif
194 if(njts==1) then
195 j=1
196 do i=max(nits-1,nids),min(nite+1,nide-1)
197 fbys(i,1)=cfield(inear(i,j),jnear(i,j))
198 enddo
199 endif
200 if(njte>=njde-1) then
201 j=njde-1
202 do i=max(nits-1,nids),min(nite+1,nide-1)
203 fbye(i,1)=cfield(inear(i,j),jnear(i,j))
204 enddo
205 endif
206 end subroutine c2b_near2d
208 subroutine n2c_near2d (&
209 cfield,nfield,ipos,jpos, &
210 cids, cide, cjds, cjde, &
211 cims, cime, cjms, cjme, &
212 cits, cite, cjts, cjte, &
213 nids, nide, njds, njde, &
214 nims, nime, njms, njme, &
215 nits, nite, njts, njte)
216 implicit none
217 integer, intent(in) :: &
218 cids, cide, cjds, cjde, &
219 cims, cime, cjms, cjme, &
220 cits, cite, cjts, cjte, &
221 nids, nide, njds, njde, &
222 nims, nime, njms, njme, &
223 nits, nite, njts, njte, &
224 ipos,jpos
225 real, intent(out) :: cfield(cims:cime,cjms:cjme)
226 real, intent(in) :: nfield(nims:nime,njms:njme)
228 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
229 integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
231 jstart=MAX(jpos+1,cjts)
232 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
234 bigj: do j=jstart,jend
235 a=mod(j,2)
236 istart=MAX(ipos+a,cits)
237 iend=MIN(ipos+(nide-nids)/nri-1,cite)
238 nj = (j-jpos)*nrj + 1
240 iloop: do i=istart,iend
241 ni = (i-ipos)*nri + 2 - a
242 cfield(i,j)=nfield(ni,nj)
243 enddo iloop
244 enddo bigj
245 end subroutine n2c_near2d
247 subroutine c2n_near2d (inear,jnear, &
248 cfield,nfield,imask, &
249 cims, cime, cjms, cjme, &
250 nids, nide, njds, njde, &
251 nims, nime, njms, njme, &
252 nits, nite, njts, njte)
253 implicit none
255 integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
257 integer, intent(in):: cims, cime, cjms, cjme, &
258 nids, nide, njds, njde, &
259 nims, nime, njms, njme, &
260 nits, nite, njts, njte
261 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
262 real, intent(in) :: cfield(cims:cime,cjms:cjme)
263 real, intent(out) :: nfield(nims:nime,njms:njme)
265 integer :: j,i,a,nx
267 nx=min(nide-1,nite)-max(nids,nits)+1
269 bigj: do j=max(njds,njts),min(njde-1,njte)
270 interploop: do i=max(nids,nits),min(nide-1,nite)
271 if(imask(i,j)/=0) cycle interploop
272 nfield(i,j)=cfield(inear(i,j),jnear(i,j))
273 enddo interploop
274 end do bigj
275 end subroutine c2n_near2d
277 subroutine c2n_near3d (inear,jnear, &
278 cfield,nfield,imask, &
279 cims, cime, cjms, cjme, ckms, ckme, &
280 nids, nide, njds, njde, nkds, nkde, &
281 nims, nime, njms, njme, nkms, nkme, &
282 nits, nite, njts, njte, nkts, nkte)
283 implicit none
285 integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
287 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
288 nids, nide, njds, njde, nkds, nkde, &
289 nims, nime, njms, njme, nkms, nkme, &
290 nits, nite, njts, njte, nkts, nkte
291 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
292 real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
293 real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
295 integer :: j,i,a,nx,k
297 nx=min(nide-1,nite)-max(nids,nits)+1
299 bigk: do k=nkds,nkde
300 bigj: do j=max(njds,njts),min(njde-1,njte)
301 interploop: do i=max(nids,nits),min(nide-1,nite)
302 if(imask(i,j)/=0) cycle interploop
303 nfield(i,j,k)=cfield(inear(i,j),jnear(i,j),k)
304 enddo interploop
305 end do bigj
306 enddo bigk
307 end subroutine c2n_near3d
309 subroutine c2n_near3dikj (inear,jnear, &
310 cfield,nfield,imask, &
311 cims, cime, cjms, cjme, ckms, ckme, &
312 nids, nide, njds, njde, nkds, nkde, &
313 nims, nime, njms, njme, nkms, nkme, &
314 nits, nite, njts, njte, nkts, nkte)
315 implicit none
317 integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
319 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
320 nids, nide, njds, njde, nkds, nkde, &
321 nims, nime, njms, njme, nkms, nkme, &
322 nits, nite, njts, njte, nkts, nkte
323 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
324 real, intent(in) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
325 real, intent(out) :: nfield(nims:nime,nkms:nkme,njms:njme)
327 integer :: j,i,a,nx,k
329 nx=min(nide-1,nite)-max(nids,nits)+1
331 bigj: do j=max(njds,njts),min(njde-1,njte)
332 bigk: do k=nkds,nkde
333 interploop: do i=max(nids,nits),min(nide-1,nite)
334 if(imask(i,j)/=0) cycle interploop
335 nfield(i,k,j)=cfield(inear(i,j),k,jnear(i,j))
336 enddo interploop
337 enddo bigk
338 end do bigj
339 end subroutine c2n_near3dikj
341 subroutine c2n_sst (inear,jnear, &
342 cfield,nfield, &
343 cims, cime, cjms, cjme, &
344 nids, nide, njds, njde, &
345 nims, nime, njms, njme, &
346 nits, nite, njts, njte)
347 implicit none
349 integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
351 integer, intent(in):: cims, cime, cjms, cjme, &
352 nids, nide, njds, njde, &
353 nims, nime, njms, njme, &
354 nits, nite, njts, njte
355 real, intent(in) :: cfield(cims:cime,cjms:cjme)
356 real, intent(out) :: nfield(nims:nime,njms:njme)
358 integer :: j,i,a,nx
360 nx=min(nide-1,nite)-max(nids,nits)+1
362 bigj: do j=max(njds,njts),min(njde-1,njte)
363 interploop: do i=max(nids,nits),min(nide-1,nite)
364 nfield(i,j)=cfield(inear(i,j),jnear(i,j))
365 enddo interploop
366 end do bigj
367 end subroutine c2n_sst
369 ! ********************************************************************
370 ! subs *_INEAR2D -- horizontally interpolate integers via nearest
371 ! neighbor method, but do not vertically interpolate. Only handles
372 ! 2D H grid integer arrays.
373 ! ********************************************************************
375 subroutine c2b_inear2d (inear,jnear,cfield, &
376 fbxs,fbxe,fbys,fbye, &
377 cims, cime, cjms, cjme, &
378 nids, nide, njds, njde, &
379 nims, nime, njms, njme, &
380 nits, nite, njts, njte)
381 implicit none
382 integer, parameter :: bdyw = 1
384 integer, intent(in):: &
385 cims, cime, cjms, cjme, &
386 nids, nide, njds, njde, &
387 nims, nime, njms, njme, &
388 nits, nite, njts, njte
390 integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
391 integer, intent(in) :: cfield(cims:cime,cjms:cjme)
393 ! Output field boundary info:
394 integer,dimension(nims:nime,bdyw) :: fbys,fbye
395 integer,dimension(njms:njme,bdyw) :: fbxs,fbxe
397 integer :: j,i,a,nx,nz,ck,k,j1,j2,add
399 nx=min(nide-1,nite)-max(nids,nits)+1
401 j1=max(njts-1,njds)
402 if(mod(j1,2)/=1) j1=j1+1
404 j2=min(njte+1,njde-1)
405 if(mod(j2,2)/=1) j2=j2-1
407 if(nits==1) then
408 i=1
409 do j=j1,j2,2
410 fbxs(j,1)=cfield(inear(i,j),jnear(i,j))
411 enddo
412 endif
413 if(nite>=nide-1) then
414 i=nide-1
415 do j=j1,j2,2
416 fbxe(j,1)=cfield(inear(i,j),jnear(i,j))
417 enddo
418 endif
419 if(njts==1) then
420 j=1
421 do i=max(nits-1,nids),min(nite+1,nide-1)
422 fbys(i,1)=cfield(inear(i,j),jnear(i,j))
423 enddo
424 endif
425 if(njte>=njde-1) then
426 j=njde-1
427 do i=max(nits-1,nids),min(nite+1,nide-1)
428 fbye(i,1)=cfield(inear(i,j),jnear(i,j))
429 enddo
430 endif
431 end subroutine c2b_inear2d
433 subroutine n2c_inear2d (&
434 cfield,nfield,ipos,jpos, &
435 cids, cide, cjds, cjde, &
436 cims, cime, cjms, cjme, &
437 cits, cite, cjts, cjte, &
438 nids, nide, njds, njde, &
439 nims, nime, njms, njme, &
440 nits, nite, njts, njte)
441 implicit none
442 integer, intent(in) :: &
443 cids, cide, cjds, cjde, &
444 cims, cime, cjms, cjme, &
445 cits, cite, cjts, cjte, &
446 nids, nide, njds, njde, &
447 nims, nime, njms, njme, &
448 nits, nite, njts, njte, &
449 ipos,jpos
450 integer, intent(out) :: cfield(cims:cime,cjms:cjme)
451 integer, intent(in) :: nfield(nims:nime,njms:njme)
453 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
454 integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
456 jstart=MAX(jpos+1,cjts)
457 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
459 bigj: do j=jstart,jend
460 a=mod(j,2)
461 istart=MAX(ipos+a,cits)
462 iend=MIN(ipos+(nide-nids)/nri-1,cite)
463 nj = (j-jpos)*nrj + 1
465 iloop: do i=istart,iend
466 ni = (i-ipos)*nri + 2 - a
467 cfield(i,j)=nfield(ni,nj)
468 enddo iloop
469 enddo bigj
470 end subroutine n2c_inear2d
472 subroutine c2n_inear2d (inear,jnear, &
473 cfield,nfield,imask, &
474 cims, cime, cjms, cjme, &
475 nids, nide, njds, njde, &
476 nims, nime, njms, njme, &
477 nits, nite, njts, njte)
478 implicit none
480 integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
482 integer, intent(in):: cims, cime, cjms, cjme, &
483 nids, nide, njds, njde, &
484 nims, nime, njms, njme, &
485 nits, nite, njts, njte
486 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
487 integer, intent(in) :: cfield(cims:cime,cjms:cjme)
488 integer, intent(out) :: nfield(nims:nime,njms:njme)
490 integer :: j,i,a,nx
492 nx=min(nide-1,nite)-max(nids,nits)+1
494 bigj: do j=max(njds,njts),min(njde-1,njte)
495 interploop: do i=max(nids,nits),min(nide-1,nite)
496 if(imask(i,j)/=0) cycle interploop
497 nfield(i,j)=cfield(inear(i,j),jnear(i,j))
498 enddo interploop
499 end do bigj
500 end subroutine c2n_inear2d
502 ! ********************************************************************
503 ! subs *_COPY2D -- horizontally interpolate but do not vertically
504 ! interpolate. Only handles 2D arrays (H or V grid).
505 ! ********************************************************************
507 subroutine c2b_copy2d (II,JJ,W1,W2,W3,W4,cfield, &
508 fbxs,fbxe,fbys,fbye, &
509 cims, cime, cjms, cjme, &
510 nids, nide, njds, njde, &
511 nims, nime, njms, njme, &
512 nits, nite, njts, njte, hgrid)
513 implicit none
514 integer, parameter :: bdyw = 1
516 integer, intent(in):: &
517 cims, cime, cjms, cjme, &
518 nids, nide, njds, njde, &
519 nims, nime, njms, njme, &
520 nits, nite, njts, njte
521 logical, intent(in) :: hgrid
522 real, intent(in), dimension(nims:nime,njms:njme) :: &
523 W1,W2,W3,W4
524 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
525 real, intent(in) :: cfield(cims:cime,cjms:cjme)
527 ! Output field boundary info:
528 real,dimension(nims:nime,bdyw) :: fbys,fbye
529 real,dimension(njms:njme,bdyw) :: fbxs,fbxe
531 integer :: j,i,a,nx,nz,ck,k,j1,j2,add
532 real :: weight
534 nx=min(nide-1,nite)-max(nids,nits)+1
536 if(hgrid) then
537 a=1
538 else
539 a=0
540 endif
542 j1=max(njts-1,njds)
543 if(mod(j1,2)/=a) j1=j1+1
545 j2=min(njte+1,njde-1)
546 if(mod(j2,2)/=a) j2=j2-1
548 if(nits==1) then
549 i=1
550 do j=j1,j2,2
551 if(hgrid) then
552 a=1-mod(JJ(i,j),2)
553 else
554 a=mod(JJ(i,j),2)
555 endif
556 fbxs(j,1)=ICOPY2D(cfield,i,j)
557 enddo
558 endif
559 if(nite>=nide-1) then
560 i=nide-1
561 do j=j1,j2,2
562 if(hgrid) then
563 a=1-mod(JJ(i,j),2)
564 else
565 a=mod(JJ(i,j),2)
566 endif
567 fbxe(j,1)=ICOPY2D(cfield,i,j)
568 enddo
569 endif
570 if(njts==1) then
571 j=1
572 do i=max(nits-1,nids),min(nite+1,nide-1)
573 if(hgrid) then
574 a=1-mod(JJ(i,j),2)
575 else
576 a=mod(JJ(i,j),2)
577 endif
578 fbys(i,1)=ICOPY2D(cfield,i,j)
579 enddo
580 endif
581 if(njte>=njde-1) then
582 j=njde-1
583 do i=max(nits-1,nids),min(nite+1,nide-1)
584 if(hgrid) then
585 a=1-mod(JJ(i,j),2)
586 else
587 a=mod(JJ(i,j),2)
588 endif
589 fbye(i,1)=ICOPY2D(cfield,i,j)
590 enddo
591 endif
592 end subroutine c2b_copy2d
593 subroutine n2c_copy2d (&
594 cfield,nfield,ipos,jpos, &
595 cids, cide, cjds, cjde, &
596 cims, cime, cjms, cjme, &
597 cits, cite, cjts, cjte, &
598 nids, nide, njds, njde, &
599 nims, nime, njms, njme, &
600 nits, nite, njts, njte, &
601 hgrid)
602 implicit none
603 integer, intent(in) :: &
604 cids, cide, cjds, cjde, &
605 cims, cime, cjms, cjme, &
606 cits, cite, cjts, cjte, &
607 nids, nide, njds, njde, &
608 nims, nime, njms, njme, &
609 nits, nite, njts, njte, &
610 ipos,jpos
611 real, intent(out) :: cfield(cims:cime,cjms:cjme)
612 real, intent(in) :: nfield(nims:nime,njms:njme)
613 logical, intent(in) :: hgrid
615 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
616 integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
618 jstart=MAX(jpos+1,cjts)
619 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
621 bigj: do j=jstart,jend
622 a=mod(j,2)
623 if(.not.hgrid) then
624 a=1-a
625 endif
626 istart=MAX(ipos+a,cits)
627 iend=MIN(ipos+(nide-nids)/nri-1,cite)
628 nj = (j-jpos)*nrj + 1
630 iloop: do i=istart,iend
631 ni = (i-ipos)*nri + 2 - a
632 cfield(i,j)=NGRAB2D(nfield,ni,nj)
633 enddo iloop
634 enddo bigj
635 end subroutine n2c_copy2d
637 subroutine n2c_max2d (&
638 cfield,nfield,ipos,jpos, &
639 cids, cide, cjds, cjde, &
640 cims, cime, cjms, cjme, &
641 cits, cite, cjts, cjte, &
642 nids, nide, njds, njde, &
643 nims, nime, njms, njme, &
644 nits, nite, njts, njte, &
645 hgrid)
646 implicit none
647 integer, intent(in) :: &
648 cids, cide, cjds, cjde, &
649 cims, cime, cjms, cjme, &
650 cits, cite, cjts, cjte, &
651 nids, nide, njds, njde, &
652 nims, nime, njms, njme, &
653 nits, nite, njts, njte, &
654 ipos,jpos
655 real, intent(out) :: cfield(cims:cime,cjms:cjme)
656 real, intent(in) :: nfield(nims:nime,njms:njme)
657 logical, intent(in) :: hgrid
659 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
660 integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
662 jstart=MAX(jpos+1,cjts)
663 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
665 bigj: do j=jstart,jend
666 a=mod(j,2)
667 if(.not.hgrid) then
668 a=1-a
669 endif
670 istart=MAX(ipos+a,cits)
671 iend=MIN(ipos+(nide-nids)/nri-1,cite)
672 nj = (j-jpos)*nrj + 1
674 iloop: do i=istart,iend
675 ni = (i-ipos)*nri + 2 - a
676 N2I_SET_MAX_IJ(cfield,nfield,i,j,ni,nj)
677 enddo iloop
678 enddo bigj
679 end subroutine n2c_max2d
681 subroutine c2n_copy2d (II,JJ,W1,W2,W3,W4, &
682 cfield,nfield,imask, &
683 cims, cime, cjms, cjme, &
684 nids, nide, njds, njde, &
685 nims, nime, njms, njme, &
686 nits, nite, njts, njte, hgrid)
687 implicit none
688 logical, intent(in) :: hgrid
689 real, intent(in), dimension(nims:nime,njms:njme) :: &
690 W1,W2,W3,W4
691 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
693 integer, intent(in):: cims, cime, cjms, cjme, &
694 nids, nide, njds, njde, &
695 nims, nime, njms, njme, &
696 nits, nite, njts, njte
697 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
698 real, intent(in) :: cfield(cims:cime,cjms:cjme)
699 real, intent(out) :: nfield(nims:nime,njms:njme)
701 integer :: j,i,a,nx
703 nx=min(nide-1,nite)-max(nids,nits)+1
705 bigj: do j=max(njds,njts),min(njde-1,njte)
706 interploop: do i=max(nids,nits),min(nide-1,nite)
707 if(imask(i,j)/=0) cycle interploop
708 if(hgrid) then
709 a=1-mod(JJ(i,j),2)
710 else
711 a=mod(JJ(i,j),2)
712 endif
713 nfield(i,j)=ICOPY2D(cfield,i,j)
714 enddo interploop
715 end do bigj
716 end subroutine c2n_copy2d
718 subroutine c2n_copy2d_nomask (II,JJ,W1,W2,W3,W4, &
719 cfield,nfield, &
720 cims, cime, cjms, cjme, &
721 nids, nide, njds, njde, &
722 nims, nime, njms, njme, &
723 nits, nite, njts, njte, hgrid)
724 implicit none
725 logical, intent(in) :: hgrid
726 real, intent(in), dimension(nims:nime,njms:njme) :: &
727 W1,W2,W3,W4
728 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
730 integer, intent(in):: cims, cime, cjms, cjme, &
731 nids, nide, njds, njde, &
732 nims, nime, njms, njme, &
733 nits, nite, njts, njte
734 real, intent(in) :: cfield(cims:cime,cjms:cjme)
735 real, intent(out) :: nfield(nims:nime,njms:njme)
737 integer :: j,i,a,nx
739 nx=min(nide-1,nite)-max(nids,nits)+1
741 bigj: do j=max(njds,njts),min(njde-1,njte)
742 interploop: do i=max(nids,nits),min(nide-1,nite)
743 if(hgrid) then
744 a=1-mod(JJ(i,j),2)
745 else
746 a=mod(JJ(i,j),2)
747 endif
748 nfield(i,j)=ICOPY2D(cfield,i,j)
749 enddo interploop
750 end do bigj
751 end subroutine c2n_copy2d_nomask
754 ! ********************************************************************
755 ! subs *_COPY3D -- horizontally interpolate but do not vertically
756 ! interpolate
757 ! ********************************************************************
759 subroutine c2b_copy3d (II,JJ,W1,W2,W3,W4,cfield, &
760 fbxs,fbxe,fbys,fbye, &
761 cims, cime, cjms, cjme, ckms, ckme, &
762 nids, nide, njds, njde, nkds, nkde, &
763 nims, nime, njms, njme, nkms, nkme, &
764 nits, nite, njts, njte, nkts, nkte, hgrid)
765 implicit none
766 integer, parameter :: bdyw = 1
767 logical, intent(in) :: hgrid
768 integer, intent(in):: &
769 cims, cime, cjms, cjme, ckms, ckme, &
770 nids, nide, njds, njde, nkds, nkde, &
771 nims, nime, njms, njme, nkms, nkme, &
772 nits, nite, njts, njte, nkts, nkte
774 real, intent(in), dimension(nims:nime,njms:njme) :: &
775 W1,W2,W3,W4
776 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
777 real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
779 ! Output field boundary info:
780 real,dimension(nims:nime,nkms:nkme,bdyw) :: fbys,fbye
781 real,dimension(njms:njme,nkms:nkme,bdyw) :: fbxs,fbxe
783 integer :: j,i,a,nx,nz,ck,k,j1,j2
784 real :: weight
786 nx=min(nide-1,nite)-max(nids,nits)+1
787 nz=nkde-nkds+1
789 j1=max(njts-1,njds)
790 if(mod(j1,2)/=0) j1=j1+1
792 j2=min(njte+1,njde-1)
793 if(mod(j2,2)/=0) j2=j2-1
795 if(nits==1) then
796 i=1
797 do j=j1,j2,2
798 if(hgrid) then
799 a=1-mod(JJ(i,j),2)
800 else
801 a=mod(JJ(i,j),2)
802 endif
803 kloop1: do k=nkds,min(nkde,nkte)
804 fbxs(j,k,1)=ICOPY(cfield,i,j,k)
805 enddo kloop1
806 enddo
807 endif
808 if(nite>=nide-1) then
809 i=nide-1
810 do j=j1,j2,2
811 if(hgrid) then
812 a=1-mod(JJ(i,j),2)
813 else
814 a=mod(JJ(i,j),2)
815 endif
816 kloop2: do k=nkds,min(nkde,nkte)
817 fbxe(j,k,1)=ICOPY(cfield,i,j,k)
818 enddo kloop2
819 enddo
820 endif
821 if(njts==1) then
822 j=1
823 do i=max(nits-1,nids),min(nite+1,nide-1)
824 if(hgrid) then
825 a=1-mod(JJ(i,j),2)
826 else
827 a=mod(JJ(i,j),2)
828 endif
829 kloop3: do k=nkts,min(nkde,nkte)
830 fbys(i,k,1)=ICOPY(cfield,i,j,k)
831 enddo kloop3
832 enddo
833 endif
834 if(njte>=njde-1) then
835 j=njde-1
836 do i=max(nits-1,nids),min(nite+1,nide-1)
837 if(hgrid) then
838 a=1-mod(JJ(i,j),2)
839 else
840 a=mod(JJ(i,j),2)
841 endif
842 kloop4: do k=nkts,min(nkde,nkte)
843 fbye(i,k,1)=ICOPY(cfield,i,j,k)
844 enddo kloop4
845 enddo
846 endif
847 end subroutine c2b_copy3d
849 subroutine n2c_copy3d (&
850 cfield,nfield,ipos,jpos, &
851 cids, cide, cjds, cjde, ckds, ckde, &
852 cims, cime, cjms, cjme, ckms, ckme, &
853 cits, cite, cjts, cjte, ckts, ckte, &
854 nids, nide, njds, njde, nkds, nkde, &
855 nims, nime, njms, njme, nkms, nkme, &
856 nits, nite, njts, njte, nkts, nkte, &
857 hgrid)
858 implicit none
859 logical, intent(in) :: hgrid
860 integer, intent(in) :: &
861 cids, cide, cjds, cjde, ckds, ckde, &
862 cims, cime, cjms, cjme, ckms, ckme, &
863 cits, cite, cjts, cjte, ckts, ckte, &
864 nids, nide, njds, njde, nkds, nkde, &
865 nims, nime, njms, njme, nkms, nkme, &
866 nits, nite, njts, njte, nkts, nkte, &
867 ipos,jpos
868 real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
869 real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
871 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
872 integer :: nx,nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
873 real :: weight
875 nx=min(nide-2,nite)-max(nids+1,nits)+1
876 nz=nkde-nkds+1
878 jstart=MAX(jpos+1,cjts)
879 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
881 bigj: do j=jstart,jend
882 if(hgrid) then
883 a=mod(j,2)
884 else
885 a=1-mod(j,2)
886 endif
887 istart=MAX(ipos+a,cits)
888 iend=MIN(ipos+(nide-nids)/nri-1,cite)
889 nj = (j-jpos)*nrj + 1
891 iloop: do i=istart,iend
892 do k=nkds,nkde
893 ni = (i-ipos)*nri + 2 - a
894 cfield(i,j,k)=NGRAB(nfield,ni,nj,k)
895 enddo
896 enddo iloop
897 enddo bigj
898 end subroutine n2c_copy3d
900 subroutine n2c_max3d (&
901 cfield,nfield,ipos,jpos, &
902 cids, cide, cjds, cjde, ckds, ckde, &
903 cims, cime, cjms, cjme, ckms, ckme, &
904 cits, cite, cjts, cjte, ckts, ckte, &
905 nids, nide, njds, njde, nkds, nkde, &
906 nims, nime, njms, njme, nkms, nkme, &
907 nits, nite, njts, njte, nkts, nkte, &
908 hgrid)
909 implicit none
910 logical, intent(in) :: hgrid
911 integer, intent(in) :: &
912 cids, cide, cjds, cjde, ckds, ckde, &
913 cims, cime, cjms, cjme, ckms, ckme, &
914 cits, cite, cjts, cjte, ckts, ckte, &
915 nids, nide, njds, njde, nkds, nkde, &
916 nims, nime, njms, njme, nkms, nkme, &
917 nits, nite, njts, njte, nkts, nkte, &
918 ipos,jpos
919 real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
920 real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
922 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
923 integer :: nx,nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
924 real :: weight
926 nx=min(nide-2,nite)-max(nids+1,nits)+1
927 nz=nkde-nkds+1
929 jstart=MAX(jpos+1,cjts)
930 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
932 bigj: do j=jstart,jend
933 if(hgrid) then
934 a=mod(j,2)
935 else
936 a=1-mod(j,2)
937 endif
938 istart=MAX(ipos+a,cits)
939 iend=MIN(ipos+(nide-nids)/nri-1,cite)
940 nj = (j-jpos)*nrj + 1
942 iloop: do i=istart,iend
943 do k=nkds,nkde
944 ni = (i-ipos)*nri + 2 - a
945 N2I_SET_MAX(cfield,nfield,i,j,k,ni,nj)
946 enddo
947 enddo iloop
948 enddo bigj
949 end subroutine n2c_max3d
951 subroutine c2n_copy3d (II,JJ,W1,W2,W3,W4, &
952 cfield,nfield,imask, &
953 cims, cime, cjms, cjme, ckms, ckme, &
954 nids, nide, njds, njde, nkds, nkde, &
955 nims, nime, njms, njme, nkms, nkme, &
956 nits, nite, njts, njte, nkts, nkte, &
957 hgrid)
958 implicit none
959 real, intent(in), dimension(nims:nime,njms:njme) :: &
960 W1,W2,W3,W4
961 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
962 logical, intent(in) :: hgrid
963 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
964 nids, nide, njds, njde, nkds, nkde, &
965 nims, nime, njms, njme, nkms, nkme, &
966 nits, nite, njts, njte, nkts, nkte
967 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
968 real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
969 real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
971 integer :: j,i,a,nx,nz,k
973 nx=min(nide-1,nite)-max(nids,nits)+1
974 nz=nkde-nkds+1
976 bigj: do j=max(njds,njts),min(njde-1,njte)
977 interploop: do i=max(nids,nits),min(nide-1,nite)
978 if(imask(i,j)/=0) cycle interploop
979 kinterploop: do k=nkds,nkde
980 if(hgrid) then
981 a=1-mod(JJ(i,j),2)
982 else
983 a=mod(JJ(i,j),2)
984 endif
985 nfield(i,j,k) = ICOPY(cfield,i,j,k)
986 enddo kinterploop
987 enddo interploop
988 end do bigj
989 end subroutine c2n_copy3d
991 ! ********************************************************************
992 ! subs *_MASS -- horizontally and vertically interpolate. Vertical
993 ! interpolation uses the iinfo and winfo arrays (and equivalent
994 ! boundary arrays) generated in the *_fulldom subroutines.
995 ! ********************************************************************
997 subroutine c2b_mass (II,JJ,W1,W2,W3,W4,cfield, &
998 ibxs,ibxe,ibys,ibye, &
999 wbxs,wbxe,wbys,wbye, &
1000 fbxs,fbxe,fbys,fbye, &
1001 emethod,evalue, &
1002 cims, cime, cjms, cjme, ckms, ckme, &
1003 nids, nide, njds, njde, nkds, nkde, &
1004 nims, nime, njms, njme, nkms, nkme, &
1005 nits, nite, njts, njte, nkts, nkte)
1006 use module_interp_store, only: kpres
1007 implicit none
1008 integer, parameter :: bdyw = 1
1010 integer, intent(in):: &
1011 cims, cime, cjms, cjme, ckms, ckme, &
1012 nids, nide, njds, njde, nkds, nkde, &
1013 nims, nime, njms, njme, nkms, nkme, &
1014 nits, nite, njts, njte, nkts, nkte
1016 integer, intent(in) :: emethod ! extrapolation method
1017 real, intent(in) :: evalue
1018 real, intent(in), dimension(nims:nime,njms:njme) :: &
1019 W1,W2,W3,W4
1020 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1021 real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
1023 ! Output field boundary info:
1024 real,dimension(nims:nime,nkms:nkme,bdyw) :: fbys,fbye
1025 real,dimension(njms:njme,nkms:nkme,bdyw) :: fbxs,fbxe
1027 ! Weights and indices:
1028 real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
1029 real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
1030 integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
1031 integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
1033 integer :: j,i,a,nx,nz,ck,k,j1,j2
1034 real :: weight
1036 nx=min(nide-1,nite)-max(nids,nits)+1
1037 nz=nkde-nkds+1
1039 j1=max(njts-1,njds)
1040 if(mod(j1,2)/=1) j1=j1+1
1042 j2=min(njte+1,njde-1)
1043 if(mod(j2,2)/=1) j2=j2-1
1045 if(nits==1) then
1046 i=1
1047 do j=j1,j2,2
1048 a=1-mod(JJ(i,j),2)
1049 kcopy1: do k=min(nkde,nkte),kpres+1,-1
1050 fbxs(j,k,1) = ICOPY(cfield,i,j,k)
1051 enddo kcopy1
1052 kloop1: do k=kpres,nkds,-1
1053 weight=wbxs(j,k,1)
1054 ck=ibxs(j,k,1)
1056 if(ck>1) then
1057 fbxs(j,k,1) = &
1058 ICOPY(cfield,i,j,ck) * weight + &
1059 ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1060 else
1061 ! Extrapolate
1062 if(emethod==EConst) then
1063 fbxs(j,k,1)=evalue
1064 else if(emethod==ECopy) then
1065 fbxs(j,k,1)=ICOPY(cfield,i,j,1)
1066 else ! assume 2: linear extrap
1067 fbxs(j,k,1)=evalue * weight + &
1068 ICOPY(cfield,i,j,1) * (1.0-weight)
1069 endif
1070 endif
1071 enddo kloop1
1072 enddo
1073 endif
1074 if(nite>=nide-1) then
1075 i=nide-1
1076 do j=j1,j2,2
1077 a=1-mod(JJ(i,j),2)
1078 kcopy2: do k=min(nkde,nkte),kpres+1,-1
1079 fbxe(j,k,1)=ICOPY(cfield,i,j,k)
1080 enddo kcopy2
1081 kloop2: do k=kpres,nkds,-1
1082 weight=wbxe(j,k,1)
1083 ck=ibxe(j,k,1)
1085 if(ck>1) then
1086 fbxe(j,k,1) = &
1087 ICOPY(cfield,i,j,ck) * weight + &
1088 ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1089 else
1090 ! Extrapolate
1091 if(emethod==EConst) then
1092 fbxe(j,k,1)=evalue
1093 else if(emethod==ECopy) then
1094 fbxe(j,k,1)=ICOPY(cfield,i,j,1)
1095 else ! assume 2: linear extrap
1096 fbxe(j,k,1)=evalue * weight + &
1097 ICOPY(cfield,i,j,1) * (1.0-weight)
1098 endif
1099 endif
1100 enddo kloop2
1101 enddo
1102 endif
1103 if(njts==1) then
1104 j=1
1105 do i=max(nits-1,nids),min(nite+1,nide-1)
1106 a=1-mod(JJ(i,j),2)
1107 kcopy3: do k=min(nkde,nkte),kpres+1,-1
1108 fbys(i,k,1) = ICOPY(cfield,i,j,k)
1109 enddo kcopy3
1110 kloop3: do k=kpres,nkds,-1
1111 weight=wbys(i,k,1)
1112 ck=ibys(i,k,1)
1114 if(ck>1) then
1115 fbys(i,k,1) = &
1116 ICOPY(cfield,i,j,ck) * weight + &
1117 ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1118 else
1119 ! Extrapolate
1120 if(emethod==EConst) then
1121 fbys(i,k,1)=evalue
1122 else if(emethod==ECopy) then
1123 fbys(i,k,1)=ICOPY(cfield,i,j,1)
1124 else ! assume 2: linear extrap
1125 fbys(i,k,1)=evalue*weight + &
1126 ICOPY(cfield,i,j,1)*(1.0-weight)
1127 endif
1128 endif
1129 enddo kloop3
1130 enddo
1131 endif
1132 if(njte>=njde-1) then
1133 j=njde-1
1134 do i=max(nits-1,nids),min(nite+1,nide-1)
1135 a=1-mod(JJ(i,j),2)
1136 kcopy4: do k=min(nkde,nkte),kpres+1,-1
1137 fbye(i,k,1) = ICOPY(cfield,i,j,k)
1138 enddo kcopy4
1139 kloop4: do k=kpres,nkds,-1
1140 weight=wbye(i,k,1)
1141 ck=ibye(i,k,1)
1143 if(ck>1) then
1144 fbye(i,k,1) = &
1145 ICOPY(cfield,i,j,ck) * weight + &
1146 ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1147 else
1148 if(emethod==EConst) then
1149 fbye(i,k,1)=evalue
1150 else if(emethod==ECopy) then
1151 fbye(i,k,1)=ICOPY(cfield,i,j,1)
1152 else ! assume 2: linear extrap
1153 fbye(i,k,1)=evalue*weight + &
1154 ICOPY(cfield,i,j,1)*(1.0-weight)
1155 endif
1156 endif
1157 enddo kloop4
1158 enddo
1159 endif
1161 end subroutine c2b_mass
1163 subroutine n2c_mass (&
1164 cfield,nfield,iinfo,winfo,ipos,jpos, &
1165 emethod, evalue, &
1166 cids, cide, cjds, cjde, ckds, ckde, &
1167 cims, cime, cjms, cjme, ckms, ckme, &
1168 cits, cite, cjts, cjte, ckts, ckte, &
1169 nids, nide, njds, njde, nkds, nkde, &
1170 nims, nime, njms, njme, nkms, nkme, &
1171 nits, nite, njts, njte, nkts, nkte)
1172 use module_interp_store, only: kpres
1173 implicit none
1174 integer, intent(in) :: &
1175 cids, cide, cjds, cjde, ckds, ckde, &
1176 cims, cime, cjms, cjme, ckms, ckme, &
1177 cits, cite, cjts, cjte, ckts, ckte, &
1178 nids, nide, njds, njde, nkds, nkde, &
1179 nims, nime, njms, njme, nkms, nkme, &
1180 nits, nite, njts, njte, nkts, nkte, &
1181 ipos,jpos, emethod
1182 real, intent(in) :: evalue
1183 real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: winfo
1184 integer, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: iinfo
1185 real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
1186 real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
1188 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1189 integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
1190 real :: weight
1192 nz=nkde-nkds+1
1194 jstart=MAX(jpos+1,cjts)
1195 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1197 kcopy: do k=ckde,kpres+1,-1
1198 jcopy: do j=jstart,jend
1199 a=mod(j,2)
1200 istart=MAX(ipos+a,cits)
1201 iend=MIN(ipos+(nide-nids)/nri-1,cite)
1202 nj = (j-jpos)*nrj + 1
1203 icopy: do i=istart,iend
1204 ni = (i-ipos)*nri + 2 - a
1205 cfield(i,j,k) = NGRAB(nfield,ni,nj,k)
1206 enddo icopy
1207 enddo jcopy
1208 enddo kcopy
1210 bigj: do j=jstart,jend
1211 a=mod(j,2)
1212 istart=MAX(ipos+a,cits)
1213 iend=MIN(ipos+(nide-nids)/nri-1,cite)
1214 nj = (j-jpos)*nrj + 1
1216 iloop: do i=istart,iend
1217 ni = (i-ipos)*nri + 2 - a
1218 kinterploop: do k=kpres,ckds,-1
1219 weight=winfo(i,j,k)
1220 nk=iinfo(i,j,k)
1222 if(nk>1) then
1223 cfield(i,j,k) = &
1224 NGRAB(nfield,ni,nj,nk) * weight + &
1225 ! pjj/cray - source line limit in Cray compiler
1226 NGRAB(nfield,ni,nj,nk-1) &
1227 * (1.0-weight)
1228 else
1229 ! Extrapolate
1230 if(emethod==EConst) then
1231 cfield(i,j,k)=evalue
1232 elseif(emethod==ECopy) then
1233 cfield(i,j,k)=NGRAB(nfield,ni,nj,1)
1234 else ! Assume 2: linear extrap
1235 cfield(i,j,k)=evalue*weight + &
1236 NGRAB(nfield,ni,nj,1)*(1.0-weight)
1237 endif
1238 endif
1239 end do kinterploop
1240 enddo iloop
1241 enddo bigj
1242 end subroutine n2c_mass
1244 subroutine n2c_massikj (&
1245 cfield,nfield,iinfo,winfo,ipos,jpos, &
1246 emethod, evalue, &
1247 cids, cide, cjds, cjde, ckds, ckde, &
1248 cims, cime, cjms, cjme, ckms, ckme, &
1249 cits, cite, cjts, cjte, ckts, ckte, &
1250 nids, nide, njds, njde, nkds, nkde, &
1251 nims, nime, njms, njme, nkms, nkme, &
1252 nits, nite, njts, njte, nkts, nkte)
1253 use module_interp_store, only: kpres
1254 implicit none
1255 integer, intent(in) :: &
1256 cids, cide, cjds, cjde, ckds, ckde, &
1257 cims, cime, cjms, cjme, ckms, ckme, &
1258 cits, cite, cjts, cjte, ckts, ckte, &
1259 nids, nide, njds, njde, nkds, nkde, &
1260 nims, nime, njms, njme, nkms, nkme, &
1261 nits, nite, njts, njte, nkts, nkte, &
1262 ipos,jpos, emethod
1263 real, intent(in) :: evalue
1264 real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: winfo
1265 integer, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: iinfo
1266 real, intent(out) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
1267 real, intent(in) :: nfield(nims:nime,nkms:nkme,njms:njme)
1269 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1270 integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
1271 real :: weight
1273 nz=nkde-nkds+1
1275 jstart=MAX(jpos+1,cjts)
1276 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1278 bigj: do j=jstart,jend
1279 a=mod(j,2)
1280 istart=MAX(ipos+a,cits)
1281 iend=MIN(ipos+(nide-nids)/nri-1,cite)
1282 nj = (j-jpos)*nrj + 1
1284 kcopyloop: do k=nkde,kpres+1,-1
1285 icopyloop: do i=istart,iend
1286 ni = (i-ipos)*nri + 2 - a
1287 cfield(i,k,j) = NGRABIKJ(nfield,ni,k,nj)
1288 enddo icopyloop
1289 enddo kcopyloop
1291 kinterploop: do k=kpres+1,nkds,-1
1292 iloop: do i=istart,iend
1293 ni = (i-ipos)*nri + 2 - a
1294 weight=winfo(i,j,k)
1295 nk=iinfo(i,j,k)
1297 if(nk>1) then
1298 cfield(i,k,j) = &
1299 NGRABIKJ(nfield,ni,nk,nj) * weight + &
1300 ! pjj/cray - source line limit in Cray compiler
1301 NGRABIKJ(nfield,ni,nk-1,nj) &
1302 * (1.0-weight)
1303 else
1304 ! Extrapolate
1305 if(emethod==EConst) then
1306 cfield(i,k,j)=evalue
1307 elseif(emethod==ECopy) then
1308 cfield(i,k,j)=NGRABIKJ(nfield,ni,1,nj)
1309 else ! Assume 2: linear extrap
1310 cfield(i,k,j)=evalue*weight + &
1311 NGRABIKJ(nfield,ni,1,nj)*(1.0-weight)
1312 endif
1313 endif
1314 enddo iloop
1315 end do kinterploop
1316 enddo bigj
1317 end subroutine n2c_massikj
1319 subroutine c2n_mass (II,JJ,W1,W2,W3,W4, &
1320 cfield,nfield,iinfo,winfo,imask, &
1321 emethod, evalue, &
1322 cims, cime, cjms, cjme, ckms, ckme, &
1323 nids, nide, njds, njde, nkds, nkde, &
1324 nims, nime, njms, njme, nkms, nkme, &
1325 nits, nite, njts, njte, nkts, nkte)
1326 use module_interp_store, only: kpres
1327 implicit none
1328 real, intent(in), dimension(nims:nime,njms:njme) :: &
1329 W1,W2,W3,W4
1330 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1332 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
1333 nids, nide, njds, njde, nkds, nkde, &
1334 nims, nime, njms, njme, nkms, nkme, &
1335 nits, nite, njts, njte, nkts, nkte, emethod
1336 real, intent(in) :: evalue
1337 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: winfo
1338 integer, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: iinfo
1339 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
1340 real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
1341 real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
1342 character*255 :: message
1343 integer :: j,i,a,nx,nz,ck,k
1344 real :: weight
1346 nx=min(nide-1,nite)-max(nids,nits)+1
1347 nz=nkde-nkds+1
1348 if(kpres<=nkds .or. kpres>=nkde) then
1349 call wrf_error_fatal('invalid kpres: outside domain bounds')
1350 end if
1351 bigj: do j=max(njds,njts),min(njde-1,njte)
1352 interploop: do i=max(nids,nits),min(nide-1,nite)
1353 if(imask(i,j)/=0) cycle interploop
1354 a=1-mod(JJ(i,j),2)
1355 kcopyloop: do k=nkde,kpres+1,-1
1356 nfield(i,j,k)=ICOPY(cfield,i,j,k)
1357 !weight=winfo(i,j,k)
1358 ck=iinfo(i,j,k)
1359 ! if(1.-weight>1e-5) then
1360 ! 2000 format(I0,", ",I0,", ",I0,": invalid weight at constant pressure level: w=",F0.5," i=",I0)
1361 ! write(message,2000) i,j,k, weight,ck
1362 ! call wrf_error_fatal(message)
1363 ! endif
1364 ! if(ck/=k) then
1365 ! 2001 format(I0,", ",I0,", ",I0,": invalid iinfo at constant pressure level: w=",F0.5," i=",I0)
1366 ! write(message,2000) i,j,k, weight,ck
1367 ! call wrf_error_fatal(message)
1368 ! endif
1369 enddo kcopyloop
1370 kinterploop: do k=kpres,nkds,-1
1371 weight=winfo(i,j,k)
1372 ck=iinfo(i,j,k)
1374 if(ck>1) then
1375 nfield(i,j,k) = &
1376 ICOPY(cfield,i,j,ck) * weight + &
1377 ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1378 else
1379 ! Extrapolate
1380 if(emethod==EConst) then
1381 nfield(i,j,k)=evalue
1382 elseif(emethod==ECopy) then
1383 nfield(i,j,k)=ICOPY(cfield,i,j,1)
1384 else ! Assume 2: linear extrap
1385 nfield(i,j,k)=evalue*weight + &
1386 ICOPY(cfield,i,j,1)*(1.0-weight)
1387 endif
1388 endif
1389 enddo kinterploop
1390 enddo interploop
1391 end do bigj
1393 end subroutine c2n_mass
1395 subroutine c2n_massikj (II,JJ,W1,W2,W3,W4, &
1396 cfield,nfield,iinfo,winfo,imask, &
1397 emethod, evalue, &
1398 cims, cime, cjms, cjme, ckms, ckme, &
1399 nids, nide, njds, njde, nkds, nkde, &
1400 nims, nime, njms, njme, nkms, nkme, &
1401 nits, nite, njts, njte, nkts, nkte)
1402 use module_interp_store, only: kpres
1403 implicit none
1404 real, intent(in), dimension(nims:nime,njms:njme) :: &
1405 W1,W2,W3,W4
1406 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1408 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
1409 nids, nide, njds, njde, nkds, nkde, &
1410 nims, nime, njms, njme, nkms, nkme, &
1411 nits, nite, njts, njte, nkts, nkte, emethod
1412 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: winfo
1413 integer, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: iinfo
1414 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
1415 real, intent(in) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
1416 real, intent(out) :: nfield(nims:nime,nkms:nkme,njms:njme)
1417 real, intent(In) :: evalue
1419 integer :: j,i,a,nx,nz,ck,k
1420 real :: weight
1422 nx=min(nide-1,nite)-max(nids,nits)+1
1423 nz=nkde-nkds+1
1425 bigj: do j=max(njds,njts),min(njde-1,njte)
1426 interploop: do i=max(nids,nits),min(nide-1,nite)
1427 if(imask(i,j)/=0) cycle interploop
1428 a=1-mod(JJ(i,j),2)
1429 kcopyloop: do k=nkde,kpres+1,-1
1430 nfield(i,k,j)=IKJCOPY(cfield,i,k,j)
1431 end do kcopyloop
1432 kinterploop: do k=kpres,nkds,-1
1433 weight=winfo(i,j,k)
1434 ck=iinfo(i,j,k)
1436 if(ck>1) then
1437 nfield(i,k,j) = &
1438 IKJCOPY(cfield,i,ck,j) * weight + &
1439 IKJCOPY(cfield,i,ck-1,j) * (1.0-weight)
1440 else
1441 ! Extrapolate
1442 if(emethod==EConst) then
1443 nfield(i,k,j)=evalue
1444 elseif(emethod==ECopy) then
1445 nfield(i,k,j)=IKJCOPY(cfield,i,1,j)
1446 else ! Assume 2: linear extrapolation
1447 nfield(i,k,j)=evalue * weight + &
1448 IKJCOPY(cfield,i,1,j) * (1.0-weight)
1449 endif
1450 endif
1451 enddo kinterploop
1452 enddo interploop
1453 end do bigj
1455 end subroutine c2n_massikj
1457 subroutine find_kpres(kpres, eta2, kds,kde, kms,kme)
1458 ! Find the level at which the pressure/sigma transition occurs.
1459 ! All levels above this are pure pressure levels, whereas kpres
1460 ! and below are sigma.
1461 integer, intent(out) :: kpres
1462 real, intent(in) :: eta2(kms:kme)
1463 integer, intent(in) :: kds,kde, kms,kme
1465 integer :: k
1467 k=kde-1
1468 do while(eta2(k) < 1e-5 .and. eta2(k-1) < 1e-5)
1469 k=k-1
1470 if(k<=kds) then
1471 call wrf_error_fatal('New NMM interpolation routines do not work in a constant pressure space.')
1472 endif
1473 enddo
1474 kpres=k
1475 end subroutine find_kpres
1477 ! ********************************************************************
1478 ! subs *_FULLDOM -- recalculates PD and PINT based on FIS if
1479 ! requested. Also, generates vertical interpolation arrays for use
1480 ! in later interpolations and interpolates T and Q while doing so.
1481 ! ********************************************************************
1483 subroutine n2c_fulldom ( &
1484 deta1,deta2, eta1,eta2, ptop,pdtop, &
1485 cfis,cpint,ct,cpd,cq, &
1486 cids, cide, cjds, cjde, ckds, ckde, &
1487 cims, cime, cjms, cjme, ckms, ckme, &
1488 cits, cite, cjts, cjte, ckts, ckte, &
1489 nfis,npint,nt,npd,nq, &
1490 ipos,jpos, &
1491 ! nri,nrj &
1492 out_iinfo,out_winfo, &
1493 nids, nide, njds, njde, nkds, nkde, &
1494 nims, nime, njms, njme, nkms, nkme, &
1495 nits, nite, njts, njte, nkts, nkte, kpres)
1497 implicit none
1499 integer, intent(in) :: ipos,jpos
1500 ! integer, intent(in) :: nri,nrj
1501 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1503 integer, intent(in):: &
1504 nids, nide, njds, njde, nkds, nkde, &
1505 nims, nime, njms, njme, nkms, nkme, &
1506 nits, nite, njts, njte, nkts, nkte, kpres
1507 integer, intent(in):: &
1508 cids, cide, cjds, cjde, ckds, ckde, &
1509 cims, cime, cjms, cjme, ckms, ckme, &
1510 cits, cite, cjts, cjte, ckts, ckte
1511 real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1512 real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1514 real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
1515 integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
1517 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
1518 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
1519 real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
1520 real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
1522 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1523 real, intent(in) :: ptop,pdtop
1525 real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
1526 real, dimension(nkde-1,nite-nits+1) :: inT0,inQ,icT,icQ, qinfo,winfo
1527 real, dimension(nkde,nite-nits+1) :: inPINT,icPINT
1528 integer, dimension(nkde-1,nite-nits+1) :: iinfo
1529 integer :: nx,nz,k,i,a,j, istart,iend,jstart,jend, ni,nj,jprint,itest,jtest
1530 character*255 :: message
1531 logical bad, warned
1533 warned=.false. ! allow one P>PSTD message
1535 nx=min(cide-2,cite)-max(cids+1,cits)+1
1536 nz=ckde-ckds+1
1538 jstart=MAX(jpos+1,cjts)
1539 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1541 bigj: do j=jstart,jend
1542 nj = (j-jpos)*nrj + 1
1544 a=mod(j,2)
1545 istart=MAX(ipos+a,cits)
1546 iend=MIN(ipos+(nide-nids)/nri-1,cite)
1547 nx=iend-istart+1
1549 ! STEP 1: Copy coarse and fine nest data into
1550 ! temporary arrays, reordering dimensions:
1551 qtloop: do k=ckts,ckte-1
1552 do i=istart,iend
1553 ni = (i-ipos)*nri + 2 - a
1555 UPCOPY(inT0,nT,i,j,k,ni,nj)
1556 UPCOPY(inQ,nQ,i,j,k,ni,nj)
1557 UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1558 enddo
1559 enddo qtloop
1561 k=nkte
1562 loop2d: do i=istart,iend
1563 ni = (i-ipos)*nri + 2 - a
1565 UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1566 UPCOPY(inPD,nPD,i,j,1,ni,nj)
1567 UPCOPY(inFIS,nFIS,i,j,1,ni,nj)
1569 icPD(1,i-istart+1)=cPD(i,j,1)
1570 ! icPD(1,i-istart+1)=use_this_pd(i,j,1)
1571 icFIS(1,i-istart+1)=cFIS(i,j,1)
1572 enddo loop2d
1574 ! Step 2: Interpolate coarse grid to fine grid in reordered
1575 ! arrays:
1576 call interp_T_PD_Q(nmm_method_linear, nmm_keep_pd, nx,nz, &
1577 deta1,deta2,eta1,eta2,ptop,pdtop, kpres, &
1578 inFIS,icFIS, inPINT,icPINT, inT0, icT, inPD,icPD, inQ,icQ, &
1579 iinfo, winfo, warned)
1581 ! Step 3: Copy back from reordered arrays to final nest arrays:
1583 qtloop2: do k=ckts,max(ckte-1,ckde-1)
1584 do i=istart,iend
1585 cT(i,j,k)=icT(k,i-istart+1)
1586 cQ(i,j,k)=icQ(k,i-istart+1)
1587 enddo
1588 enddo qtloop2
1590 izloop: do k=ckts,max(ckte-1,ckde-1)
1591 ixloop: do i=istart,iend
1592 out_iinfo(i,j,k)=iinfo(k,i-istart+1)
1593 out_winfo(i,j,k)=winfo(k,i-istart+1)
1594 enddo ixloop
1595 enddo izloop
1597 iendloop: do i=istart,iend
1598 out_iinfo(i,j,nkde)=nkde
1599 out_winfo(i,j,nkde)=1.0
1600 end do iendloop
1602 k=nkte+1
1603 loop2d2: do i=istart,iend
1604 cPD(i,j,1)=icPD(1,i-istart+1)
1605 enddo loop2d2
1606 end do bigj
1607 end subroutine n2c_fulldom
1609 subroutine n2c_fulldom_new ( &
1610 deta1,deta2, eta1,eta2, ptop,pdtop, &
1611 cfis,cpint,ct,cpd,cq, &
1612 cids, cide, cjds, cjde, ckds, ckde, &
1613 cims, cime, cjms, cjme, ckms, ckme, &
1614 cits, cite, cjts, cjte, ckts, ckte, &
1615 nfis,npint,nt,npd,nq, &
1616 ipos,jpos, &
1617 ! nri,nrj &
1618 out_iinfo,out_winfo, &
1619 nids, nide, njds, njde, nkds, nkde, &
1620 nims, nime, njms, njme, nkms, nkme, &
1621 nits, nite, njts, njte, nkts, nkte, kpres)
1623 implicit none
1625 integer, intent(in) :: ipos,jpos
1626 ! integer, intent(in) :: nri,nrj
1627 integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1629 integer, intent(in):: &
1630 nids, nide, njds, njde, nkds, nkde, &
1631 nims, nime, njms, njme, nkms, nkme, &
1632 nits, nite, njts, njte, nkts, nkte, kpres
1633 integer, intent(in):: &
1634 cids, cide, cjds, cjde, ckds, ckde, &
1635 cims, cime, cjms, cjme, ckms, ckme, &
1636 cits, cite, cjts, cjte, ckts, ckte
1637 real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1638 real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1640 real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
1641 integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
1643 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
1644 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
1645 real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
1646 real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
1648 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1649 real, intent(in) :: ptop,pdtop
1651 real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
1652 real, dimension(kpres+1,nite-nits+1) :: inT0,inQ,icT,icQ, qinfo,winfo
1653 real, dimension(kpres+2,nite-nits+1) :: inPINT,icPINT
1654 integer, dimension(kpres+1,nite-nits+1) :: iinfo
1655 integer :: nx,nz,k,i,a,j, istart,iend,jstart,jend, ni,nj,jprint,itest,jtest
1656 character*255 :: message
1657 logical bad, warned
1659 warned=.false. ! Allow one P>PSTD message
1661 nx=min(cide-2,cite)-max(cids+1,cits)+1
1662 nz=ckde-ckds+1
1664 jstart=MAX(jpos+1,cjts)
1665 jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1667 bigj: do j=jstart,jend
1668 nj = (j-jpos)*nrj + 1
1670 a=mod(j,2)
1671 istart=MAX(ipos+a,cits)
1672 iend=MIN(ipos+(nide-nids)/nri-1,cite)
1673 nx=iend-istart+1
1675 ! STEP 1: Copy coarse and fine nest data into
1676 ! temporary arrays, reordering dimensions:
1677 qtloop: do k=ckts,kpres+1
1678 do i=istart,iend
1679 ni = (i-ipos)*nri + 2 - a
1681 UPCOPY(inT0,nT,i,j,k,ni,nj)
1682 UPCOPY(inQ,nQ,i,j,k,ni,nj)
1683 UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1684 enddo
1685 enddo qtloop
1687 k=kpres+1
1688 loop2d: do i=istart,iend
1689 ni = (i-ipos)*nri + 2 - a
1691 UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1692 UPCOPY(inPD,nPD,i,j,1,ni,nj)
1693 UPCOPY(inFIS,nFIS,i,j,1,ni,nj)
1695 icPD(1,i-istart+1)=cPD(i,j,1)
1696 ! icPD(1,i-istart+1)=use_this_pd(i,j,1)
1697 icFIS(1,i-istart+1)=cFIS(i,j,1)
1698 enddo loop2d
1700 ! Step 2: Interpolate coarse grid to fine grid in reordered
1701 ! arrays:
1702 call interp_T_PD_Q_kpres(nmm_method_linear, nmm_keep_pd, nx,nz, &
1703 deta1,deta2,eta1,eta2,ptop,pdtop, kpres, kpres+2, &
1704 inFIS,icFIS, inPINT,icPINT, inT0, icT, inPD,icPD, inQ,icQ, &
1705 iinfo, winfo, warned)
1707 ! Step 3: Copy back from reordered arrays to final nest arrays:
1709 qtloop2: do k=ckts,kpres+1
1710 do i=istart,iend
1711 cT(i,j,k)=icT(k,i-istart+1)
1712 cQ(i,j,k)=icQ(k,i-istart+1)
1713 enddo
1714 enddo qtloop2
1716 izloop: do k=ckts,kpres+1
1717 ixloop: do i=istart,iend
1718 out_iinfo(i,j,k)=iinfo(k,i-istart+1)
1719 out_winfo(i,j,k)=winfo(k,i-istart+1)
1720 enddo ixloop
1721 enddo izloop
1723 k=nkte+1
1724 loop2d2: do i=istart,iend
1725 cPD(i,j,1)=icPD(1,i-istart+1)
1726 enddo loop2d2
1727 end do bigj
1729 kcopy: do k=kpres+1,ckde-1
1730 jcopy: do j=jstart,jend
1731 nj = (j-jpos)*nrj + 1
1733 a=mod(j,2)
1734 istart=MAX(ipos+a,cits)
1735 iend=MIN(ipos+(nide-nids)/nri-1,cite)
1736 icopy: do i=istart,iend
1737 ni = (i-ipos)*nri + 2 - a
1738 cT(i,j,k)=NGRAB(nT,ni,nj,k)
1739 cQ(i,j,k)=NGRAB(nQ,ni,nj,k)
1740 out_iinfo(i,j,k)=k
1741 out_winfo(i,j,k)=1.0
1742 enddo icopy
1743 if(k==ckde-1) then
1744 icopy2: do i=istart,iend
1745 out_iinfo(i,j,nkde)=nkde
1746 out_winfo(i,j,nkde)=1.0
1747 enddo icopy2
1748 endif
1749 end do jcopy
1750 end do kcopy
1751 end subroutine n2c_fulldom_new
1753 subroutine c2b_fulldom (II,JJ,W1,W2,W3,W4,&
1754 deta1,deta2, eta1,eta2, ptop,pdtop, &
1755 cfis,cpint,ct,cpd,cq, nfis, &
1756 cims, cime, cjms, cjme, ckms, ckme, &
1757 nids, nide, njds, njde, nkds, nkde, &
1758 nims, nime, njms, njme, nkms, nkme, &
1759 nits, nite, njts, njte, nkts, nkte, &
1760 kpres, &
1761 ibxs, ibxe, ibys, ibye, &
1762 wbxs, wbxe, wbys, wbye, &
1763 pdbxs, pdbxe, pdbys, pdbye, &
1764 tbxs, tbxe, tbys, tbye, &
1765 qbxs, qbxe, qbys, qbye)
1766 implicit none
1767 integer, intent(in):: &
1768 cims, cime, cjms, cjme, ckms, ckme, &
1769 nids, nide, njds, njde, nkds, nkde, &
1770 nims, nime, njms, njme, nkms, nkme, &
1771 nits, nite, njts, njte, nkts, nkte
1773 integer, parameter :: bdyw=1
1775 ! Domain info:
1776 real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
1777 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1778 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1779 real, intent(in) :: ptop,pdtop
1780 integer, intent(in) :: kpres
1782 ! Parent fields:
1783 real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1784 real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1786 ! Nest terrain info:
1787 real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
1789 ! T, Q, PINT, PD boundary info:
1790 real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
1791 real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
1792 real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
1793 real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
1795 ! Weights and indices:
1796 real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
1797 real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
1798 integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
1799 integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
1801 integer :: i,j,k,a,b,nx,nz,used,j1,j2,used1
1803 real, dimension(1,2*(nite-nits+5)+2*(njte-njts+5)) :: inFIS,inPD,icFIS,icPD
1804 real, dimension(nkde-1,2*(nite-nits+5)+2*(njte-njts+5)) :: inT,inQ,icT,icQ,winfo
1805 integer, dimension(nkde-1,2*(nite-nits+5)+2*(njte-njts+5)) :: iinfo
1806 real, dimension(nkde,2*(nite-nits+5)+2*(njte-njts+5)) :: inPINT,icPINT
1807 logical :: warned
1809 warned=.false. ! Allow one P>PSTD message
1811 nx=min(nide-1,nite)-max(nids,nits)+1
1812 nz=nkde-nkds+1
1814 j1=max(njts-1,njds)
1815 if(mod(j1,2)/=1) j1=j1+1
1817 j2=min(njte+1,njde-1)
1818 if(mod(j2,2)/=1) j2=j2-1
1820 used=0
1821 bdyloop: do b=1,4
1822 if_xbdy: if(b==1 .or. b==2) then
1823 if(b==1) then
1824 if(nits/=1) cycle bdyloop
1825 i=1
1826 endif
1827 if(b==2) then
1828 if(nite<nide-1) cycle bdyloop
1829 i=nide-1
1830 endif
1831 do j=j1,j2,2
1832 a=1-mod(JJ(i,j),2)
1833 used=used+1
1834 do k=nkts,nkte-1
1835 uFCOPY(icT,cT,i,j,k)
1836 uFCOPY(icQ,cQ,i,j,k)
1837 uFCOPY(icPINT,cPINT,i,j,k)
1838 enddo
1840 k=nkte
1841 a=1-mod(JJ(i,j),2)
1843 uFCOPY(icPINT,cPINT,i,j,k)
1844 uFCOPY(icPD,cPD,i,j,1)
1845 uFCOPY(icFIS,cFIS,i,j,1)
1847 inFIS(1,used)=nFIS(i,j,1)
1849 enddo
1850 endif if_xbdy
1851 if_ybdy: if(b==3 .or. b==4) then
1852 if(b==3) then
1853 if(njts/=1) cycle bdyloop
1854 j=1
1855 endif
1856 if(b==4) then
1857 if(njte<njde-1) cycle bdyloop
1858 j=njde-1
1859 endif
1860 do i=max(nits-1,nids),min(nite+1,nide-1)
1861 used=used+1
1862 do k=nkts,nkte-1
1863 a=1-mod(JJ(i,j),2)
1864 uFCOPY(icT,cT,i,j,k)
1865 uFCOPY(icQ,cQ,i,j,k)
1866 uFCOPY(icPINT,cPINT,i,j,k)
1867 enddo
1869 k=nkte
1870 a=1-mod(JJ(i,j),2)
1872 uFCOPY(icPINT,cPINT,i,j,k)
1873 uFCOPY(icPD,cPD,i,j,1)
1874 uFCOPY(icFIS,cFIS,i,j,1)
1876 inFIS(1,used)=nFIS(i,j,1)
1877 enddo
1878 endif if_ybdy
1879 enddo bdyloop
1881 if(used==0) then
1882 ! No boundary points on this tile so we are done
1883 return
1884 endif
1886 call interp_T_PD_Q_kpres(nmm_method_linear, nmm_interp_pd, used,nz, &
1887 deta1,deta2,eta1,eta2,ptop,pdtop, kpres, nz, &
1888 icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
1889 iinfo, winfo, warned)
1891 used1=used
1892 used=0
1894 if_bxs: if(nits==1) then
1895 i=1
1896 do j=j1,j2,2
1897 used=used+1
1898 do k=nkts,nkte-1
1899 tbxs(j,k,1)=inT(k,used)
1900 qbxs(j,k,1)=inQ(k,used)
1901 ibxs(j,k,1)=iinfo(k,used)
1902 wbxs(j,k,1)=winfo(k,used)
1903 enddo
1904 ibxs(j,nkde,1)=nkde
1905 wbxs(j,nkde,1)=1.0
1906 pdbxs(j,1,1)=inPD(1,used)
1907 enddo
1908 endif if_bxs
1910 if_bxe: if(nite>=nide-1) then
1911 i=nide-1
1912 do j=j1,j2,2
1913 used=used+1
1914 do k=nkts,nkte-1
1915 tbxe(j,k,1)=inT(k,used)
1916 qbxe(j,k,1)=inQ(k,used)
1917 ibxe(j,k,1)=iinfo(k,used)
1918 wbxe(j,k,1)=winfo(k,used)
1919 enddo
1920 ibxe(j,nkde,1)=nkde
1921 wbxe(j,nkde,1)=1.0
1922 pdbxe(j,1,1)=inPD(1,used)
1923 enddo
1924 endif if_bxe
1926 if_bys: if(njts==1) then
1927 j=1
1928 do i=max(nits-1,nids),min(nite+1,nide-1)
1929 used=used+1
1930 do k=nkts,nkte-1
1931 tbys(i,k,1)=inT(k,used)
1932 qbys(i,k,1)=inQ(k,used)
1933 ibys(i,k,1)=iinfo(k,used)
1934 wbys(i,k,1)=winfo(k,used)
1935 enddo
1936 ibys(i,nkde,1)=nkde
1937 wbys(i,nkde,1)=1.0
1938 pdbys(i,1,1)=inPD(1,used)
1939 enddo
1940 endif if_bys
1942 if_bye: if(njte>=njde-1) then
1943 j=njde-1
1944 do i=max(nits-1,nids),min(nite+1,nide-1)
1945 used=used+1
1946 do k=nkts,nkte-1
1947 tbye(i,k,1)=inT(k,used)
1948 qbye(i,k,1)=inQ(k,used)
1949 ibye(i,k,1)=iinfo(k,used)
1950 wbye(i,k,1)=winfo(k,used)
1951 enddo
1952 ibye(i,nkde,1)=nkde
1953 wbye(i,nkde,1)=1.0
1954 pdbye(i,1,1)=inPD(1,used)
1955 enddo
1956 endif if_bye
1958 if(used/=used1) then
1959 call wrf_error_fatal('Number of input and output points does not match.')
1960 endif
1962 end subroutine c2b_fulldom
1964 ! ####################################################################
1966 subroutine c2b_fulldom_new (II,JJ,W1,W2,W3,W4,&
1967 deta1,deta2, eta1,eta2, ptop,pdtop, &
1968 cfis,cpint,ct,cpd,cq, nfis, &
1969 cims, cime, cjms, cjme, ckms, ckme, &
1970 nids, nide, njds, njde, nkds, nkde, &
1971 nims, nime, njms, njme, nkms, nkme, &
1972 nits, nite, njts, njte, nkts, nkte, &
1973 kpres, &
1974 ibxs, ibxe, ibys, ibye, &
1975 wbxs, wbxe, wbys, wbye, &
1976 pdbxs, pdbxe, pdbys, pdbye, &
1977 tbxs, tbxe, tbys, tbye, &
1978 qbxs, qbxe, qbys, qbye)
1979 implicit none
1980 integer, intent(in):: &
1981 cims, cime, cjms, cjme, ckms, ckme, &
1982 nids, nide, njds, njde, nkds, nkde, &
1983 nims, nime, njms, njme, nkms, nkme, &
1984 nits, nite, njts, njte, nkts, nkte
1986 integer, parameter :: bdyw=1
1988 ! Domain info:
1989 real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
1990 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1991 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1992 real, intent(in) :: ptop,pdtop
1993 integer, intent(in) :: kpres
1995 ! Parent fields:
1996 real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1997 real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1999 ! Nest terrain info:
2000 real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
2002 ! T, Q, PINT, PD boundary info:
2003 real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
2004 real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
2005 real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
2006 real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
2008 ! Weights and indices:
2009 real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
2010 real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
2011 integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
2012 integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
2014 integer :: i,j,k,a,b,nx,nz,used,j1,j2,used1
2016 real, dimension(1,2*(nite-nits+5)+2*(njte-njts+5)) :: inFIS,inPD,icFIS,icPD
2017 real, dimension(kpres+1,2*(nite-nits+5)+2*(njte-njts+5)) :: inT,inQ,icT,icQ,winfo
2018 integer, dimension(kpres+1,2*(nite-nits+5)+2*(njte-njts+5)) :: iinfo
2019 real, dimension(kpres+2,2*(nite-nits+5)+2*(njte-njts+5)) :: inPINT,icPINT
2020 logical :: warned
2022 warned=.false. ! Allow one P>PSTD message
2024 nx=min(nide-1,nite)-max(nids,nits)+1
2025 nz=nkde-nkds+1
2027 j1=max(njts-1,njds)
2028 if(mod(j1,2)/=1) j1=j1+1
2030 j2=min(njte+1,njde-1)
2031 if(mod(j2,2)/=1) j2=j2-1
2033 used=0
2034 bdyloop: do b=1,4
2035 if_xbdy: if(b==1 .or. b==2) then
2036 if(b==1) then
2037 if(nits/=1) cycle bdyloop
2038 i=1
2039 endif
2040 if(b==2) then
2041 if(nite<nide-1) cycle bdyloop
2042 i=nide-1
2043 endif
2044 do j=j1,j2,2
2045 used=used+1
2046 do k=nkts,kpres+1
2047 a=1-mod(JJ(i,j),2)
2048 uFCOPY(icT,cT,i,j,k)
2049 uFCOPY(icQ,cQ,i,j,k)
2050 uFCOPY(icPINT,cPINT,i,j,k)
2051 enddo
2053 k=kpres+2
2054 a=1-mod(JJ(i,j),2)
2056 uFCOPY(icPINT,cPINT,i,j,k)
2057 uFCOPY(icPD,cPD,i,j,1)
2058 uFCOPY(icFIS,cFIS,i,j,1)
2060 inFIS(1,used)=nFIS(i,j,1)
2062 enddo
2063 endif if_xbdy
2064 if_ybdy: if(b==3 .or. b==4) then
2065 if(b==3) then
2066 if(njts/=1) cycle bdyloop
2067 j=1
2068 endif
2069 if(b==4) then
2070 if(njte<njde-1) cycle bdyloop
2071 j=njde-1
2072 endif
2073 do i=max(nits-1,nids),min(nite+1,nide-1)
2074 used=used+1
2075 do k=nkts,kpres+1
2076 a=1-mod(JJ(i,j),2)
2077 uFCOPY(icT,cT,i,j,k)
2078 uFCOPY(icQ,cQ,i,j,k)
2079 uFCOPY(icPINT,cPINT,i,j,k)
2080 enddo
2082 k=kpres+2
2083 a=1-mod(JJ(i,j),2)
2085 uFCOPY(icPINT,cPINT,i,j,k)
2086 uFCOPY(icPD,cPD,i,j,1)
2087 uFCOPY(icFIS,cFIS,i,j,1)
2089 inFIS(1,used)=nFIS(i,j,1)
2090 enddo
2091 endif if_ybdy
2092 enddo bdyloop
2094 if(used==0) then
2095 ! No boundary points on this tile so we are done
2096 return
2097 endif
2099 call interp_T_PD_Q_kpres(nmm_method_linear, nmm_interp_pd, used,nz, &
2100 deta1,deta2,eta1,eta2,ptop,pdtop, kpres, kpres+2, &
2101 icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
2102 iinfo, winfo, warned)
2104 used1=used
2105 used=0
2107 if_bxs: if(nits==1) then
2108 i=1
2109 do j=j1,j2,2
2110 used=used+1
2111 do k=nkts,kpres
2112 tbxs(j,k,1)=inT(k,used)
2113 qbxs(j,k,1)=inQ(k,used)
2114 ibxs(j,k,1)=iinfo(k,used)
2115 wbxs(j,k,1)=winfo(k,used)
2116 enddo
2117 ibxs(j,nkde,1)=nkde
2118 wbxs(j,nkde,1)=1.0
2119 pdbxs(j,1,1)=inPD(1,used)
2120 enddo
2121 do k=kpres+1,nkde-1
2122 do j=j1,j2,2
2123 a=1-mod(JJ(i,j),2)
2124 tbxs(j,k,1)=ICOPY(cT,i,j,k)
2125 qbxs(j,k,1)=ICOPY(cQ,i,j,k)
2126 ibxs(j,k,1)=k
2127 wbxs(j,k,1)=1.0
2128 enddo
2129 enddo
2130 endif if_bxs
2132 if_bxe: if(nite>=nide-1) then
2133 i=nide-1
2134 do j=j1,j2,2
2135 used=used+1
2136 do k=nkts,kpres
2137 tbxe(j,k,1)=inT(k,used)
2138 qbxe(j,k,1)=inQ(k,used)
2139 ibxe(j,k,1)=iinfo(k,used)
2140 wbxe(j,k,1)=winfo(k,used)
2141 enddo
2142 ibxe(j,nkde,1)=nkde
2143 wbxe(j,nkde,1)=1.0
2144 pdbxe(j,1,1)=inPD(1,used)
2145 enddo
2146 do k=kpres+1,nkde-1
2147 do j=j1,j2,2
2148 a=1-mod(JJ(i,j),2)
2149 tbxe(j,k,1)=ICOPY(cT,i,j,k)
2150 qbxe(j,k,1)=ICOPY(cQ,i,j,k)
2151 ibxe(j,k,1)=k
2152 wbxe(j,k,1)=1.0
2153 enddo
2154 enddo
2155 endif if_bxe
2157 if_bys: if(njts==1) then
2158 j=1
2159 do i=max(nits-1,nids),min(nite+1,nide-1)
2160 used=used+1
2161 do k=nkts,kpres
2162 tbys(i,k,1)=inT(k,used)
2163 qbys(i,k,1)=inQ(k,used)
2164 ibys(i,k,1)=iinfo(k,used)
2165 wbys(i,k,1)=winfo(k,used)
2166 enddo
2167 ibys(i,nkde,1)=nkde
2168 wbys(i,nkde,1)=1.0
2169 pdbys(i,1,1)=inPD(1,used)
2170 enddo
2171 do k=kpres+1,nkde-1
2172 do i=max(nits-1,nids),min(nite+1,nide-1)
2173 a=1-mod(JJ(i,j),2)
2174 tbys(i,k,1)=ICOPY(cT,i,j,k)
2175 qbys(i,k,1)=ICOPY(cQ,i,j,k)
2176 ibys(i,k,1)=k
2177 wbys(i,k,1)=1.0
2178 enddo
2179 enddo
2180 endif if_bys
2182 if_bye: if(njte>=njde-1) then
2183 j=njde-1
2184 do i=max(nits-1,nids),min(nite+1,nide-1)
2185 used=used+1
2186 do k=nkts,kpres
2187 tbye(i,k,1)=inT(k,used)
2188 qbye(i,k,1)=inQ(k,used)
2189 ibye(i,k,1)=iinfo(k,used)
2190 wbye(i,k,1)=winfo(k,used)
2191 enddo
2192 ibye(i,nkde,1)=nkde
2193 wbye(i,nkde,1)=1.0
2194 pdbye(i,1,1)=inPD(1,used)
2195 enddo
2196 do k=kpres+1,nkde-1
2197 do i=max(nits-1,nids),min(nite+1,nide-1)
2198 a=1-mod(JJ(i,j),2)
2199 tbye(i,k,1)=ICOPY(cT,i,j,k)
2200 qbye(i,k,1)=ICOPY(cQ,i,j,k)
2201 ibye(i,k,1)=k
2202 wbye(i,k,1)=1.0
2203 enddo
2204 enddo
2205 endif if_bye
2207 if(used/=used1) then
2208 call wrf_error_fatal('Number of input and output points does not match.')
2209 endif
2211 end subroutine c2b_fulldom_new
2213 ! ####################################################################
2215 subroutine c2n_fulldom (II,JJ,W1,W2,W3,W4,&
2216 deta1,deta2, eta1,eta2, ptop,pdtop, &
2217 cfis,cpint,ct,cpd,cq, &
2218 cims, cime, cjms, cjme, ckms, ckme, &
2219 nfis,npint,nt,npd,nq, &
2220 out_iinfo,out_winfo,imask,&
2221 nids, nide, njds, njde, nkds, nkde, &
2222 nims, nime, njms, njme, nkms, nkme, &
2223 nits, nite, njts, njte, nkts, nkte, kpres)
2224 implicit none
2225 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
2226 nids, nide, njds, njde, nkds, nkde, &
2227 nims, nime, njms, njme, nkms, nkme, &
2228 nits, nite, njts, njte, nkts, nkte, kpres
2229 real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
2230 real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
2232 real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_winfo
2233 integer, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_iinfo
2235 real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
2236 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
2237 real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
2238 real, intent(inout), dimension(nims:nime,njms:njme,1) :: nPD
2239 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
2241 real, intent(in), dimension(nims:nime,njms:njme) :: &
2242 W1,W2,W3,W4
2243 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
2245 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
2246 real, intent(in) :: ptop,pdtop
2248 real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
2249 real, dimension(nkde-1,nite-nits+1) :: inT,inQ,icT,icQ,winfo
2250 integer, dimension(nkde-1,nite-nits+1) :: iinfo
2251 real, dimension(nkde,nite-nits+1) :: inPINT,icPINT
2253 real :: pdcheck
2254 integer :: i,j,k,a, nx,nz,used
2255 logical :: badbad, warned
2257 warned=.false. ! Allow one P>PSTD message
2259 nx=min(nide-1,nite)-max(nids,nits)+1
2260 nz=nkde-nkds+1
2262 bigj: do j=max(njds,njts),min(njde-1,njte)
2263 ! STEP 1: Copy coarse and fine nest data into
2264 ! temporary arrays, reordering dimensions:
2265 used=0
2266 iloop1: do i=max(nids,nits),min(nide-1,nite)
2267 if(imask(i,j)/=0) cycle iloop1
2268 used=used+1
2270 qtloop: do k=nkts,nkte-1
2271 a=1-mod(JJ(i,j),2)
2272 uFCOPY(icT,cT,i,j,k)
2273 uFCOPY(icQ,cQ,i,j,k)
2274 uFCOPY(icPINT,cPINT,i,j,k)
2275 enddo qtloop
2277 k=nkte
2278 a=1-mod(JJ(i,j),2)
2280 uFCOPY(icPINT,cPINT,i,j,k)
2281 uFCOPY(icPD,cPD,i,j,1)
2282 uFCOPY(icFIS,cFIS,i,j,1)
2284 ! nPD(i,j,1)=use_this_pd(i,j,1) ! FIXME: remove this line
2285 ! if(nPD(i,j,1)> 1.3e5 .or. nPD(i,j,1)<2e4) then
2286 ! 3131 format('Invalid input nest PD in C2N interpolation: PD(',I0,',',I0,') = ',F0.7)
2287 ! write(0,3131) i,j,nPD(i,j,1)
2288 ! call wrf_error_fatal('Invalid input nest PD')
2289 ! endif
2290 inPD(1,used)=nPD(i,j,1)
2291 inFIS(1,used)=nFIS(i,j,1)
2292 enddo iloop1
2294 ! Step 2: Interpolate coarse grid to fine grid in reordered
2295 ! arrays:
2297 if(used==0) then
2298 ! No points in this row require interpolation, so we're done
2299 ! with this row:
2300 cycle bigj
2301 else
2302 call interp_T_PD_Q(nmm_method_linear, nmm_interp_pd, used,nz, &
2303 deta1,deta2,eta1,eta2,ptop,pdtop,kpres, &
2304 icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
2305 iinfo, winfo, warned)
2306 endif
2308 ! Step 3: Copy back from reordered arrays to final nest arrays:
2309 used=0
2310 iloop2: do i=max(nids,nits),min(nide-1,nite)
2311 if(imask(i,j)/=0) cycle iloop2
2312 used=used+1
2314 qtloop2: do k=nkts,nkte-1
2315 nT(i,j,k)=inT(k,used)
2316 nQ(i,j,k)=inQ(k,used)
2317 nQ(i,j,k)=inQ(k,used)
2318 enddo qtloop2
2320 izloop: do k=nkts,nkte-1
2321 out_iinfo(i,j,k)=iinfo(k,used)
2322 out_winfo(i,j,k)=winfo(k,used)
2323 enddo izloop
2324 out_iinfo(i,j,nkde)=nkde
2325 out_winfo(i,j,nkde)=1.0
2327 k=nkte+1
2328 a=1-mod(JJ(i,j),2)
2329 nPD(i,j,1)=inPD(1,used)
2330 pdcheck=npd(i,j,1)+ptop+pdtop
2331 ! if(pdcheck<50000.0 .or. pdcheck>105000.0) then
2332 ! 3131 format('Invalid output nest PD in C2N interpolation: PD(',I0,',',I0,') = ',F0.7,' which is ',F0.7,' Pa')
2333 ! write(0,3131) i,j,nPD(i,j,1),pdcheck
2334 ! call wrf_error_fatal('Invalid output nest PD')
2335 ! endif
2336 enddo iloop2
2337 enddo bigj
2339 end subroutine c2n_fulldom
2341 ! ********************************************************************
2342 ! INTERP_T_PD_Q -- this is the actual vertical interpolation
2343 ! subroutine, the implementation of the *_fulldom subroutines. It
2344 ! takes arrays of atmospheric columns, with no knowledge of the
2345 ! horizontal layout of those columns. This routine will recalculate
2346 ! PD if requested, then generate interpolation information and
2347 ! interpolate T and Q, within each atmospheric column.
2348 !
2349 ! Other variables are handled by the various other c2n, n2c and c2b
2350 ! routines in this module. PD, T and Q are handled specially since
2351 ! they are critical to the stability of the hydrostatic solver, and
2352 ! hence have more expensive extrapolation and mass rebalancing
2353 ! methods.
2354 !
2355 ! Since this routine has no knowledge of the horizontal layout of
2356 ! the atmospheric columns it is given, the columns do not have to
2357 ! lie on an I- or J-row. The boundary interpolation routines take
2358 ! advantage of that and provide all boundary information to this
2359 ! routine in a single call. Of course, the caller must handle all
2360 ! horizontal interpolation.
2361 ! ********************************************************************
2363 subroutine interp_T_PD_Q(method, pd_interp, nx, nz, &
2364 deta1,deta2, eta1,eta2, ptop,pdtop, kpres, &
2365 fisA,fisB, pintA,pintB, tA,tB, pdA,pdB, qA,qB, &
2366 iinfo, winfo, warned)
2367 implicit none
2369 integer, intent(in) :: pd_interp,method, kpres
2370 ! real, intent(in) :: dtA,dtB
2372 ! Coordinate system definitions must be the same for all domains:
2373 real, intent(in) :: deta1(nz),deta2(nz),eta1(nz),eta2(nz),ptop,pdtop
2374 integer, intent(in) :: nx,nz
2376 ! Surface height and mass field information for source (A):
2377 real, intent(in) :: fisA(nx),tA(nz-1,nx),pdA(nx),qA(nz-1,nx),pintA(nz,nx)
2379 ! Surface height and mass field information for target (B):
2380 real, intent(inout) :: fisB(nx),tB(nz-1,nx),pdB(nx),qB(nz-1,nx),pintB(nz,nx)
2382 ! Interpolation or extrapolation information for use in other
2383 ! calls later on:
2384 real, intent(out) :: winfo(nz-1,nx)
2385 integer, intent(out) :: iinfo(nz-1,nx)
2386 logical, intent(inout) :: warned
2388 ! ==================== Local variables ====================
2390 character*255 :: message
2391 integer :: ix,iz,izA,izB,xpr
2392 real :: zA,zB,znext,apelp,rtopp,dz,weight,A,B,pstd1,z,pstd2,pstd12
2393 real :: tsfc(nx), slp(nx), zmslp(nx), z0mid(nx), zbelow(nx), &
2394 tbelow(nx), RHbelow(nx)
2395 real :: pb,pb1,pb2,pa,pa1,pa2,pnext,pa3, wnum,wdenom, QC, P0
2397 ! Constants from UPP params.f for RH calculation (all mks):
2398 real, parameter :: PQ0=379.90516
2399 real, parameter :: A2=17.2693882
2400 real, parameter :: A3=273.16
2401 real, parameter :: A4=35.86
2402 real, parameter :: RHmin=1.0E-6 ! minimal RH bound
2404 if(method/=nmm_method_linear) then
2405 call wrf_error_fatal('only linear interpolation is supported')
2406 endif
2408 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2409 !! Step 1: calculate near-surface values !!!!!!!!!!!!!!!!!!!!!
2410 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2412 pstd1=p_ref ! pstd(1) from base_state_parent
2413 ! pstd(2) from base_state_parent:
2414 pstd2=eta1(2)*pdtop + eta2(2)*(p_ref-pdtop-ptop) + ptop
2415 pstd12=exp((alog(pstd1)+alog(pstd2))*0.5)
2417 do ix=1,nx
2418 ! These calculations are from base_state_parent:
2419 APELP = (pintA(2,ix)+pintA(1,ix))
2420 RTOPP = TRG*tA(1,ix)*(1.0+qA(1,ix)*P608)/APELP
2421 DZ = RTOPP*(DETA1(1)*PDTOP+DETA2(1)*pdA(ix))
2423 Z0MID(ix) = fisA(ix)/g + dz/2
2425 zA=fisA(ix)/g
2427 TSFC(ix) = TA(1,ix)*(1.+D608*QA(1,ix)) &
2428 + LAPSR*(zA + zA+dz)*0.5
2430 A = LAPSR*zA/TSFC(ix)
2431 SLP(ix) = PINTA(1,ix)*(1-A)**COEF2 ! sea level pressure
2432 B = (pstd1/SLP(ix))**COEF3
2433 ZMSLP(ix)= TSFC(ix)*LAPSI*(1.0 - B) ! Height at 1030. mb level
2435 TBELOW(ix) = TA(1,ix) + LAPSR*(Z0MID(ix)-ZMSLP(ix))
2437 ! Calculate lowest level RH. This calculation is from UPP
2438 ! CALRH.f:
2439 P0=pdA(ix)+ptop+pdtop ! Use hydrostatic lowest model level P
2440 QC=PQ0/P0*EXP(A2*(TA(1,ix)-A3)/(TA(1,ix)-A4))
2441 RHbelow(ix)=max(RHmin,min(1.0,QA(1,ix)/QC))
2442 enddo
2446 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2447 !! Step 2: figure out the new surface pressures !!!!!!!!!!!!!!
2448 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2452 if_pd_interp: if(pd_interp==nmm_keep_pd) then
2453 ! PD interpolation is turned off, so we use the original PD.
2454 elseif(pd_interp==nmm_interp_pd) then
2455 ! PD interpolation is requested. As with the old base_state_parent,
2456 ! determine PD by interpolating or extrapolating the non-hydrostatic
2457 ! pressure field (PINT) in source grid to the surface height of the
2458 ! target grid.
2459 xloop: do ix=1,nx
2460 if(pintA(1,ix)>p_ref) then
2461 ! Cannot extrapolate PD below ground because pressure is
2462 ! unrealistically high. Follow base_state_parent method:
2463 ! when this happens, assign input pressure to output
2464 ! pressure.
2465 if(.not.warned) then
2466 call wrf_message2('WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD')
2467 write(message,*) 'PINT(1),PD(1),PSTD(1)',pintA(1,ix),pdA(ix),p_ref
2468 call wrf_message2(message)
2469 warned=.true.
2470 endif
2471 pdB(ix)=pdA(ix)
2472 cycle xloop
2473 endif
2475 zA=fisA(ix)/g
2476 zB=fisB(ix)/g
2478 if(zB<zA) then
2479 ! Need to extrapolate PD below ground.
2480 ! Follow base_state_parent method:
2481 ! For pressures between pint(1) and 1030 mbar, interpolate
2482 iz=1
2483 weight=abs((zB-zmslp(ix))/(zA-zmslp(ix)))
2484 pdB(ix) = exp(weight*log(pdA(ix)+pdtop+ptop) + &
2485 (1.0-weight)*log(P_REF)) &
2486 - pdtop - ptop
2487 ! write(0,*) 'WARNING: extrapolating below ground at ix=',ix
2488 !1302 format('Extrap at ix=',I0,' zA=',F0.3,' zB=',F0.3,' zmslp=',F0.3,' weight=',F0.3,' pdB+pdtop+ptop=',F0.3,' pintA=',F0.3,' pdA(ix)+pdtop+ptop=',F0.3,' P_REF=',F0.3)
2489 ! write(0,1302) ix,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2490 cycle xloop
2491 endif
2493 ! PD in grid B lies above ground in grid A. Find out where by
2494 ! walking up from the ground until we find the two levels the
2495 ! pressure lies between.
2496 znext=-9e9
2497 z=zA
2498 do iz=1,nz-1
2499 ! Calculations from base_state_parent:
2500 APELP = (pintA(iz+1,ix)+pintA(iz,ix))
2501 RTOPP = TRG*tA(iz,ix)*(1.0+qA(iz,ix)*P608)/APELP
2502 DZ = RTOPP*(DETA1(iz)*PDTOP+DETA2(iz)*pdA(ix))
2503 znext=z+dz
2505 if(znext>zB) then
2506 ! Interpolate between these two levels
2507 weight=(zB-z)/dz
2508 pdB(ix) = exp(weight*log(pintA(iz+1,ix)) + &
2509 (1.0-weight)*log(pintA(iz,ix))) &
2510 - pdtop - ptop
2511 ! if(iz>1) then
2512 ! write(0,*) 'WARNING: interpolate above ground at ix=',ix
2513 !1303 format('Interp at ix=',I0,' with iz=',I0,' zA=',F0.3,' zB=',F0.3,' zmslp=',F0.3,' weight=',F0.3,' pdB+pdtop+ptop=',F0.3,' pintA=',F0.3,' pdA(ix)+pdtop+ptop=',F0.3,' P_REF=',F0.3)
2514 ! write(0,1303) ix,iz,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2515 ! endif
2516 cycle xloop
2517 else
2518 z=znext
2519 endif
2521 enddo
2523 201 format('interp_T_PD_Q: Target domain surface height ',F0.4,'m is higher than the model top ',F0.4,'m of the source domain. ABORTING')
2524 write(message,201) zB,znext
2525 call wrf_error_fatal(message)
2526 enddo xloop
2527 else
2528 202 format('Invalid value ',I0,' for pd_interp in interp_T_PD_Q')
2529 write(message,202) pd_interp
2530 call wrf_error_fatal(message)
2531 endif if_pd_interp
2535 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2536 !! Step 3: interpolate T and Q in pressure space !!!!!!!!!!!!!
2537 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2541 outer: do ix=1,nx
2542 ! For constant pressure levels, do a straight level-by-level copy:
2543 iz=nz-1
2544 copyloop: do while(iz>kpres)
2545 tB(iz,ix)=tA(iz,ix)
2546 qB(iz,ix)=qA(iz,ix)
2548 ! Save interpolation information
2549 iinfo(iz,ix)=iz
2550 winfo(iz,ix)=1.0
2552 iz=iz-1
2553 enddo copyloop
2555 ! For sigma levels where target domain lies within source,
2556 ! interpolate from top down, stopping when we go below ground
2557 ! in the source domain.
2558 izA=iz
2559 izB=iz
2561 ! FIXME: REMOVE THIS CHECK
2562 if(iz>nz) then
2563 ! Make sure the entire vertical extent isn't sigma levels:
2564 call wrf_error_fatal('ERROR: WRF-NMM does not support pure sigma levels (only sigma-pressure hybrid)')
2565 endif
2567 pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2568 pB1=log(eta1(izB+1)*pdtop + eta2(izB+1)*pdB(ix) + ptop)
2569 pB=(pB2+pB1)*0.5
2571 pA2=log(eta1(izA)*pdtop + eta2(izA)*pdA(ix) + ptop)
2572 pA1=log(eta1(izA+1)*pdtop + eta2(izA+1)*pdA(ix) + ptop)
2573 pA=(pA2+pA1)*0.5
2575 ! Find the lowest mass level izA that is above pB
2576 interpinit: do while(izA>1)
2577 pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2578 pnext=(pA2+pA3)*0.5
2579 wdenom=pnext-pA
2580 wnum=pnext-pb
2581 if(pA<=pB .and. wnum>1e-5) then
2582 exit interpinit
2583 else
2584 pA1=pA2
2585 pA2=pa3
2586 izA=izA-1
2587 pA=pnext
2588 endif
2589 enddo interpinit
2591 ! Loop over all remaining B points, interpolating to them.
2592 interploop: do while(izB>0 .and. izA>1)
2593 ! Find the source domain levels that this target domain level lies between:
2594 zinterp: do while(izA>1)
2595 if(pnext>=pB) then
2596 ! We found the two levels, so interpolate and move on to
2597 ! the next target level:
2598 weight=max(0.,wnum/wdenom)
2599 tB(izB,ix)=weight*tA(izA,ix) + (1.0-weight)*tA(izA-1,ix)
2600 qB(izB,ix)=weight*qA(izA,ix) + (1.0-weight)*qA(izA-1,ix)
2602 ! Save interpolation info
2603 iinfo(izB,ix)=izA ! upper level
2604 winfo(izB,ix)=weight ! linear interpolation weight
2606 ! We interpolated to a B point, so go to the next one:
2607 pB1=pB2
2608 izB=izB-1
2609 if(izB>0) then
2610 pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2611 pB=(pb1+pb2)*0.5
2612 else
2613 exit interploop
2614 endif
2615 wnum=pnext-pb
2617 exit zinterp
2618 else
2619 izA=izA-1
2620 pA1=pA2
2621 pA2=pa3
2622 pA=pnext
2623 if(izA>1) then
2624 pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2625 pnext=(pA2+pA3)*0.5
2626 wdenom=pnext-pa
2627 wnum=pnext-pb
2628 else
2629 exit interploop
2630 endif
2631 endif
2632 enddo zinterp
2633 enddo interploop
2635 ! Follow what base_state_parent did for temperature:
2636 ! interpolate between the second to last level and P_REF using
2637 ! a lapse rate atmosphere. For Q, we use constant RH below
2638 ! ground, as suggested by Greg Thompson.
2639 extraploop: do while(izB>=1)
2640 ! Decide extrapolation weight:
2641 weight=(pB-pA)/(pstd12-pA)
2643 ! Extrapolate Q by copying lowest sigma layer value:
2644 !qB(izB,ix)=qA(1,ix)
2646 ! Extrapolate Q by linearly interpolating between 0 and surface value:
2647 !qB(izB,ix)=(1.0-weight)*qA(1,ix)
2649 ! Extrapolate T using a lapse rate atmosphere
2650 tB(izB,ix)=weight*tbelow(ix) + (1.0-weight)*tA(1,ix)
2652 ! Extrapolate Q using constant RH below ground, as suggested
2653 ! by Greg Thompson. This is the inversion of the RH
2654 ! calculation used earlier in this function:
2655 P0=eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop
2656 QC=PQ0/P0*EXP(A2*(TB(izB,ix)-A3)/(TB(izB,ix)-A4))
2657 QB(izB,ix)=QC*RHbelow(ix)
2659 ! Save extrapolation information
2660 iinfo(izB,ix)=0
2661 winfo(izB,ix)=weight
2663 ! Move to the next B level
2664 izB=izB-1
2665 if(izB>0) then
2666 pB1=pB2
2667 pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2668 pB=(pb1+pb2)*0.5
2669 endif
2670 enddo extraploop
2671 enddo outer
2672 end subroutine interp_T_PD_Q
2674 subroutine interp_T_PD_Q_kpres(method, pd_interp, nx, nz, &
2675 deta1,deta2, eta1,eta2, ptop,pdtop, kpres, nz2, &
2676 fisA,fisB, pintA,pintB, tA,tB, pdA,pdB, qA,qB, &
2677 iinfo, winfo, warned)
2678 implicit none
2680 integer, intent(in) :: pd_interp,method, kpres, nz2
2681 ! real, intent(in) :: dtA,dtB
2683 ! Coordinate system definitions must be the same for all domains:
2684 real, intent(in) :: deta1(nz),deta2(nz),eta1(nz),eta2(nz),ptop,pdtop
2685 integer, intent(in) :: nx,nz
2687 ! Surface height and mass field information for source (A):
2688 real, intent(in) :: fisA(nx),tA(nz2-1,nx),pdA(nx),qA(nz2-1,nx),pintA(nz2,nx)
2690 ! Surface height and mass field information for target (B):
2691 real, intent(inout) :: fisB(nx),tB(nz2-1,nx),pdB(nx),qB(nz2-1,nx),pintB(nz2,nx)
2693 ! Interpolation or extrapolation information for use in other
2694 ! calls later on:
2695 real, intent(out) :: winfo(nz2-1,nx)
2696 integer, intent(out) :: iinfo(nz2-1,nx)
2698 logical, intent(inout) :: warned
2700 ! ==================== Local variables ====================
2702 character*255 :: message
2703 integer :: ix,iz,izA,izB,xpr
2704 real :: zA,zB,znext,apelp,rtopp,dz,weight,A,B,pstd1,z,pstd2,pstd12
2705 real :: tsfc(nx), slp(nx), zmslp(nx), z0mid(nx), zbelow(nx), &
2706 tbelow(nx), RHbelow(nx)
2707 real :: pb,pb1,pb2,pa,pa1,pa2,pnext,pa3, wnum,wdenom, QC, P0
2709 ! Constants from UPP params.f for RH calculation (all mks):
2710 real, parameter :: PQ0=379.90516
2711 real, parameter :: A2=17.2693882
2712 real, parameter :: A3=273.16
2713 real, parameter :: A4=35.86
2714 real, parameter :: RHmin=1.0E-6 ! minimal RH bound
2716 if(method/=nmm_method_linear) then
2717 call wrf_error_fatal('only linear interpolation is supported')
2718 endif
2720 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2721 !! Step 1: calculate near-surface values !!!!!!!!!!!!!!!!!!!!!
2722 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2726 pstd1=p_ref ! pstd(1) from base_state_parent
2727 ! pstd(2) from base_state_parent:
2728 pstd2=eta1(2)*pdtop + eta2(2)*(p_ref-pdtop-ptop) + ptop
2729 pstd12=exp((alog(pstd1)+alog(pstd2))*0.5)
2731 do ix=1,nx
2732 ! These calculations are from base_state_parent:
2733 APELP = (pintA(2,ix)+pintA(1,ix))
2734 RTOPP = TRG*tA(1,ix)*(1.0+qA(1,ix)*P608)/APELP
2735 DZ = RTOPP*(DETA1(1)*PDTOP+DETA2(1)*pdA(ix))
2737 Z0MID(ix) = fisA(ix)/g + dz/2
2739 zA=fisA(ix)/g
2741 TSFC(ix) = TA(1,ix)*(1.+D608*QA(1,ix)) &
2742 + LAPSR*(zA + zA+dz)*0.5
2744 A = LAPSR*zA/TSFC(ix)
2745 SLP(ix) = PINTA(1,ix)*(1-A)**COEF2 ! sea level pressure
2746 B = (pstd1/SLP(ix))**COEF3
2747 ZMSLP(ix)= TSFC(ix)*LAPSI*(1.0 - B) ! Height at 1030. mb level
2749 TBELOW(ix) = TA(1,ix) + LAPSR*(Z0MID(ix)-ZMSLP(ix))
2751 ! Calculate lowest level RH. This calculation is from UPP
2752 ! CALRH.f:
2753 P0=pdA(ix)+ptop+pdtop ! Use hydrostatic lowest model level P
2754 QC=PQ0/P0*EXP(A2*(TA(1,ix)-A3)/(TA(1,ix)-A4))
2755 RHbelow(ix)=max(RHmin,min(1.0,QA(1,ix)/QC))
2756 enddo
2760 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2761 !! Step 2: figure out the new surface pressures !!!!!!!!!!!!!!
2762 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2766 if_pd_interp: if(pd_interp==nmm_keep_pd) then
2767 ! PD interpolation is turned off, so we use the original PD.
2768 elseif(pd_interp==nmm_interp_pd) then
2769 ! PD interpolation is requested. As with the old base_state_parent,
2770 ! determine PD by interpolating or extrapolating the non-hydrostatic
2771 ! pressure field (PINT) in source grid to the surface height of the
2772 ! target grid.
2773 xloop: do ix=1,nx
2774 if(pintA(1,ix)>p_ref) then
2775 ! Cannot extrapolate PD below ground because pressure is
2776 ! unrealistically high. Follow base_state_parent method:
2777 ! when this happens, assign input pressure to output
2778 ! pressure.
2779 if(.not.warned) then
2780 call wrf_message2('WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD')
2781 WRITE(message,*) 'PINT(1),PD(1),PSTD(1)',pintA(1,ix),pdA(ix),p_ref
2782 call wrf_message2(message)
2783 warned=.true.
2784 endif
2785 pdB(ix)=pdA(ix)
2786 cycle xloop
2787 endif
2789 zA=fisA(ix)/g
2790 zB=fisB(ix)/g
2792 if(zB<zA) then
2793 ! Need to extrapolate PD below ground.
2794 ! Follow base_state_parent method:
2795 ! For pressures between pint(1) and 1030 mbar, interpolate
2796 iz=1
2797 weight=abs((zB-zmslp(ix))/(zA-zmslp(ix)))
2798 pdB(ix) = exp(weight*log(pdA(ix)+pdtop+ptop) + &
2799 (1.0-weight)*log(P_REF)) &
2800 - pdtop - ptop
2801 ! write(0,*) 'WARNING: extrapolating below ground at ix=',ix
2802 !1302 format('Extrap at ix=',I0,' zA=',F0.3,' zB=',F0.3,' zmslp=',F0.3,' weight=',F0.3,' pdB+pdtop+ptop=',F0.3,' pintA=',F0.3,' pdA(ix)+pdtop+ptop=',F0.3,' P_REF=',F0.3)
2803 ! write(0,1302) ix,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2804 cycle xloop
2805 endif
2807 ! PD in grid B lies above ground in grid A. Find out where by
2808 ! walking up from the ground until we find the two levels the
2809 ! pressure lies between.
2810 znext=-9e9
2811 z=zA
2812 do iz=1,nz2-1
2813 ! Calculations from base_state_parent:
2814 APELP = (pintA(iz+1,ix)+pintA(iz,ix))
2815 RTOPP = TRG*tA(iz,ix)*(1.0+qA(iz,ix)*P608)/APELP
2816 DZ = RTOPP*(DETA1(iz)*PDTOP+DETA2(iz)*pdA(ix))
2817 znext=z+dz
2819 if(znext>zB) then
2820 ! Interpolate between these two levels
2821 weight=(zB-z)/dz
2822 pdB(ix) = exp(weight*log(pintA(iz+1,ix)) + &
2823 (1.0-weight)*log(pintA(iz,ix))) &
2824 - pdtop - ptop
2825 ! if(iz>1) then
2826 ! write(0,*) 'WARNING: interpolate above ground at ix=',ix
2827 !1303 format('Interp at ix=',I0,' with iz=',I0,' zA=',F0.3,' zB=',F0.3,' zmslp=',F0.3,' weight=',F0.3,' pdB+pdtop+ptop=',F0.3,' pintA=',F0.3,' pdA(ix)+pdtop+ptop=',F0.3,' P_REF=',F0.3)
2828 ! write(0,1303) ix,iz,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2829 ! endif
2830 cycle xloop
2831 else
2832 z=znext
2833 endif
2835 enddo
2837 201 format('interp_T_PD_Q_kpres: Target domain surface height ',F0.4,'m is higher than the model top ',F0.4,'m of the source domain. ABORTING')
2838 write(message,201) zB,znext
2839 call wrf_error_fatal(message)
2840 enddo xloop
2841 else
2842 202 format('Invalid value ',I0,' for pd_interp in interp_T_PD_Q')
2843 write(message,202) pd_interp
2844 call wrf_error_fatal(message)
2845 endif if_pd_interp
2849 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2850 !! Step 3: interpolate T and Q in pressure space !!!!!!!!!!!!!
2851 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2855 outer: do ix=1,nx
2856 ! For constant pressure levels, do a straight level-by-level copy:
2857 iz=nz2-1
2858 copyloop: do while(iz>kpres)
2859 tB(iz,ix)=tA(iz,ix)
2860 qB(iz,ix)=qA(iz,ix)
2862 ! Save interpolation information
2863 iinfo(iz,ix)=iz
2864 winfo(iz,ix)=1.0
2866 iz=iz-1
2867 enddo copyloop
2869 ! For sigma levels where target domain lies within source,
2870 ! interpolate from top down, stopping when we go below ground
2871 ! in the source domain.
2872 izA=iz
2873 izB=iz
2875 ! FIXME: REMOVE THIS CHECK
2876 if(iz>nz2) then
2877 ! Make sure the entire vertical extent isn't sigma levels:
2878 call wrf_error_fatal('ERROR: WRF-NMM does not support pure sigma levels (only sigma-pressure hybrid)')
2879 endif
2881 pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2882 pB1=log(eta1(izB+1)*pdtop + eta2(izB+1)*pdB(ix) + ptop)
2883 pB=(pB2+pB1)*0.5
2885 pA2=log(eta1(izA)*pdtop + eta2(izA)*pdA(ix) + ptop)
2886 pA1=log(eta1(izA+1)*pdtop + eta2(izA+1)*pdA(ix) + ptop)
2887 pA=(pA2+pA1)*0.5
2889 ! Find the lowest mass level izA that is above pB
2890 interpinit: do while(izA>1)
2891 pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2892 pnext=(pA2+pA3)*0.5
2893 wdenom=pnext-pA
2894 wnum=pnext-pb
2895 if(pA<=pB .and. wnum>1e-5) then
2896 exit interpinit
2897 else
2898 pA1=pA2
2899 pA2=pa3
2900 izA=izA-1
2901 pA=pnext
2902 endif
2903 enddo interpinit
2905 ! Loop over all remaining B points, interpolating to them.
2906 interploop: do while(izB>0 .and. izA>1)
2907 ! Find the source domain levels that this target domain level lies between:
2908 zinterp: do while(izA>1)
2909 if(pnext>=pB) then
2910 ! We found the two levels, so interpolate and move on to
2911 ! the next target level:
2912 weight=max(0.,wnum/wdenom)
2913 tB(izB,ix)=weight*tA(izA,ix) + (1.0-weight)*tA(izA-1,ix)
2914 qB(izB,ix)=weight*qA(izA,ix) + (1.0-weight)*qA(izA-1,ix)
2916 ! Save interpolation info
2917 iinfo(izB,ix)=izA ! upper level
2918 winfo(izB,ix)=weight ! linear interpolation weight
2920 ! We interpolated to a B point, so go to the next one:
2921 pB1=pB2
2922 izB=izB-1
2923 if(izB>0) then
2924 pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2925 pB=(pb1+pb2)*0.5
2926 else
2927 exit interploop
2928 endif
2929 wnum=pnext-pb
2931 exit zinterp
2932 else
2933 izA=izA-1
2934 pA1=pA2
2935 pA2=pa3
2936 pA=pnext
2937 if(izA>1) then
2938 pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2939 pnext=(pA2+pA3)*0.5
2940 wdenom=pnext-pa
2941 wnum=pnext-pb
2942 else
2943 exit interploop
2944 endif
2945 endif
2946 enddo zinterp
2947 enddo interploop
2949 ! Follow what base_state_parent did for temperature:
2950 ! interpolate between the second to last level and P_REF using
2951 ! a lapse rate atmosphere. For Q, we use constant RH below
2952 ! ground, as suggested by Greg Thompson.
2953 extraploop: do while(izB>=1)
2954 ! Decide extrapolation weight:
2955 weight=(pB-pA)/(pstd12-pA)
2957 ! Extrapolate Q by copying lowest sigma layer value:
2958 !qB(izB,ix)=qA(1,ix)
2960 ! Extrapolate Q by linearly interpolating between 0 and surface value:
2961 !qB(izB,ix)=(1.0-weight)*qA(1,ix)
2963 ! Extrapolate T using a lapse rate atmosphere
2964 tB(izB,ix)=weight*tbelow(ix) + (1.0-weight)*tA(1,ix)
2966 ! Extrapolate Q using constant RH below ground, as suggested
2967 ! by Greg Thompson. This is the inversion of the RH
2968 ! calculation used earlier in this function:
2969 P0=eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop
2970 QC=PQ0/P0*EXP(A2*(TB(izB,ix)-A3)/(TB(izB,ix)-A4))
2971 QB(izB,ix)=QC*RHbelow(ix)
2973 ! Save extrapolation information
2974 iinfo(izB,ix)=0
2975 winfo(izB,ix)=weight
2977 ! Move to the next B level
2978 izB=izB-1
2979 if(izB>0) then
2980 pB1=pB2
2981 pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2982 pB=(pb1+pb2)*0.5
2983 endif
2984 enddo extraploop
2985 enddo outer
2986 end subroutine interp_T_PD_Q_kpres
2987 end module module_interp_nmm
2989 ! ********************************************************************
2990 ! subs EXT_*_FULLDOM -- simple wrappers around module_interp_nmm's
2991 ! *_FULLDOM routines. These exist solely to allow module_dm to
2992 ! interface with this module -- this is needed because the WRF build
2993 ! scripts compile module_dm before module_interp_nmm.
2994 ! ********************************************************************
2996 subroutine ext_c2n_fulldom (II,JJ,W1,W2,W3,W4,&
2997 deta1,deta2, eta1,eta2, ptop,pdtop, &
2998 cpint,ct,cpd,cq, &
2999 cims, cime, cjms, cjme, ckms, ckme, &
3000 npint,nt,npd,nq, &
3001 out_iinfo,out_winfo,imask,&
3002 nids, nide, njds, njde, nkds, nkde, &
3003 nims, nime, njms, njme, nkms, nkme, &
3004 nits, nite, njts, njte, nkts, nkte)
3005 ! This subroutine is an alias for c2n_fulldom. It exists to allow
3006 ! the routine to be called from module_dm
3007 use module_interp_store, only: parent_fis, nest_fis, kpres
3008 use module_interp_nmm, only: c2n_fulldom, find_kpres
3009 implicit none
3010 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
3011 nids, nide, njds, njde, nkds, nkde, &
3012 nims, nime, njms, njme, nkms, nkme, &
3013 nits, nite, njts, njte, nkts, nkte
3014 real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
3015 real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD
3017 real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_winfo
3018 integer, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_iinfo
3020 real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
3021 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
3022 real, intent(inout), dimension(nims:nime,njms:njme,1) :: nPD
3023 integer, intent(in), dimension(nims:nime,njms:njme) :: imask
3025 real, intent(in), dimension(nims:nime,njms:njme) :: &
3026 W1,W2,W3,W4
3027 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
3029 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
3030 real, intent(in) :: ptop,pdtop
3031 integer :: i,j,k
3032 logical :: badbad
3034 call find_kpres(kpres,eta2,nkds,nkde,nkms,nkme)
3035 call c2n_fulldom(II,JJ,W1,W2,W3,W4,&
3036 deta1,deta2, eta1,eta2, ptop,pdtop, &
3037 parent_fis,cpint,ct,cpd,cq, &
3038 cims, cime, cjms, cjme, ckms, ckme, &
3039 nest_fis,npint,nt,npd,nq, &
3040 out_iinfo,out_winfo,imask, &
3041 nids, nide, njds, njde, nkds, nkde, &
3042 nims, nime, njms, njme, nkms, nkme, &
3043 nits, nite, njts, njte, nkts, nkte, &
3044 kpres)
3046 end subroutine ext_c2n_fulldom
3048 subroutine ext_n2c_fulldom ( &
3049 deta1,deta2, eta1,eta2, ptop,pdtop, &
3050 cpint,ct,cpd,cq, &
3051 cids, cide, cjds, cjde, ckds, ckde, &
3052 cims, cime, cjms, cjme, ckms, ckme, &
3053 cits, cite, cjts, cjte, ckts, ckte, &
3054 npint,nt,npd,nq, &
3055 ipos,jpos, &
3056 out_iinfo,out_winfo, &
3057 nids, nide, njds, njde, nkds, nkde, &
3058 nims, nime, njms, njme, nkms, nkme, &
3059 nits, nite, njts, njte, nkts, nkte)
3060 ! This subroutine is an alias for n2c_fulldom. It exists to allow
3061 ! the routine to be called from module_dm
3062 use module_interp_store, only: parent_fis, nest_fis, kpres
3063 use module_interp_nmm, only: n2c_fulldom, find_kpres, n2c_fulldom_new
3064 implicit none
3065 integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
3066 cids, cide, cjds, cjde, ckds, ckde, &
3067 cits, cite, cjts, cjte, ckts, ckte, &
3068 nids, nide, njds, njde, nkds, nkde, &
3069 nims, nime, njms, njme, nkms, nkme, &
3070 nits, nite, njts, njte, nkts, nkte, &
3071 ipos, jpos ! nest start loctions within parent
3072 real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
3073 real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD
3075 real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
3076 integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
3078 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
3079 real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
3080 real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
3082 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
3083 real, intent(in) :: ptop,pdtop
3085 call find_kpres(kpres,eta2,nkds,nkde,nkms,nkme)
3086 call n2c_fulldom_new( &
3087 deta1,deta2, eta1,eta2, ptop,pdtop, &
3088 parent_fis,cpint,ct,cpd,cq, &
3089 cids, cide, cjds, cjde, ckds, ckde, &
3090 cims, cime, cjms, cjme, ckms, ckme, &
3091 cits, cite, cjts, cjte, ckts, ckte, &
3092 nest_fis,npint,nt,npd,nq, &
3093 ipos,jpos, &
3094 out_iinfo,out_winfo, &
3095 nids, nide, njds, njde, nkds, nkde, &
3096 nims, nime, njms, njme, nkms, nkme, &
3097 nits, nite, njts, njte, nkts, nkte, kpres)
3098 end subroutine ext_n2c_fulldom
3100 subroutine ext_c2b_fulldom (II,JJ,W1,W2,W3,W4,&
3101 deta1,deta2, eta1,eta2, ptop,pdtop, &
3102 cpint,ct,cpd,cq, &
3103 cims, cime, cjms, cjme, ckms, ckme, &
3104 nids, nide, njds, njde, nkds, nkde, &
3105 nims, nime, njms, njme, nkms, nkme, &
3106 nits, nite, njts, njte, nkts, nkte, &
3107 ibxs, ibxe, ibys, ibye, &
3108 wbxs, wbxe, wbys, wbye, &
3109 pdbxs, pdbxe, pdbys, pdbye, &
3110 tbxs, tbxe, tbys, tbye, &
3111 qbxs, qbxe, qbys, qbye)
3112 use module_interp_nmm, only: c2b_fulldom, find_kpres, c2b_fulldom_new
3113 use module_interp_store, only: parent_fis, nest_fis, kpres
3114 implicit none
3115 integer, intent(in):: &
3116 cims, cime, cjms, cjme, ckms, ckme, &
3117 nids, nide, njds, njde, nkds, nkde, &
3118 nims, nime, njms, njme, nkms, nkme, &
3119 nits, nite, njts, njte, nkts, nkte
3120 integer, parameter :: bdyw = 1
3121 ! Domain info:
3122 real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
3123 integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
3124 real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
3125 real, intent(in) :: ptop,pdtop
3127 ! Parent fields:
3128 real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
3129 real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD
3131 ! T, Q, PINT, PD boundary info:
3132 real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
3133 real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
3134 real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
3135 real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
3137 ! Weights and indices:
3138 real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
3139 real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
3140 integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
3141 integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
3143 call find_kpres(kpres,eta2,nkds,nkde,nkms,nkme)
3144 call c2b_fulldom_new (II,JJ,W1,W2,W3,W4,&
3145 deta1,deta2, eta1,eta2, ptop,pdtop, &
3146 parent_fis,cpint,ct,cpd,cq,nest_fis, &
3147 cims, cime, cjms, cjme, ckms, ckme, &
3148 nids, nide, njds, njde, nkds, nkde, &
3149 nims, nime, njms, njme, nkms, nkme, &
3150 nits, nite, njts, njte, nkts, nkte, &
3151 kpres, &
3152 ibxs, ibxe, ibys, ibye, &
3153 wbxs, wbxe, wbys, wbye, &
3154 pdbxs, pdbxe, pdbys, pdbye, &
3155 tbxs, tbxe, tbys, tbye, &
3156 qbxs, qbxe, qbys, qbye)
3157 end subroutine ext_c2b_fulldom