diff --git a/parm/post_avblflds.xml b/parm/post_avblflds.xml
index ac518b3aa..19261b3ca 100755
--- a/parm/post_avblflds.xml
+++ b/parm/post_avblflds.xml
@@ -8479,5 +8479,21 @@
3.0
+
+ 1021
+ VPOT_ON_ISOBARIC_SFC_FROM_WIND_FLD
+ VPOT
+ isobaric_sfc
+ 3.0
+
+
+
+ 1022
+ STRM_ON_ISOBARIC_SFC_FROM_WIND_FLD
+ STRM
+ isobaric_sfc
+ 3.0
+
+
diff --git a/parm/sfs/postcntrl_sfs.xml b/parm/sfs/postcntrl_sfs.xml
index 00dd6acb8..a57aa8599 100644
--- a/parm/sfs/postcntrl_sfs.xml
+++ b/parm/sfs/postcntrl_sfs.xml
@@ -68,6 +68,18 @@
5.0
+
+ VPOT_ON_ISOBARIC_SFC_FROM_WIND_FLD
+ 20000.
+ 3.0
+
+
+
+ STRM_ON_ISOBARIC_SFC_FROM_WIND_FLD
+ 20000.
+ 3.0
+
+
MSLET_ON_MEAN_SEA_LVL
6.0
diff --git a/parm/sfs/postxconfig-NT-sfs.txt b/parm/sfs/postxconfig-NT-sfs.txt
index 9dd920d84..1a74deeb1 100644
--- a/parm/sfs/postxconfig-NT-sfs.txt
+++ b/parm/sfs/postxconfig-NT-sfs.txt
@@ -1,5 +1,5 @@
1
-114
+116
GFSPRS
0
ncep_nco
@@ -352,6 +352,90 @@ isobaric_sfc
?
?
?
+1021
+VPOT_ON_ISOBARIC_SFC_FROM_WIND_FLD
+?
+1
+tmpl4_0
+VPOT
+?
+?
+isobaric_sfc
+0
+?
+1
+20000.
+?
+0
+?
+0
+?
+?
+?
+?
+0
+0.0
+0
+0.0
+?
+0
+0.0
+0
+0.0
+0
+0.0
+0
+0.0
+1
+3.0
+0
+0
+0
+?
+?
+?
+1022
+STRM_ON_ISOBARIC_SFC_FROM_WIND_FLD
+?
+1
+tmpl4_0
+STRM
+?
+?
+isobaric_sfc
+0
+?
+1
+20000.
+?
+0
+?
+0
+?
+?
+?
+?
+0
+0.0
+0
+0.0
+?
+0
+0.0
+0
+0.0
+0
+0.0
+0
+0.0
+1
+3.0
+0
+0
+0
+?
+?
+?
23
MSLET_ON_MEAN_SEA_LVL
?
diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f
index 7f4f62f7b..746ea3541 100644
--- a/sorc/ncep_post.fd/MDL2P.f
+++ b/sorc/ncep_post.fd/MDL2P.f
@@ -39,6 +39,7 @@
!> 2023-08-24 | Y Mao | Add gtg_on option for GTG interpolation
!> 2023-09-12 | J Kenyon | Prevent spurious supercooled rain and cloud water
!> 2024-04-23 | E James | Adding smoke emissions (ebb) from RRFS
+!> 2024-09-23 | K Asmar | Add velocity potential and streamfunction from wind vectors
!>
!> @author T Black W/NP2 @date 1999-09-23
!--------------------------------------------------------------------------------------
@@ -75,7 +76,7 @@ SUBROUTINE MDL2P(iostatusD3D)
IEND_2U, slrutah_on, gtg_on
use rqstfld_mod, only: IGET, LVLS, ID, IAVBLFLD, LVLSXML
use gridspec_mod, only: GRIDTYPE, MAPTYPE, DXVAL
- use upp_physics, only: FPVSNEW, CALRH, CALVOR, CALSLR_ROEBBER, CALSLR_UUTAH
+ use upp_physics, only: FPVSNEW, CALRH, CALVOR, CALSLR_ROEBBER, CALSLR_UUTAH, CALCHIPSI
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
@@ -107,6 +108,7 @@ SUBROUTINE MDL2P(iostatusD3D)
INTEGER, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: NL1X, NL1XF
real, dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LSM) :: TPRS, QPRS, FPRS
real, dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LSM) :: RHPRS
+ real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: CHI, PSI
!
INTEGER K, NSMOOTH
!
@@ -228,6 +230,7 @@ SUBROUTINE MDL2P(iostatusD3D)
(IGET(257) > 0) .OR. (IGET(258) > 0) .OR. &
(IGET(294) > 0) .OR. (IGET(268) > 0) .OR. &
(IGET(331) > 0) .OR. (IGET(326) > 0) .OR. &
+ (IGET(1021) > 0) .OR. (IGET(1022) > 0) .OR. &
! add D3D fields
(IGET(354) > 0) .OR. (IGET(355) > 0) .OR. &
(IGET(356) > 0) .OR. (IGET(357) > 0) .OR. &
@@ -1816,6 +1819,64 @@ SUBROUTINE MDL2P(iostatusD3D)
ENDIF
ENDIF
!
+!*** STREAMFUNCTION (PSI) AND VELOCITY POTENTIAL (CHI)
+!
+ IF ( (IGET(1021) > 0 .or. IGET(1022) > 0) .and. MODELNAME == 'GFS' ) THEN
+ IF (LVLS(LP,IGET(1021)) > 0 .or. LVLS(LP,IGET(1022)) > 0) THEN
+ CALL CALCHIPSI(USL,VSL,CHI,PSI)
+! print *,'me=',me,'EGRID1=',EGRID1(1:10,JSTA)
+!
+!*** CHI
+!
+ IF (LVLS(LP,IGET(1021)) > 0) THEN
+!$omp parallel do private(i,j)
+ DO J=JSTA,JEND
+ DO I=ISTA,IEND
+ GRID1(I,J) = CHI(I,J)
+ ENDDO
+ ENDDO
+ if(grib == 'grib2')then
+ cfld = cfld + 1
+ fld_info(cfld)%ifld=IAVBLFLD(IGET(1021))
+ fld_info(cfld)%lvl=LVLSXML(LP,IGET(1021))
+!$omp parallel do private(i,j,ii,jj)
+ do j=1,jend-jsta+1
+ jj = jsta+j-1
+ do i=1,iend-ista+1
+ ii=ista+i-1
+ datapd(i,j,cfld) = GRID1(ii,jj)
+ enddo
+ enddo
+ endif
+ ENDIF !CHI
+!
+!*** PSI
+!
+ IF (LVLS(LP,IGET(1022)) > 0) THEN
+!$omp parallel do private(i,j)
+ DO J=JSTA,JEND
+ DO I=ISTA,IEND
+ GRID1(I,J) = PSI(I,J)
+ ENDDO
+ ENDDO
+ if(grib == 'grib2')then
+ cfld = cfld + 1
+ fld_info(cfld)%ifld=IAVBLFLD(IGET(1022))
+ fld_info(cfld)%lvl=LVLSXML(LP,IGET(1022))
+!$omp parallel do private(i,j,ii,jj)
+ do j=1,jend-jsta+1
+ jj = jsta+j-1
+ do i=1,iend-ista+1
+ ii=ista+i-1
+ datapd(i,j,cfld) = GRID1(ii,jj)
+ enddo
+ enddo
+ endif
+ ENDIF !PSI
+ ENDIF !LVLS(CHIPSI)
+ ENDIF !CHIPSI
+!
+!
! GEOSTROPHIC STREAMFUNCTION.
IF (IGET(086) > 0) THEN
IF (LVLS(LP,IGET(086)) > 0) THEN
diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f
index b2d3f95d5..297fd76a0 100644
--- a/sorc/ncep_post.fd/UPP_PHYSICS.f
+++ b/sorc/ncep_post.fd/UPP_PHYSICS.f
@@ -26,6 +26,8 @@
!>
!> tvirtual() computes virtual temperature.
!>
+!> calchipsi() computes streamfunction and velocity potential.
+!>
!> ### Program history log:
!> Date | Programmer | Comments
!> -----|------------|---------
@@ -33,6 +35,7 @@
!> 2022-07-11 | Jesse Meng | CALSLR_ROEBBER
!> 2023-02-14 | Jesse Meng | CALSLR_UUTAH
!> 2023-03-22 | Sam Trahan | Fix out-of-bounds access by not calling BOUND
+!> 2024-11-21 | K. Asmar, J. Meng, G. Vandenberghe | CALCHIPSI
!>
!> @author Jesse Meng @date 2020-05-20
module upp_physics
@@ -49,7 +52,7 @@ module upp_physics
public :: CALRH_GFS, CALRH_GSD, CALRH_NAM
public :: CALRH_PW
public :: CALSLR_ROEBBER, CALSLR_UUTAH
- public :: CALVOR
+ public :: CALVOR, CALCHIPSI
public :: FPVSNEW
public :: TVIRTUAL
@@ -4498,5 +4501,472 @@ END SUBROUTINE CALSLR_UUTAH
!
!-------------------------------------------------------------------------------------
!
- end module upp_physics
+!> Computes streamfunction and velocity potential from absolute vorticity
+!> and divergence (computed as in calvor subroutine).
+!>
+!> Applies a poisson solver with 300,000 iterations plus a convergence
+!> condition to exit the loop when the error is below 50.
+!>
+!> @param[in] UWND U-wind (m/s) at mass-points
+!> @param[in] VWND V-wind (m/s) at mass-points
+!> @param[out] CHI velocity potential (m^2/s) at mass-points
+!> @param[out] PSI streamfunction (m^2/s) at mass-points
+!>
+!> ### Program history log:
+!> Date | Programmer | Comments
+!> -----|------------|---------
+!> 2024-10-28 | K. Asmar and J. Meng | Initial
+!> 2024-11-21 | George Vandenberghe | Add convergence condition
+!>
+!> @author(s) K. Asmar, J. Meng, G. Vandenberghe @date 2024-11-21
+ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI)
+!
+ use vrbls2d, only: f
+ use masks, only: gdlat, gdlon, dx, dy
+ use params_mod, only: d00, dtr, small, erad
+ use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, global, &
+ jsta, jend, im, jm, jsta_m, jend_m, gdsdegr,&
+ ista, iend, ista_m, iend_m, ista_2l, iend_2u, me, num_procs
+ use gridspec_mod, only: gridtype, dyval
+ use upp_math, only: DVDXDUDY, DDVDX, DDUDY, UUAVG
+ use mpi
+!
+ implicit none
+!
+! DECLARE VARIABLES.
+!
+ REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(in) :: UWND, VWND
+ REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: ABSV, DIV
+ REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(inout) :: CHI, PSI
+ REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: PTMP, ATMP, DTMP
+ REAL, dimension(IM,2) :: GLATPOLES, COSLPOLES, UPOLES, VPOLES, AVPOLES
+ REAL, dimension(IM,JSTA:JEND) :: COSLTEMP, AVTEMP
+!
+ real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), cosl(:,:)
+ INTEGER, allocatable :: IHE(:),IHW(:), IE(:),IW(:)
+!
+ integer, parameter :: npass2=2, npass3=3
+ integer I,J,ip1,im1,ii,iir,iil,jj,JMT2,imb2, npass, nn, jtem
+ real R2DX,R2DY,DVDX,DUDY,UAVG,TPH1,TPHI, tx1(im+2), tx2(im+2)
+ real rtmp, rerr, err,pval,errmax,errmin,edif
+ integer ier,jjk, mype
+!
+!***************************************************************************
+! START CALCHIPSI HERE.
+!
+! LOOP TO COMPUTE ABSOLUTE VORTICITY FROM WINDS.
+!
+!$omp parallel do private(i,j)
+ DO J=JSTA_2L,JEND_2U
+ DO I=ISTA_2L,IEND_2U
+ ABSV(I,J) = SPVAL
+ DIV(I,J) = SPVAL
+ CHI(I,J) = SPVAL
+ PSI(I,J) = SPVAL
+ ENDDO
+ ENDDO
+!
+ CALL EXCH(UWND)
+ CALL EXCH(VWND)
+!
+ CALL EXCH(GDLAT(ISTA_2L,JSTA_2L))
+ CALL EXCH(GDLON(ISTA_2L,JSTA_2L))
+!
+ allocate (wrk1(ista:iend,jsta:jend), wrk2(ista:iend,jsta:jend), &
+ & wrk3(ista:iend,jsta:jend), cosl(ista_2l:iend_2u,jsta_2l:jend_2u))
+ allocate(iw(im),ie(im))
+ imb2 = im/2
+!$omp parallel do private(i)
+ do i=ista,iend
+ ie(i) = i+1
+ iw(i) = i-1
+ enddo
+! iw(1) = im
+! ie(im) = 1
+!
+! if(1>=jsta .and. 1<=jend)then
+! if(cos(gdlat(1,1)*dtr)= SMALL) then
+ wrk1(i,j) = 1.0 / (ERAD*cosl(i,j))
+ else
+ wrk1(i,j) = 0.
+ end if
+ if(i == im .or. i == 1) then
+ wrk2(i,j) = 1.0 / ((360.+GDLON(ip1,J)-GDLON(im1,J))*DTR) !1/dlam
+ else
+ wrk2(i,j) = 1.0 / ((GDLON(ip1,J)-GDLON(im1,J))*DTR) !1/dlam
+ end if
+ enddo
+ enddo
+! CALL EXCH(cosl(1,JSTA_2L))
+ CALL EXCH(cosl)
+!
+ call fullpole( cosl(ista_2l:iend_2u,jsta_2l:jend_2u),coslpoles)
+ call fullpole(gdlat(ista_2l:iend_2u,jsta_2l:jend_2u),glatpoles)
+!
+!$omp parallel do private(i,j,ii)
+ DO J=JSTA,JEND
+ if (j == 1) then
+ if(gdlat(ista,j) > 0.) then ! count from north to south
+ do i=ista,iend
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ ! wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J+1)-GDLAT(II,J))*DTR) !1/dphi
+ wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J+1)-GLATPOLES(ii,1))*DTR) !1/dphi
+ enddo
+ else ! count from south to north
+ do i=ista,iend
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ ! wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J+1)+GDLAT(II,J))*DTR) !1/dphi
+ wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J+1)+GLATPOLES(ii,1))*DTR) !1/dphi
+!
+ enddo
+ end if
+ elseif (j == JM) then
+ if(gdlat(ista,j) < 0.) then ! count from north to south
+ do i=ista,iend
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ ! wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J-1)+GDLAT(II,J))*DTR)
+ wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J-1)+GLATPOLES(ii,2))*DTR)
+ enddo
+ else ! count from south to north
+ do i=ista,iend
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ ! wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J-1)-GDLAT(II,J))*DTR)
+ wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J-1)-GLATPOLES(ii,2))*DTR)
+ enddo
+ end if
+ else
+ do i=ista,iend
+ wrk3(i,j) = 1.0 / ((GDLAT(I,J-1)-GDLAT(I,J+1))*DTR) !1/dphi
+ enddo
+ endif
+ enddo
+!
+ npass = 0
+!
+ jtem = jm / 18 + 1
+!
+ call fullpole(uwnd(ista_2l:iend_2u,jsta_2l:jend_2u),upoles)
+ call fullpole(vwnd(ista_2l:iend_2u,jsta_2l:jend_2u),vpoles)
+!
+!$omp parallel do private(i,j,ip1,im1,ii,jj,tx1,tx2)
+ DO J=JSTA,JEND
+ IF(J == 1) then ! Near North or South pole
+ if(gdlat(ista,j) > 0.) then ! count from north to south
+ IF(cosl(ista,j) >= SMALL) THEN !not a pole point
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ if(VWND(ip1,J)==SPVAL .or. VWND(im1,J)==SPVAL .or. &
+! UWND(II,J)==SPVAL .or. UWND(I,J+1)==SPVAL) cycle
+ UPOLES(II,1)==SPVAL .or. UWND(I,J+1)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,J)-VWND(im1,J))*wrk2(i,j) &
+ & + (upoles(II,1)*coslpoles(II,1) &
+ & + UWND(I,J+1)*COSL(I,J+1))*wrk3(i,j)) * wrk1(i,j) &
+ & + F(I,J)
+ DIV(I,J) = ((UWND(ip1,J)-UWND(im1,J))*wrk2(i,j) &
+ & - (VPOLES(II,1)*COSLPOLEs(II,1) &
+ & + VWND(I,J+1)*COSL(I,J+1))*wrk3(i,j)) * wrk1(i,j)
+ enddo
+ ELSE !pole point, compute at j=2
+ jj = 2
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ if(VWND(ip1,JJ)==SPVAL .or. VWND(im1,JJ)==SPVAL .or. &
+ UWND(I,J)==SPVAL .or. UWND(I,jj+1)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,JJ)-VWND(im1,JJ))*wrk2(i,jj) &
+ & - (UWND(I,J)*COSL(I,J) &
+ - UWND(I,jj+1)*COSL(I,Jj+1))*wrk3(i,jj)) * wrk1(i,jj) &
+ & + F(I,Jj)
+ DIV(I,J) = ((UWND(ip1,JJ)-UWND(im1,JJ))*wrk2(i,jj) &
+ & + (VWND(I,J)*COSL(I,J) &
+ - VWND(I,jj+1)*COSL(I,jj+1))*wrk3(i,jj)) * wrk1(i,jj)
+ enddo
+ ENDIF
+ else
+ IF(cosl(ista,j) >= SMALL) THEN !not a pole point
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ if(VWND(ip1,J)==SPVAL .or. VWND(im1,J)==SPVAL .or. &
+! UWND(II,J)==SPVAL .or. UWND(I,J+1)==SPVAL) cycle
+ UPOLES(II,1)==SPVAL .or. UWND(I,J+1)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,J)-VWND(im1,J))*wrk2(i,j) &
+ & - (upoles(II,1)*coslpoles(II,1) &
+ & + UWND(I,J+1)*COSL(I,J+1))*wrk3(i,j)) * wrk1(i,j) &
+ & + F(I,J)
+ DIV(I,J) = ((UWND(ip1,J)-UWND(im1,J))*wrk2(i,j) &
+ & + (vpoles(II,1)*coslpoles(II,1) &
+ & + VWND(I,J+1)*COSL(I,J+1))*wrk3(i,j)) * wrk1(i,j)
+ enddo
+ ELSE !pole point, compute at j=2
+ jj = 2
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ if(VWND(ip1,JJ)==SPVAL .or. VWND(im1,JJ)==SPVAL .or. &
+ UWND(I,J)==SPVAL .or. UWND(I,jj+1)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,JJ)-VWND(im1,JJ))*wrk2(i,jj) &
+ & + (UWND(I,J)*COSL(I,J) &
+ - UWND(I,jj+1)*COSL(I,Jj+1))*wrk3(i,jj)) * wrk1(i,jj) &
+ & + F(I,Jj)
+ DIV(I,J) = ((UWND(ip1,JJ)-UWND(im1,JJ))*wrk2(i,jj) &
+ & - (VWND(I,J)*COSL(I,J) &
+ - VWND(I,jj+1)*COSL(I,Jj+1))*wrk3(i,jj)) * wrk1(i,jj)
+ enddo
+ ENDIF
+ endif
+ ELSE IF(J == JM) THEN ! Near North or South Pole
+ if(gdlat(ista,j) < 0.) then ! count from north to south
+ IF(cosl(ista,j) >= SMALL) THEN !not a pole point
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ if(VWND(ip1,J)==SPVAL .or. VWND(im1,J)==SPVAL .or. &
+! UWND(I,J-1)==SPVAL .or. UWND(II,J)==SPVAL) cycle
+ UWND(I,J-1)==SPVAL .or. UPOLES(II,2)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,J)-VWND(im1,J))*wrk2(i,j) &
+ & - (UWND(I,J-1)*COSL(I,J-1) &
+ & + upoles(II,2)*coslpoles(II,2))*wrk3(i,j)) * wrk1(i,j) &
+ & + F(I,J)
+ DIV(I,J) = ((UWND(ip1,J)-UWND(im1,J))*wrk2(i,j) &
+ & + (VWND(I,J-1)*COSL(I,J-1) &
+ & + vpoles(II,2)*coslpoles(II,2))*wrk3(i,j)) * wrk1(i,j)
+ enddo
+ ELSE !pole point,compute at jm-1
+ jj = jm-1
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ if(VWND(ip1,JJ)==SPVAL .or. VWND(im1,JJ)==SPVAL .or. &
+ UWND(I,jj-1)==SPVAL .or. UWND(I,J)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,JJ)-VWND(im1,JJ))*wrk2(i,jj) &
+ & - (UWND(I,jj-1)*COSL(I,Jj-1) &
+ & - UWND(I,J)*COSL(I,J))*wrk3(i,jj)) * wrk1(i,jj) &
+ & + F(I,Jj)
+ DIV(I,J) = ((UWND(ip1,JJ)-UWND(im1,JJ))*wrk2(i,jj) &
+ & + (VWND(I,jj-1)*COSL(I,Jj-1) &
+ & - VWND(I,J)*COSL(I,J))*wrk3(i,jj)) * wrk1(i,jj)
+ enddo
+ ENDIF
+ else
+ IF(cosl(ista,j) >= SMALL) THEN !not a pole point
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ ii = i + imb2
+ if (ii > im) ii = ii - im
+ if(VWND(ip1,J)==SPVAL .or. VWND(im1,J)==SPVAL .or. &
+! UWND(I,J-1)==SPVAL .or. UWND(II,J)==SPVAL) cycle
+ UWND(I,J-1)==SPVAL .or. UPOLES(II,2)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,J)-VWND(im1,J))*wrk2(i,j) &
+ & + (UWND(I,J-1)*COSL(I,J-1) &
+ & + upoles(II,2)*coslpoles(II,2))*wrk3(i,j)) * wrk1(i,j) &
+ & + F(I,J)
+ DIV(I,J) = ((UWND(ip1,J)-UWND(im1,J))*wrk2(i,j) &
+ & - (VWND(I,J-1)*COSL(I,J-1) &
+ & + vpoles(II,2)*coslpoles(II,2))*wrk3(i,j)) * wrk1(i,j)
+ enddo
+ ELSE !pole point,compute at jm-1
+ jj = jm-1
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ if(VWND(ip1,JJ)==SPVAL .or. VWND(im1,JJ)==SPVAL .or. &
+ UWND(I,jj-1)==SPVAL .or. UWND(I,J)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,JJ)-VWND(im1,JJ))*wrk2(i,jj) &
+ & + (UWND(I,jj-1)*COSL(I,Jj-1) &
+ & - UWND(I,J)*COSL(I,J))*wrk3(i,jj)) * wrk1(i,jj) &
+ & + F(I,Jj)
+ DIV(I,J) = ((UWND(ip1,JJ)-UWND(im1,JJ))*wrk2(i,jj) &
+ & - (VWND(I,jj-1)*COSL(I,Jj-1) &
+ & - VWND(I,J)*COSL(I,J))*wrk3(i,jj)) * wrk1(i,jj)
+ enddo
+ ENDIF
+ endif
+ ELSE
+ DO I=ISTA,IEND
+ ip1 = ie(i)
+ im1 = iw(i)
+ if(VWND(ip1,J)==SPVAL .or. VWND(im1,J)==SPVAL .or. &
+ UWND(I,J-1)==SPVAL .or. UWND(I,J+1)==SPVAL) cycle
+ ABSV(I,J) = ((VWND(ip1,J)-VWND(im1,J))*wrk2(i,j) &
+ & - (UWND(I,J-1)*COSL(I,J-1) &
+ - UWND(I,J+1)*COSL(I,J+1))*wrk3(i,j)) * wrk1(i,j) &
+ + F(I,J)
+ DIV(I,J) = ((UWND(ip1,J)-UWND(im1,J))*wrk2(i,j) &
+ & + (VWND(I,J-1)*COSL(I,J-1) &
+ - VWND(I,J+1)*COSL(I,J+1))*wrk3(i,j)) * wrk1(i,j)
+ ENDDO
+ END IF
+ if (npass > 0) then
+ do i=ista,iend
+ tx1(i) = absv(i,j)
+ enddo
+ do nn=1,npass
+ do i=ista,iend
+ tx2(i+1) = tx1(i)
+ enddo
+ tx2(1) = tx2(im+1)
+ tx2(im+2) = tx2(2)
+ do i=2,im+1
+ tx1(i-1) = 0.25 * (tx2(i-1) + tx2(i+1)) + 0.5*tx2(i)
+ enddo
+ enddo
+ do i=ista,iend
+ absv(i,j) = tx1(i)
+ enddo
+ endif
+ END DO ! end of J loop
+
+! deallocate (wrk1, wrk2, wrk3, cosl)
+! GFS use lon avg as one scaler value for pole point
+!
+ ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1,jsta),SPVAL,ABSV(1,jsta))
+!
+ call exch(absv(ista_2l:iend_2u,jsta_2l:jend_2u))
+ call fullpole(absv(ista_2l:iend_2u,jsta_2l:jend_2u),avpoles)
+!
+ cosltemp=spval
+ if(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
+ if(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
+ avtemp=spval
+ if(jsta== 1) avtemp(1:im, 1)=avpoles(1:im,1)
+ if(jend==jm) avtemp(1:im,jm)=avpoles(1:im,2)
+!
+ call poleavg(IM,JM,JSTA,JEND,SMALL,cosltemp(1,jsta),SPVAL,avtemp(1,jsta))
+!
+ if(jsta== 1) absv(ista:iend, 1)=avtemp(ista:iend, 1)
+ if(jend==jm) absv(ista:iend,jm)=avtemp(ista:iend,jm)
+!
+! deallocate (wrk1, wrk11, wrk2, wrk3, cosl, iw, ie)
+!
+ call exch(absv(ista_2l:iend_2u,jsta_2l:jend_2u))
+ call exch(div(ista_2l:iend_2u,jsta_2l:jend_2u))
+!
+! store absv and div factors before poisson loops
+!$omp parallel do private(i,j)
+ DO J=JSTA,JEND
+ DO I=ISTA,IEND
+ ATMP(I,J)=0.25*(ABSV(I,J)-F(I,J))/(wrk2(i,j)*wrk1(i,j)*wrk3(i,j)*wrk1(i,j)*COSL(i,j)*4.)
+ DTMP(I,J)=0.25*DIV(I,J)/(wrk2(i,j)*wrk1(i,j)*wrk3(i,j)*wrk1(i,j)*COSL(i,j)*4.)
+ ENDDO
+ ENDDO
+!
+! poisson solver for psi and chi
+ PSI=0.
+ do jjk=1,1000
+ DO jj=1,300
+ call exch(psi(ista_2l:iend_2u,jsta_2l:jend_2u))
+ PTMP=PSI
+ err=0
+ DO J=JSTA,JEND
+ DO I=ISTA,IEND
+ IF(J>1 .and. J1 .and. J