diff --git a/src/ecwam/CMakeLists.txt b/src/ecwam/CMakeLists.txt index 9f97b48f..9514524c 100644 --- a/src/ecwam/CMakeLists.txt +++ b/src/ecwam/CMakeLists.txt @@ -291,6 +291,7 @@ list( APPEND ecwam_srcs wgrib_edition.F90 wgribencode.F90 wgribencode_model.F90 + wgribencode_values.F90 wgribenout.F90 wgribout.F90 wnfluxes.F90 diff --git a/src/ecwam/wgribencode.F90 b/src/ecwam/wgribencode.F90 index fc6353cd..4030a36c 100644 --- a/src/ecwam/wgribencode.F90 +++ b/src/ecwam/wgribencode.F90 @@ -1,5 +1,5 @@ ! (C) Copyright 1989- ECMWF. -! +! ! This software is licensed under the terms of the Apache Licence Version 2.0 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. ! In applying this licence, ECMWF does not waive the privileges and immunities @@ -82,6 +82,7 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & IMPLICIT NONE #include "abort1.intfb.h" #include "difdate.intfb.h" +#include "wgribencode_values.intfb.h" INTEGER(KIND=JWIM), INTENT(IN) :: IU06, ITEST, I1, I2 INTEGER(KIND=JWIM), INTENT(IN) :: ITABLE, IPARAM, KLEV, ITMIN, ITMAX, IK, IM, IFCST @@ -137,6 +138,7 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & INTEGER(KIND=JWIM) :: NWINOFF INTEGER(KIND=JWIM) :: NPROMA, MTHREADS, JC, JCS, JCL, JJ, ITHRS + REAL(KIND=JWRB) :: TEMP REAL(KIND=JWRB) :: ZMINSPEC, PPMAX, PPMIN, DELTAPP, ABSPPREC REAL(KIND=JPHOOK) :: ZHOOK_HANDLE @@ -180,114 +182,15 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & !* 0. PUT FIELD INTO GLOBAL MATRIX VALUES. ! ----------------------------------- - - IF (IRGG == 1 .OR. CLDOMAIN == 'm' ) THEN - ICOUNT=1 - ELSEIF (CLDOMAIN == 's' ) THEN - ICOUNT=1 - ELSE - ICOUNT = (NINT((90._JWRB - AMONOP ) / XDELLA))*I1 + 1 - VALUES(:)=ZMISS - ENDIF - -! PAD THE POLES IF INCLUDED IN THE GRID - IF (LPADPOLES) THEN - IF ((NINT((90._JWRB - AMONOP ) / XDELLA)) == 0) THEN - TEMP=0._JWRB - NN=0 - J=2 - JSN=I2-J+1 - DO I=1,NLONRGG(JSN) - IF (FIELD(I,J) /= ZMISS) THEN - TEMP=TEMP+FIELD(I,J) - NN=NN+1 - ENDIF - ENDDO - IF (NN > 0) THEN - TEMP=TEMP/NN - ELSE - TEMP=ZMISS - ENDIF - J=1 - JSN=I2-J+1 - DO I=1,NLONRGG(JSN) - FIELD(I,J)=TEMP - ENDDO - ENDIF - IF ((NINT((-90._JWRB - AMOSOP ) / XDELLA)) == 0) THEN - TEMP=0._JWRB - NN=0 - J=I2-1 - JSN=I2-J+1 - DO I=1,NLONRGG(JSN) - IF (FIELD(I,J) /= ZMISS) THEN - TEMP=TEMP+FIELD(I,J) - NN=NN+1 - ENDIF - ENDDO - IF (NN > 0) THEN - TEMP=TEMP/NN - ELSE - TEMP=ZMISS - ENDIF - J=I2 - JSN=I2-J+1 - DO I=1,NLONRGG(JSN) - FIELD(I,J)=TEMP - ENDDO - ENDIF - ENDIF - -! FILL THE ENCODING ARRAY -! ----------------------- - DO J = 1, I2 - JSN = I2-J+1 - VALUES(ICOUNT:ICOUNT+NLONRGG(JSN)-1)=FIELD(1:NLONRGG(JSN), J) - ICOUNT=ICOUNT+NLONRGG(JSN) - ENDDO - - IF (ITABPAR == 140251 .AND. .NOT. LLSPECNOT251 ) THEN - -! spectra are encoded on a log10 scale and small values below a certain threshold are set to missing - DELTAPP=(2**NGRBRESS-1)*PPRESOL - ABSPPREC=ABS(PPREC) - ZMINSPEC = LOG10(PPEPS)+ABSPPREC - - ALLOCATE(VALM(MTHREADS)) -!$OMP PARALLEL DO SCHEDULE(STATIC) PRIVATE(JC, JCS, JCL, JJ, ITHRS) - DO JC= 1, NTENCODE, NPROMA - JCS = JC - JCL = MIN(JCS+NPROMA-1, NTENCODE) - ITHRS = JC/NPROMA + 1 - VALM(ITHRS) = ZMINSPEC - DO JJ = JCS, JCL - IF(VALUES(JJ) /= ZMISS ) THEN - VALUES(JJ) = LOG10(VALUES(JJ)+PPEPS)+ABSPPREC - VALM(ITHRS) = MAX(VALM(ITHRS), VALUES(JJ)) - ELSE - VALUES(JJ) = ZMINSPEC - ENDIF - ENDDO - ENDDO -!$OMP END PARALLEL DO - - PPMAX=MAXVAL(VALM(:)) - DEALLOCATE(VALM) - - PPMIN=MIN(PPMISS, PPMAX-DELTAPP) - PPMIN=MIN(PPMIN, PPMIN_RESET) - -!$OMP PARALLEL DO SCHEDULE(STATIC) PRIVATE(JC, JCS, JCL, JJ) - DO JC= 1, NTENCODE, NPROMA - JCS = JC - JCL = MIN(JCS+NPROMA-1, NTENCODE) - DO JJ = JCS, JCL - IF ( VALUES(JJ) <= PPMIN ) VALUES(JJ)=ZMISS - ENDDO - ENDDO -!$OMP END PARALLEL DO - - ENDIF + CALL WGRIBENCODE_VALUES ( I1, I2, & +& FIELD, & +& ITABPAR, LLSPECNOT251, & +& PPMISS, PPEPS, PPREC, PPRESOL, PPMIN_RESET, NTENCODE, & +& NGRBRESS, LPADPOLES, & +& NLONRGG_SIZE, NLONRGG, IRGG, & +& AMONOP, AMOSOP, XDELLA, CLDOMAIN, & +& ZMISS, & +& VALUES) !* 1. FIX PARAMETERS AND PACK DATA. ! ----------------------------- @@ -298,8 +201,8 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & IF ( ITMIN /= 0 .OR. ITMAX /= 0 ) THEN ! NEED TO CHANGE TO SPECIFIC TEMPLATE FOR ENCODING ITMIN AND ITMAX CALL IGRIB_GET_VALUE(IGRIB_HANDLE,'numberOfForecastsInEnsemble',IDUM,KRET=IRET ) - IF ( IRET == 0 ) THEN - IF ( IDUM > 0 ) THEN + IF ( IRET == 0 ) THEN + IF ( IDUM > 0 ) THEN CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'productDefinitionTemplateNumber', NTRG2TMPP) ELSE CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'productDefinitionTemplateNumber', NTRG2TMPD) @@ -316,7 +219,7 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'scaleFactorOfUpperWavePeriodLimit', 0) CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'scaledValueOfUpperWavePeriodLimit', ITMAX) ELSEIF ( ITMIN /= 0 ) THEN -! [ ITMIN +! [ ITMIN CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'typeOfWavePeriodInterval', 3) CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'scaleFactorOfLowerWavePeriodLimit', 0) CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'scaledValueOfLowerWavePeriodLimit', ITMIN) @@ -350,24 +253,24 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'paramId',ITABPAR) WRITE(NULERR,*) ' THE PARAMETER SHOULD BE ADDED TO THE LIST OF' - WRITE(NULERR,*) ' PARAMETERS KNOWN BY ECCODES !!!' - WRITE(NULERR,*) ' IN THE MEAN TIME, THE PROGRAM WILL CONTINUE.' - WRITE(NULERR,*) ' USING EXPERIMENTAL PARAMETER TABLE 212' + WRITE(NULERR,*) ' PARAMETERS KNOWN BY ECCODES !!!' + WRITE(NULERR,*) ' IN THE MEAN TIME, THE PROGRAM WILL CONTINUE.' + WRITE(NULERR,*) ' USING EXPERIMENTAL PARAMETER TABLE 212' WRITE(NULERR,*) ' WITH paramId= ', ITABPAR WRITE(NULERR,*) ' *********************************************' WRITE(IU06,*) ' THE PARAMETER SHOULD BE ADDED TO THE LIST OF' - WRITE(IU06,*) ' PARAMETERS KNOWN BY ECCODES !!!' - WRITE(IU06,*) ' USING EXPERIMENTAL PARAMETER TABLE 212' + WRITE(IU06,*) ' PARAMETERS KNOWN BY ECCODES !!!' + WRITE(IU06,*) ' USING EXPERIMENTAL PARAMETER TABLE 212' WRITE(IU06,*) ' WITH paramId= ', ITABPAR WRITE(IU06,*) ' *********************************************' ELSE WRITE(NULERR,*) ' THE PARAMETER SHOULD BE ADDED TO THE LIST OF' - WRITE(NULERR,*) ' PARAMETERS KNOWN BY ECCODES !!!' - WRITE(NULERR,*) ' ' + WRITE(NULERR,*) ' PARAMETERS KNOWN BY ECCODES !!!' + WRITE(NULERR,*) ' ' WRITE(NULERR,*) ' IN THE MEANTIME, YOU MIGHT WANT TO USE EXPERIMENTAL PARAMETER TABLE 212' - WRITE(NULERR,*) ' SET LLRSTGRIBPARAM TO TRUE IN THE INPUT NAMELIST AND RERUN' - WRITE(NULERR,*) ' ADAPT THE ARCHIVING ACCORDINGLY.' + WRITE(NULERR,*) ' SET LLRSTGRIBPARAM TO TRUE IN THE INPUT NAMELIST AND RERUN' + WRITE(NULERR,*) ' ADAPT THE ARCHIVING ACCORDINGLY.' WRITE(NULERR,*) ' ************************************************************************' CALL ABORT1 ENDIF @@ -387,7 +290,7 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & CALL IGRIB_GET_VALUE(IGRIB_HANDLE,'level',IDUM,KRET=IRET) IF ( IRET == JPGRIB_SUCCESS ) CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'level',KLEV) - IF (.NOT. LGRHDIFS .OR. & + IF (.NOT. LGRHDIFS .OR. & & (MARSTYPE == 'an' .AND. IFCST == 0) .OR. & & (MARSTYPE == '4v') .OR. & (MARSTYPE == 'fg' .AND. IFCST == 0) ) THEN @@ -403,7 +306,7 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & !!! for compatibility with previous coding, impose: CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'timeRangeIndicator',10) CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'endStep',IFCST) - + IF ( MARSTYPE == 'fg' .AND. IFCST == 0 ) THEN ICLASS = 1 ELSEIF ( MARSTYPE == '4v' ) THEN @@ -442,7 +345,7 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & ! this only works for 12 hour analysis windows !!!1 IF (IH2+3 == 12 ) THEN NWINOFF=0 - ELSE + ELSE NWINOFF=12-MOD(IH2+3,12) ENDIF IF ( ITABPAR == 140251 .AND. IGRIB_VERSION == 1 ) THEN @@ -452,11 +355,11 @@ SUBROUTINE WGRIBENCODE ( IU06, ITEST, & ! in hours CALL IGRIB_GET_VALUE(IGRIB_HANDLE,'offsetToEndOf4DvarWindow',IDUM,KRET=IRET) IF ( IRET == 0 ) CALL IGRIB_SET_VALUE(IGRIB_HANDLE,'offsetToEndOf4DvarWindow',NWINOFF) - ENDIF + ENDIF ELSEIF ( MARSTYPE == 'cf' ) THEN - ICLASS = 10 + ICLASS = 10 ELSEIF ( MARSTYPE == 'pf' ) THEN - ICLASS = 11 + ICLASS = 11 ELSE ICLASS = 9 ENDIF diff --git a/src/ecwam/wgribencode_values.F90 b/src/ecwam/wgribencode_values.F90 new file mode 100644 index 00000000..e62650c9 --- /dev/null +++ b/src/ecwam/wgribencode_values.F90 @@ -0,0 +1,203 @@ +! (C) Copyright 1989- ECMWF. +! +! This software is licensed under the terms of the Apache Licence Version 2.0 +! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +! In applying this licence, ECMWF does not waive the privileges and immunities +! granted to it by virtue of its status as an intergovernmental organisation +! nor does it submit to any jurisdiction. +! + +SUBROUTINE WGRIBENCODE_VALUES ( I1, I2, & +& FIELD, & +& ITABPAR, LLSPECNOT251, & +& PPMISS, PPEPS, PPREC, PPRESOL, PPMIN_RESET, NTENCODE, & +& NGRBRESS, LPADPOLES, & +& NLONRGG_SIZE, NLONRGG, IRGG, & +& AMONOP, AMOSOP, XDELLA, CLDOMAIN, & +& ZMISS, & +& VALUES ) + + +! ---------------------------------------------------------------------- + USE PARKIND_WAVE, ONLY : JWIM, JWRB, JWRU + + USE YOWGRIBHD, ONLY : NTRG2TMPD, NTRG2TMPP, LLRSTGRIBPARAM + + USE YOWGRIB , ONLY : IGRIB_GET_VALUE, IGRIB_SET_VALUE, JPGRIB_SUCCESS + USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK + USE EC_LUN , ONLY : NULERR + USE OML_MOD , ONLY : OML_GET_MAX_THREADS + +! ---------------------------------------------------------------------- + IMPLICIT NONE +#include "abort1.intfb.h" + + INTEGER(KIND=JWIM), INTENT(IN) :: I1, I2 + INTEGER(KIND=JWIM), INTENT(IN) :: ITABPAR + LOGICAL, INTENT(IN) :: LLSPECNOT251 + + REAL(KIND=JWRB), INTENT(INOUT) :: FIELD(I1,I2) + + ! From yowgribhd + REAL(KIND=JWRB), INTENT(IN) :: PPMISS ! every spectral values less or equal ppmiss are replaced by the missing data indicator + REAL(KIND=JWRB), INTENT(IN) :: PPEPS ! Small number used in spectral packing of 140251 + REAL(KIND=JWRB), INTENT(IN) :: PPREC ! Reference value for spectral packing of 140251 + REAL(KIND=JWRB), INTENT(IN) :: PPRESOL ! Maximun resolution possible when encoding spectra (parameter 140251). + REAL(KIND=JWRB), INTENT(IN) :: PPMIN_RESET ! Can be used to set the minimum of ppmin in wgribout to a lower value. + INTEGER(KIND=JWIM), INTENT(IN) :: NTENCODE ! Total number of grid points for encoding + INTEGER(KIND=JWIM), INTENT(IN) :: NGRBRESS ! Number of bits used to encode spectra + LOGICAL, INTENT(IN) :: LPADPOLES ! True if poles are padded when savind to grib. + + ! From yowgrid + INTEGER(KIND=JWIM), INTENT(IN) :: NLONRGG_SIZE + INTEGER(KIND=JWIM), INTENT(IN) :: NLONRGG(NLONRGG_SIZE) + + ! From yowmap + INTEGER(KIND=JWIM), INTENT(IN) :: IRGG ! Grid code: 0 = regular, 1 = irregular. + REAL(KIND=JWRB), INTENT(IN) :: AMONOP ! Most Northern latitude in grid (degree). + REAL(KIND=JWRB), INTENT(IN) :: AMOSOP ! Most Southern latitude in grid (degree). + REAL(KIND=JWRB), INTENT(IN) :: XDELLA ! Grid increment for latitude (degree). + + ! From yowparam + CHARACTER(LEN=1), INTENT(IN) :: CLDOMAIN ! Defines the domain of the model (for the fdb and for selection of some variables) + + ! From yowpcons + REAL(KIND=JWRB), INTENT(IN) :: ZMISS ! Missing data indicator (set in chief or via the ifs). + + REAL(KIND=JWRB), DIMENSION(:), INTENT(INOUT) :: VALUES + + +! ---------------------------------------------------------------------- + INTEGER(KIND=JWIM) :: ICOUNT, NN, I, J, JSN + INTEGER(KIND=JWIM) :: NPROMA, MTHREADS, JC, JCS, JCL, JJ, ITHRS +!$ INTEGER,EXTERNAL :: OMP_GET_MAX_THREADS + + REAL(KIND=JWRB) :: TEMP + REAL(KIND=JWRB) :: ZMINSPEC, PPMAX, PPMIN, DELTAPP, ABSPPREC + + + + REAL(KIND=JPHOOK) :: ZHOOK_HANDLE + REAL(KIND=JWRB), ALLOCATABLE :: VALM(:) + +! ---------------------------------------------------------------------- + IF (LHOOK) CALL DR_HOOK('WGRIBENCODE_VALUES',0,ZHOOK_HANDLE) + + MTHREADS=OMP_GET_MAX_THREADS() + NPROMA=NTENCODE/MTHREADS + 1 + + +!* 0. PUT FIELD INTO GLOBAL MATRIX VALUES. +! ----------------------------------- + + IF (IRGG == 1 .OR. CLDOMAIN == 'm' ) THEN + ICOUNT=1 + ELSEIF (CLDOMAIN == 's' ) THEN + ICOUNT=1 + ELSE + ICOUNT = (NINT((90._JWRB - AMONOP ) / XDELLA))*I1 + 1 + VALUES(:)=ZMISS + ENDIF + +! PAD THE POLES IF INCLUDED IN THE GRID + IF (LPADPOLES) THEN + IF ((NINT((90._JWRB - AMONOP ) / XDELLA)) == 0) THEN + TEMP=0._JWRB + NN=0 + J=2 + JSN=I2-J+1 + DO I=1,NLONRGG(JSN) + IF (FIELD(I,J) /= ZMISS) THEN + TEMP=TEMP+FIELD(I,J) + NN=NN+1 + ENDIF + ENDDO + IF (NN > 0) THEN + TEMP=TEMP/NN + ELSE + TEMP=ZMISS + ENDIF + J=1 + JSN=I2-J+1 + DO I=1,NLONRGG(JSN) + FIELD(I,J)=TEMP + ENDDO + ENDIF + IF ((NINT((-90._JWRB - AMOSOP ) / XDELLA)) == 0) THEN + TEMP=0._JWRB + NN=0 + J=I2-1 + JSN=I2-J+1 + DO I=1,NLONRGG(JSN) + IF (FIELD(I,J) /= ZMISS) THEN + TEMP=TEMP+FIELD(I,J) + NN=NN+1 + ENDIF + ENDDO + IF (NN > 0) THEN + TEMP=TEMP/NN + ELSE + TEMP=ZMISS + ENDIF + J=I2 + JSN=I2-J+1 + DO I=1,NLONRGG(JSN) + FIELD(I,J)=TEMP + ENDDO + ENDIF + ENDIF + +! FILL THE ENCODING ARRAY +! ----------------------- + DO J = 1, I2 + JSN = I2-J+1 + VALUES(ICOUNT:ICOUNT+NLONRGG(JSN)-1)=FIELD(1:NLONRGG(JSN), J) + ICOUNT=ICOUNT+NLONRGG(JSN) + ENDDO + + IF (ITABPAR == 140251 .AND. .NOT. LLSPECNOT251 ) THEN + +! spectra are encoded on a log10 scale and small values below a certain threshold are set to missing + DELTAPP=(2**NGRBRESS-1)*PPRESOL + ABSPPREC=ABS(PPREC) + ZMINSPEC = LOG10(PPEPS)+ABSPPREC + + ALLOCATE(VALM(MTHREADS)) +!$OMP PARALLEL DO SCHEDULE(STATIC) PRIVATE(JC, JCS, JCL, JJ, ITHRS) + DO JC= 1, NTENCODE, NPROMA + JCS = JC + JCL = MIN(JCS+NPROMA-1, NTENCODE) + ITHRS = JC/NPROMA + 1 + VALM(ITHRS) = ZMINSPEC + DO JJ = JCS, JCL + IF(VALUES(JJ) /= ZMISS ) THEN + VALUES(JJ) = LOG10(VALUES(JJ)+PPEPS)+ABSPPREC + VALM(ITHRS) = MAX(VALM(ITHRS), VALUES(JJ)) + ELSE + VALUES(JJ) = ZMINSPEC + ENDIF + ENDDO + ENDDO +!$OMP END PARALLEL DO + + PPMAX=MAXVAL(VALM(:)) + DEALLOCATE(VALM) + + PPMIN=MIN(PPMISS, PPMAX-DELTAPP) + PPMIN=MIN(PPMIN, PPMIN_RESET) + +!$OMP PARALLEL DO SCHEDULE(STATIC) PRIVATE(JC, JCS, JCL, JJ) + DO JC= 1, NTENCODE, NPROMA + JCS = JC + JCL = MIN(JCS+NPROMA-1, NTENCODE) + DO JJ = JCS, JCL + IF ( VALUES(JJ) <= PPMIN ) VALUES(JJ)=ZMISS + ENDDO + ENDDO +!$OMP END PARALLEL DO + + ENDIF + + IF (LHOOK) CALL DR_HOOK('WGRIBENCODE_VALUES',1,ZHOOK_HANDLE) + +END SUBROUTINE WGRIBENCODE_VALUES