UPP  001
 All Data Structures Files Functions Pages
CALVOR.f
Go to the documentation of this file.
1 
26  SUBROUTINE calvor(UWND,VWND,ABSV)
27 
28 !
29 !
30  use vrbls2d, only: f
31  use masks, only: gdlat, gdlon, dx, dy
32  use params_mod, only: d00, dtr, small, erad
33  use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, global, &
34  jsta, jend, im, jm, jsta_m, jend_m, gdsdegr
35  use gridspec_mod, only: gridtype, dyval
36  use upp_math, only: dvdxdudy, ddvdx, ddudy, uuavg
37 
38  implicit none
39 !
40 ! DECLARE VARIABLES.
41 !
42  REAL, dimension(im,jsta_2l:jend_2u), intent(in) :: uwnd, vwnd
43  REAL, dimension(im,jsta_2l:jend_2u), intent(inout) :: absv
44 !
45  real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), cosl(:,:)
46  INTEGER, allocatable :: ihe(:),ihw(:), ie(:),iw(:)
47 !
48  integer, parameter :: npass2=2, npass3=3
49  integer i,j,ip1,im1,ii,iir,iil,jj,jmt2,imb2, npass, nn, jtem
50  real r2dx,r2dy,dvdx,dudy,uavg,tph1,tphi, tx1(im+2), tx2(im+2)
51 !
52 !***************************************************************************
53 ! START CALVOR HERE.
54 !
55 ! LOOP TO COMPUTE ABSOLUTE VORTICITY FROM WINDS.
56 !
57  IF(modelname == 'RAPR') then
58 !$omp parallel do private(i,j)
59  DO j=jsta_2l,jend_2u
60  DO i=1,im
61  absv(i,j) = d00
62  ENDDO
63  ENDDO
64  else
65 !$omp parallel do private(i,j)
66  DO j=jsta_2l,jend_2u
67  DO i=1,im
68  absv(i,j) = spval
69  ENDDO
70  ENDDO
71  endif
72 
73 ! print*,'dyval in CALVOR= ',DYVAL
74 
75  CALL exch_f(uwnd)
76 !
77  IF (modelname == 'GFS' .or. global) THEN
78  CALL exch(gdlat(1,jsta_2l))
79 
80  allocate (wrk1(im,jsta:jend), wrk2(im,jsta:jend), &
81  & wrk3(im,jsta:jend), cosl(im,jsta_2l:jend_2u))
82  allocate(iw(im),ie(im))
83 
84  imb2 = im/2
85 !$omp parallel do private(i)
86  do i=1,im
87  ie(i) = i+1
88  iw(i) = i-1
89  enddo
90  iw(1) = im
91  ie(im) = 1
92 
93 ! if(1>=jsta .and. 1<=jend)then
94 ! if(cos(gdlat(1,1)*dtr)<small)poleflag=.T.
95 ! end if
96 ! call mpi_bcast(poleflag,1,MPI_LOGICAL,0,mpi_comm_comp,iret)
97 
98 !$omp parallel do private(i,j,ip1,im1)
99  DO j=jsta,jend
100  do i=1,im
101  ip1 = ie(i)
102  im1 = iw(i)
103  cosl(i,j) = cos(gdlat(i,j)*dtr)
104  IF(cosl(i,j) >= small) then
105  wrk1(i,j) = 1.0 / (erad*cosl(i,j))
106  else
107  wrk1(i,j) = 0.
108  end if
109  if(i == im .or. i == 1) then
110  wrk2(i,j) = 1.0 / ((360.+gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
111  else
112  wrk2(i,j) = 1.0 / ((gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
113  end if
114  enddo
115  enddo
116 ! CALL EXCH(cosl(1,JSTA_2L))
117  CALL exch(cosl)
118 
119 !$omp parallel do private(i,j,ii)
120  DO j=jsta,jend
121  if (j == 1) then
122  if(gdlat(1,j) > 0.) then ! count from north to south
123  do i=1,im
124  ii = i + imb2
125  if (ii > im) ii = ii - im
126  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j+1)-gdlat(ii,j))*dtr) !1/dphi
127  enddo
128  else ! count from south to north
129  do i=1,im
130  ii = i + imb2
131  if (ii > im) ii = ii - im
132  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j+1)+gdlat(ii,j))*dtr) !1/dphi
133 !
134  enddo
135  end if
136  elseif (j == jm) then
137  if(gdlat(1,j) < 0.) then ! count from north to south
138  do i=1,im
139  ii = i + imb2
140  if (ii > im) ii = ii - im
141  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j-1)+gdlat(ii,j))*dtr)
142  enddo
143  else ! count from south to north
144  do i=1,im
145  ii = i + imb2
146  if (ii > im) ii = ii - im
147  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j-1)-gdlat(ii,j))*dtr)
148  enddo
149  end if
150  else
151  do i=1,im
152  wrk3(i,j) = 1.0 / ((gdlat(i,j-1)-gdlat(i,j+1))*dtr) !1/dphi
153  enddo
154  endif
155  enddo
156 
157  npass = 0
158 
159  jtem = jm / 18 + 1
160 !$omp parallel do private(i,j,ip1,im1,ii,jj,tx1,tx2)
161  DO j=jsta,jend
162 ! npass = npass2
163 ! if (j > jm-jtem+1 .or. j < jtem) npass = npass3
164  IF(j == 1) then ! Near North or South pole
165  if(gdlat(1,j) > 0.) then ! count from north to south
166  IF(cosl(1,j) >= small) THEN !not a pole point
167  DO i=1,im
168  ip1 = ie(i)
169  im1 = iw(i)
170  ii = i + imb2
171  if (ii > im) ii = ii - im
172  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
173  uwnd(ii,j)==spval .or. uwnd(i,j+1)==spval) cycle
174  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
175  & + (uwnd(ii,j)*cosl(ii,j) &
176  & + uwnd(i,j+1)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j) &
177  & + f(i,j)
178  enddo
179  ELSE !pole point, compute at j=2
180  jj = 2
181  DO i=1,im
182  ip1 = ie(i)
183  im1 = iw(i)
184  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
185  uwnd(i,j)==spval .or. uwnd(i,jj+1)==spval) cycle
186  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
187  & - (uwnd(i,j)*cosl(i,j) &
188  - uwnd(i,jj+1)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj) &
189  & + f(i,jj)
190  enddo
191  ENDIF
192  else
193  IF(cosl(1,j) >= small) THEN !not a pole point
194  DO i=1,im
195  ip1 = ie(i)
196  im1 = iw(i)
197  ii = i + imb2
198  if (ii > im) ii = ii - im
199  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
200  uwnd(ii,j)==spval .or. uwnd(i,j+1)==spval) cycle
201  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
202  & - (uwnd(ii,j)*cosl(ii,j) &
203  & + uwnd(i,j+1)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j) &
204  & + f(i,j)
205  enddo
206  ELSE !pole point, compute at j=2
207  jj = 2
208  DO i=1,im
209  ip1 = ie(i)
210  im1 = iw(i)
211  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
212  uwnd(i,j)==spval .or. uwnd(i,jj+1)==spval) cycle
213  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
214  & + (uwnd(i,j)*cosl(i,j) &
215  - uwnd(i,jj+1)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj) &
216  & + f(i,jj)
217  enddo
218  ENDIF
219  endif
220  ELSE IF(j == jm) THEN ! Near North or South Pole
221  if(gdlat(1,j) < 0.) then ! count from north to south
222  IF(cosl(1,j) >= small) THEN !not a pole point
223  DO i=1,im
224  ip1 = ie(i)
225  im1 = iw(i)
226  ii = i + imb2
227  if (ii > im) ii = ii - im
228  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
229  uwnd(i,j-1)==spval .or. uwnd(ii,j)==spval) cycle
230  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
231  & - (uwnd(i,j-1)*cosl(i,j-1) &
232  & + uwnd(ii,j)*cosl(ii,j))*wrk3(i,j)) * wrk1(i,j) &
233  & + f(i,j)
234  enddo
235  ELSE !pole point,compute at jm-1
236  jj = jm-1
237  DO i=1,im
238  ip1 = ie(i)
239  im1 = iw(i)
240  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
241  uwnd(i,jj-1)==spval .or. uwnd(i,j)==spval) cycle
242  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
243  & - (uwnd(i,jj-1)*cosl(i,jj-1) &
244  & - uwnd(i,j)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj) &
245  & + f(i,jj)
246  enddo
247  ENDIF
248  else
249  IF(cosl(1,j) >= small) THEN !not a pole point
250  DO i=1,im
251  ip1 = ie(i)
252  im1 = iw(i)
253  ii = i + imb2
254  if (ii > im) ii = ii - im
255  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
256  uwnd(i,j-1)==spval .or. uwnd(ii,j)==spval) cycle
257  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
258  & + (uwnd(i,j-1)*cosl(i,j-1) &
259  & + uwnd(ii,j)*cosl(ii,j))*wrk3(i,j)) * wrk1(i,j) &
260  & + f(i,j)
261  enddo
262  ELSE !pole point,compute at jm-1
263  jj = jm-1
264  DO i=1,im
265  ip1 = ie(i)
266  im1 = iw(i)
267  if(vwnd(ip1,jj)==spval .or. vwnd(im1,jj)==spval .or. &
268  uwnd(i,jj-1)==spval .or. uwnd(i,j)==spval) cycle
269  absv(i,j) = ((vwnd(ip1,jj)-vwnd(im1,jj))*wrk2(i,jj) &
270  & + (uwnd(i,jj-1)*cosl(i,jj-1) &
271  & - uwnd(i,j)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj) &
272  & + f(i,jj)
273  enddo
274  ENDIF
275  endif
276  ELSE
277  DO i=1,im
278  ip1 = ie(i)
279  im1 = iw(i)
280  if(vwnd(ip1,j)==spval .or. vwnd(im1,j)==spval .or. &
281  uwnd(i,j-1)==spval .or. uwnd(i,j+1)==spval) cycle
282  absv(i,j) = ((vwnd(ip1,j)-vwnd(im1,j))*wrk2(i,j) &
283  & - (uwnd(i,j-1)*cosl(i,j-1) &
284  - uwnd(i,j+1)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j) &
285  + f(i,j)
286  ENDDO
287  END IF
288 ! if(ABSV(I,J)>1.0)print*,'Debug CALVOR',i,j,VWND(ip1,J),VWND(im1,J), &
289 ! wrk2(i,j),UWND(I,J-1),COSL(I,J-1),UWND(I,J+1),COSL(I,J+1),wrk3(i,j),cosl(i,j),F(I,J),ABSV(I,J)
290  if (npass > 0) then
291  do i=1,im
292  tx1(i) = absv(i,j)
293  enddo
294  do nn=1,npass
295  do i=1,im
296  tx2(i+1) = tx1(i)
297  enddo
298  tx2(1) = tx2(im+1)
299  tx2(im+2) = tx2(2)
300  do i=2,im+1
301  tx1(i-1) = 0.25 * (tx2(i-1) + tx2(i+1)) + 0.5*tx2(i)
302  enddo
303  enddo
304  do i=1,im
305  absv(i,j) = tx1(i)
306  enddo
307  endif
308  END DO ! end of J loop
309 
310 ! deallocate (wrk1, wrk2, wrk3, cosl)
311 ! GFS use lon avg as one scaler value for pole point
312 
313  call poleavg(im,jm,jsta,jend,small,cosl(1,jsta),spval,absv(1,jsta))
314  deallocate (wrk1, wrk2, wrk3, cosl, iw, ie)
315 
316  ELSE !(MODELNAME == 'GFS' .or. global)
317 
318  IF (gridtype == 'B')THEN
319  CALL exch_f(vwnd)
320  ENDIF
321 
322  CALL dvdxdudy(uwnd,vwnd)
323 
324  IF(gridtype == 'A')THEN
325 !$omp parallel do private(i,j,jmt2,tphi,r2dx,r2dy,dvdx,dudy,uavg)
326  DO j=jsta_m,jend_m
327  jmt2 = jm/2+1
328  tphi = (j-jmt2)*(dyval/gdsdegr)*dtr
329  DO i=2,im-1
330  IF(vwnd(i+1,j)<spval.AND.vwnd(i-1,j)<spval.AND. &
331  & uwnd(i,j+1)<spval.AND.uwnd(i,j-1)<spval) THEN
332  dvdx = ddvdx(i,j)
333  dudy = ddudy(i,j)
334  uavg = uuavg(i,j)
335 ! is there a (f+tan(phi)/erad)*u term?
336  IF(modelname == 'RAPR') then
337  absv(i,j) = dvdx - dudy + f(i,j) ! for run RAP over north pole
338  else
339  absv(i,j) = dvdx - dudy + f(i,j) + uavg*tan(gdlat(i,j)*dtr)/erad ! not sure about this???
340  endif
341  END IF
342  END DO
343  END DO
344 
345  ELSE IF (gridtype == 'E')THEN
346  allocate(ihw(jsta_2l:jend_2u), ihe(jsta_2l:jend_2u))
347 !$omp parallel do private(j)
348  DO j=jsta_2l,jend_2u
349  ihw(j) = -mod(j,2)
350  ihe(j) = ihw(j)+1
351  ENDDO
352 !$omp parallel do private(i,j,jmt2,tphi,r2dx,r2dy,dvdx,dudy,uavg)
353  DO j=jsta_m,jend_m
354  jmt2 = jm/2+1
355  tphi = (j-jmt2)*(dyval/1000.)*dtr
356  tphi = (j-jmt2)*(dyval/gdsdegr)*dtr
357  DO i=2,im-1
358  IF(vwnd(i+ihe(j),j) < spval.AND.vwnd(i+ihw(j),j) < spval .AND. &
359  & uwnd(i,j+1) < spval .AND.uwnd(i,j-1) < spval) THEN
360  dvdx = ddvdx(i,j)
361  dudy = ddudy(i,j)
362  uavg = uuavg(i,j)
363 ! is there a (f+tan(phi)/erad)*u term?
364  absv(i,j) = dvdx - dudy + f(i,j) + uavg*tan(tphi)/erad
365  END IF
366  END DO
367  END DO
368  deallocate(ihw, ihe)
369  ELSE IF (gridtype == 'B')THEN
370 ! CALL EXCH_F(VWND) !done before dvdxdudy() Jesse 20200520
371  DO j=jsta_m,jend_m
372  jmt2 = jm/2+1
373  tphi = (j-jmt2)*(dyval/gdsdegr)*dtr
374  DO i=2,im-1
375  if(vwnd(i, j)==spval .or. vwnd(i, j-1)==spval .or. &
376  vwnd(i-1,j)==spval .or. vwnd(i-1,j-1)==spval .or. &
377  uwnd(i, j)==spval .or. uwnd(i-1,j)==spval .or. &
378  uwnd(i,j-1)==spval .or. uwnd(i-1,j-1)==spval) cycle
379  dvdx = ddvdx(i,j)
380  dudy = ddudy(i,j)
381  uavg = uuavg(i,j)
382 ! is there a (f+tan(phi)/erad)*u term?
383  absv(i,j) = dvdx - dudy + f(i,j) + uavg*tan(tphi)/erad
384  END DO
385  END DO
386  END IF
387  END IF
388 !
389 ! END OF ROUTINE.
390 !
391  RETURN
392  END
393 
410  SUBROUTINE caldiv(UWND,VWND,DIV)
411  use masks, only: gdlat, gdlon
412  use params_mod, only: d00, dtr, small, erad
413  use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, global, &
414  jsta, jend, im, jm, jsta_m, jend_m, lm
415  use gridspec_mod, only: gridtype
416 
417  implicit none
418 !
419 ! DECLARE VARIABLES.
420 !
421  REAL, dimension(im,jsta_2l:jend_2u,lm), intent(in) :: uwnd,vwnd
422  REAL, dimension(im,jsta:jend,lm), intent(inout) :: div
423 !
424  real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), cosl(:,:)
425  INTEGER, allocatable :: ihe(:),ihw(:), ie(:),iw(:)
426 !
427  real :: dnpole, dspole, tem
428  integer i,j,ip1,im1,ii,iir,iil,jj,imb2, l
429 !
430 !***************************************************************************
431 ! START CALDIV HERE.
432 !
433 ! LOOP TO COMPUTE DIVERGENCE FROM WINDS.
434 !
435  CALL exch(gdlat(1,jsta_2l))
436 
437  allocate (wrk1(im,jsta:jend), wrk2(im,jsta:jend), &
438  & wrk3(im,jsta:jend), cosl(im,jsta_2l:jend_2u))
439  allocate(iw(im),ie(im))
440 
441  imb2 = im/2
442 !$omp parallel do private(i)
443  do i=1,im
444  ie(i) = i+1
445  iw(i) = i-1
446  enddo
447  iw(1) = im
448  ie(im) = 1
449 
450 
451 !$omp parallel do private(i,j,ip1,im1)
452  DO j=jsta,jend
453  do i=1,im
454  ip1 = ie(i)
455  im1 = iw(i)
456  cosl(i,j) = cos(gdlat(i,j)*dtr)
457  IF(cosl(i,j) >= small) then
458  wrk1(i,j) = 1.0 / (erad*cosl(i,j))
459  else
460  wrk1(i,j) = 0.
461  end if
462  if(i == im .or. i == 1) then
463  wrk2(i,j) = 1.0 / ((360.+gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
464  else
465  wrk2(i,j) = 1.0 / ((gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
466  end if
467  enddo
468  ENDDO
469 
470  CALL exch(cosl)
471 
472 !$omp parallel do private(i,j,ii)
473  DO j=jsta,jend
474  if (j == 1) then
475  if(gdlat(1,j) > 0.) then ! count from north to south
476  do i=1,im
477  ii = i + imb2
478  if (ii > im) ii = ii - im
479  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j+1)-gdlat(ii,j))*dtr) !1/dphi
480  enddo
481  else ! count from south to north
482  do i=1,im
483  ii = i + imb2
484  if (ii > im) ii = ii - im
485  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j+1)+gdlat(ii,j))*dtr) !1/dphi
486  enddo
487  end if
488  elseif (j == jm) then
489  if(gdlat(1,j) < 0.) then ! count from north to south
490  do i=1,im
491  ii = i + imb2
492  if (ii > im) ii = ii - im
493  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j-1)+gdlat(ii,j))*dtr)
494  enddo
495  else ! count from south to north
496  do i=1,im
497  ii = i + imb2
498  if (ii > im) ii = ii - im
499  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j-1)-gdlat(ii,j))*dtr)
500  enddo
501  end if
502  else
503  do i=1,im
504  wrk3(i,j) = 1.0 / ((gdlat(i,j-1)-gdlat(i,j+1))*dtr) !1/dphi
505  enddo
506  endif
507  enddo
508 
509  do l=1,lm
510 !$omp parallel do private(i,j)
511  DO j=jsta,jend
512  DO i=1,im
513  div(i,j,l) = spval
514  ENDDO
515  ENDDO
516 
517  CALL exch_f(vwnd(1,jsta_2l,l))
518 
519 !$omp parallel do private(i,j,ip1,im1,ii,jj)
520  DO j=jsta,jend
521  IF(j == 1) then ! Near North pole
522  if(gdlat(1,j) > 0.) then ! count from north to south
523  IF(cosl(1,j) >= small) THEN !not a pole point
524  DO i=1,im
525  ip1 = ie(i)
526  im1 = iw(i)
527  ii = i + imb2
528  if (ii > im) ii = ii - im
529  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
530  & - (vwnd(ii,j,l)*cosl(ii,j) &
531  & + vwnd(i,j+1,l)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j)
532  enddo
533 !--
534  ELSE !North pole point, compute at j=2
535  jj = 2
536  do i=1,im
537  ip1 = ie(i)
538  im1 = iw(i)
539  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
540  & + (vwnd(i,j,l)*cosl(i,j) &
541  - vwnd(i,jj+1,l)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj)
542  enddo
543 !--
544  ENDIF
545  else
546  IF(cosl(1,j) >= small) THEN !not a pole point
547  DO i=1,im
548  ip1 = ie(i)
549  im1 = iw(i)
550  ii = i + imb2
551  if (ii > im) ii = ii - im
552  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
553  & + (vwnd(ii,j,l)*cosl(ii,j) &
554  & + vwnd(i,j+1,l)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j)
555  enddo
556 !--
557  ELSE !North pole point, compute at j=2
558  jj = 2
559  do i=1,im
560  ip1 = ie(i)
561  im1 = iw(i)
562  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
563  & - (vwnd(i,j,l)*cosl(i,j) &
564  - vwnd(i,jj+1,l)*cosl(i,jj+1))*wrk3(i,jj)) * wrk1(i,jj)
565  enddo
566  ENDIF
567  endif
568  ELSE IF(j == jm) THEN ! Near South pole
569  if(gdlat(1,j) < 0.) then ! count from north to south
570  IF(cosl(1,j) >= small) THEN !not a pole point
571  DO i=1,im
572  ip1 = ie(i)
573  im1 = iw(i)
574  ii = i + imb2
575  if (ii > im) ii = ii - im
576  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
577  & + (vwnd(i,j-1,l)*cosl(i,j-1) &
578  & + vwnd(ii,j,l)*cosl(ii,j))*wrk3(i,j)) * wrk1(i,j)
579  enddo
580 !--
581  ELSE !South pole point,compute at jm-1
582  jj = jm-1
583  do i=1,im
584  ip1 = ie(i)
585  im1 = iw(i)
586  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
587  & + (vwnd(i,jj-1,l)*cosl(i,jj-1) &
588  & - vwnd(i,j,l)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj)
589 
590  enddo
591  ENDIF
592  else
593  IF(cosl(1,j) >= small) THEN !not a pole point
594  DO i=1,im
595  ip1 = ie(i)
596  im1 = iw(i)
597  ii = i + imb2
598  if (ii > im) ii = ii - im
599  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
600  & - (vwnd(i,j-1,l)*cosl(i,j-1) &
601  & + vwnd(ii,j,l)*cosl(ii,j))*wrk3(i,j)) * wrk1(i,j)
602  enddo
603 !--
604  ELSE !South pole point,compute at jm-1
605  jj = jm-1
606  do i=1,im
607  ip1 = ie(i)
608  im1 = iw(i)
609  div(i,j,l) = ((uwnd(ip1,jj,l)-uwnd(im1,jj,l))*wrk2(i,jj) &
610  & - (vwnd(i,jj-1,l)*cosl(i,jj-1) &
611  & - vwnd(i,j,l)*cosl(i,j))*wrk3(i,jj)) * wrk1(i,jj)
612 
613  enddo
614  ENDIF
615  endif
616  ELSE
617  DO i=1,im
618  ip1 = ie(i)
619  im1 = iw(i)
620  div(i,j,l) = ((uwnd(ip1,j,l)-uwnd(im1,j,l))*wrk2(i,j) &
621  & + (vwnd(i,j-1,l)*cosl(i,j-1) &
622  - vwnd(i,j+1,l)*cosl(i,j+1))*wrk3(i,j)) * wrk1(i,j)
623 !sk06132016
624  if(div(i,j,l)>1.0)print*,'Debug in CALDIV',i,j,uwnd(ip1,j,l),uwnd(im1,j,l), &
625  & wrk2(i,j),vwnd(i,j-1,l),cosl(i,j-1),vwnd(i,j+1,l),cosl(i,j+1), &
626  & wrk3(i,j),wrk1(i,j),div(i,j,l)
627 !--
628  ENDDO
629  ENDIF
630  ENDDO ! end of J loop
631 
632 ! GFS use lon avg as one scaler value for pole point
633  call poleavg(im,jm,jsta,jend,small,cosl(1,jsta),spval,div(1,jsta,l))
634 !sk06142016e
635  if(div(1,jsta,l)>1.0)print*,'Debug in CALDIV',jsta,div(1,jsta,l)
636 ! print*,'Debug in CALDIV',' jsta= ',jsta,DIV(1,jsta,l)
637 
638  enddo ! end of l looop
639 !--
640  deallocate (wrk1, wrk2, wrk3, cosl, iw, ie)
641 
642 
643  END SUBROUTINE caldiv
644 
645  SUBROUTINE calgradps(PS,PSX,PSY)
646 
661  use masks, only: gdlat, gdlon
662  use params_mod, only: dtr, d00, small, erad
663  use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, global, &
664  jsta, jend, im, jm, jsta_m, jend_m
665  use gridspec_mod, only: gridtype
666 
667  implicit none
668 !
669 ! DECLARE VARIABLES.
670 !
671  REAL, dimension(im,jsta_2l:jend_2u), intent(in) :: ps
672  REAL, dimension(im,jsta_2l:jend_2u), intent(inout) :: psx,psy
673 !
674  real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), cosl(:,:)
675  INTEGER, allocatable :: ihe(:),ihw(:), ie(:),iw(:)
676 !
677  integer i,j,ip1,im1,ii,iir,iil,jj,imb2
678 !
679 !***************************************************************************
680 ! START CALGRADPS HERE.
681 !
682 ! LOOP TO COMPUTE ZONAL AND MERIDIONAL GRADIENTS OF PS OR LNPS
683 !
684 !sk06162016 DO J=JSTA_2L,JEND_2U
685 !$omp parallel do private(i,j)
686  DO j=jsta,jend
687  DO i=1,im
688  psx(i,j) = spval
689  psy(i,j) = spval
690 !sk PSX(I,J) = D00
691 !sk PSY(I,J) = D00
692  ENDDO
693  ENDDO
694 
695  CALL exch_f(ps)
696 
697 ! IF (MODELNAME == 'GFS' .or. global) THEN
698  CALL exch(gdlat(1,jsta_2l))
699 
700  allocate (wrk1(im,jsta:jend), wrk2(im,jsta:jend), &
701  & wrk3(im,jsta:jend), cosl(im,jsta_2l:jend_2u))
702  allocate(iw(im),ie(im))
703 
704  imb2 = im/2
705 !$omp parallel do private(i)
706  do i=1,im
707  ie(i) = i+1
708  iw(i) = i-1
709  enddo
710  iw(1) = im
711  ie(im) = 1
712 
713 
714 !$omp parallel do private(i,j,ip1,im1)
715  DO j=jsta,jend
716  do i=1,im
717  ip1 = ie(i)
718  im1 = iw(i)
719  cosl(i,j) = cos(gdlat(i,j)*dtr)
720  if(cosl(i,j) >= small) then
721  wrk1(i,j) = 1.0 / (erad*cosl(i,j))
722  else
723  wrk1(i,j) = 0.
724  end if
725  if(i == im .or. i == 1) then
726  wrk2(i,j) = 1.0 / ((360.+gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
727  else
728  wrk2(i,j) = 1.0 / ((gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
729  end if
730  enddo
731  ENDDO
732 
733  CALL exch(cosl)
734 
735 !$omp parallel do private(i,j,ii)
736  DO j=jsta,jend
737  if (j == 1) then
738  if(gdlat(1,j) > 0.) then ! count from north to south
739  do i=1,im
740  ii = i + imb2
741  if (ii > im) ii = ii - im
742  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j+1)-gdlat(ii,j))*dtr) !1/dphi
743  enddo
744  else ! count from south to north
745  do i=1,im
746  ii = i + imb2
747  if (ii > im) ii = ii - im
748  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j+1)+gdlat(ii,j))*dtr) !1/dphi
749  enddo
750  end if
751  elseif (j == jm) then
752  if(gdlat(1,j) < 0.) then ! count from north to south
753  do i=1,im
754  ii = i + imb2
755  if (ii > im) ii = ii - im
756  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j-1)+gdlat(ii,j))*dtr)
757  enddo
758  else ! count from south to north
759  do i=1,im
760  ii = i + imb2
761  if (ii > im) ii = ii - im
762  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j-1)-gdlat(ii,j))*dtr)
763  enddo
764  end if
765  else
766  do i=1,im
767  wrk3(i,j) = 1.0 / ((gdlat(i,j-1)-gdlat(i,j+1))*dtr) !1/dphi
768  enddo
769  endif
770  ENDDO
771 
772 !$omp parallel do private(i,j,ip1,im1,ii,jj)
773  DO j=jsta,jend
774  IF(j == 1) then ! Near North pole
775  if(gdlat(1,j) > 0.) then ! count from north to south
776  IF(cosl(1,j) >= small) THEN !not a pole point
777  DO i=1,im
778  ip1 = ie(i)
779  im1 = iw(i)
780  ii = i + imb2
781  if (ii > im) ii = ii - im
782  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
783  psy(i,j) = (ps(ii,j)-ps(i,j+1))*wrk3(i,j)/erad
784  enddo
785  ELSE !North pole point, compute at j=2
786  jj = 2
787  DO i=1,im
788  ip1 = ie(i)
789  im1 = iw(i)
790  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
791  psy(i,j) = (ps(i,j)-ps(i,jj+1))*wrk3(i,jj)/erad
792  enddo
793  ENDIF
794  else
795  IF(cosl(1,j) >= small) THEN !not a pole point
796  DO i=1,im
797  ip1 = ie(i)
798  im1 = iw(i)
799  ii = i + imb2
800  if (ii > im) ii = ii - im
801  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
802  psy(i,j) = - (ps(ii,j)-ps(i,j+1))*wrk3(i,j)/erad
803  enddo
804  ELSE !North pole point, compute at j=2
805  jj = 2
806  DO i=1,im
807  ip1 = ie(i)
808  im1 = iw(i)
809  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
810  psy(i,j) = - (ps(i,j)-ps(i,jj+1))*wrk3(i,jj)/erad
811  enddo
812  ENDIF
813  endif
814  ELSE IF(j == jm) THEN ! Near South pole
815  if(gdlat(1,j) < 0.) then ! count from north to south
816  IF(cosl(1,j) >= small) THEN !not a pole point
817  DO i=1,im
818  ip1 = ie(i)
819  im1 = iw(i)
820  ii = i + imb2
821  if (ii > im) ii = ii - im
822  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
823  psy(i,j) = (ps(i,j-1)-ps(ii,j))*wrk3(i,j)/erad
824  enddo
825  ELSE !South pole point,compute at jm-1
826  jj = jm-1
827  DO i=1,im
828  ip1 = ie(i)
829  im1 = iw(i)
830  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
831  psy(i,j) = (ps(i,jj-1)-ps(i,j))*wrk3(i,jj)/erad
832  enddo
833  ENDIF
834  else
835  IF(cosl(1,j) >= small) THEN !not a pole point
836  DO i=1,im
837  ip1 = ie(i)
838  im1 = iw(i)
839  ii = i + imb2
840  if (ii > im) ii = ii - im
841  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
842  psy(i,j) = - (ps(i,j-1)-ps(ii,j))*wrk3(i,j)/erad
843  enddo
844  ELSE !South pole point,compute at jm-1
845  jj = jm-1
846  DO i=1,im
847  ip1 = ie(i)
848  im1 = iw(i)
849  psx(i,j) = (ps(ip1,jj)-ps(im1,jj))*wrk2(i,jj)*wrk1(i,jj)
850  psy(i,j) = - (ps(i,jj-1)-ps(i,j))*wrk3(i,jj)/erad
851  enddo
852  ENDIF
853  endif
854  ELSE
855  DO i=1,im
856  ip1 = ie(i)
857  im1 = iw(i)
858  psx(i,j) = (ps(ip1,j)-ps(im1,j))*wrk2(i,j)*wrk1(i,j)
859  psy(i,j) = (ps(i,j-1)-ps(i,j+1))*wrk3(i,j)/erad
860 !sk06142016A
861  if(psx(i,j)>100.0)print*,'Debug in CALGRADPS: PSX',i,j,ps(ip1,j),ps(im1,j), &
862 ! print*,'Debug in CALGRADPS',i,j,PS(ip1,J),PS(im1,J), &
863  & wrk2(i,j),wrk1(i,j),psx(i,j)
864  if(psy(i,j)>100.0)print*,'Debug in CALGRADPS: PSY',i,j,ps(i,j-1),ps(i,j+1), &
865 ! print*,'Debug in CALGRADPS',i,j,PS(i,J-1),PS(i,J+1), &
866  & wrk3(i,j),erad,psy(i,j)
867 !--
868  ENDDO
869  END IF
870 !
871  ENDDO ! end of J loop
872 
873  deallocate (wrk1, wrk2, wrk3, cosl, iw, ie)
874 
875 ! END IF
876 
877  END SUBROUTINE calgradps
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1
upp_math
dvdxdudy() computes dudy, dvdx, uwnd
Definition: UPP_MATH.f:17
exch_f
subroutine exch_f(a)
Definition: EXCH.f:65