From 26aaeb0bd8ba8956d9710c6ec808f3ad6e7f9cfd Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA Date: Mon, 21 Oct 2024 20:19:10 +0000 Subject: [PATCH 01/17] add spectral strm and vpot at 200mb for SFS --- parm/post_avblflds.xml | 16 +++ parm/sfs/postcntrl_sfs.xml | 12 ++ parm/sfs/postxconfig-NT-sfs.txt | 86 +++++++++++++- sorc/ncep_post.fd/CALCHIPSI.f | 176 ++++++++++++++++++++++++++++ sorc/ncep_post.fd/CMakeLists.txt | 1 + sorc/ncep_post.fd/GRIDSPEC.f | 1 + sorc/ncep_post.fd/INITPOST_NETCDF.f | 5 +- sorc/ncep_post.fd/MDL2P.f | 85 ++++++++++++++ 8 files changed, 379 insertions(+), 3 deletions(-) create mode 100644 sorc/ncep_post.fd/CALCHIPSI.f 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/CALCHIPSI.f b/sorc/ncep_post.fd/CALCHIPSI.f new file mode 100644 index 000000000..4da03a3d0 --- /dev/null +++ b/sorc/ncep_post.fd/CALCHIPSI.f @@ -0,0 +1,176 @@ +!> @file +!> @brief Subroutine that computes the velocity potential and +!> streamfunction from isobaric winds. +!> +!>
+!> This routine is based on the CFS genpsiandchi program that
+!> computes velocity potential and streamfunction from the
+!> isobaric wind components. The program was authored and provided
+!> by Saha and H. Chuang.
+!> Given the U-V wind components at P-points, this routine
+!> collects the winds in the full IM,JM,LSM domain,
+!> transforms them back to spectrum space and computes divergence,
+!> vorticity, streamfunction and potential. The routine returns: 
+!> 	  PSI: the streamfunction in global domain
+!> 	  CHI: the velocity potential in global domain
+!>
+!> +!> @param[in] UISO real U-wind component (m/s) at all P-points. +!> @param[in] VISO real V-wind component (m/s) at all P-points. +!> @param[out] CHI real velocity potential (m^2/s) in full grid domain at all P-points. +!> @param[out] PSI real streamfunction (m^2/s) in full grid domain at all P-points +!> +!> ### Program history log: +!> Date | Programmer | Comments +!> -----------|--------------|--------- +!> 2024-07-17 | Karina Asmar | Initial +!> 2024-07-25 | Jesse Meng | Add MPI scatterv +!> +!> @author Karina Asmar EMC/VPPPG @date 2024-07-17 +!----------------------------------------------------------------------- +!> @brief Subroutine that computes velocity potential and streamfunction +!> from isobaric winds. +!> +!> @param[in] UISO real U-wind component (m/s) at all P-points. +!> @param[in] VISO real V-wind component (m/s) at all P-points. +!> @param[out] CHI real velocity potential (m^2/s) in full grid domain at P-points. +!> @param[out] PSI real streamfunction (m^2/s) in full grid domain at P-points +!----------------------------------------------------------------------- + SUBROUTINE CALCHIPSI(UISO,VISO,CHI,PSI) +! +! INCLUDE ETA GRID DIMENSIONS. SET/DERIVE OTHER PARAMETERS. +! + use gridspec_mod, only: IDRT + use ctlblk_mod, only: ISTA, IEND, JSTA, JEND, IM, JM, LSM, ME, SPVAL, MPI_COMM_COMP,& + num_procs, icnt, idsp, isxa, iexa, jsxa, jexa + use rqstfld_mod, only: IGET, LVLS +! +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none + + include 'mpif.h' +! +! DECLARE VARIABLES. +! + integer :: JCAP, I, J, L, IERR + REAL, dimension(ISTA:IEND,JSTA:JEND,LSM), intent(in) :: UISO, VISO + REAL, dimension(IM,JM,LSM), intent(out) :: CHI, PSI + + integer k, m + real, allocatable :: CHI1(:),CHISUB(:),PSI1(:),PSISUB(:),COL_UWIND(:,:),COL_VWIND(:,:), & + IN_UWIND(:,:),IN_VWIND(:,:),OUT_UWIND(:,:),OUT_VWIND(:,:), & + DIV(:,:),ZO(:,:),CHI_OUT(:,:),PSI_OUT(:,:) + +! +!*************************************************************************** +! START CALCHIPSI HERE. +! +! SAVE ALL P LEVELS OF U/V WINDS AT GLOBAL GRID + + ALLOCATE(COL_UWIND(IM,JM)) + ALLOCATE(COL_VWIND(IM,JM)) + + ALLOCATE(IN_UWIND(IM,JM)) + ALLOCATE(IN_VWIND(IM,JM)) + ALLOCATE(OUT_UWIND(IM,JM)) + ALLOCATE(OUT_VWIND(IM,JM)) + ALLOCATE(DIV(IM,JM)) + ALLOCATE(ZO(IM,JM)) + ALLOCATE(CHI_OUT(IM,JM)) + ALLOCATE(PSI_OUT(IM,JM)) + + ALLOCATE(CHI1(im*jm)) + ALLOCATE(CHISUB(icnt(me))) + ALLOCATE(PSI1(im*jm)) + ALLOCATE(PSISUB(icnt(me))) + + CHI = SPVAL + PSI = SPVAL + + DO L=1,LSM + IF(LVLS(L,IGET(1021)) > 0)THEN + + CALL COLLECT_ALL(UISO(ISTA:IEND,JSTA:JEND,L),COL_UWIND) + CALL COLLECT_ALL(VISO(ISTA:IEND,JSTA:JEND,L),COL_VWIND) +!$omp parallel do private(i,j) + DO J=1,JM + DO I=1,IM + IN_UWIND(I,J)=COL_UWIND(I,J) + IN_VWIND(I,J)=COL_VWIND(I,J) + ENDDO + ENDDO + + IF (ME==0) THEN + + ! SET MAX WAVELENGTH FOR SPECTRAL TRUNCATION + IF(IDRT == 0)THEN + JCAP = (JM-3)/2 + ELSE + JCAP = JM-1 + ENDIF + + ! COMPUTE CHI/PSI FROM WIND VECTORS IN SPECTRAL SPACE + CALL SPTRUNV(0,JCAP,IDRT,IM, & + JM,IDRT,IM,JM,1, & + 0,0,0,0, & + 0,0,0,0, & + IN_UWIND(1,1),IN_VWIND(1,1), & + .FALSE.,OUT_UWIND(1,1),OUT_VWIND(1,1), & + .FALSE.,DIV,ZO, & + .TRUE.,CHI_OUT(1,1),PSI_OUT(1,1)) + + ENDIF ! END OF ME=0 BLOCK + + CALL MPI_BARRIER(MPI_COMM_COMP, IERR) + + IF (ME==0) THEN + k=0 + DO m=0,num_procs-1 + DO J=jsxa(m),jexa(m) + DO I=isxa(m),iexa(m) + k=k+1 + CHI1(k)=CHI_OUT(I,J) + PSI1(k)=PSI_OUT(I,J) + ENDDO + ENDDO + ENDDO + ENDIF + + CALL MPI_SCATTERV(CHI1,icnt,idsp,MPI_REAL, & + CHISUB,icnt(me),MPI_REAL,0,MPI_COMM_WORLD,IERR) + CALL MPI_SCATTERV(PSI1,icnt,idsp,MPI_REAL, & + PSISUB,icnt(me),MPI_REAL,0,MPI_COMM_WORLD,IERR) + + k=0 + DO J=JSTA,JEND + DO I=ISTA,IEND + k=k+1 + CHI(I,J,L)=CHISUB(k) + PSI(I,J,L)=PSISUB(k) + ENDDO + ENDDO + + ENDIF + ENDDO + + DEALLOCATE(CHI1) + DEALLOCATE(CHISUB) + DEALLOCATE(PSI1) + DEALLOCATE(PSISUB) + + DEALLOCATE(IN_UWIND) + DEALLOCATE(IN_VWIND) + DEALLOCATE(OUT_UWIND) + DEALLOCATE(OUT_VWIND) + DEALLOCATE(DIV) + DEALLOCATE(ZO) + DEALLOCATE(CHI_OUT) + DEALLOCATE(PSI_OUT) + + DEALLOCATE(COL_UWIND) + DEALLOCATE(COL_VWIND) +! +! +! END OF ROUTINE. + RETURN + END diff --git a/sorc/ncep_post.fd/CMakeLists.txt b/sorc/ncep_post.fd/CMakeLists.txt index 5dad96ae1..3fe624c6b 100644 --- a/sorc/ncep_post.fd/CMakeLists.txt +++ b/sorc/ncep_post.fd/CMakeLists.txt @@ -4,6 +4,7 @@ list(APPEND LIB_SRC AVIATION.f BNDLYR.f BOUND.f + CALCHIPSI.f CALDRG.f CALDWP.f CALGUST.f diff --git a/sorc/ncep_post.fd/GRIDSPEC.f b/sorc/ncep_post.fd/GRIDSPEC.f index 3f5936d9f..6fb94f966 100644 --- a/sorc/ncep_post.fd/GRIDSPEC.f +++ b/sorc/ncep_post.fd/GRIDSPEC.f @@ -35,6 +35,7 @@ module GRIDSPEC_mod integer cenlonv !< center longitude of grid integer latlastv !< latitude of last grid point (upper right corner latitude) integer lonlastv !< longitude of last grid point (upper right corner longitude) + integer idrt !< grid identifier real PSMAPF !< map scale factor character(len=1) gridtype !< type of grid staggering as in Arakawa grids (Arakawa-A through Arakawa-E) ! diff --git a/sorc/ncep_post.fd/INITPOST_NETCDF.f b/sorc/ncep_post.fd/INITPOST_NETCDF.f index 6e1a9ecd0..4d9a4885c 100644 --- a/sorc/ncep_post.fd/INITPOST_NETCDF.f +++ b/sorc/ncep_post.fd/INITPOST_NETCDF.f @@ -58,6 +58,7 @@ !> 2024-06-25 | Wen Meng | Add capability to read fhzero as either an integer or float !> 2024-08-26 | Karina Asmar | Add temporal u/v, speed max wind components at 10m agl !> 2024-10-11 | Sam Trahan | Fixed an incorrect array length in read_netcdf_3d_para +!> 2024-10-21 | Karina Asmar | Read in and store idrt in gridspec_mod !> !> @author Hui-Ya Chuang @date 2016-03-04 !---------------------------------------------------------------------- @@ -127,7 +128,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, & dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & latstartv, latlastv,cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r, STANDLON, & - latse,lonse,latnw,lonnw + latse,lonse,latnw,lonnw,idrt use upp_physics, only: fpvsnew !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none @@ -195,7 +196,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) !jw integer ii,jj,js,je,iyear,imn,iday,itmp,ioutcount,istatus, & I,J,L,ll,k,k1,kf,irtn,igdout,n,Index,nframe, & - nframed2,iunitd3d,ierr,idum,iret,nrec,idrt + nframed2,iunitd3d,ierr,idum,iret,nrec integer ncid3d,ncid2d,varid,nhcas,varid_bl,iret_bl real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & tvll,pmll,tv, tx1, tx2, zpbltop diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index 7f4f62f7b..524cc9d1b 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -1,3 +1,4 @@ + !> @file !> @brief mdl2p() computes vertical interpolation of model levels to pressure. !> @@ -39,6 +40,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-07-17 | K Asmar | Add velocity potential and streamfunction from wind vectors !> !> @author T Black W/NP2 @date 1999-09-23 !-------------------------------------------------------------------------------------- @@ -107,6 +109,9 @@ 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:IEND,JSTA:JEND,LSM) :: USLP, VSLP + real, dimension(IM,JM,LSM) :: CHI, PSI ! INTEGER K, NSMOOTH ! @@ -228,6 +233,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. & @@ -1769,6 +1775,19 @@ SUBROUTINE MDL2P(iostatusD3D) endif ENDIF ENDIF +! +! *** K. ASMAR - SAVE ALL P-LEVELS OF U/V WINDS FOR VELOCITY POTENTIAL AND STREAMFUNCTION +! + IF (IGET(1021)>0 .OR. IGET(1022)>0) THEN +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=ISTA,IEND + USLP(I,J,LP)=USL(I,J) + VSLP(I,J,LP)=VSL(I,J) + ENDDO + ENDDO + ENDIF + ! !*** ABSOLUTE VORTICITY ! @@ -4319,6 +4338,72 @@ SUBROUTINE MDL2P(iostatusD3D) endif ENDIF ! +! *** K. ASMAR - COMPUTE VELOCITY POTENTIAL AND STREAMFUNCTION +! + IF (IGET(1021) > 0 .OR. IGET(1022) > 0) THEN + CALL CALCHIPSI(USLP,VSLP,CHI,PSI) + ! SEPARATE VERTICAL LOOP TO STORE VELOCITY POTENTIAL AND STREAMFUNCTION + DO LP=1,LSM + ! VELOCITY POTENTIAL + IF(IGET(1021) > 0) THEN + IF(LVLS(LP,IGET(1021)) > 0)THEN +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=ISTA,IEND + IF(CHI(I,J,LP) < SPVAL) THEN + GRID1(I,J) = CHI(I,J,LP) + ELSE + GRID1(I,J) = SPVAL + ENDIF + 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 + ENDIF + + !STREAMFUNCTION + IF(IGET(1022) > 0) THEN + IF(LVLS(LP,IGET(1022)) > 0)THEN +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=ISTA,IEND + IF(PSI(I,J,LP) < SPVAL) THEN + GRID1(I,J) = PSI(I,J,LP) + ELSE + GRID1(I,J) = SPVAL + ENDIF + 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 + ENDIF + ENDDO ! END OF VERTICAL LOOP LW=1,LSM FOR VPOT AND STRM + ENDIF ! END OF IF BLOCK FOR CALVPOTSTRM +! if(allocated(d3dsl)) deallocate(d3dsl) if(allocated(smokesl)) deallocate(smokesl) if(allocated(fv3dustsl)) deallocate(fv3dustsl) From 21c43db377039591f340d72ad85ebb9152276656 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA Date: Thu, 21 Nov 2024 19:30:57 +0000 Subject: [PATCH 02/17] replace full domain spectral with parallelized numerical computation for strm and vpot --- sorc/ncep_post.fd/CALCHIPSI.f | 176 ---------- sorc/ncep_post.fd/CMakeLists.txt | 1 - sorc/ncep_post.fd/COLLECT_LOC.f | 13 +- sorc/ncep_post.fd/GRIDSPEC.f | 1 - sorc/ncep_post.fd/INITPOST_NETCDF.f | 5 +- sorc/ncep_post.fd/MDL2P.f | 171 +++++----- sorc/ncep_post.fd/UPP_PHYSICS.f | 486 +++++++++++++++++++++++++++- 7 files changed, 593 insertions(+), 260 deletions(-) delete mode 100644 sorc/ncep_post.fd/CALCHIPSI.f diff --git a/sorc/ncep_post.fd/CALCHIPSI.f b/sorc/ncep_post.fd/CALCHIPSI.f deleted file mode 100644 index 4da03a3d0..000000000 --- a/sorc/ncep_post.fd/CALCHIPSI.f +++ /dev/null @@ -1,176 +0,0 @@ -!> @file -!> @brief Subroutine that computes the velocity potential and -!> streamfunction from isobaric winds. -!> -!>
-!> This routine is based on the CFS genpsiandchi program that
-!> computes velocity potential and streamfunction from the
-!> isobaric wind components. The program was authored and provided
-!> by Saha and H. Chuang.
-!> Given the U-V wind components at P-points, this routine
-!> collects the winds in the full IM,JM,LSM domain,
-!> transforms them back to spectrum space and computes divergence,
-!> vorticity, streamfunction and potential. The routine returns: 
-!> 	  PSI: the streamfunction in global domain
-!> 	  CHI: the velocity potential in global domain
-!>
-!> -!> @param[in] UISO real U-wind component (m/s) at all P-points. -!> @param[in] VISO real V-wind component (m/s) at all P-points. -!> @param[out] CHI real velocity potential (m^2/s) in full grid domain at all P-points. -!> @param[out] PSI real streamfunction (m^2/s) in full grid domain at all P-points -!> -!> ### Program history log: -!> Date | Programmer | Comments -!> -----------|--------------|--------- -!> 2024-07-17 | Karina Asmar | Initial -!> 2024-07-25 | Jesse Meng | Add MPI scatterv -!> -!> @author Karina Asmar EMC/VPPPG @date 2024-07-17 -!----------------------------------------------------------------------- -!> @brief Subroutine that computes velocity potential and streamfunction -!> from isobaric winds. -!> -!> @param[in] UISO real U-wind component (m/s) at all P-points. -!> @param[in] VISO real V-wind component (m/s) at all P-points. -!> @param[out] CHI real velocity potential (m^2/s) in full grid domain at P-points. -!> @param[out] PSI real streamfunction (m^2/s) in full grid domain at P-points -!----------------------------------------------------------------------- - SUBROUTINE CALCHIPSI(UISO,VISO,CHI,PSI) -! -! INCLUDE ETA GRID DIMENSIONS. SET/DERIVE OTHER PARAMETERS. -! - use gridspec_mod, only: IDRT - use ctlblk_mod, only: ISTA, IEND, JSTA, JEND, IM, JM, LSM, ME, SPVAL, MPI_COMM_COMP,& - num_procs, icnt, idsp, isxa, iexa, jsxa, jexa - use rqstfld_mod, only: IGET, LVLS -! -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none - - include 'mpif.h' -! -! DECLARE VARIABLES. -! - integer :: JCAP, I, J, L, IERR - REAL, dimension(ISTA:IEND,JSTA:JEND,LSM), intent(in) :: UISO, VISO - REAL, dimension(IM,JM,LSM), intent(out) :: CHI, PSI - - integer k, m - real, allocatable :: CHI1(:),CHISUB(:),PSI1(:),PSISUB(:),COL_UWIND(:,:),COL_VWIND(:,:), & - IN_UWIND(:,:),IN_VWIND(:,:),OUT_UWIND(:,:),OUT_VWIND(:,:), & - DIV(:,:),ZO(:,:),CHI_OUT(:,:),PSI_OUT(:,:) - -! -!*************************************************************************** -! START CALCHIPSI HERE. -! -! SAVE ALL P LEVELS OF U/V WINDS AT GLOBAL GRID - - ALLOCATE(COL_UWIND(IM,JM)) - ALLOCATE(COL_VWIND(IM,JM)) - - ALLOCATE(IN_UWIND(IM,JM)) - ALLOCATE(IN_VWIND(IM,JM)) - ALLOCATE(OUT_UWIND(IM,JM)) - ALLOCATE(OUT_VWIND(IM,JM)) - ALLOCATE(DIV(IM,JM)) - ALLOCATE(ZO(IM,JM)) - ALLOCATE(CHI_OUT(IM,JM)) - ALLOCATE(PSI_OUT(IM,JM)) - - ALLOCATE(CHI1(im*jm)) - ALLOCATE(CHISUB(icnt(me))) - ALLOCATE(PSI1(im*jm)) - ALLOCATE(PSISUB(icnt(me))) - - CHI = SPVAL - PSI = SPVAL - - DO L=1,LSM - IF(LVLS(L,IGET(1021)) > 0)THEN - - CALL COLLECT_ALL(UISO(ISTA:IEND,JSTA:JEND,L),COL_UWIND) - CALL COLLECT_ALL(VISO(ISTA:IEND,JSTA:JEND,L),COL_VWIND) -!$omp parallel do private(i,j) - DO J=1,JM - DO I=1,IM - IN_UWIND(I,J)=COL_UWIND(I,J) - IN_VWIND(I,J)=COL_VWIND(I,J) - ENDDO - ENDDO - - IF (ME==0) THEN - - ! SET MAX WAVELENGTH FOR SPECTRAL TRUNCATION - IF(IDRT == 0)THEN - JCAP = (JM-3)/2 - ELSE - JCAP = JM-1 - ENDIF - - ! COMPUTE CHI/PSI FROM WIND VECTORS IN SPECTRAL SPACE - CALL SPTRUNV(0,JCAP,IDRT,IM, & - JM,IDRT,IM,JM,1, & - 0,0,0,0, & - 0,0,0,0, & - IN_UWIND(1,1),IN_VWIND(1,1), & - .FALSE.,OUT_UWIND(1,1),OUT_VWIND(1,1), & - .FALSE.,DIV,ZO, & - .TRUE.,CHI_OUT(1,1),PSI_OUT(1,1)) - - ENDIF ! END OF ME=0 BLOCK - - CALL MPI_BARRIER(MPI_COMM_COMP, IERR) - - IF (ME==0) THEN - k=0 - DO m=0,num_procs-1 - DO J=jsxa(m),jexa(m) - DO I=isxa(m),iexa(m) - k=k+1 - CHI1(k)=CHI_OUT(I,J) - PSI1(k)=PSI_OUT(I,J) - ENDDO - ENDDO - ENDDO - ENDIF - - CALL MPI_SCATTERV(CHI1,icnt,idsp,MPI_REAL, & - CHISUB,icnt(me),MPI_REAL,0,MPI_COMM_WORLD,IERR) - CALL MPI_SCATTERV(PSI1,icnt,idsp,MPI_REAL, & - PSISUB,icnt(me),MPI_REAL,0,MPI_COMM_WORLD,IERR) - - k=0 - DO J=JSTA,JEND - DO I=ISTA,IEND - k=k+1 - CHI(I,J,L)=CHISUB(k) - PSI(I,J,L)=PSISUB(k) - ENDDO - ENDDO - - ENDIF - ENDDO - - DEALLOCATE(CHI1) - DEALLOCATE(CHISUB) - DEALLOCATE(PSI1) - DEALLOCATE(PSISUB) - - DEALLOCATE(IN_UWIND) - DEALLOCATE(IN_VWIND) - DEALLOCATE(OUT_UWIND) - DEALLOCATE(OUT_VWIND) - DEALLOCATE(DIV) - DEALLOCATE(ZO) - DEALLOCATE(CHI_OUT) - DEALLOCATE(PSI_OUT) - - DEALLOCATE(COL_UWIND) - DEALLOCATE(COL_VWIND) -! -! -! END OF ROUTINE. - RETURN - END diff --git a/sorc/ncep_post.fd/CMakeLists.txt b/sorc/ncep_post.fd/CMakeLists.txt index bad8c172c..42bdece4e 100644 --- a/sorc/ncep_post.fd/CMakeLists.txt +++ b/sorc/ncep_post.fd/CMakeLists.txt @@ -4,7 +4,6 @@ list(APPEND LIB_SRC AVIATION.f BNDLYR.f BOUND.f - CALCHIPSI.f CALDRG.f CALDWP.f CALGUST.f diff --git a/sorc/ncep_post.fd/COLLECT_LOC.f b/sorc/ncep_post.fd/COLLECT_LOC.f index 66cf7b18e..63823ed33 100644 --- a/sorc/ncep_post.fd/COLLECT_LOC.f +++ b/sorc/ncep_post.fd/COLLECT_LOC.f @@ -10,7 +10,8 @@ !> Date | Programmer | Comments !> -----------|---------------------|---------- !> 2000-01-06 | Jim Tuccillo | Initial -!> 2021-06-01 | George Vandenberghe | 2D Decomposition +!> 2021-06-01 | George Vandenberghe | 2D Decomposition +!> 2024-11-19 | George Vandenberghe | Add timers !> !> @author Jim Tuccillo IBM @date 2000-01-06 !-------------------------------------------------------------------------------- @@ -33,10 +34,12 @@ SUBROUTINE COLLECT_LOC ( A, B ) real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(in) :: a real, dimension(im,jm), intent(out) :: b integer ierr,n + real*8 ta,tb,tc,td,te real, allocatable :: rbufs(:) allocate(buff(im*jm)) jj=( jexa(me)-jsxa(me)+1) * (iexa(me)-isxa(me)+1) allocate( rbufs(( jexa(me)-jsxa(me)+1) * (iexa(me)-isxa(me)+1)) ) + ta=mpi_wtime() ! if ( num_procs <= 1 ) then b = a @@ -81,6 +84,9 @@ SUBROUTINE COLLECT_LOC ( A, B ) deallocate(buff) deallocate(rbufs) + tb=mpi_wtime() + if(me .eq. 0) print 109,' GWVX COLLECT TIME ',im,jm,tb-ta + 109 format(a,2i10,f20.10) end ! !----------------------------------------------------------------------- @@ -104,6 +110,8 @@ SUBROUTINE COLLECT_ALL ( A, B ) real, dimension(im,jm), intent(out) :: b integer ierr,n real, allocatable :: rbufs(:) + real*8 tb,ta + ta=mpi_wtime() allocate(buff(im*jm)) jj=( jexa(me)-jsxa(me)+1) * (iexa(me)-isxa(me)+1) allocate( rbufs(( jexa(me)-jsxa(me)+1) * (iexa(me)-isxa(me)+1)) ) @@ -146,6 +154,9 @@ SUBROUTINE COLLECT_ALL ( A, B ) deallocate(buff) deallocate(rbufs) + tb=mpi_wtime() + if(me .eq. 0) print 109,' GWVX COLLECT_ALL',tb-ta + 109 format(a,f20.10) end diff --git a/sorc/ncep_post.fd/GRIDSPEC.f b/sorc/ncep_post.fd/GRIDSPEC.f index 6fb94f966..3f5936d9f 100644 --- a/sorc/ncep_post.fd/GRIDSPEC.f +++ b/sorc/ncep_post.fd/GRIDSPEC.f @@ -35,7 +35,6 @@ module GRIDSPEC_mod integer cenlonv !< center longitude of grid integer latlastv !< latitude of last grid point (upper right corner latitude) integer lonlastv !< longitude of last grid point (upper right corner longitude) - integer idrt !< grid identifier real PSMAPF !< map scale factor character(len=1) gridtype !< type of grid staggering as in Arakawa grids (Arakawa-A through Arakawa-E) ! diff --git a/sorc/ncep_post.fd/INITPOST_NETCDF.f b/sorc/ncep_post.fd/INITPOST_NETCDF.f index 42f142e45..f52ca19de 100644 --- a/sorc/ncep_post.fd/INITPOST_NETCDF.f +++ b/sorc/ncep_post.fd/INITPOST_NETCDF.f @@ -58,7 +58,6 @@ !> 2024-06-25 | Wen Meng | Add capability to read fhzero as either an integer or float !> 2024-08-26 | Karina Asmar | Add temporal u/v, speed max wind components at 10m agl !> 2024-10-11 | Sam Trahan | Fixed an incorrect array length in read_netcdf_3d_para -!> 2024-10-21 | Karina Asmar | Read in and store idrt in gridspec_mod !> !> @author Hui-Ya Chuang @date 2016-03-04 !---------------------------------------------------------------------- @@ -128,7 +127,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, & dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & latstartv, latlastv,cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r, STANDLON, & - latse,lonse,latnw,lonnw,idrt + latse,lonse,latnw,lonnw use upp_physics, only: fpvsnew !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none @@ -196,7 +195,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) !jw integer ii,jj,js,je,iyear,imn,iday,itmp,ioutcount,istatus, & I,J,L,ll,k,k1,kf,irtn,igdout,n,Index,nframe, & - nframed2,iunitd3d,ierr,idum,iret,nrec + nframed2,iunitd3d,ierr,idum,iret,nrec,idrt integer ncid3d,ncid2d,varid,nhcas,varid_bl,iret_bl real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & tvll,pmll,tv, tx1, tx2, zpbltop diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index 524cc9d1b..707acbf6f 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -1,4 +1,3 @@ - !> @file !> @brief mdl2p() computes vertical interpolation of model levels to pressure. !> @@ -40,7 +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-07-17 | K Asmar | Add velocity potential and streamfunction from wind vectors +!> 2024-09-23 | K Asmar | Add velocity potential and streamfunction from wind vectors !> !> @author T Black W/NP2 @date 1999-09-23 !-------------------------------------------------------------------------------------- @@ -77,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 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -109,9 +108,8 @@ 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:IEND,JSTA:JEND,LSM) :: USLP, VSLP - real, dimension(IM,JM,LSM) :: CHI, PSI + real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: CHI, PSI ! INTEGER K, NSMOOTH ! @@ -1779,7 +1777,6 @@ SUBROUTINE MDL2P(iostatusD3D) ! *** K. ASMAR - SAVE ALL P-LEVELS OF U/V WINDS FOR VELOCITY POTENTIAL AND STREAMFUNCTION ! IF (IGET(1021)>0 .OR. IGET(1022)>0) THEN -!$omp parallel do private(i,j) DO J=JSTA,JEND DO I=ISTA,IEND USLP(I,J,LP)=USL(I,J) @@ -1787,7 +1784,6 @@ SUBROUTINE MDL2P(iostatusD3D) ENDDO ENDDO ENDIF - ! !*** ABSOLUTE VORTICITY ! @@ -1835,6 +1831,101 @@ SUBROUTINE MDL2P(iostatusD3D) ENDIF ENDIF ! +!*** CHIPSI +! + IF (IGET(1021) > 0 .or. IGET(1022) > 0) 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 (SMFLAG .or. ioform == 'binarympiio' ) THEN + call AllGETHERV(GRID1) + if (ioform == 'binarympiio') then +! nsmooth = max(2, min(30,nint(jm/94.0))) +! do k=1,5 + CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) + CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) +! enddo + else + NSMOOTH = nint(4.*(13500./dxm)) +! endif + do k=1,NSMOOTH + CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) + end do + endif + ENDIF + + 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 (SMFLAG .or. ioform == 'binarympiio' ) THEN + call AllGETHERV(GRID1) + if (ioform == 'binarympiio') then +! nsmooth = max(2, min(30,nint(jm/94.0))) +! do k=1,5 + CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) + CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) +! enddo + else + NSMOOTH = nint(4.*(13500./dxm)) +! endif + do k=1,NSMOOTH + CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) + end do + endif + ENDIF + + 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 @@ -4338,72 +4429,6 @@ SUBROUTINE MDL2P(iostatusD3D) endif ENDIF ! -! *** K. ASMAR - COMPUTE VELOCITY POTENTIAL AND STREAMFUNCTION -! - IF (IGET(1021) > 0 .OR. IGET(1022) > 0) THEN - CALL CALCHIPSI(USLP,VSLP,CHI,PSI) - ! SEPARATE VERTICAL LOOP TO STORE VELOCITY POTENTIAL AND STREAMFUNCTION - DO LP=1,LSM - ! VELOCITY POTENTIAL - IF(IGET(1021) > 0) THEN - IF(LVLS(LP,IGET(1021)) > 0)THEN -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=ISTA,IEND - IF(CHI(I,J,LP) < SPVAL) THEN - GRID1(I,J) = CHI(I,J,LP) - ELSE - GRID1(I,J) = SPVAL - ENDIF - 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 - ENDIF - - !STREAMFUNCTION - IF(IGET(1022) > 0) THEN - IF(LVLS(LP,IGET(1022)) > 0)THEN -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=ISTA,IEND - IF(PSI(I,J,LP) < SPVAL) THEN - GRID1(I,J) = PSI(I,J,LP) - ELSE - GRID1(I,J) = SPVAL - ENDIF - 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 - ENDIF - ENDDO ! END OF VERTICAL LOOP LW=1,LSM FOR VPOT AND STRM - ENDIF ! END OF IF BLOCK FOR CALVPOTSTRM -! if(allocated(d3dsl)) deallocate(d3dsl) if(allocated(smokesl)) deallocate(smokesl) if(allocated(fv3dustsl)) deallocate(fv3dustsl) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index b2d3f95d5..97ceafc98 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 | George Vandenbergue | 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 @@ -1737,7 +1740,6 @@ end function TVIRTUAL !> 2019-10-17 | Y Mao | Skip calculation when U/V is SPVAL !> 2020-11-06 | J Meng | Use UPP_MATH Module !> 2022-05-26 | H Chuang | Use GSL approach for FV3R -!> 2024-10-16 | J Kenyon | Initialize ABSV as SPVAL for MPAS applications !> !> @author Russ Treadon W/NP2 @date 1992-12-22 @@ -1747,7 +1749,7 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) 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, submodelname, global, & + 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 @@ -1774,7 +1776,7 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) ! ! LOOP TO COMPUTE ABSOLUTE VORTICITY FROM WINDS. ! - IF(MODELNAME == 'RAPR' .AND. SUBMODELNAME /= 'MPAS') then ! for RAP / HRRR only + IF(MODELNAME == 'RAPR') then !$omp parallel do private(i,j) DO J=JSTA_2L,JEND_2U DO I=ISTA_2L,IEND_2U @@ -2063,7 +2065,7 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) if(jsta== 1) absv(ista:iend, 1)=avtemp(ista:iend, 1) if(jend==jm) absv(ista:iend,jm)=avtemp(ista:iend,jm) - deallocate (wrk1, wrk2, wrk3, cosl, iw, ie) +! deallocate (wrk1, wrk2, wrk3, cosl, iw, ie) ELSE !(MODELNAME == 'GFS' .or. global) @@ -2139,6 +2141,8 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) END DO END IF END IF + + deallocate (wrk1, wrk2, wrk3, cosl, iw, ie) ! ! END OF ROUTINE. ! @@ -4497,6 +4501,478 @@ SUBROUTINE CALSLR_UUTAH(SLR) END SUBROUTINE CALSLR_UUTAH ! !------------------------------------------------------------------------------------- +! +!> Computes streamfunction and velocity potential from absolute vorticity +!> and divergence (computed as in CALVOR subroutine). +!> +!> Applies a Poisson Solver with 300,000 iterations (developed by J. Meng) +!> plus a convergence condition to exit the loop when the error is below 50. +!> +!> @param[in] UWND U-wind (m/s) at P-points +!> @param[in] VWND V-wind (m/s) at P-points +!> @param[out] CHI velocity potential (m^2/s) at P-points +!> @param[out] PSI streamfunction (m^2/s) at P-points +!> +!> ### Program history log: +!> Date | Programmer | Comments +!> -----|------------|--------- +!> 2024-11-21 | George Vandenberghe | Initial +!> +!> @author(s) George 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 + 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 + real*8 ta,tb,tc,td,de,tf + 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 +! npass = npass2 +! if (j > jm-jtem+1 .or. j < jtem) npass = npass3 + 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)) + + PSI=0. + ta=mpi_wtime() + 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 Date: Fri, 22 Nov 2024 09:56:27 -0500 Subject: [PATCH 03/17] add modelname==GFS condition and update author documentation --- sorc/ncep_post.fd/UPP_PHYSICS.f | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index 97ceafc98..b5cf29591 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -4505,8 +4505,8 @@ END SUBROUTINE CALSLR_UUTAH !> Computes streamfunction and velocity potential from absolute vorticity !> and divergence (computed as in CALVOR subroutine). !> -!> Applies a Poisson Solver with 300,000 iterations (developed by J. Meng) -!> plus a convergence condition to exit the loop when the error is below 50. +!> 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 P-points !> @param[in] VWND V-wind (m/s) at P-points @@ -4516,9 +4516,10 @@ END SUBROUTINE CALSLR_UUTAH !> ### Program history log: !> Date | Programmer | Comments !> -----|------------|--------- -!> 2024-11-21 | George Vandenberghe | Initial +!> 2024-10-28 | K. Asmar and J. Meng | Initial +!> 2024-11-21 | George Vandenberghe | Add convergence condition !> -!> @author(s) George Vandenberghe @date 2024-11-21 +!> @author(s) K. Asmar, J. Meng, G. Vandenberghe @date 2024-11-21 SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) ! @@ -4571,6 +4572,8 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) CALL EXCH(UWND) CALL EXCH(VWND) ! + + IF (MODELNAME == 'GFS' .or. global) THEN CALL EXCH(GDLAT(ISTA_2L,JSTA_2L)) CALL EXCH(GDLON(ISTA_2L,JSTA_2L)) From 31cc868ff473942c31192136076704d49d26f6e2 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Fri, 22 Nov 2024 09:58:07 -0500 Subject: [PATCH 04/17] add authors calchipsi --- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index b5cf29591..d07e8f8d9 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -35,7 +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 | George Vandenbergue | CALCHIPSI +!> 2024-11-21 | K. Asmar, J. Meng, G. Vandenberghe | CALCHIPSI !> !> @author Jesse Meng @date 2020-05-20 module upp_physics From a5d8a777f1b8744e424be9a5acf47033896a24ba Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Fri, 22 Nov 2024 10:06:11 -0500 Subject: [PATCH 05/17] add endif for GFS model block --- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index d07e8f8d9..6ee6e27f0 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -4869,6 +4869,8 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) ! deallocate (wrk1, wrk11, wrk2, wrk3, cosl, iw, ie) + ENDIF ! END of MODELNAME=='GFS' BLOCK + call exch(absv(ista_2l:iend_2u,jsta_2l:jend_2u)) call exch(div(ista_2l:iend_2u,jsta_2l:jend_2u)) From 5134baee629bb9c44e6ada7fbd0b4f42195a1c86 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 06:14:19 -0500 Subject: [PATCH 06/17] remove uslp/vslp --- sorc/ncep_post.fd/MDL2P.f | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index 707acbf6f..512c07dd0 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -108,7 +108,6 @@ 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:IEND,JSTA:JEND,LSM) :: USLP, VSLP real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: CHI, PSI ! INTEGER K, NSMOOTH @@ -1773,17 +1772,6 @@ SUBROUTINE MDL2P(iostatusD3D) endif ENDIF ENDIF -! -! *** K. ASMAR - SAVE ALL P-LEVELS OF U/V WINDS FOR VELOCITY POTENTIAL AND STREAMFUNCTION -! - IF (IGET(1021)>0 .OR. IGET(1022)>0) THEN - DO J=JSTA,JEND - DO I=ISTA,IEND - USLP(I,J,LP)=USL(I,J) - VSLP(I,J,LP)=VSL(I,J) - ENDDO - ENDDO - ENDIF ! !*** ABSOLUTE VORTICITY ! From 98d0ac122c283f3133938746bc89cd41fb3a72bd Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 06:17:29 -0500 Subject: [PATCH 07/17] add j kenyon log --- sorc/ncep_post.fd/UPP_PHYSICS.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index 6ee6e27f0..e0cb4f194 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -1740,6 +1740,7 @@ end function TVIRTUAL !> 2019-10-17 | Y Mao | Skip calculation when U/V is SPVAL !> 2020-11-06 | J Meng | Use UPP_MATH Module !> 2022-05-26 | H Chuang | Use GSL approach for FV3R +!> 2024-10-16 | J Kenyon | Initialize ABSV as SPVAL for MPAS applications !> !> @author Russ Treadon W/NP2 @date 1992-12-22 @@ -1776,7 +1777,7 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) ! ! LOOP TO COMPUTE ABSOLUTE VORTICITY FROM WINDS. ! - IF(MODELNAME == 'RAPR') then + IF(MODELNAME == 'RAPR' .AND. SUBMODELNAME /= 'MPAS') then ! for RAP / HRRR only !$omp parallel do private(i,j) DO J=JSTA_2L,JEND_2U DO I=ISTA_2L,IEND_2U From fd8abf279f7eef33c03e05563c74e3fdb6c98881 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 06:23:07 -0500 Subject: [PATCH 08/17] fix: add submodelname to calvor parameters --- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index e0cb4f194..c76c62b30 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -1750,7 +1750,7 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) 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, & + use ctlblk_mod, only: jsta_2l, jend_2u, spval, modelname, submodelname, 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 From 4627df4a1cdb05a6aafe3730688370735b099476 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 06:25:19 -0500 Subject: [PATCH 09/17] remove deallocate from calvor --- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 -- 1 file changed, 2 deletions(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index c76c62b30..73e8d64ca 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -2142,8 +2142,6 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) END DO END IF END IF - - deallocate (wrk1, wrk2, wrk3, cosl, iw, ie) ! ! END OF ROUTINE. ! From f88bf42f14691f6a855380c86d82ea033f445bcc Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 06:28:02 -0500 Subject: [PATCH 10/17] uncomment deallocate in calvor --- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index 73e8d64ca..dfd6af188 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -2066,7 +2066,7 @@ SUBROUTINE CALVOR(UWND,VWND,ABSV) if(jsta== 1) absv(ista:iend, 1)=avtemp(ista:iend, 1) if(jend==jm) absv(ista:iend,jm)=avtemp(ista:iend,jm) -! deallocate (wrk1, wrk2, wrk3, cosl, iw, ie) + deallocate (wrk1, wrk2, wrk3, cosl, iw, ie) ELSE !(MODELNAME == 'GFS' .or. global) From 3a2370488a8998eb6ac7369f96911a092871aeb3 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 09:28:07 -0500 Subject: [PATCH 11/17] comment out smoothing of strm/vpot --- sorc/ncep_post.fd/MDL2P.f | 64 +++++++++++++++++++-------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index 512c07dd0..d10036722 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -1836,22 +1836,22 @@ SUBROUTINE MDL2P(iostatusD3D) ENDDO ENDDO - IF (SMFLAG .or. ioform == 'binarympiio' ) THEN - call AllGETHERV(GRID1) - if (ioform == 'binarympiio') then -! nsmooth = max(2, min(30,nint(jm/94.0))) -! do k=1,5 - CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) - CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) -! enddo - else - NSMOOTH = nint(4.*(13500./dxm)) -! endif - do k=1,NSMOOTH - CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) - end do - endif - ENDIF +! IF (SMFLAG .or. ioform == 'binarympiio' ) THEN +! call AllGETHERV(GRID1) +! if (ioform == 'binarympiio') then +!! nsmooth = max(2, min(30,nint(jm/94.0))) +!! do k=1,5 +! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) +! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) +!! enddo +! else +! NSMOOTH = nint(4.*(13500./dxm)) +!! endif +! do k=1,NSMOOTH +! CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) +! end do +! endif +! ENDIF if(grib == 'grib2')then cfld = cfld + 1 @@ -1878,22 +1878,22 @@ SUBROUTINE MDL2P(iostatusD3D) ENDDO ENDDO - IF (SMFLAG .or. ioform == 'binarympiio' ) THEN - call AllGETHERV(GRID1) - if (ioform == 'binarympiio') then -! nsmooth = max(2, min(30,nint(jm/94.0))) -! do k=1,5 - CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) - CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) -! enddo - else - NSMOOTH = nint(4.*(13500./dxm)) -! endif - do k=1,NSMOOTH - CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) - end do - endif - ENDIF +! IF (SMFLAG .or. ioform == 'binarympiio' ) THEN +! call AllGETHERV(GRID1) +! if (ioform == 'binarympiio') then +!! nsmooth = max(2, min(30,nint(jm/94.0))) +!! do k=1,5 +! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) +! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) +!! enddo +! else +! NSMOOTH = nint(4.*(13500./dxm)) +!! endif +! do k=1,NSMOOTH +! CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) +! end do +! endif +! ENDIF if(grib == 'grib2')then cfld = cfld + 1 From 36631297fa4c2a14e0e128c9def94556bff40e56 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 11:45:41 -0500 Subject: [PATCH 12/17] add modelname=GFS for 1021/1022 variables --- sorc/ncep_post.fd/MDL2P.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index d10036722..fb2e4a597 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -1821,7 +1821,7 @@ SUBROUTINE MDL2P(iostatusD3D) ! !*** CHIPSI ! - IF (IGET(1021) > 0 .or. IGET(1022) > 0) THEN +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) From 740cba92f878c2f32e11be0414f44cefef433068 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA Date: Wed, 4 Dec 2024 16:53:59 +0000 Subject: [PATCH 13/17] fix space typo in . and. --- sorc/ncep_post.fd/MDL2P.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index fb2e4a597..f57c505e1 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -1821,7 +1821,7 @@ SUBROUTINE MDL2P(iostatusD3D) ! !*** CHIPSI ! -IF ( (IGET(1021) > 0 .or. IGET(1022) > 0). and. MODELNAME == 'GFS' ) THEN +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) From 49bb1a2989f75211d6e03d23813160c8597bf511 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA Date: Wed, 4 Dec 2024 16:59:14 +0000 Subject: [PATCH 14/17] remove smoothing function for strm/vpot --- sorc/ncep_post.fd/MDL2P.f | 36 ------------------------------------ 1 file changed, 36 deletions(-) diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index f57c505e1..ce48547a4 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -1835,24 +1835,6 @@ SUBROUTINE MDL2P(iostatusD3D) GRID1(I,J) = CHI(I,J) ENDDO ENDDO - -! IF (SMFLAG .or. ioform == 'binarympiio' ) THEN -! call AllGETHERV(GRID1) -! if (ioform == 'binarympiio') then -!! nsmooth = max(2, min(30,nint(jm/94.0))) -!! do k=1,5 -! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) -! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) -!! enddo -! else -! NSMOOTH = nint(4.*(13500./dxm)) -!! endif -! do k=1,NSMOOTH -! CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) -! end do -! endif -! ENDIF - if(grib == 'grib2')then cfld = cfld + 1 fld_info(cfld)%ifld=IAVBLFLD(IGET(1021)) @@ -1877,24 +1859,6 @@ SUBROUTINE MDL2P(iostatusD3D) GRID1(I,J) = PSI(I,J) ENDDO ENDDO - -! IF (SMFLAG .or. ioform == 'binarympiio' ) THEN -! call AllGETHERV(GRID1) -! if (ioform == 'binarympiio') then -!! nsmooth = max(2, min(30,nint(jm/94.0))) -!! do k=1,5 -! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,0.5) -! CALL SMOOTHC(GRID1,SDUMMY,IM,JM,-0.5) -!! enddo -! else -! NSMOOTH = nint(4.*(13500./dxm)) -!! endif -! do k=1,NSMOOTH -! CALL SMOOTH(GRID1,SDUMMY,IM,JM,0.5) -! end do -! endif -! ENDIF - if(grib == 'grib2')then cfld = cfld + 1 fld_info(cfld)%ifld=IAVBLFLD(IGET(1022)) From c62f225929b59c40505063ec733cf5ae026eb9fe Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA Date: Wed, 4 Dec 2024 17:03:45 +0000 Subject: [PATCH 15/17] remove modelname==GFS if block from calchipsi --- sorc/ncep_post.fd/UPP_PHYSICS.f | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index dfd6af188..f5761c044 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -4571,16 +4571,14 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) CALL EXCH(UWND) CALL EXCH(VWND) ! + CALL EXCH(GDLAT(ISTA_2L,JSTA_2L)) + CALL EXCH(GDLON(ISTA_2L,JSTA_2L)) - IF (MODELNAME == 'GFS' .or. global) THEN - 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), & + 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)) + allocate(iw(im),ie(im)) - imb2 = im/2 + imb2 = im/2 !$omp parallel do private(i) do i=ista,iend ie(i) = i+1 @@ -4867,9 +4865,7 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) if(jend==jm) absv(ista:iend,jm)=avtemp(ista:iend,jm) ! deallocate (wrk1, wrk11, wrk2, wrk3, cosl, iw, ie) - - ENDIF ! END of MODELNAME=='GFS' BLOCK - +! call exch(absv(ista_2l:iend_2u,jsta_2l:jend_2u)) call exch(div(ista_2l:iend_2u,jsta_2l:jend_2u)) From 6a8aa871de58546c2f5af09c93d62a0d2a3587e1 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 12:06:12 -0500 Subject: [PATCH 16/17] remove commented lines for npass and j loop --- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 -- 1 file changed, 2 deletions(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index f5761c044..07a45ba8a 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -4667,8 +4667,6 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) !$omp parallel do private(i,j,ip1,im1,ii,jj,tx1,tx2) DO J=JSTA,JEND -! npass = npass2 -! if (j > jm-jtem+1 .or. j < jtem) npass = npass3 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 From d4889ebb2bdcddd350649ac975f94ac5005f50c4 Mon Sep 17 00:00:00 2001 From: KarinaAsmar-NOAA <148993962+KarinaAsmar-NOAA@users.noreply.github.com> Date: Wed, 4 Dec 2024 12:47:41 -0500 Subject: [PATCH 17/17] fix comments, indentations for standardized format --- sorc/ncep_post.fd/UPP_PHYSICS.f | 664 ++++++++++++++++---------------- 1 file changed, 331 insertions(+), 333 deletions(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index 07a45ba8a..54ccb10e3 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -4502,15 +4502,15 @@ END SUBROUTINE CALSLR_UUTAH !------------------------------------------------------------------------------------- ! !> Computes streamfunction and velocity potential from absolute vorticity -!> and divergence (computed as in CALVOR subroutine). +!> and divergence (computed as in calvor subroutine). !> -!> Applies a Poisson Solver with 300,000 iterations plus a convergence +!> 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 P-points -!> @param[in] VWND V-wind (m/s) at P-points -!> @param[out] CHI velocity potential (m^2/s) at P-points -!> @param[out] PSI streamfunction (m^2/s) at P-points +!> @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 @@ -4520,7 +4520,6 @@ END SUBROUTINE CALSLR_UUTAH !> !> @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 @@ -4531,7 +4530,7 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) use gridspec_mod, only: gridtype, dyval use upp_math, only: DVDXDUDY, DDVDX, DDUDY, UUAVG use mpi - +! implicit none ! ! DECLARE VARIABLES. @@ -4550,8 +4549,8 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) 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 - real*8 ta,tb,tc,td,de,tf - integer ier,jjk, mype + real*8 ta,tb,tc,td,de,tf + integer ier,jjk, mype ! !*************************************************************************** ! START CALCHIPSI HERE. @@ -4559,26 +4558,26 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) ! 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 - + 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)) + & wrk3(ista:iend,jsta:jend), cosl(ista_2l:iend_2u,jsta_2l:jend_2u)) + allocate(iw(im),ie(im)) - imb2 = im/2 + imb2 = im/2 !$omp parallel do private(i) do i=ista,iend ie(i) = i+1 @@ -4586,384 +4585,383 @@ SUBROUTINE CALCHIPSI (UWND,VWND,CHI,PSI) 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 + DO J=JSTA,JEND + do i=ista,iend + ip1 = ie(i) + im1 = iw(i) + cosl(i,j) = cos(gdlat(i,j)*dtr) + IF(cosl(i,j) >= 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 -! 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) - + 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 + 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 + 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 + 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) - + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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) & + 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 + 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 - 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 + tx2(i+1) = tx1(i) enddo - do i=ista,iend - absv(i,j) = tx1(i) + 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 - endif - END DO ! end of J loop + enddo + do i=ista,iend + absv(i,j) = tx1(i) + enddo + endif + END DO ! end of J loop -! deallocate (wrk1, wrk2, wrk3, cosl) +! 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) - +! + 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)) - +! +! poisson solver for psi and chi PSI=0. - ta=mpi_wtime() - do jjk=1,1000 + ta=mpi_wtime() + 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. J1 .and. J1 .and. J