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)