#include "wwm_functions.h"
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_WIND_INPUT
#ifdef NCDF
USE NETCDF
#endif
USE DATAPOOL
IMPLICIT NONE
INTEGER :: IT, IFILE
REAL(rkind) :: WDIRT
REAL(rkind) :: cf_w1, cf_w2
WINDXY(:,:) = 0.0
#ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
MNP_WIND=MNP
allocate(XP_WIND(MNP_WIND), YP_WIND(MNP_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 1')
XP_WIND=XP
YP_WIND=YP
ELSE
MNP_WIND=np_total
IF (myrank .eq. 0) THEN
allocate(XP_WIND(MNP_WIND), YP_WIND(MNP_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 1')
XP_WIND=XPtotal
YP_WIND=YPtotal
END IF
END IF
#else
MNP_WIND=MNP
allocate(XP_WIND(MNP_WIND), YP_WIND(MNP_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 1')
XP_WIND=XP
YP_WIND=YP
#endif
IF (LSTWD) THEN
IF (LCWIN) THEN
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'HOMOGENOUS STEADY WIND FIELD IS USED'
WRITE(WINDBG%FHNDL,'("+TRACE...",A,I10)') 'WIND IS COMING FROM WWM - WINDFORMAT', IWINDFORMAT, LWDIR
FLUSH(WINDBG%FHNDL)
END IF
IF (LWDIR) THEN
CALL DEG2NAUT(WDIR, WDIRT, LNAUTIN)
WINDXY(:,1) = WVEL * COS(WDIRT * DEGRAD)
WINDXY(:,2) = WVEL * SIN(WDIRT * DEGRAD)
ELSE
WINDXY(:,1) = CWINDX
WINDXY(:,2) = CWINDY
END IF
ELSE ! LCWIN
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("+TRACE...",A,I10)') 'WIND IS COMING FROM WWM - WINDFORMAT', IWINDFORMAT
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'SPATIAL VARIABLE WIND FIELD IS USED'
FLUSH(WINDBG%FHNDL)
END IF
IF (IWINDFORMAT == 1) THEN
CALL CSEVAL( WIN%FHNDL, TRIM(WIN%FNAME), .FALSE., 2, WINDXY, MULTIPLE_IN_WIND)
#ifdef NCDF
ELSE IF (IWINDFORMAT == 2) THEN ! NETCDF created using ncl_convert2nc using DWD grib
CALL INIT_NETCDF_DWD
CALL FIND_WIND_NEAREST_LOWER_IDX(SEWI%BMJD, idxWind)
IFILE=WIND_TIME_IFILE(idxWind)
IT=WIND_TIME_IT(idxWind)
CALL READ_NETCDF_DWD(IFILE, IT, WINDXY)
ELSE IF (IWINDFORMAT == 3) THEN ! NETCDF created using cdo -f nc copy file.grb file.nc this is CFRS
CALL INIT_NETCDF_CRFS
CALL FIND_WIND_NEAREST_LOWER_IDX(SEWI%BMJD, idxWind)
IFILE=WIND_TIME_IFILE(idxWind)
IT=WIND_TIME_IT(idxWind)
CALL READ_NETCDF_CRFS(IFILE, IT, WINDXY)
ELSE IF (IWINDFORMAT == 4) THEN ! NETCDF NARR downloaded from NOMAD
CALL INIT_NETCDF_NARR
CALL FIND_WIND_NEAREST_LOWER_IDX(SEWI%BMJD, idxWind)
IFILE=WIND_TIME_IFILE(idxWind)
IT=WIND_TIME_IT(idxWind)
CALL READ_NETCDF_NARR(IFILE, IT, WINDXY)
ELSE IF (IWINDFORMAT == 5) THEN ! NETCDF CF_COMPLIANT STATIONARY FIELD
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'COMPUTING CF INTERPOLATION COEFS AND LOADING WIND_TIME_MJD'
FLUSH(WINDBG%FHNDL)
END IF
CALL INIT_NETCDF_CF_WWM_WIND(eVAR_WIND)
ALLOCATE(tmp_wind1(MNP,2),tmp_wind2(MNP,2), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 1')
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
CALL READ_INTERP_NETCDF_CF_WWM_WIND(REC1_wind_new,tmp_wind1)
IF (cf_w1.NE.1) THEN
CALL READ_INTERP_NETCDF_CF_WWM_WIND(REC2_wind_new,tmp_wind2)
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
ELSE IF (IWINDFORMAT == 6) THEN ! DIRECT WWM forcing (no interp)
CALL INIT_DIRECT_NETCDF_CF(eVAR_WIND, MULTIPLE_IN_WIND, WIN%FNAME, "Uwind")
ALLOCATE(tmp_wind1(MNP,2),tmp_wind2(MNP,2), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 1')
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
CALL READ_DIRECT_NETCDF_CF(eVAR_wind, REC1_wind_new,tmp_wind1)
IF (cf_w1.NE.1) THEN
CALL READ_DIRECT_NETCDF_CF(eVAR_wind, REC2_wind_new,tmp_wind2)
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
#endif
#ifdef GRIB_API_ECMWF
ELSE IF (IWINDFORMAT == 7) THEN ! GRIB forcing from ecmwf
ALLOCATE(tmp_wind1(MNP,2),tmp_wind2(MNP,2), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 1')
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
CALL READ_GRIB_WIND(REC1_wind_new,tmp_wind1)
IF (cf_w1.NE.1) THEN
CALL READ_GRIB_WIND(REC2_wind_new,tmp_wind2)
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
#endif
ELSE
CALL wwm_abort('Wrong choice of IWINDFORMAT (maybe need NETCDF or GRIB)')
END IF
ENDIF
ELSE IF (LSEWD) THEN
IF (LCWIN) THEN
CALL wwm_abort('LSEWD + LCWIN NOT READY')
! CALL READ_WIND_TIME_SERIES(IT) ! set time according to wwminput.nml and get initial time step
! CALL SET_INITIAL_WIND(IT) !
ELSE
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("+TRACE...",A,I10)') 'WIND IS COMING FROM WWM - WINDFORMAT', IWINDFORMAT
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'NONSTATIONARY WIND FIELD IS USED '
FLUSH(WINDBG%FHNDL)
END IF
SEWI%TOTL = (SEWI%EMJD - SEWI%BMJD) * DAY2SEC
SEWI%ISTP = NINT( SEWI%TOTL / SEWI%DELT ) + 1
SEWI%TMJD = SEWI%BMJD
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) SEWI%BEGT, SEWI%ENDT, SEWI%ISTP, SEWI%TOTL/3600.0, SEWI%DELT
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'SPATIAL VARIABLE WIND FIELD IS USED'
WRITE(WINDBG%FHNDL,*) 'IWINDFORMAT=', IWINDFORMAT
FLUSH(WINDBG%FHNDL)
END IF
IF (IWINDFORMAT == 1) THEN
OPEN(WIN%FHNDL, FILE = TRIM(WIN%FNAME), STATUS = 'OLD')
CALL CSEVAL( WIN%FHNDL, TRIM(WIN%FNAME), .TRUE., 2, WINDXY, MULTIPLE_IN_WIND)
#ifdef NCDF
ELSE IF (IWINDFORMAT == 2) THEN ! NETCDF created using ncl_convert2nc using DWD grib
CALL INIT_NETCDF_DWD
CALL FIND_WIND_NEAREST_LOWER_IDX(MAIN%TMJD, idxWind)
IFILE=WIND_TIME_IFILE(idxWind)
IT=WIND_TIME_IT(idxWind)
CALL READ_NETCDF_DWD(IFILE, IT, WINDXY)
ELSE IF (IWINDFORMAT == 3) THEN ! NETCDF created using cdo -f nc copy file.grb file.nc
CALL INIT_NETCDF_CRFS
CALL FIND_WIND_NEAREST_LOWER_IDX(MAIN%TMJD, idxWind)
IFILE=WIND_TIME_IFILE(idxWind)
IT=WIND_TIME_IT(idxWind)
CALL READ_NETCDF_CRFS(IFILE, IT, WINDXY)
ELSE IF (IWINDFORMAT == 4) THEN ! NETCDF created using cdo -f nc copy file.grb file.nc
CALL INIT_NETCDF_NARR
CALL FIND_WIND_NEAREST_LOWER_IDX(MAIN%TMJD, idxWind)
IFILE=WIND_TIME_IFILE(idxWind)
IT=WIND_TIME_IT(idxWind)
CALL READ_NETCDF_NARR(IFILE, IT, WINDXY)
ELSE IF (IWINDFORMAT == 5) THEN
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'SPATIAL/TEMPORAL VARIABLE WIND FIELD IS USED CF NETCDF'
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'COMPUTING CF INTERPOLATION COEFS AND LOADING WIND_TIME_MJD'
FLUSH(WINDBG%FHNDL)
END IF
CALL INIT_NETCDF_CF_WWM_WIND(eVAR_WIND)
ALLOCATE(tmp_wind1(MNP,2), tmp_wind2(MNP,2), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 2')
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
CALL READ_INTERP_NETCDF_CF_WWM_WIND(REC1_wind_new,tmp_wind1)
IF (cf_w1.NE.1) THEN
CALL READ_INTERP_NETCDF_CF_WWM_WIND(REC2_wind_new,tmp_wind2)
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
ELSE IF (IWINDFORMAT == 6) THEN
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'SPATIAL/TEMPORAL VARIABLE WIND FIELD IS USED CF NETCDF'
WRITE(WINDBG%FHNDL,'("+TRACE...",A)') 'COMPUTING CF INTERPOLATION COEFS AND LOADING WIND_TIME_MJD'
FLUSH(WINDBG%FHNDL)
END IF
CALL INIT_DIRECT_NETCDF_CF(eVAR_WIND, MULTIPLE_IN_WIND, WIN%FNAME, "Uwind")
ALLOCATE(tmp_wind1(MNP,2), tmp_wind2(MNP,2), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 2')
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
CALL READ_DIRECT_NETCDF_CF(eVAR_wind, REC1_wind_new,tmp_wind1)
IF (cf_w1.NE.1) THEN
CALL READ_DIRECT_NETCDF_CF(eVAR_wind, REC2_wind_new,tmp_wind2)
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
#endif
#ifdef GRIB_API_ECMWF
ELSE IF (IWINDFORMAT == 7) THEN
CALL INIT_GRIB_WIND
ALLOCATE(tmp_wind1(MNP,2), tmp_wind2(MNP,2), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 2')
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
CALL READ_GRIB_WIND(REC1_wind_new,tmp_wind1)
IF (cf_w1.NE.1) THEN
CALL READ_GRIB_WIND(REC2_wind_new,tmp_wind2)
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
#endif
ENDIF
ENDIF
ENDIF
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("+TRACE... Done with CF init, Uwind ",F7.2,2x,F7.2)')minval(WINDXY(:,1)),maxval(WINDXY(:,1))
WRITE(WINDBG%FHNDL,'("+TRACE... Done with CF init, Vwind ",F7.2,2x,F7.2)')minval(WINDXY(:,2)),maxval(WINDXY(:,2))
FLUSH(WINDBG%FHNDL)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE UPDATE_WIND(K)
#ifdef NCDF
USE NETCDF
#endif
USE DATAPOOL
IMPLICIT NONE
REAL(rkind) :: TMP(MNP,2)
#if defined NCDF || defined GRIB_API_ECMWF
REAL(rkind) :: cf_w1, cf_w2
INTEGER :: IT, IFILE
#endif
INTEGER, intent(in) :: K
!AR: All crap ... defining K without using means that nobody has ever checked the results or anything else, so why coding at all?
!AR: Mathieu can you please fix this !!!
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'MAIN%TMJD=', MAIN%TMJD
WRITE(WINDBG%FHNDL,*) 'SEWI(TMJD,EMJD)=', SEWI%TMJD, SEWI%EMJD
END IF
IF ( LSEWD .AND. (MAIN%TMJD .ge. SEWI%TMJD-1.E-8) .AND. (MAIN%TMJD .le. SEWI%EMJD+1.e-8) ) THEN
IF (IWINDFORMAT == 1) THEN
!NDM: Need to add the facility for LINTERWD
CALL CSEVAL( WIN%FHNDL, WIN%FNAME, .TRUE., 2, TMP, MULTIPLE_IN_WIND)
DVWIND = (TMP-WINDXY)/SEWI%DELT*MAIN%DELT
#ifdef NCDF
ELSE IF (IWINDFORMAT == 2) THEN ! DWD_NETCDF
CALL MOVE_BY_ONE_INDEX(IFILE, IT)
CALL READ_NETCDF_DWD(IFILE, IT, TMP)
DVWIND = (TMP-WINDXY)/SEWI%DELT*MAIN%DELT
ELSE IF (IWINDFORMAT == 3) THEN ! NOAA CFRS ... the 1st step is analysis and then we have 5 + 1 forecasts, which give one the option to use either only the 6 forecast's after the analysis or use the analysis with 5 forecast's
CALL MOVE_BY_ONE_INDEX(IFILE, IT)
CALL READ_NETCDF_CRFS(IFILE, IT, TMP)
DVWIND = (TMP-WINDXY)/SEWI%DELT*MAIN%DELT
ELSE IF (IWINDFORMAT == 4) THEN ! NOAA NARR ...
CALL MOVE_BY_ONE_INDEX(IFILE, IT)
CALL READ_NETCDF_NARR(IFILE, IT, TMP)
DVWIND = (TMP-WINDXY)/SEWI%DELT*MAIN%DELT
ELSE IF (IWINDFORMAT == 5) THEN
IF (K.EQ.1) THEN
REC1_wind_old = 0
REC2_wind_old = 0
END IF
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
IF (REC1_wind_new.NE.REC1_wind_old) THEN
CALL READ_INTERP_NETCDF_CF_WWM_WIND(REC1_wind_new,tmp_wind1)
END IF
IF (REC2_wind_new.NE.REC2_wind_old) THEN
CALL READ_INTERP_NETCDF_CF_WWM_WIND(REC2_wind_new,tmp_wind2)
END IF
IF (cf_w1.NE.1) THEN
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
REC1_wind_old = REC1_wind_new
REC2_wind_old = REC2_wind_new
ELSE IF (IWINDFORMAT == 6) THEN
IF (K.EQ.1) THEN
REC1_wind_old = 0
REC2_wind_old = 0
END IF
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
IF (REC1_wind_new.NE.REC1_wind_old) THEN
CALL READ_DIRECT_NETCDF_CF(eVAR_wind, REC1_wind_new,tmp_wind1)
END IF
IF (REC2_wind_new.NE.REC2_wind_old) THEN
CALL READ_DIRECT_NETCDF_CF(eVAR_wind, REC2_wind_new,tmp_wind2)
END IF
IF (cf_w1.NE.1) THEN
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
REC1_wind_old = REC1_wind_new
REC2_wind_old = REC2_wind_new
#endif
#ifdef GRIB_API_ECMWF
ELSE IF (IWINDFORMAT == 7) THEN
IF (K.EQ.1) THEN
REC1_wind_old = 0
REC2_wind_old = 0
END IF
CALL GET_CF_TIME_INDEX(eVAR_WIND, REC1_wind_new,REC2_wind_new,cf_w1,cf_w2)
IF (REC1_wind_new.NE.REC1_wind_old) THEN
CALL READ_GRIB_WIND(REC1_wind_new,tmp_wind1)
END IF
IF (REC2_wind_new.NE.REC2_wind_old) THEN
CALL READ_GRIB_WIND(REC2_wind_new,tmp_wind2)
END IF
IF (cf_w1.NE.1) THEN
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)+cf_w2*tmp_wind2(:,:)
ELSE
WINDXY(:,:) = cf_w1*tmp_wind1(:,:)
END IF
REC1_wind_old = REC1_wind_new
REC2_wind_old = REC2_wind_new
#endif
END IF
SEWI%TMJD = SEWI%TMJD + SEWI%DELT*SEC2DAY
END IF
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,'("max WINDXY:",2F7.2)')maxval(WINDXY(:,1)),maxval(WINDXY(:,2))
WRITE(WINDBG%FHNDL,'("min WINDXY:",2F7.2)')minval(WINDXY(:,1)),minval(WINDXY(:,2))
FLUSH(WINDBG%FHNDL)
END IF
IF (LWINDSWAN) THEN
WRITE(3333,*) SEWI%TMJD
WRITE(3333,*) WINDXY(:,1)
WRITE(3333,*) WINDXY(:,2)
END IF
IF (LSEWD.AND.(IWINDFORMAT.NE.5).AND.(IWINDFORMAT.NE.6) ) THEN
WINDXY = WINDXY + DVWIND
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE MOVE_BY_ONE_INDEX(IFILE, IT)
USE DATAPOOL
implicit none
integer, intent(out) :: IFILE, IT
idxWind =idxWind+1
IF (idxWind .gt. NDT_WIND_ALL_FILES) THEN
CALL WWM_ABORT('Need wind after the time')
END IF
IFILE=WIND_TIME_IFILE(idxWind)
IT=WIND_TIME_IT(idxWind)
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE FIND_WIND_NEAREST_LOWER_IDX(eTime, idx)
USE DATAPOOL
implicit none
real(rkind), intent(in) :: eTime
integer, intent(out) :: idx
CHARACTER(LEN=15) :: eTimeStr
integer eIdxF, eIdx
eIdxF=-1
DO eIdx=1,NDT_WIND_ALL_FILES
IF (WIND_TIME_ALL_FILES(eIdx) .le. eTime + THR8) THEN
eIdxF=eIdx
ENDIF
END DO
IF (eIdxF .eq. -1) THEN
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'NDT_WIND_ALL_FILES=', NDT_WIND_ALL_FILES
DO eIdx=1,NDT_WIND_ALL_FILES
CALL MJD2CT(WIND_TIME_ALL_FILES(eIdx),eTimeStr)
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) ' eIdx=', eIdx
WRITE(WINDBG%FHNDL,*) ' eTime=', WIND_TIME_ALL_FILES(eIdx)
WRITE(WINDBG%FHNDL,*) ' eTimeStr=', eTimeStr
END IF
END DO
IF (WRITEWINDBGFLAG == 1) CALL FLUSH(WINDBG%FHNDL)
CALL WWM_ABORT('We failed to find the wind index')
END IF
idx=eIdxF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE KERNEL_INTERP_UV_WINDFD(outwind)
USE DATAPOOL
IMPLICIT NONE
INTEGER I, J
REAL(rkind), INTENT(out) :: outwind(MNP_WIND,2)
REAL(rkind) :: Uw, Vw
INTEGER IX, IY
REAL(rkind) :: cf_scale_factor, cf_add_offset
cf_scale_factor = eVAR_WIND % cf_scale_factor
cf_add_offset = eVAR_WIND % cf_add_offset
DO I = 1, MNP_WIND
Uw=ZERO
Vw=ZERO
IX=CF_IX(I)
IY=CF_IY(I)
DO J=1,4
Uw=Uw + CF_COEFF(J,I)*UWIND_FD(IX+SHIFTXY(J,1),IY+SHIFTXY(J,2))
Vw=Vw + CF_COEFF(J,I)*VWIND_FD(IX+SHIFTXY(J,1),IY+SHIFTXY(J,2))
END DO
outwind(I,1)=Uw*cf_scale_factor + cf_add_offset
outwind(I,2)=Vw*cf_scale_factor + cf_add_offset
END DO
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'KERNEL_INTERP_UV_WINDFD'
WRITE(WINDBG%FHNDL,*) 'UWIND_FD, min/max=', minval(UWIND_FD), maxval(UWIND_FD)
WRITE(WINDBG%FHNDL,*) 'VWIND_FD, min/max=', minval(VWIND_FD), maxval(VWIND_FD)
WRITE(WINDBG%FHNDL,*) 'UWIND_FE, min/max=', minval(outwind(:,1)), maxval(outwind(:,1))
WRITE(WINDBG%FHNDL,*) 'VWIND_FE, min/max=', minval(outwind(:,2)), maxval(outwind(:,2))
WRITE(WINDBG%FHNDL,*) 'max(CF_COEFF)=', maxval(abs(CF_COEFF))
WRITE(WINDBG%FHNDL,*) 'cf_scale_factor=', cf_scale_factor
WRITE(WINDBG%FHNDL,*) 'cf_add_offset=', cf_add_offset
FLUSH(WINDBG%FHNDL)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE LOAD_INTERP_ARRAY(FileSave, success)
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
logical, intent(out) :: success
character(len=256), intent(in) :: FileSave
character (len = *), parameter :: CallFct = "LOAD_INTERP_ARRAY"
integer, allocatable :: CF_IX_GLOBAL(:), CF_IY_GLOBAL(:)
real(rkind), allocatable :: CF_COEFF_GLOBAL(:,:)
integer, allocatable :: ListFirstMNP(:)
integer, allocatable :: CF_IX_loc(:), CF_IY_loc(:)
real(rkind), allocatable :: CF_COEFF_loc(:,:)
integer iret, ncid, varid
integer IP, IPglob, iPROC, NPloc, IPloc
INQUIRE(FILE=TRIM(FileSave), EXIST=LPRECOMP_EXIST)
IF (LPRECOMP_EXIST .eqv. .FALSE.) THEN
success=.FALSE.
RETURN
END IF
success=.TRUE.
#ifdef MPI_PARALL_GRID
IF (myrank .eq. 0) THEN
#endif
allocate(CF_IX_GLOBAL(np_total), CF_IY_GLOBAL(np_total), CF_COEFF_GLOBAL(4,np_total), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
!
iret=nf90_open(TRIM(FileSave), NF90_NOWRITE, ncid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, iret)
!
iret=nf90_inq_varid(ncid, "CF_IX", varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
iret=NF90_GET_VAR(ncid, varid, CF_IX_GLOBAL)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
!
iret=nf90_inq_varid(ncid, "CF_IY", varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
iret=NF90_GET_VAR(ncid, varid, CF_IY_GLOBAL)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
!
iret=nf90_inq_varid(ncid, "CF_COEFF", varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
iret=NF90_GET_VAR(ncid, varid, CF_COEFF_GLOBAL)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
!
iret=nf90_close(ncid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 27, iret)
!
#ifdef MPI_PARALL_GRID
END IF
#endif
!
#ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .eqv. .FALSE.) THEN
CF_IX=CF_IX_GLOBAL
CF_IY=CF_IY_GLOBAL
CF_COEFF=CF_COEFF_GLOBAL
deallocate(CF_IX_GLOBAL, CF_IY_GLOBAL, CF_COEFF_GLOBAL)
ELSE
IF (myrank .eq. 0) THEN
allocate(ListFirstMNP(nproc), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
ListFirstMNP=0
DO iProc=2,nproc
ListFirstMNP(iProc)=ListFirstMNP(iProc-1) + ListMNP(iProc-1)
END DO
DO IP=1,MNP
IPglob=iplg(IP)
CF_IX(IP)=CF_IX_GLOBAL(IPglob)
CF_IY(IP)=CF_IY_GLOBAL(IPglob)
CF_COEFF(:,IP)=CF_COEFF_GLOBAL(:,IPglob)
END DO
!
DO iPROC=2,nproc
NPloc=ListMNP(iPROC)
allocate(CF_IX_loc(NPloc), CF_IY_loc(NPloc), CF_COEFF_loc(4,NPloc), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
DO IPloc=1,NPloc
IPglob=ListIPLG(IPloc + ListFirstMNP(iPROC))
CF_IX_loc(IPloc)=CF_IX_GLOBAL(IPglob)
CF_IY_loc(IPloc)=CF_IY_GLOBAL(IPglob)
CF_COEFF_loc(:, IPloc)=CF_COEFF_GLOBAL(:, IPglob)
END DO
CALL MPI_SEND(CF_IX_loc, NPloc, itype, iPROC-1, 711, comm, ierr)
CALL MPI_SEND(CF_IY_loc, NPloc, itype, iPROC-1, 712, comm, ierr)
CALL MPI_SEND(CF_COEFF_loc, 4*NPloc, rtype, iPROC-1, 713, comm, ierr)
deallocate(CF_IX_loc, CF_IY_loc, CF_COEFF_loc)
END DO
deallocate(ListFirstMNP)
deallocate(CF_IX_GLOBAL, CF_IY_GLOBAL, CF_COEFF_GLOBAL)
ELSE
CALL MPI_RECV(CF_IX, MNP, itype, 0, 711, comm, istatus, ierr)
CALL MPI_RECV(CF_IY, MNP, itype, 0, 712, comm, istatus, ierr)
CALL MPI_RECV(CF_COEFF, 4*MNP, rtype, 0, 713, comm, istatus, ierr)
END IF
END IF
#else
CF_IX=CF_IX_GLOGAL
CF_IY=CF_IY_GLOBAL
CF_COEFF=CF_COEFF_GLOBAL
deallocate(CF_IX_GLOBAL, CF_IY_GLOBAL, CF_COEFF_GLOBAL)
#endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE SAVE_INTERP_ARRAY(FileSave)
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
character(len=256), intent(in) :: FileSave
character (len = *), parameter :: CallFct = "SAVE_INTERP_ARRAY"
integer, allocatable :: CF_IX_GLOBAL(:), CF_IY_GLOBAL(:)
real(rkind), allocatable :: CF_COEFF_GLOBAL(:,:)
integer, allocatable :: CF_IX_loc(:), CF_IY_loc(:)
real(rkind), allocatable :: CF_COEFF_loc(:,:)
integer, allocatable :: ListFirstMNP(:)
integer ncid, iret, var_id
integer mnp_dims, four_dims
integer IP, IPglob, iPROC, NP_RESloc, IPloc
IF (WRITESTATFLAG == 1) THEN
WRITE(STAT%FHNDL,*) 'minval(CF_IX)=', minval(CF_IX)
WRITE(STAT%FHNDL,*) 'minval(CF_IY)=', minval(CF_IY)
WRITE(STAT%FHNDL,*) 'minval(CF_COEFF)=', minval(CF_COEFF)
END IF
#ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .eqv. .TRUE.) THEN
IF (myrank .eq. 0) THEN
allocate(ListFirstMNP(nproc), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
ListFirstMNP=0
DO iProc=2,nproc
ListFirstMNP(iProc)=ListFirstMNP(iProc-1) + ListMNP(iProc-1)
END DO
!
allocate(CF_IX_GLOBAL(np_total), CF_IY_GLOBAL(np_total), CF_COEFF_GLOBAL(4,np_total), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
!
DO IP=1,NP_RES
IPglob=iplg(IP)
CF_IX_GLOBAL(IPglob)=CF_IX(IP)
CF_IY_GLOBAL(IPglob)=CF_IY(IP)
CF_COEFF_GLOBAL(:, IPglob)=CF_COEFF(:,IP)
END DO
DO iPROC=2,nproc
NP_RESloc=ListNP_RES(iPROC)
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) ' iPROC=', iPROC, ' NP_RES_loc=', NP_RESloc
allocate(CF_IX_loc(NP_RESloc), CF_IY_loc(NP_RESloc), CF_COEFF_loc(4,NP_RESloc), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) 'Step 1'
CALL MPI_RECV(CF_IX_loc, NP_RESloc, itype, iProc-1, 611, comm, istatus, ierr)
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) 'Step 2'
CALL MPI_RECV(CF_IY_loc, NP_RESloc, itype, iProc-1, 612, comm, istatus, ierr)
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) 'Step 3'
CALL MPI_RECV(CF_COEFF_loc, 4*NP_RESloc, rtype, iProc-1, 613, comm, istatus, ierr)
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) 'Step 4'
DO IPloc=1,NP_RESloc
IPglob=ListIPLG(IPloc + ListFirstMNP(iProc))
CF_IX_GLOBAL(IPglob)=CF_IX_loc(IPloc)
CF_IY_GLOBAL(IPglob)=CF_IY_loc(IPloc)
CF_COEFF_GLOBAL(:, IPglob)=CF_COEFF_loc(:, IPloc)
END DO
deallocate(CF_IX_loc, CF_IY_loc, CF_COEFF_loc)
END DO
deallocate(ListFirstMNP)
ELSE
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) ' NP_RES=', NP_RES
CALL MPI_SEND(CF_IX, NP_RES, itype, 0, 611, comm, ierr)
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) 'Step 1'
CALL MPI_SEND(CF_IY, NP_RES, itype, 0, 612, comm, ierr)
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) 'Step 2'
CALL MPI_SEND(CF_COEFF, 4*NP_RES, rtype, 0, 613, comm, ierr)
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL,*) 'Step 3'
END IF
ELSE
allocate(CF_IX_GLOBAL(np_total), CF_IY_GLOBAL(np_total), CF_COEFF_GLOBAL(4,np_total), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
CF_IX_GLOBAL=CF_IX
CF_IY_GLOBAL=CF_IY
CF_COEFF_GLOBAL=CF_COEFF
END IF
#else
allocate(CF_IX_GLOBAL(np_total), CF_IY_GLOBAL(np_total), CF_COEFF_GLOBAL(4,np_total), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
CF_IX_GLOBAL=CF_IX
CF_IY_GLOBAL=CF_IY
CF_COEFF_GLOBAL=CF_COEFF
#endif
!
! Now writing up
!
#ifdef MPI_PARALL_GRID
IF (myrank .eq. 0) THEN
#endif
iret=nf90_create(TRIM(FileSave), nf90_CLOBBER, ncid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, iret)
!
iret = nf90_def_dim(ncid, 'mnp', np_total, mnp_dims)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, iret)
!
iret = nf90_def_dim(ncid, 'four', 4, four_dims)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, iret)
!
iret=nf90_def_var(ncid,'CF_IX',NF90_INT,(/ mnp_dims/),var_id)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, iret)
!
iret=nf90_def_var(ncid,'CF_IY',NF90_INT,(/ mnp_dims/),var_id)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, iret)
!
iret=nf90_def_var(ncid,'CF_COEFF',NF90_RUNTYPE,(/ four_dims, mnp_dims/),var_id)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, iret)
!
iret=nf90_close(ncid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, iret)
!
! Now writing the data
!
IF (WRITESTATFLAG == 1) THEN
WRITE(STAT%FHNDL,*) 'minval(CF_IX_GLOBAL)=', minval(CF_IX_GLOBAL)
WRITE(STAT%FHNDL,*) 'minval(CF_IY_GLOBAL)=', minval(CF_IY_GLOBAL)
WRITE(STAT%FHNDL,*) 'minval(CF_COEFF_GLOBAL)=', minval(CF_COEFF_GLOBAL)
END IF
!
iret=nf90_open(TRIM(FileSave), NF90_WRITE, ncid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, iret)
!
iret=nf90_inq_varid(ncid,'CF_IX',var_id)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, iret)
!
iret=nf90_put_var(ncid,var_id,CF_IX_GLOBAL,start = (/ 1 /), count=(/np_total/))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, iret)
!
iret=nf90_inq_varid(ncid,'CF_IY',var_id)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, iret)
!
iret=nf90_put_var(ncid,var_id,CF_IY_GLOBAL,start = (/ 1 /), count=(/np_total/))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, iret)
!
iret=nf90_inq_varid(ncid,'CF_COEFF',var_id)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, iret)
!
iret=nf90_put_var(ncid,var_id,CF_COEFF_GLOBAL,start = (/ 1, 1 /), count=(/4, np_total/))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, iret)
!
iret=nf90_close(ncid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, iret)
!
deallocate(CF_IX_GLOBAL, CF_IY_GLOBAL, CF_COEFF_GLOBAL)
#ifdef MPI_PARALL_GRID
END IF
#endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE COMPUTE_SINGLE_INTERPOLATION_INFO(TheInfo, EXTRAPO_IN, eX, eY, eCF_IX, eCF_IY, eCF_COEFF, EXTRAPO_OUT)
USE DATAPOOL
IMPLICIT NONE
type(FD_FORCING_GRID), intent(in) :: TheInfo
logical, intent(in) :: EXTRAPO_IN
real(rkind), intent(in) :: eX, eY
integer, intent(out) :: eCF_IX, eCF_IY
real(rkind), intent(out) :: eCF_COEFF(4)
logical, intent(out) :: EXTRAPO_OUT
!
integer IX, IY
integer IXs, IYs
integer IXmin, IYmin, IXmax, IYmax
integer nx, ny
integer aShift
REAL(rkind) :: WI(3), X(3), Y(3), a, b
REAL(rkind) :: MinDist, eDist
nx = TheInfo % nx_dim
ny = TheInfo % ny_dim
MinDist=LARGE
! WRITE(WINDBG%FHNDL,*) 'Start finding closest coordinate'
! FLUSH(WINDBG%FHNDL)
DO IX=1,nx-1
DO IY=1,ny-1
eDist=(eX-TheInfo % LON(IX,IY))**2 + (eY-TheInfo % LAT(IX,IY))**2
IF (eDist .lt. MinDist) THEN
MinDist=eDist
IXs=IX
IYs=IY
END IF
END DO
END DO
aShift=1
! WRITE(WINDBG%FHNDL,*) 'lon(eX)=', eX
! WRITE(WINDBG%FHNDL,*) 'lat(eY)=', eY
! WRITE(WINDBG%FHNDL,*) 'LON(IXs+-1,IYs)=', TheInfo % LON(IXs-1,IYs), TheInfo % LON(IXs,IYs),TheInfo % LON(IXs+1,IYs)
! WRITE(WINDBG%FHNDL,*) 'LAT(IXs,IYs+-1)=', TheInfo % LAT(IXs,IYs-1), TheInfo % LAT(IXs,IYs),TheInfo % LAT(IXs,IYs+1)
! WRITE(WINDBG%FHNDL,*) 'MinDist(km)=', SQRT(MinDist)*110
! FLUSH(WINDBG%FHNDL)
EXTRAPO_OUT=.TRUE.
DO
IXmin=max(1, IXs - aShift)
IYmin=max(1, IYs - aShift)
IXmax=min(nx-1, IXs+aShift)
IYmax=min(ny-1, IYs+aShift)
DO IX=IXmin,IXmax
DO IY=IYmin,IYmax
!
! First triangle
!
X(1)=TheInfo % LON(IX, IY)
X(2)=TheInfo % LON(IX+1, IY)
X(3)=TheInfo % LON(IX, IY+1)
Y(1)=TheInfo % LAT(IX, IY)
Y(2)=TheInfo % LAT(IX+1, IY)
Y(3)=TheInfo % LAT(IX, IY+1)
CALL INTELEMENT_COEF(X,Y,eX,eY,WI)
IF (minval(WI) .ge. -THR) THEN
EXTRAPO_OUT=.FALSE.
eCF_IX=IX
eCF_IY=IY
a=WI(2)
b=WI(3)
eCF_COEFF(1)=(1-a)*(1-b)
eCF_COEFF(2)=a*(1-b)
eCF_COEFF(3)=(1-a)*b
eCF_COEFF(4)=a*b
RETURN
END IF
!
! Second triangle
!
X(1)=TheInfo % LON(IX+1, IY+1)
X(2)=TheInfo % LON(IX+1, IY)
X(3)=TheInfo % LON(IX, IY+1)
Y(1)=TheInfo % LAT(IX+1, IY+1)
Y(2)=TheInfo % LAT(IX+1, IY)
Y(3)=TheInfo % LAT(IX, IY+1)
CALL INTELEMENT_COEF(X,Y,eX,eY,WI)
IF (minval(WI) .ge. -THR) THEN
EXTRAPO_OUT=.FALSE.
eCF_IX=IX
eCF_IY=IY
a=1 - WI(3)
b=1 - WI(2)
eCF_COEFF(1)=(1-a)*(1-b)
eCF_COEFF(2)=a*(1-b)
eCF_COEFF(3)=(1-a)*b
eCF_COEFF(4)=a*b
END IF
END DO
END DO
IF ((IXmin .eq. 1).and.(IYmin .eq. 1).and.(IXmax .eq. nx-1).and.(IYmax .eq. ny-1)) THEN
EXIT
END IF
aShift=aShift + 1
END DO
IF ( (EXTRAPO_IN .eqv. .FALSE.) .and. (EXTRAPO_OUT .eqv. .TRUE.) ) THEN
WRITE(STAT % FHNDL,*) 'aShift=', aShift
WRITE(STAT % FHNDL,*) 'eX=', eX, 'eY=', eY
FLUSH(STAT % FHNDL)
CALL WWM_ABORT('We find a model point outside of the available forcing grid')
ELSE
eCF_IX = IXs
eCF_IY = IYs
eCF_COEFF(1)=1
eCF_COEFF(2)=0
eCF_COEFF(3)=0
eCF_COEFF(4)=0
! EXTRAPO_OUT=.TRUE.
! WRITE(STAT % FHNDL,*) 'Point ', eX, '/', eY, ' outside grid'
! WRITE(STAT % FHNDL,*) 'MinDist=', MinDist
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE COMPUTE_CF_COEFFICIENTS(TheInfo)
USE DATAPOOL
IMPLICIT NONE
type(FD_FORCING_GRID), intent(in) :: TheInfo
integer I, IX, IY, IXs, IYs, IXmin, IYmin, IXmax, IYmax
integer aShift, WeFind
real(rkind) eDist, MinDist
real(rkind), allocatable :: dist(:,:)
real(rkind) closest_r(2)
integer closest(2)
real(rkind) d_lon, d_lat
integer i11, j11, i12, j12, i21, j21
integer eCF_IX, eCF_IY
real(rkind) eCF_COEFF(4)
integer :: nbExtrapolation = 0
real(rkind) :: MaxMinDist = 0
character(len=256) :: FileSave = "wwm_filesave_interp_array.nc"
logical success
logical EXTRAPO_OUT
real(rkind) eX, eY
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'Starting node loop for calcs of coefs'
allocate(CF_IX(MNP_WIND), CF_IY(MNP_WIND), SHIFTXY(4,2), CF_COEFF(4,MNP_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 52')
SHIFTXY(1,1)=0
SHIFTXY(1,2)=0
SHIFTXY(2,1)=1
SHIFTXY(2,2)=0
SHIFTXY(3,1)=0
SHIFTXY(3,2)=1
SHIFTXY(4,1)=1
SHIFTXY(4,2)=1
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'LSAVE_INTERP_ARRAY=', LSAVE_INTERP_ARRAY
IF (LSAVE_INTERP_ARRAY) THEN
CALL LOAD_INTERP_ARRAY(FileSave, success)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'success=', success
IF (success .eqv. .TRUE.) RETURN
END IF
CF_IX=0
CF_IY=0
CF_COEFF=0
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'min(lon)=', minval(TheInfo % LON)
WRITE(WINDBG%FHNDL,*) 'max(lon)=', maxval(TheInfo % LON)
WRITE(WINDBG%FHNDL,*) 'min(lat)=', minval(TheInfo % LAT)
WRITE(WINDBG%FHNDL,*) 'max(lat)=', maxval(TheInfo % LAT)
WRITE(WINDBG%FHNDL,*) 'MNP_WIND=', MNP_WIND
WRITE(WINDBG%FHNDL,*) 'START BIG SLOW LOOP'
FLUSH(WINDBG%FHNDL)
END IF
DO I = 1, MNP_WIND
IF (I .eq. 1) THEN
IXs=1
IYs=1
ELSE
IXs=CF_IX(I-1)
IYs=CF_IX(I-1)
END IF
eX=XP_WIND(I)
eY=YP_WIND(I)
CALL COMPUTE_SINGLE_INTERPOLATION_INFO(TheInfo, EXTRAPOLATION_ALLOWED_WIND, eX, eY, eCF_IX, eCF_IY, eCF_COEFF, EXTRAPO_OUT)
CF_IX(I) = eCF_IX
CF_IY(I) = eCF_IY
CF_COEFF(:,I) = eCF_COEFF
IF (EXTRAPO_OUT .eqv. .TRUE.) THEN
nbExtrapolation = nbExtrapolation + 1
END IF
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'I(', MNP_WIND,')=', I
FLUSH(WINDBG%FHNDL)
END IF
END DO
IF (LSAVE_INTERP_ARRAY) THEN
CALL SAVE_INTERP_ARRAY(FileSave)
END IF
IF (EXTRAPOLATION_ALLOWED_WIND .eqv. .TRUE. .AND. WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) ' nbExtrapolation=', nbExtrapolation
WRITE(WINDBG%FHNDL,*) ' MaxMinDist=', sqrt(MaxMinDist)
END IF
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) ' done interp calcs'
FLUSH(WINDBG%FHNDL)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE GET_CF_TIME_INDEX(eVAR, REC1, REC2, w1, w2)
! For given wwm_time and wind_time return records to get and weights for time
! interpolation F(wwm_time)=F(rec1)*w1 + F(rec2)*w2
!
USE DATAPOOL
IMPLICIT NONE
TYPE(VAR_NETCDF_CF), intent(in) :: eVAR
REAL(rkind), INTENT(OUT) :: w1, w2
INTEGER, INTENT(OUT) :: REC1, REC2
REAL(rkind) :: eTime1, eTime2
INTEGER :: iTime
DO iTime=2,eVAR % nbTime
eTime1=eVAR % ListTime(iTime-1)
eTime2=eVAR % ListTime(iTime)
IF ((eTime1 .le. MAIN%TMJD).and.(MAIN%TMJD .le. eTime2)) THEN
REC2=iTime
REC1=iTime-1
w2=(MAIN % TMJD - eTime1)/(eTime2-eTime1)
w1=(eTime2 - MAIN % TMJD)/(eTime2-eTime1)
RETURN
END IF
END DO
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'Time error in wind for CF'
WRITE(WINDBG%FHNDL,*) 'MAIN % TMJD=', MAIN%TMJD
WRITE(WINDBG%FHNDL,*) 'min(wind_time_mjd)=', minval(eVAR % ListTime)
WRITE(WINDBG%FHNDL,*) 'max(wind_time_mjd)=', maxval(eVAR % ListTime)
FLUSH(WINDBG%FHNDL)
END IF
CALL WWM_ABORT('Error in CF wind forcing time setup')
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE SYNCHRONIZE_WIND_TIME_IFILE_IT
USE DATAPOOL
IMPLICIT NONE
integer IPROC, eInt(1)
#ifdef MPI_PARALL_GRID
IF (.NOT. MULTIPLE_IN_WIND) THEN
IF (myrank .eq. 0) THEN
eInt(1)=NDT_WIND_ALL_FILES
DO IPROC=2,nproc
CALL MPI_SEND(eInt,1,itype, iProc-1, 811, comm, ierr)
END DO
DO IPROC=2,nproc
CALL MPI_SEND(WIND_TIME_ALL_FILES,NDT_WIND_ALL_FILES,rtype, iProc-1, 812, comm, ierr)
CALL MPI_SEND(WIND_TIME_IFILE,NDT_WIND_ALL_FILES,itype, iProc-1, 813, comm, ierr)
CALL MPI_SEND(WIND_TIME_IT,NDT_WIND_ALL_FILES,itype, iProc-1, 814, comm, ierr)
END DO
ELSE
CALL MPI_RECV(eInt,1,itype, 0, 811, comm, istatus, ierr)
NDT_WIND_ALL_FILES=eInt(1)
ALLOCATE(WIND_TIME_ALL_FILES(NDT_WIND_ALL_FILES), WIND_TIME_IFILE(NDT_WIND_ALL_FILES), WIND_TIME_IT(NDT_WIND_ALL_FILES), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 3')
CALL MPI_RECV(WIND_TIME_ALL_FILES,NDT_WIND_ALL_FILES,rtype, 0, 812, comm, istatus, ierr)
CALL MPI_RECV(WIND_TIME_IFILE,NDT_WIND_ALL_FILES,rtype, 0, 813, comm, istatus, ierr)
CALL MPI_RECV(WIND_TIME_IT,NDT_WIND_ALL_FILES,rtype, 0, 814, comm, istatus, ierr)
END IF
END IF
#endif
SEWI%DELT = ( WIND_TIME_ALL_FILES(2) - WIND_TIME_ALL_FILES(1) ) * DAY2SEC
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
#ifdef NCDF
SUBROUTINE INIT_NETCDF_DWD
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
INTEGER :: IT, IFILE
INTEGER :: ILON_ID, ILAT_ID, ITIME_ID, I, J, COUNTER
REAL(rkind) :: DTMP
REAL(rkind), ALLOCATABLE :: WIND_TIME(:)
character ( len = 20 ) chrtmp
character ( len = 15 ) chrdate
integer, dimension(nf90_max_var_dims) :: dimIDs
character (len = *), parameter :: CallFct="INIT_NETCDF_DWD"
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
OPEN(WIN%FHNDL,FILE=WIN%FNAME,STATUS='OLD')
!
! count number of netcdf files in list ...
!
NUM_NETCDF_FILES = 0
DO
READ( WIN%FHNDL, *, IOSTAT = ISTAT )
IF ( ISTAT /= 0 ) EXIT
NUM_NETCDF_FILES = NUM_NETCDF_FILES + 1
END DO
REWIND (WIN%FHNDL)
ALLOCATE(NETCDF_FILE_NAMES(NUM_NETCDF_FILES), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 3')
DO IT = 1, NUM_NETCDF_FILES
READ( WIN%FHNDL, *) NETCDF_FILE_NAMES(IT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) IT, NETCDF_FILE_NAMES(IT)
END DO
CLOSE (WIN%FHNDL)
!
! check number of time steps in netcdf file ... it is assumed that all files have the same ammount of time steps ...
!
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(1)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(1), NF90_NOWRITE, WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, 'initial_time0_encoded', ITIME_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, ITIME_ID, dimids = dimids)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = NDT_WIND_FILE)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
!
! check dimensions in the netcdf ... again it is assumed that this is not changing for all files ...
!
ISTAT = nf90_inq_varid(WIND_NCID, 'g0_lon_2', ILON_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, ILON_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = NDX_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, 'g0_lat_1', ILAT_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, ILAT_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = NDY_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ALLOCATE (COORD_WIND_X(NDX_WIND), COORD_WIND_Y(NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 4')
!
! read coordinates from files ....
!
ISTAT = NF90_GET_VAR(WIND_NCID, ILON_ID, COORD_WIND_X)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, ILAT_ID, COORD_WIND_Y)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT)
!
! estimate offset ...
!
OFFSET_X_WIND = MINVAL(COORD_WIND_X)
OFFSET_Y_WIND = MINVAL(COORD_WIND_Y)
!
! resolution ...
!
DX_WIND = ABS(MAXVAL(COORD_WIND_X)-MINVAL(COORD_WIND_X))/(NDX_WIND-1)
DY_WIND = ABS(MAXVAL(COORD_WIND_Y)-MINVAL(COORD_WIND_Y))/(NDY_WIND-1)
!
! close netcdf file ...
!
ISTAT = NF90_CLOSE(WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, ISTAT)
!
! total number of time steps ... in all files
!
NDT_WIND_ALL_FILES = NDT_WIND_FILE * NUM_NETCDF_FILES
ALLOCATE (WIND_TIME(NDT_WIND_FILE), WIND_TIME_ALL_FILES(NDT_WIND_ALL_FILES), WIND_TIME_IFILE(NDT_WIND_ALL_FILES), WIND_TIME_IT(NDT_WIND_ALL_FILES), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 5')
!
! read all time steps in the proper format and transform in wwm time line
!
DO IFILE = 1, NUM_NETCDF_FILES
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(IFILE)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, ITIME_ID, WIND_TIME)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, ISTAT)
DO IT = 1, NDT_WIND_FILE
WRITE (3001, *) WIND_TIME(IT)
END DO
REWIND(3001)
DO IT = 1, NDT_WIND_FILE
READ (3001, '(A20)') CHRTMP
! YYYYMMDD
CHRDATE(1:8) = CHRTMP(4:12)
CHRDATE(9:9) = '.'
! HH
CHRDATE(10:11) = CHRTMP(12:13)
! MMSS
CHRDATE(12:15) = '0000' ! Construct propper format YYYYMMDDHHMMSS character .... len = 15 ... :)
CALL CT2MJD(CHRDATE,DTMP) !
WIND_TIME(IT) = DTMP ! Double time with respect to 19000101.000000
END DO
CLOSE(3001)
DO IT = 1, NDT_WIND_FILE
WIND_TIME_ALL_FILES(IT+(IFILE-1)*NDT_WIND_FILE) = WIND_TIME(IT)
WIND_TIME_IFILE(IT+(IFILE-1)*NDT_WIND_FILE) = IFILE
WIND_TIME_IT (IT+(IFILE-1)*NDT_WIND_FILE) = IT
END DO
END DO ! IFILE
IF (LWRITE_ORIG_WIND) THEN
WRITE (3010, '(I10)') 0
WRITE (3010, '(I10)') NDX_WIND * NDY_WIND
COUNTER = 0
DO I = 1, NDY_WIND
DO J = 1, NDX_WIND
WRITE (3010, '(I10,3F15.4)') COUNTER, OFFSET_X_WIND+(J-1)*DX_WIND ,OFFSET_Y_WIND+(I-1)*DY_WIND , 0.0
COUNTER = COUNTER + 1
END DO
END DO
WRITE (3010, *) (NDX_WIND-1)*(NDY_WIND-1)*2
DO J = 0, NDY_WIND-2
DO I = 0, NDX_WIND-2
WRITE (3010, '(5I10)') I+J*NDX_WIND , NDX_WIND+I+J* NDX_WIND, NDX_WIND+I+1+J*NDX_WIND, 0, 0
WRITE (3010, '(5I10)') NDX_WIND+I+1+J*NDX_WIND, I+1+J*NDX_WIND , I+J*NDX_WIND , 0, 0
END DO
END DO
OPEN(3011, FILE = 'ergwindorig.bin', FORM = 'UNFORMATTED')
END IF
ALLOCATE (WIND_X(NDX_WIND,NDY_WIND), WIND_Y(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 9')
# ifdef MPI_PARALL_GRID
END IF
# endif
CALL SYNCHRONIZE_WIND_TIME_IFILE_IT
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE READ_NETCDF_DWD(IFILE, IT, eField)
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
!
! READS WIND_Y, WIND_X and PRESSURE from a given NCID within one DWD file
!
INTEGER, INTENT(IN) :: IFILE, IT
REAL(rkind), intent(inout) :: eField(MNP,2)
REAL(rkind) :: Vtotal1(MNP_WIND)
REAL(rkind) :: Vtotal2(MNP_WIND)
REAL(rkind) :: Vlocal(MNP)
character (len = *), parameter :: CallFct="READ_NETCDF_DWD"
INTEGER :: DWIND_X_ID, DWIND_Y_ID
INTEGER :: numLons, numLats, numTime, iy, counter, ip, i, j
REAL(rkind), ALLOCATABLE :: TMP(:,:)
REAL(rkind), ALLOCATABLE :: U(:), V(:), H(:)
REAL(rkind), SAVE :: TIME
INTEGER, DIMENSION (nf90_max_var_dims) :: dimIDs
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
IF (IFILE .GT. NUM_NETCDF_FILES) CALL WWM_ABORT('SOMETHING IS WRONG WE RUN OUT OF WIND TIME')
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(IFILE)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, 'V_GDS0_HTGL_13', DWIND_X_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, DWIND_X_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = numLons)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(2), len = numLats)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(3), len = numTime)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, 'U_GDS0_HTGL_13', DWIND_Y_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, DWIND_Y_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = numLons)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(2), len = numLats)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(3), len = numTime)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, DWIND_X_ID, WIND_X, start = (/ 1, 1, IT /), count = (/ NDX_WIND, NDY_WIND, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, DWIND_Y_ID, WIND_Y, start = (/ 1, 1, IT /), count = (/ NDX_WIND, NDY_WIND, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, ISTAT)
IF (LINVERTY) THEN
ALLOCATE(TMP(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 11')
DO IY = 1, NDY_WIND
TMP(:,NDY_WIND-(IY-1)) = wind_x(:,IY)
END DO
wind_x = TMP
DO IY = 1, NDY_WIND
TMP(:,NDY_WIND-(IY-1)) = wind_y(:,IY)
END DO
wind_y = TMP
DO IY = 1, NDY_WIND
TMP(:,NDY_WIND-(IY-1)) = atmo_press(:,IY)
END DO
atmo_press = TMP
DEALLOCATE(TMP)
END IF
IF (LWRITE_ORIG_WIND) THEN
ALLOCATE(U(NDX_WIND*NDY_WIND), V(NDX_WIND*NDY_WIND), H(NDX_WIND*NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 15')
COUNTER = 1
DO J = 1, NDY_WIND
DO I = 1, NDX_WIND
U(COUNTER) = WIND_X(I,J)
V(COUNTER) = WIND_Y(I,J)
IF (ABS(U(COUNTER)) .GT. 1000.) U(COUNTER) = 0.
IF (ABS(V(COUNTER)) .GT. 1000.) V(COUNTER) = 0.
H(COUNTER) = SQRT((U(COUNTER)**2.+V(COUNTER)**2.))
COUNTER = COUNTER + 1
END DO
END DO
TIME = TIME + 1.
WRITE(3011) TIME
WRITE(3011) (U(IP), V(IP), H(IP), IP = 1, numLons*numLats)
DEALLOCATE(U,V,H)
END IF
ISTAT = NF90_CLOSE(WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, ISTAT)
CALL INTER_STRUCT_DATA(NDX_WIND,NDY_WIND,DX_WIND,DY_WIND,OFFSET_X_WIND,OFFSET_Y_WIND,WIND_X,Vtotal1)
CALL INTER_STRUCT_DATA(NDX_WIND,NDY_WIND,DX_WIND,DY_WIND,OFFSET_X_WIND,OFFSET_Y_WIND,WIND_Y,Vtotal2)
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
eField(:,1)=Vtotal1
eField(:,2)=Vtotal2
ELSE
CALL SCATTER_ONED_ARRAY(Vtotal1, Vlocal)
eField(:,1)=Vlocal
CALL SCATTER_ONED_ARRAY(Vtotal2, Vlocal)
eField(:,2)=Vlocal
END IF
# else
eField(:,1)=Vtotal1
eField(:,2)=Vtotal2
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_NETCDF_CRFS
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
!for data description consult ftp://nomads.ncdc.noaa.gov/CFSR/HP_time_series/200307/wnd10m.l.gdas.200307.grb2.inv
INTEGER :: IT, IFILE
INTEGER :: ILON_ID, ILAT_ID, ITIME_ID, I, J, COUNTER
REAL(rkind) :: START_TIME
REAL(rkind) , ALLOCATABLE :: WIND_TIME(:), WIND_TIME_NETCDF(:)
character ( len = 15 ) chrdate
character ( len = 40 ) beginn_time
CHARACTER(LEN=15) :: eTimeStr
character (len = *), parameter :: CallFct="INIT_NETCDF_CRFS"
integer, dimension(nf90_max_var_dims) :: dimIDs
logical PREF_ANALYZED
integer idx
REAL(rkind) :: ePresTime, eNewTime
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'Begin INIT_NETCDF_CRFS'
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
OPEN(WIN%FHNDL,FILE=WIN%FNAME,STATUS='OLD',IOSTAT = ISTAT)
!
! count number of netcdf files in list ...
! CRFS has analyzed at time 0 and forecast at times +1, +2, +3, +4, +5 +6
! Ordering in the file is
! --analyzed 0
! --fcst +1
! --fcst +2
! --fcst +3
! --fcst +4
! --fcst +5
! --fcst +6
! So, it goes by blocks of 7.
! PREF_ANALYZED = TRUE. For prefering analyzed fields when available
! FALSE For using the analyzed field only at first step.
PREF_ANALYZED=.FALSE.
NUM_NETCDF_FILES = 0
DO
READ( WIN%FHNDL, *, IOSTAT = ISTAT )
IF ( ISTAT /= 0 ) EXIT
NUM_NETCDF_FILES = NUM_NETCDF_FILES + 1
END DO
REWIND (WIN%FHNDL)
ALLOCATE(NETCDF_FILE_NAMES(NUM_NETCDF_FILES), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 18')
DO IT = 1, NUM_NETCDF_FILES
READ( WIN%FHNDL, *) NETCDF_FILE_NAMES(IT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) IT, NETCDF_FILE_NAMES(IT)
END DO
CLOSE (WIN%FHNDL)
!
! check number of time steps in netcdf file ... it is assumed that all files have the same ammount of time steps ...
!
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(1)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(1), NF90_NOWRITE, WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, 'time', ITIME_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, ITIME_ID, dimids = dimids)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = NDT_WIND_FILE)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
ISTAT = nf90_get_att(WIND_NCID, ITIME_ID, 'units', beginn_time)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
CHRDATE(1:4) = beginn_time(13:17)
CHRDATE(5:6) = beginn_time(18:19)
CHRDATE(7:8) = beginn_time(21:22)
CHRDATE(9:9) = '.'
CHRDATE(10:15)= '000000'
CALL CT2MJD(CHRDATE,START_TIME)
!
! check dimensions in the netcdf ... again it is assumed that this is not changing for all files ...
!
ISTAT = nf90_inq_varid(WIND_NCID, 'lon', ILON_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, ILON_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = NDX_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, 'lat', ILAT_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, ILAT_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = NDY_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
ALLOCATE (DCOORD_WIND_X(NDX_WIND), DCOORD_WIND_Y(NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 19')
!
! read cooridantes from files ....
!
ISTAT = NF90_GET_VAR(WIND_NCID, ILON_ID, DCOORD_WIND_X)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, ILAT_ID, DCOORD_WIND_Y)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, ISTAT)
DCOORD_WIND_Y = DCOORD_WIND_Y !+ 90.0_rkind
DO I = 1, NDX_WIND
DCOORD_WIND_X(I) = DCOORD_WIND_X(I) - 180.0_rkind
END DO
!
! estimate offset ...
!
OFFSET_X_WIND = MINVAL(DCOORD_WIND_X)
OFFSET_Y_WIND = MINVAL(DCOORD_WIND_Y)
!
! resolution ...
!
DX_WIND = ABS(MAXVAL(DCOORD_WIND_X)-MINVAL(DCOORD_WIND_X))/(NDX_WIND-1)
DY_WIND = ABS(MAXVAL(DCOORD_WIND_Y)-MINVAL(DCOORD_WIND_Y))/(NDY_WIND-1)
!
! close netcdf file ...
!
ISTAT = NF90_CLOSE(WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, ISTAT)
!
! total number of time steps ... in all files
!
ALLOCATE (WIND_X(NDX_WIND,NDY_WIND), WIND_Y(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 35')
ALLOCATE (WIND_TIME(NDT_WIND_FILE),WIND_TIME_NETCDF(NDT_WIND_FILE), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 20')
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'NDT_WIND_FILE=', NDT_WIND_FILE
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'START_TIME=', START_TIME
!
! read all time steps in the proper format and transform in wwm time line
!
NDT_WIND_ALL_FILES=0
ePresTime=-100000000
DO IFILE = 1, NUM_NETCDF_FILES
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(IFILE)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, ITIME_ID, WIND_TIME_NETCDF)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 16, ISTAT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'IFILE=', IFILE
DO IT = 1, NDT_WIND_FILE
eNewTime=START_TIME+WIND_TIME_NETCDF(IT)*3600. * SEC2DAY
IF (eNewTime .gt. ePresTime + THR8) THEN
NDT_WIND_ALL_FILES=NDT_WIND_ALL_FILES + 1
ePresTime=eNewTime
END IF
END DO
END DO
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'NDT_WIND_ALL_FILES=', NDT_WIND_ALL_FILES
ALLOCATE (WIND_TIME_ALL_FILES(NDT_WIND_ALL_FILES), WIND_TIME_IFILE(NDT_WIND_ALL_FILES), WIND_TIME_IT(NDT_WIND_ALL_FILES), &
& stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 21')
idx=0
ePresTime=-100000000
DO IFILE = 1, NUM_NETCDF_FILES
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(IFILE)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 17, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, ITIME_ID, WIND_TIME_NETCDF)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 18, ISTAT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'IFILE=', IFILE
DO IT = 1, NDT_WIND_FILE
eNewTime=START_TIME+WIND_TIME_NETCDF(IT)*3600. * SEC2DAY
CALL MJD2CT(eNewTime,eTimeStr)
IF (PREF_ANALYZED) THEN
IF (eNewTime .gt. ePresTime + THR8) THEN
ePresTime=eNewTime
idx=idx+1
WIND_TIME_ALL_FILES(idx) = eNewTime
WIND_TIME_IFILE(idx) = IFILE
WIND_TIME_IT(idx) = IT
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,110) IT, idx, WIND_TIME_NETCDF(IT) * 3600. * SEC2DAY, &
& eNewTime, eTimeStr
ENDIF
ELSE
IF (eNewTime .gt. ePresTime + THR8) THEN
ePresTime=eNewTime
idx=idx+1
ENDIF
WIND_TIME_ALL_FILES(idx) = eNewTime
WIND_TIME_IFILE(idx) = IFILE
WIND_TIME_IT(idx) = IT
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,110) IT, idx, WIND_TIME_NETCDF(IT) * 3600. * SEC2DAY, eNewTime, eTimeStr
ENDIF
END DO
END DO ! IFILE
110 FORMAT (I4, ' ', I4, ' ', F15.3, ' ', F15.3, ' ', a15)
IF (LWRITE_ORIG_WIND) THEN
WRITE (3010, '(I10)') 0
WRITE (3010, '(I10)') NDX_WIND * NDY_WIND
COUNTER = 0
DO I = 1, NDY_WIND
DO J = 1, NDX_WIND
WRITE (3010, '(I10,3F15.4)') COUNTER, OFFSET_X_WIND+(J-1)*DX_WIND ,OFFSET_Y_WIND+(I-1)*DY_WIND , 0.0
COUNTER = COUNTER + 1
END DO
END DO
WRITE (3010, *) (NDX_WIND-1)*(NDY_WIND-1)*2
DO J = 0, NDY_WIND-2
DO I = 0, NDX_WIND-2
WRITE (3010, '(5I10)') I+J*NDX_WIND , NDX_WIND+I+J* NDX_WIND, NDX_WIND+I+1+J*NDX_WIND, 0, 0
WRITE (3010, '(5I10)') NDX_WIND+I+1+J*NDX_WIND, I+1+J*NDX_WIND , I+J*NDX_WIND , 0, 0
END DO
END DO
OPEN(3011, FILE = 'ergwindorig.bin', FORM = 'UNFORMATTED')
END IF
# ifdef MPI_PARALL_GRID
END IF
# endif
CALL SYNCHRONIZE_WIND_TIME_IFILE_IT
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'End INIT_NETCDF_CRFS'
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE INIT_NETCDF_NARR
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
!
!2do update ...
!for data description consult ftp://nomads.ncdc.noaa.gov/CFSR/HP_time_series/200307/wnd10m.l.gdas.200307.grb2.inv
!
INTEGER :: IT, IFILE, II, IP
INTEGER :: ILON_ID, ILAT_ID, ITIME_ID, I, J, COUNTER
REAL(rkind) :: START_TIME, OFFSET_TIME
REAL(rkind), ALLOCATABLE :: WIND_TIME_NETCDF(:)
character ( len = 4 ) ch4
character ( len = 15 ) chrdate
character (len = *), parameter :: CallFct="INIT_NETCDF_NARR"
integer, dimension(nf90_max_var_dims) :: dimIDs
integer, allocatable :: COUNTERMAT(:,:)
integer, allocatable :: IMAT(:), JMAT(:)
integer :: NbPoint, nbFail
real(rkind) :: Wi(3), XPW(3), YPW(3)
INTEGER NI(3)
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'Begin INIT_NETCDF_NARR'
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'MULTIPLE_IN_WIND=', MULTIPLE_IN_WIND
!
! I make the assumption that the year when the dataset beginns at the year indicated in the bouc section
!
CHRDATE(1:4) = SEBO%BEGT(1:4)
CHRDATE(5:6) = '01'
CHRDATE(7:8) = '01'
CHRDATE(9:9) = '.'
CHRDATE(10:15)= '000000'
CALL CT2MJD(CHRDATE,START_TIME)
OPEN(WIN%FHNDL,FILE=WIN%FNAME,STATUS='OLD',IOSTAT = ISTAT)
!
! count number of netcdf files in list ...
!
NUM_NETCDF_FILES = 0
DO
READ( WIN%FHNDL, *, IOSTAT = ISTAT )
IF ( ISTAT /= 0 ) EXIT
NUM_NETCDF_FILES = NUM_NETCDF_FILES + 1
END DO
REWIND (WIN%FHNDL)
ALLOCATE(NETCDF_FILE_NAMES(NUM_NETCDF_FILES), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 22')
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'NUM_NETCDF_FILES=', NUM_NETCDF_FILES
DO IT = 1, NUM_NETCDF_FILES
READ( WIN%FHNDL, *) NETCDF_FILE_NAMES(IT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'IT=', IT, 'file=', NETCDF_FILE_NAMES(IT)
END DO
CLOSE (WIN%FHNDL)
!
! check number of time steps in netcdf file ... it is assumed that all files have the same ammount of time steps ...
!
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(1)))
ISTAT = NF90_OPEN(TRIM(NETCDF_FILE_NAMES(1)), NF90_NOWRITE, WINDX_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = nf90_inq_varid(WINDX_NCID, 'time', ITIME_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WINDX_NCID, ITIME_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = nf90_inquire_dimension(WINDX_NCID, dimIDs(1), len = NDT_WIND_FILE)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'NDT_WIND_FILE=', NDT_WIND_FILE
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(2)))
ISTAT = NF90_OPEN(TRIM(NETCDF_FILE_NAMES(2)), NF90_NOWRITE, WINDY_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
!
! check dimensions in the netcdf ... again it is assumed that this is not changing for all files ...
!
ISTAT = nf90_inq_varid(WINDX_NCID, 'lon', ILON_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WINDX_NCID, ILON_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
ISTAT = nf90_inquire_dimension(WINDX_NCID, dimIDs(1), len = NDX_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = nf90_inquire_dimension(WINDX_NCID, dimIDs(2), len = NDY_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
ISTAT = nf90_inq_varid(WINDY_NCID, 'lat', ILAT_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WINDY_NCID, ILAT_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
ISTAT = nf90_inquire_dimension(WINDY_NCID, dimIDs(1), len = NDX_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT)
ISTAT = nf90_inquire_dimension(WINDY_NCID, dimIDs(2), len = NDY_WIND)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, ISTAT)
ALLOCATE (DCOORD_WIND_X2(NDX_WIND,NDY_WIND), DCOORD_WIND_Y2(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 23')
!
! read coordinates from files ....
!
ISTAT = NF90_GET_VAR(WINDX_NCID, ILON_ID, DCOORD_WIND_X2)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, ISTAT)
ISTAT = NF90_GET_VAR(WINDY_NCID, ILAT_ID, DCOORD_WIND_Y2)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, ISTAT)
!
! estimate offset ...
!
OFFSET_X_WIND = MINVAL(DCOORD_WIND_X2)
OFFSET_Y_WIND = MINVAL(DCOORD_WIND_Y2)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'OFFSET_X_WIND=', OFFSET_X_WIND
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'OFFSET_Y_WIND=', OFFSET_Y_WIND
!
! close netcdf file ...
!
ISTAT = NF90_CLOSE(WINDX_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 16, ISTAT)
ISTAT = NF90_CLOSE(WINDY_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 17, ISTAT)
!
! total number of time steps ... in all files
!
NDT_WIND_ALL_FILES = NDT_WIND_FILE * NUM_NETCDF_FILES/2
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'NDT_WIND_ALL_FILES=', NDT_WIND_ALL_FILES
ALLOCATE (WIND_TIME_NETCDF(NDT_WIND_FILE), WIND_TIME_ALL_FILES(NDT_WIND_ALL_FILES), WIND_TIME_IFILE(NDT_WIND_ALL_FILES), WIND_TIME_IT(NDT_WIND_ALL_FILES), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 24')
!
! read all time steps in the proper format and transform in wwm time line
!
DO IFILE = 1, NUM_NETCDF_FILES, 2
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(IFILE)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, WINDX_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 18, ISTAT)
ISTAT = NF90_GET_VAR(WINDX_NCID, ITIME_ID, WIND_TIME_NETCDF)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 19, ISTAT)
ISTAT = NF90_CLOSE(WINDX_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 20, ISTAT)
WIND_TIME_NETCDF = WIND_TIME_NETCDF/24 ! Transform to days ...
OFFSET_TIME = MINVAL(WIND_TIME_NETCDF) - START_TIME ! in this dataset the time start at 1800 in WWM it start at 1900
WIND_TIME_NETCDF = WIND_TIME_NETCDF - OFFSET_TIME ! Now in WWM timeline ...
DO IT = 1, NDT_WIND_FILE
WIND_TIME_ALL_FILES(IT+(IFILE-1)*NDT_WIND_FILE) = WIND_TIME_NETCDF(IT)
WIND_TIME_IFILE(IT+(IFILE-1)*NDT_WIND_FILE) = IFILE
WIND_TIME_IT(IT+(IFILE-1)*NDT_WIND_FILE) = IT
END DO
CH4 = CHRDATE(1:4)
READ (CH4,'(i4)') I
I = I + 1
WRITE (11111,'(i4)') I
REWIND(11111)
READ(11111,*) ch4
CHRDATE(1:4) = ch4
CHRDATE(5:6) = '01'
CHRDATE(7:8) = '01'
CHRDATE(9:9) = '.'
CHRDATE(10:15)= '000000'
CALL CT2MJD(CHRDATE,START_TIME)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) CHRDATE, START_TIME
END DO ! IFILE
!
! Now the geographic interpolation
!
NE_WIND = (NDX_WIND-1)*(NDY_WIND-1)*2
NP_WIND = NDX_WIND*NDY_WIND
IF (LWRITE_ORIG_WIND) THEN
WRITE (3010, '(I10)') 0
WRITE (3010, '(I10)') NP_WIND
END IF
COUNTER = 0
NbPoint=NDX_WIND*NDY_WIND
ALLOCATE(XYPWIND(2,NbPoint), IMAT(NbPoint), JMAT(NbPoint), COUNTERMAT(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 28')
DO I = 1, NDX_WIND
DO J = 1, NDY_WIND
IF (DCOORD_WIND_X2(I,J) .GT. 0.) THEN
! Transformed to a unified domain extending below -180.
DCOORD_WIND_X2(I,J) = -1 * DCOORD_WIND_X2(I,J) - (180.-DCOORD_WIND_X2(I,J)) * 2
END IF
COUNTER = COUNTER + 1
XYPWIND(1,COUNTER) = DCOORD_WIND_X2(I,J)
XYPWIND(2,COUNTER) = DCOORD_WIND_Y2(I,J)
IMAT(COUNTER)=I
JMAT(COUNTER)=J
COUNTERMAT(I,J)=COUNTER
IF (LWRITE_ORIG_WIND) WRITE (3010, '(I10,3F15.4)') COUNTER, XYPWIND(1,COUNTER), XYPWIND(2,COUNTER), 0.0
END DO
END DO
IF (LWRITE_ORIG_WIND) WRITE (3010, *) NE_WIND
NE_WIND = (NDX_WIND-1)*(NDY_WIND-1)*2
ALLOCATE(INE_WIND(3,NE_WIND), UWND_NARR(NDX_WIND*NDY_WIND), VWND_NARR(NDX_WIND*NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 30')
INE_WIND = 0
UWND_NARR = 0.
VWND_NARR = 0.
II = 0
DO I = 1, NDX_WIND-1
DO J = 1, NDY_WIND-1
II = II + 1
INE_WIND(1,II) = COUNTERMAT(I,J)
INE_WIND(2,II) = COUNTERMAT(I+1,J)
INE_WIND(3,II) = COUNTERMAT(I,J+1)
II = II + 1
INE_WIND(1,II) = COUNTERMAT(I+1,J+1)
INE_WIND(2,II) = COUNTERMAT(I,J+1)
INE_WIND(3,II) = COUNTERMAT(I+1,J)
END DO
END DO
IF (LWRITE_ORIG_WIND) OPEN(3011, FILE = 'ergwindorig.bin', FORM = 'UNFORMATTED')
ALLOCATE(WIND_ELE(MNP_WIND), WI_NARR(MNP_WIND, 3), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 32')
WIND_ELE = 0
WI_NARR = 0
nbFail=0
DO IP = 1, MNP_WIND
CALL FIND_ELE_WIND( NE_WIND, NP_WIND, INE_WIND, XYPWIND, XP_WIND(IP), YP_WIND(IP), WIND_ELE(IP))
IF (WIND_ELE(IP) .eq. 0) THEN
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'POINT OF THE MESH IS OUT OF THE WIND FIELD', IP, XP(IP), YP(IP)
nbFail=nbFail+1
ELSE
NI=INE_WIND(:,WIND_ELE(IP))
XPW=XYPWIND(1,NI)
YPW=XYPWIND(2,NI)
CALL INTELEMENT_COEF(XPW, YPW,XP_WIND(IP),YP_WIND(IP),Wi)
WI_NARR(IP,:)=Wi
! WRITE(WINDBG%FHNDL,*) 'IP=', MNP, ' sumWi=', sum(Wi)
! WRITE(WINDBG%FHNDL,*) 'IP=', MNP, ' minW=', minval(Wi), ' maxW=', maxval(Wi)
ENDIF
END DO
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'MNP_WIND=', MNP_WIND, ' nbFail=', nbFail
WRITE(WINDBG%FHNDL,*) 'NDX_WIND=', NDX_WIND
WRITE(WINDBG%FHNDL,*) 'NDY_WIND=', NDY_WIND
FLUSH(WINDBG%FHNDL)
END IF
DEALLOCATE(IMAT, JMAT, COUNTERMAT)
ALLOCATE (WIND_X4(NDX_WIND,NDY_WIND), WIND_Y4(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 41')
CALL SYNCHRONIZE_WIND_TIME_IFILE_IT
# ifdef MPI_PARALL_GRID
END IF
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE READ_NETCDF_CRFS(IFILE, IT, eField)
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
!
! READS WIND_Y, WIND_X and PRESSURE from a given NCID within one CRFS file
!
INTEGER, INTENT(IN) :: IFILE, IT
REAL(rkind), intent(inout) :: eField(MNP,2)
character (len = *), parameter :: CallFct="READ_NETCDF_CRFS"
INTEGER :: DWIND_X_ID, DWIND_Y_ID
INTEGER :: numLons, numLats, numTime, numHeights, iy, counter, ip, i, j, ix
REAL(rkind), ALLOCATABLE :: TMP(:,:)
REAL(rkind), ALLOCATABLE :: U(:), V(:)
REAL(rkind), SAVE :: TIME
REAL(rkind) :: Vtotal1(MNP_WIND)
REAL(rkind) :: Vtotal2(MNP_WIND)
REAL(rkind) :: Vlocal(MNP)
INTEGER, DIMENSION (nf90_max_var_dims) :: dimIDs
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'READ_NETCDF_CRFS IFILE=', IFILE, ' IT=', IT
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
IF (IFILE .GT. NUM_NETCDF_FILES) CALL WWM_ABORT('SOMETHING IS WRONG WE RUN OUT OF WIND TIME')
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(IFILE)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(IFILE), NF90_NOWRITE, WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, '10u', DWIND_X_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, DWIND_X_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = numLons)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(2), len = numLats)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(3), len = numHeights)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(4), len = numTime)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
ISTAT = nf90_inq_varid(WIND_NCID, '10v', DWIND_Y_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WIND_NCID, DWIND_Y_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(1), len = numLons)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(2), len = numLats)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(3), len = numHeights)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT)
ISTAT = nf90_inquire_dimension(WIND_NCID, dimIDs(4), len = numTime)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 13, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, DWIND_X_ID, WIND_X, start = (/ 1, 1, 1, IT /), count = (/ NDX_WIND, NDY_WIND,1,1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, ISTAT)
ISTAT = NF90_GET_VAR(WIND_NCID, DWIND_Y_ID, WIND_Y, start = (/ 1, 1, 1, IT /), count = (/ NDX_WIND, NDY_WIND,1,1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, ISTAT)
ISTAT = NF90_CLOSE(WIND_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 16, ISTAT)
CALL INTER_STRUCT_DATA(NDX_WIND,NDY_WIND,DX_WIND,DY_WIND,OFFSET_X_WIND,OFFSET_Y_WIND,WIND_X,Vtotal1)
CALL INTER_STRUCT_DATA(NDX_WIND,NDY_WIND,DX_WIND,DY_WIND,OFFSET_X_WIND,OFFSET_Y_WIND,WIND_Y,Vtotal2)
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
eField(:,1)=Vtotal1
eField(:,2)=Vtotal2
ELSE
CALL SCATTER_ONED_ARRAY(Vtotal1, Vlocal)
eField(:,1)=Vlocal
CALL SCATTER_ONED_ARRAY(Vtotal2, Vlocal)
eField(:,2)=Vlocal
END IF
# else
eField(:,1)=Vtotal1
eField(:,2)=Vtotal2
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE READ_NETCDF_NARR(IFILE, IT, eField)
USE DATAPOOL
USE NETCDF
IMPLICIT NONE
!
! READS WIND_Y, WIND_X and PRESSURE from a given NCID within one NARR file
!
INTEGER, INTENT(IN) :: IFILE, IT
REAL(rkind), intent(out) :: eField(MNP,2)
INTEGER :: DWIND_X_ID, DWIND_Y_ID
INTEGER :: numLons, numLats, counter, ip, i, j, ix
character (len = *), parameter :: CallFct="READ_NETCDF_NARR"
REAL(rkind), ALLOCATABLE :: TMP(:,:)
REAL(rkind),SAVE :: TIME
REAL(rkind) :: scale_factor
REAL(rkind) :: sumWi, eF1, eF2
REAL(rkind) :: Vtotal1(MNP_WIND)
REAL(rkind) :: Vtotal2(MNP_WIND)
REAL(rkind) :: Vlocal(MNP)
INTEGER, DIMENSION (nf90_max_var_dims) :: dimIDs
real(rkind) :: ErrorCoord, XPinterp, YPinterp
INTEGER IEwind
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
IF (2*IFILE .GT. NUM_NETCDF_FILES) THEN
CALL WWM_ABORT('NARR ERROR: Not enough files')
END IF
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(2*IFILE-1)))
ISTAT = NF90_OPEN(TRIM(NETCDF_FILE_NAMES(2*IFILE-1)), NF90_NOWRITE, WINDX_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = nf90_inq_varid(WINDX_NCID, 'uwnd', DWIND_X_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WINDX_NCID, DWIND_X_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = nf90_inquire_dimension(WINDX_NCID, dimIDs(1), len = numLons)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
ISTAT = nf90_inquire_dimension(WINDX_NCID, dimIDs(2), len = numLats)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
ISTAT = nf90_get_att(WINDX_NCID, DWIND_X_ID, 'scale_factor', scale_factor)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(NETCDF_FILE_NAMES(2*IFILE)))
ISTAT = NF90_OPEN(NETCDF_FILE_NAMES(2*IFILE), NF90_NOWRITE, WINDY_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = nf90_inq_varid(WINDY_NCID, 'vwnd', DWIND_Y_ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
ISTAT = NF90_INQUIRE_VARIABLE(WINDY_NCID, DWIND_Y_ID, dimids = dimIDs)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ISTAT = nf90_inquire_dimension(WINDY_NCID, dimIDs(1), len = numLons)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
ISTAT = nf90_inquire_dimension(WINDY_NCID, dimIDs(2), len = numLats)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT)
ISTAT = NF90_GET_VAR(WINDX_NCID, DWIND_X_ID, WIND_X4, start = (/ 1, 1, IT /), count = (/ NDX_WIND, NDY_WIND, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 14, ISTAT)
ISTAT = NF90_GET_VAR(WINDY_NCID, DWIND_Y_ID, WIND_Y4, start = (/ 1, 1, IT /), count = (/ NDX_WIND, NDY_WIND, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 15, ISTAT)
ISTAT = NF90_CLOSE(WINDX_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 16, ISTAT)
ISTAT = NF90_CLOSE(WINDY_NCID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 17, ISTAT)
IF (.FALSE.) THEN
ALLOCATE(TMP(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 43')
DO IX = 1, NDX_WIND
tmp(NDX_WIND-(IX-1),:) = wind_x4(IX,:)
END DO
wind_x4 = tmp
DO IX = 1, NDX_WIND
tmp(NDX_WIND-(IX-1),:) = wind_y4(IX,:)
END DO
wind_y4 = tmp
! DO IY = 1, NDY_WIND
! tmp(:,NDY_WIND-(IY-1)) = atmo_press(:,IY)
! END DO
! atmo_press = tmp
DEALLOCATE(TMP)
END IF
IF (.FALSE.) THEN
ALLOCATE(TMP(NDX_WIND,NDY_WIND), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 44')
DO IX = 1, NDX_WIND
IF (IX .GT. NDX_WIND/2) THEN
tmp(IX-NDX_WIND/2,:) = wind_x(IX,:)
ELSE
tmp(IX+NDX_WIND/2,:) = wind_x(IX,:)
END IF
END DO
wind_x = tmp
DO IX = 1, NDX_WIND
IF (IX .GT. NDX_WIND/2) THEN
tmp(IX-NDX_WIND/2,:) = wind_y(IX,:)
ELSE
tmp(IX+NDX_WIND/2,:) = wind_y(IX,:)
END IF
END DO
wind_y = tmp
DEALLOCATE(TMP)
END IF
COUNTER = 0
DO I = 1, NDX_WIND
DO J = 1, NDY_WIND
COUNTER = COUNTER + 1
IF (WIND_X4(I,J) .GT. 32766 .OR. WIND_Y4(I,J) .GT. 32766) THEN
UWND_NARR(COUNTER) = 0.
VWND_NARR(COUNTER) = 0.
ELSE
UWND_NARR(COUNTER) = WIND_X4(I,J) * scale_factor
VWND_NARR(COUNTER) = WIND_Y4(I,J) * scale_factor
END IF
END DO
END DO
IF (LWRITE_ORIG_WIND) THEN
TIME = TIME + 1.
WRITE(3011) TIME
WRITE(3011) (UWND_NARR(IP), VWND_NARR(IP), SQRT(UWND_NARR(IP)**2.+VWND_NARR(IP)**2.), IP = 1, NDX_WIND*NDY_WIND)
END IF
ErrorCoord=0
DO IP = 1, MNP_WIND
IEwind=WIND_ELE(IP)
IF (IEwind .gt. 0) then
XPinterp=0
YPinterp=0
eF1=0
eF2=0
sumWi=0
DO I=1,3
XPinterp=XPinterp + WI_NARR(IP,I)*XYPWIND(1,INE_WIND(I,IEwind))
YPinterp=YPinterp + WI_NARR(IP,I)*XYPWIND(2,INE_WIND(I,IEwind))
eF1=eF1+WI_NARR(IP,I)*UWND_NARR(INE_WIND(I,IEwind))
eF2=eF2+WI_NARR(IP,I)*VWND_NARR(INE_WIND(I,IEwind))
sumWi=sumWi + WI_NARR(IP,I)
END DO
Vtotal1(IP)=eF1
Vtotal2(IP)=eF2
ErrorCoord=ErrorCoord + abs(XPinterp - XP_WIND(IP)) + abs(YPinterp - YP_WIND(IP)) + abs(sumWi - 1)
ELSE
Vtotal1(IP)=ZERO
Vtotal2(IP)=ZERO
ENDIF
END DO
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'ErrorCoord=', ErrorCoord
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
eField(:,1)=Vtotal1
eField(:,2)=Vtotal2
ELSE
CALL SCATTER_ONED_ARRAY(Vtotal1, Vlocal)
eField(:,1)=Vlocal
CALL SCATTER_ONED_ARRAY(Vtotal2, Vlocal)
eField(:,2)=Vlocal
END IF
# else
eField(:,1)=Vtotal1
eField(:,2)=Vtotal2
# endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE CHECK_WIND_TIME(nbtime_mjd, WIND_TIME_MJD)
USE DATAPOOL, only : SEWI, WINDBG, rkind, THR, wwmerr, WRITEWINDBGFLAG
IMPLICIT NONE
integer, intent(in) :: nbtime_mjd
real(rkind), intent(in) :: WIND_TIME_MJD(nbtime_mjd)
CHARACTER(LEN=15) :: eTimeStr
IF (SEWI%BMJD .LT. minval(WIND_TIME_MJD) - THR) THEN
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'END OF RUN'
WRITE(WINDBG%FHNDL,*) 'WIND START TIME is outside CF wind_time range!'
END IF
CALL MJD2CT(SEWI%BMJD,eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'SEWI%BMJD=', SEWI%BMJD, ' date=', eTimeStr
CALL MJD2CT(SEWI%EMJD,eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'SEWI%EMJD=', SEWI%EMJD, ' date=', eTimeStr
CALL MJD2CT(minval(WIND_TIME_MJD),eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'min(WIND_TIME_MJD)=', minval(WIND_TIME_MJD), ' date=', eTimeStr
CALL MJD2CT(maxval(WIND_TIME_MJD),eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'max(WIND_TIME_MJD)=', maxval(WIND_TIME_MJD), ' date=', eTimeStr
IF (WRITEWINDBGFLAG == 1) FLUSH(WINDBG%FHNDL)
IF (WRITEWINDBGFLAG == 1) WRITE(wwmerr, *) 'Error in WIND_TIME_MJD 1, read ', TRIM(WINDBG%FNAME)
CALL WWM_ABORT(wwmerr)
END IF
IF (SEWI%EMJD .GT. maxval(WIND_TIME_MJD) + THR) THEN
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'END OF RUN'
WRITE(WINDBG%FHNDL,*) 'WIND END TIME is outside CF wind_time range!'
END IF
CALL MJD2CT(SEWI%BMJD,eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'SEWI%BMJD=', SEWI%BMJD, ' date=', eTimeStr
CALL MJD2CT(SEWI%EMJD,eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'SEWI%EMJD=', SEWI%EMJD, ' date=', eTimeStr
CALL MJD2CT(minval(WIND_TIME_MJD),eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'min(WIND_TIME_MJD)=', minval(WIND_TIME_MJD), ' date=', eTimeStr
CALL MJD2CT(maxval(WIND_TIME_MJD),eTimeStr)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'max(WIND_TIME_MJD)=', maxval(WIND_TIME_MJD), ' date=', eTimeStr
IF (WRITEWINDBGFLAG == 1) FLUSH(WINDBG%FHNDL)
IF (WRITEWINDBGFLAG == 1) WRITE(wwmerr, *) 'Error in WIND_TIME_MJD 2, read ', TRIM(WINDBG%FNAME)
CALL WWM_ABORT(wwmerr)
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE READ_INTERP_NETCDF_CF_WWM_WIND(RECORD_IN, outwind)
USE NETCDF
USE DATAPOOL
IMPLICIT NONE
INTEGER, INTENT(in) :: RECORD_IN
REAL(rkind), INTENT(out) :: outwind(MNP,2)
REAL(rkind) :: varTotal(MNP_WIND,2), Vtotal(MNP_WIND), Vlocal(MNP)
character (len = *), parameter :: CallFct="READ_INTERP_NETCDF_CF_WWM_WIND"
INTEGER :: FID, ID
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
CALL TEST_FILE_EXIST_DIE("Missing wind file : ", TRIM(WIN%FNAME))
ISTAT = NF90_OPEN(WIN%FNAME, NF90_NOWRITE, FID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = NF90_inq_varid(FID, 'Uwind', ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_GET_VAR(FID, ID, UWIND_FD, start = (/ 1, 1, RECORD_IN /), count = (/ NDX_WIND_FD, NDY_WIND_FD, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = NF90_inq_varid(FID, 'Vwind', ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
ISTAT = NF90_GET_VAR(FID, ID, VWIND_FD, start = (/ 1, 1, RECORD_IN /), count = (/ NDX_WIND_FD, NDY_WIND_FD, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
ISTAT = NF90_CLOSE(FID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
CALL KERNEL_INTERP_UV_WINDFD(varTotal)
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
outwind=varTotal
ELSE
Vtotal=varTotal(:,1)
CALL SCATTER_ONED_ARRAY(Vtotal, Vlocal)
outwind(:,1)=Vlocal
!
Vtotal=varTotal(:,2)
CALL SCATTER_ONED_ARRAY(Vtotal, Vlocal)
outwind(:,2)=Vlocal
END IF
# else
outwind=varTotal
# endif
END SUBROUTINE
!****************************************************************************
!* CF_COMPLIANT WIND *
!* This is the standard way to write netcdf data. *
!* See *
!* http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.pdf *
!* for details *
!****************************************************************************
SUBROUTINE INIT_NETCDF_CF_WWM_WIND(eVAR)
USE NETCDF
USE DATAPOOL
IMPLICIT NONE
INTEGER :: fid, varid, dimids(2), dimidsB(3)
integer nbChar
TYPE(VAR_NETCDF_CF), intent(inout) :: eVAR
REAL(rkind), ALLOCATABLE :: CF_LON(:,:), CF_LAT(:,:)
character (len = *), parameter :: CallFct="INIT_NETCDF_CF"
character (len=200) :: CoordString
character (len=100) :: Xname, Yname, eStrUnitTime
character (len=20) :: WindTimeStr
real(rkind) :: ConvertToDay
real(rkind) :: eTimeStart
character(len=100) :: CHRERR
integer posBlank, alen
type(FD_FORCING_GRID) TheInfo
integer IX, IY
CALL INIT_DIRECT_NETCDF_CF(eVAR, MULTIPLE_IN_WIND, WIN%FNAME, "Uwind")
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
ISTAT = nf90_open(WIN%FNAME, nf90_nowrite, fid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = nf90_inq_varid(fid, "Uwind", varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = nf90_get_att(fid, varid, "coordinates", CoordString)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
alen=LEN_TRIM(CoordString)
posBlank=INDEX(CoordString(1:alen), ' ')
Xname=CoordString(1:posBlank-1)
Yname=CoordString(posBlank+1:alen)
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'Xname=', TRIM(Xname)
WRITE(WINDBG%FHNDL,*) 'Yname=', TRIM(Yname)
FLUSH(WINDBG%FHNDL)
END IF
! Reading lontitude/latitude array
ISTAT = nf90_inq_varid(fid, Xname, varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
ISTAT = nf90_inquire_variable(fid, varid, dimids=dimids)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
ISTAT = nf90_inquire_dimension(fid, dimids(1), len=NDX_WIND_FD)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = nf90_inquire_dimension(fid, dimids(2), len=NDY_WIND_FD)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'NDX_WIND_FD=', NDX_WIND_FD
WRITE(WINDBG%FHNDL,*) 'NYX_WIND_FD=', NDY_WIND_FD
FLUSH(WINDBG%FHNDL)
END IF
allocate(CF_LON(NDX_WIND_FD, NDY_WIND_FD), CF_LAT(NDX_WIND_FD, NDY_WIND_FD), UWIND_FD(NDX_WIND_FD, NDY_WIND_FD), VWIND_FD(NDX_WIND_FD, NDY_WIND_FD), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 47')
ISTAT = nf90_get_var(fid, varid, CF_LON)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ISTAT = nf90_inq_varid(fid, Yname, varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
ISTAT = nf90_get_var(fid, varid, CF_LAT)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 12, ISTAT)
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'MIN,MAX CF_LON=', MINVAL(CF_LON), MAXVAL(CF_LON)
WRITE(WINDBG%FHNDL,*) 'MIN,MAX CF_LAT=', MINVAL(CF_LAT), MAXVAL(CF_LAT)
FLUSH(WINDBG%FHNDL)
END IF
TheInfo % nx_dim = NDX_WIND_FD
TheInfo % ny_dim = NDY_WIND_FD
allocate(TheInfo % LON(NDX_WIND_FD, NDY_WIND_FD), TheInfo % LAT(NDX_WIND_FD, NDY_WIND_FD), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 47')
DO IX=1,NDX_WIND_FD
DO IY=1,NDY_WIND_FD
TheInfo % LON(IX,IY) = CF_LON(IX,IY)
TheInfo % LAT(IX,IY) = CF_LAT(IX,IY)
END DO
END DO
DEALLOCATE(CF_LON, CF_LAT)
CALL COMPUTE_CF_COEFFICIENTS(TheInfo)
Deallocate(TheInfo % LON, TheInfo % LAT)
# ifdef MPI_PARALL_GRID
END IF
# endif
END SUBROUTINE
!**********************************************************************
!* This is for direct to elements forcing in netcdf *
!* wind_time as usual *
!* float Uwind(wind_time, mnp) *
!* float Vwind(wind_time, mnp) *
!* should be better than text. *
!**********************************************************************
SUBROUTINE READ_DIRECT_NETCDF_CF(eVAR, RECORD_IN, outwind)
USE NETCDF
USE DATAPOOL
IMPLICIT NONE
TYPE(VAR_NETCDF_CF) :: eVAR
INTEGER, INTENT(in) :: RECORD_IN
REAL(rkind), INTENT(out) :: outwind(MNP,2)
character (len = *), parameter :: CallFct="READ_DIRECT_NETCDF_CF"
INTEGER :: FID, ID
real(rkind) :: U_tot(np_total), V_tot(np_total)
real(rkind) :: Vtotal1(np_total), Vtotal2(np_total)
real(rkind) :: Vlocal(MNP)
real(rkind) :: cf_scale_factor, cf_add_offset
# ifdef MPI_PARALL_GRID
integer IP_glob, IP
# endif
character(len=10) :: eStrU, eStrV
cf_scale_factor = eVAR % cf_scale_factor
cf_add_offset = eVAR % cf_add_offset
IF (eVAR % idVar == 1) THEN
eStrU='Uwind'
eStrV='Vwind'
ELSE
eStrU='UsurfCurr'
eStrV='VsurfCurr'
END IF
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
CALL TEST_FILE_EXIST_DIE("Missing file : ", TRIM(eVAR % eFileName))
ISTAT = NF90_OPEN(TRIM(eVAR % eFileName), NF90_NOWRITE, FID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = NF90_inq_varid(FID, TRIM(eStrU), ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_GET_VAR(FID, ID, U_tot, start = (/ 1, RECORD_IN /), count = (/ np_total, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = NF90_inq_varid(FID, TRIM(eStrV), ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
ISTAT = NF90_GET_VAR(FID, ID, V_tot, start = (/ 1, RECORD_IN /), count = (/ np_total, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
ISTAT = NF90_CLOSE(FID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
Vtotal1 = cf_add_offset + cf_scale_factor*U_tot
Vtotal2 = cf_add_offset + cf_scale_factor*V_tot
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
DO IP=1,MNP
IP_glob=iplg(IP)
outwind(IP,1)=Vtotal1(IP_glob)
outwind(IP,2)=Vtotal2(IP_glob)
END DO
ELSE
CALL SCATTER_ONED_ARRAY(Vtotal1, Vlocal)
outwind(:,1)=Vlocal
CALL SCATTER_ONED_ARRAY(Vtotal2, Vlocal)
outwind(:,2)=Vlocal
END IF
# else
outwind(:,1)=Vtotal1
outwind(:,2)=Vtotal2
# endif
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'READ_DIRECT_NETCDF_CF'
WRITE(WINDBG%FHNDL,*) 'RECORD_IN=', RECORD_IN
WRITE(WINDBG%FHNDL,*) 'UWIND_FE, min/max=', minval(outwind(:,1)), maxval(outwind(:,1))
WRITE(WINDBG%FHNDL,*) 'VWIND_FE, min/max=', minval(outwind(:,2)), maxval(outwind(:,2))
FLUSH(WINDBG%FHNDL)
END IF
END SUBROUTINE
!****************************************************************************
!* *
!****************************************************************************
SUBROUTINE READ_DIRECT_NETCDF_CF1(eVAR, RECORD_IN, outvar)
USE NETCDF
USE DATAPOOL
IMPLICIT NONE
TYPE(VAR_NETCDF_CF) :: eVAR
INTEGER, INTENT(in) :: RECORD_IN
REAL(rkind), INTENT(out) :: outvar(MNP)
character (len = *), parameter :: CallFct="READ_DIRECT_NETCDF_CF"
INTEGER :: FID, ID
real(rkind) :: VAR_tot(np_total)
real(rkind) :: Vtotal(np_total)
real(rkind) :: Vlocal(MNP)
real(rkind) :: cf_scale_factor, cf_add_offset
# ifdef MPI_PARALL_GRID
integer IP_glob, IP
# endif
character(len=10) :: eStrU, eStrV
cf_scale_factor = eVAR % cf_scale_factor
cf_add_offset = eVAR % cf_add_offset
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
CALL TEST_FILE_EXIST_DIE("Missing file : ", TRIM(eVAR % eFileName))
ISTAT = NF90_OPEN(TRIM(eVAR % eFileName), NF90_NOWRITE, FID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
ISTAT = NF90_inq_varid(FID, TRIM(eVAR % eString), ID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = NF90_GET_VAR(FID, ID, VAR_tot, start = (/ 1, RECORD_IN /), count = (/ np_total, 1 /))
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = NF90_CLOSE(FID)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
Vtotal = cf_add_offset + cf_scale_factor*VAR_tot
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
DO IP=1,MNP
IP_glob=iplg(IP)
outvar(IP)=Vtotal(IP_glob)
END DO
ELSE
CALL SCATTER_ONED_ARRAY(Vtotal, outvar)
END IF
# else
outvar=Vtotal
# endif
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'READ_DIRECT_NETCDF_CF1'
WRITE(WINDBG%FHNDL,*) 'RECORD_IN=', RECORD_IN
FLUSH(WINDBG%FHNDL)
END IF
END SUBROUTINE
!****************************************************************************
!* *
!****************************************************************************
SUBROUTINE INIT_DIRECT_NETCDF_CF(eVAR, MULTIPLE_IN, eFileName, eString)
USE NETCDF
USE DATAPOOL
IMPLICIT NONE
TYPE(VAR_NETCDF_CF), intent(inout) :: eVAR
logical, intent(in) :: MULTIPLE_IN
character(len=100), intent(in) :: eFileName
character(len=*), intent(in) :: eString
!
INTEGER :: fid, varid
INTEGER :: dimidsB(3), dimids(3)
integer nbChar
character (len=20) :: WindTimeStr
character(len=100) :: CHRERR
character (len = *), parameter :: CallFct="INIT_DIRECT_NETCDF_CF"
character (len=100) :: eStrUnitTime
real(rkind) :: ConvertToDay
real(rkind) :: eTimeStart
real(rkind) :: cf_scale_factor, cf_add_offset
integer nbtime_mjd
real(rkind), allocatable :: wind_time_mjd(:)
integer eInt(1)
real(rkind) :: eReal(2)
integer IPROC
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN .or. (myrank .eq. 0)) THEN
# endif
ISTAT = nf90_open(TRIM(eFileName), nf90_nowrite, fid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 1, ISTAT)
! Reading wind attributes
! Print *, 'eString=', TRIM(eString)
! Print *, 'FNAME=', TRIM(eFileName)
ISTAT = nf90_inq_varid(fid, TRIM(eString), varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 2, ISTAT)
ISTAT = nf90_inquire_variable(fid, varid, dimids=dimidsB)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 3, ISTAT)
ISTAT = nf90_inquire_dimension(fid, dimidsB(3), name=WindTimeStr)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 4, ISTAT)
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'variable used for time=', TRIM(WindTimeStr)
FLUSH(WINDBG%FHNDL)
END IF
ISTAT = nf90_get_att(fid, varid, "scale_factor", cf_scale_factor)
IF (ISTAT /= 0) THEN
CHRERR = nf90_strerror(ISTAT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'CHRERR=', TRIM(CHRERR)
cf_scale_factor=ONE
ENDIF
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'cf_scale_factor=', cf_scale_factor
FLUSH(WINDBG%FHNDL)
END IF
ISTAT = nf90_get_att(fid, varid, "add_offset", cf_add_offset)
IF (ISTAT /= 0) THEN
CHRERR = nf90_strerror(ISTAT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'CHRERR=', TRIM(CHRERR)
cf_add_offset=ZERO
ENDIF
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'cf_add_offset=', cf_add_offset
FLUSH(WINDBG%FHNDL)
END IF
! Reading time
ISTAT = nf90_inq_varid(fid, WindTimeStr, varid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 5, ISTAT)
ISTAT = nf90_inquire_attribute(fid, varid, "units", len=nbChar)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 6, ISTAT)
ISTAT = nf90_get_att(fid, varid, "units", eStrUnitTime)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 7, ISTAT)
CALL CF_EXTRACT_TIME(eStrUnitTime, ConvertToDay, eTimeStart)
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL,*) 'eTimeStart=', eTimeStart
FLUSH(WINDBG%FHNDL)
END IF
ISTAT = nf90_inquire_variable(fid, varid, dimids=dimids)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 8, ISTAT)
ISTAT = nf90_inquire_dimension(fid, dimids(1), len=nbtime_mjd)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 9, ISTAT)
allocate(wind_time_mjd(nbtime_mjd), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 48')
ISTAT = nf90_get_var(fid, varid, wind_time_mjd)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 10, ISTAT)
ISTAT = nf90_close(fid)
CALL GENERIC_NETCDF_ERROR_WWM(CallFct, 11, ISTAT)
wind_time_mjd(:) = wind_time_mjd(:)*ConvertToDay + eTimeStart
CALL CHECK_WIND_TIME(nbtime_mjd, WIND_TIME_MJD)
# ifdef MPI_PARALL_GRID
END IF
# endif
#ifdef MPI_PARALL_GRID
IF (.NOT. MULTIPLE_IN_WIND) THEN
IF (myrank .eq. 0) THEN
eInt(1)=nbtime_mjd
DO IPROC=2,nproc
CALL MPI_SEND(eInt,1,itype, iProc-1, 811, comm, ierr)
END DO
eReal(1)=cf_scale_factor
eReal(2)=cf_add_offset
DO IPROC=2,nproc
CALL MPI_SEND(eReal,2,rtype, iProc-1, 813, comm, ierr)
END DO
DO IPROC=2,nproc
CALL MPI_SEND(WIND_TIME_MJD,nbtime_mjd,rtype, iProc-1, 812, comm, ierr)
END DO
ELSE
CALL MPI_RECV(eInt,1,itype, 0, 811, comm, istatus, ierr)
nbtime_mjd=eInt(1)
CALL MPI_RECV(eReal,2,rtype, 0, 813, comm, istatus, ierr)
cf_scale_factor=eReal(1)
cf_add_offset=eReal(2)
ALLOCATE(WIND_TIME_MJD(nbtime_mjd), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 3')
CALL MPI_RECV(WIND_TIME_MJD,nbtime_mjd,rtype, 0, 812, comm, istatus, ierr)
END IF
END IF
#endif
eVAR % cf_scale_factor = cf_scale_factor
eVAR % cf_add_offset = cf_add_offset
eVAR % nbTime = nbtime_mjd
allocate(eVAR % ListTime(nbtime_mjd))
eVAR % ListTime = WIND_TIME_MJD
eVAR % eFileName = eFileName
eVAR % eString = eString
IF (TRIM(eString) == 'Uwind') THEN
eVAR % idVar = 1
ELSE
eVAR % idVar = 2
END IF
END SUBROUTINE
#endif
#ifdef GRIB_API_ECMWF
!****************************************************************************
!* Raw reading of time entry for GRIB *
!****************************************************************************
SUBROUTINE RAW_READ_TIME_OF_GRIB_FILE(ifile, eGrib, STEPRANGE_IN, eTimeOut)
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
integer, intent(in) :: ifile
integer, intent(in) :: eGrib
LOGICAL, intent(in) :: STEPRANGE_IN
real(rkind), intent(out) :: eTimeOut
!
LOGICAL :: USE_DATATIME = .TRUE.
integer eYear, eMonth, eDay, resYear, resMonth
integer eHour, eMin, eSec, resHour, resMin
integer dataDate, stepRange, dataTime
character (len=15) :: eStrTime
REAL(rkind) :: eTimeBase
call grib_get(eGrib, 'dataDate', dataDate)
eYear=(dataDate - mod(dataDate,10000))/10000
resYear=dataDate - 10000*eYear
eMonth=(resYear - mod(resYear,100))/100
resMonth=resYear - 100*eMonth;
eDay=resMonth
IF (STEPRANGE_IN) THEN
call grib_get(eGrib, 'stepRange', stepRange)
ELSE
stepRange=0
END IF
IF (USE_DATATIME) THEN
call grib_get(eGrib, 'dataTime', dataTime)
eHour=(dataTime - mod(dataTime,100))/100
eMin=dataTime - 100*eHour
eSec=0
ELSE
eHour=0
eMin=0
eSec=0
END IF
WRITE(eStrTime,10) eYear, eMonth, eDay, eHour, eMin, eSec
10 FORMAT(i4.4,i2.2,i2.2,'.',i2.2,i2.2,i2.2)
CALL CT2MJD(eStrTime, eTimeBase)
eTimeOut=eTimeBase + MyREAL(stepRange)/24.0_rkind
END SUBROUTINE
!****************************************************************************
!* Reading time from a GRIB file *
!****************************************************************************
SUBROUTINE READ_TIME_OF_GRIB_FILE(eTimeOut, eFile, STEPRANGE_IN)
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
REAL(rkind), intent(out) :: eTimeOut
CHARACTER(len=140), intent(in) :: eFile
LOGICAL, intent(in) :: STEPRANGE_IN
integer ifile, i, n
integer, allocatable :: igrib(:)
CALL TEST_FILE_EXIST_DIE("Missing grib file: ", TRIM(eFile))
CALL GRIB_OPEN_FILE(ifile, TRIM(eFile), 'r')
call grib_count_in_file(ifile,n)
allocate(igrib(n))
i=1
call grib_new_from_file(ifile, igrib(i))
CALL RAW_READ_TIME_OF_GRIB_FILE(ifile, igrib(i), STEPRANGE_IN, eTimeOut)
CALL GRIB_CLOSE_FILE(ifile)
deallocate(igrib)
END SUBROUTINE
!****************************************************************************
!* Reading grid information from a GRIB file (case 1) *
!****************************************************************************
SUBROUTINE READ_GRID_INFO_FROM_GRIB_TYPE1(TheInfo, eGrib)
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
type(FD_FORCING_GRID), intent(out) :: TheInfo
integer, intent(in) :: eGrib
!
REAL(rkind) ::longitudeOfFirstPointInDegrees, latitudeOfFirstPointInDegrees, longitudeOfLastPointInDegrees, latitudeOfLastPointInDegrees
REAL(rkind) :: deltaLAT, deltaLON
REAL(rkind) :: iDirectionIncrement, jDirectionIncrement
integer nx_dim, ny_dim, iX, iY
call grib_get(eGrib,"numberOfPointsAlongAParallel", nx_dim)
call grib_get(eGrib,"numberOfPointsAlongAMeridian", ny_dim)
IF (WRITESTATFLAG == 1) THEN
WRITE(STAT%FHNDL, *) 'nx_dim=', nx_dim
WRITE(STAT%FHNDL, *) 'ny_dim=', ny_dim
END IF
TheInfo % nx_dim = nx_dim
TheInfo % ny_dim = ny_dim
allocate(TheInfo % LON(nx_dim, ny_dim), TheInfo % LAT(nx_dim, ny_dim), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 47')
call grib_get(eGrib, 'longitudeOfFirstGridPointInDegrees', longitudeOfFirstPointInDegrees)
call grib_get(eGrib, 'latitudeOfFirstGridPointInDegrees', latitudeOfFirstPointInDegrees)
call grib_get(eGrib, 'longitudeOfLastGridPointInDegrees', longitudeOfLastPointInDegrees)
call grib_get(eGrib, 'latitudeOfLastGridPointInDegrees', latitudeOfLastPointInDegrees)
call grib_get(eGrib, 'iDirectionIncrementInDegrees', iDirectionIncrement)
call grib_get(eGrib, 'jDirectionIncrementInDegrees', jDirectionIncrement)
IF (WRITESTATFLAG == 1) THEN
WRITE(STAT%FHNDL, *) 'LONGITUDE'
WRITE(STAT%FHNDL, *) 'longitudeOfFirstGridPointInDegrees=', longitudeOfFirstPointInDegrees
WRITE(STAT%FHNDL, *) 'longitudeOfLastGridPointInDegrees=', longitudeOfLastPointInDegrees
WRITE(STAT%FHNDL, *) 'LATITUDE'
WRITE(STAT%FHNDL, *) 'latitudeOfFirstGridPointInDegrees=', latitudeOfFirstPointInDegrees
WRITE(STAT%FHNDL, *) 'latitudeOfLastGridPointInDegrees=', latitudeOfLastPointInDegrees
WRITE(STAT%FHNDL, *) 'iDirectionIncrement=', iDirectionIncrement
WRITE(STAT%FHNDL, *) 'jDirectionIncrement=', jDirectionIncrement
END IF
deltaLON=(longitudeOfLastPointInDegrees - longitudeOfFirstPointInDegrees)/(nx_dim - 1)
deltaLAT=(latitudeOfLastPointInDegrees - latitudeOfFirstPointInDegrees)/(ny_dim - 1)
DO iX=1,nx_dim
DO iY=1,ny_dim
TheInfo % LON(iX,iY)=longitudeOfFirstPointInDegrees + (iX-1)*deltaLON
TheInfo % LAT(iX,iY)=latitudeOfFirstPointInDegrees + (iY-1)*deltaLAT
END DO
END DO
END SUBROUTINE
!****************************************************************************
!* Reading grid information from a GRIB file (case 2) *
!****************************************************************************
SUBROUTINE phirot2phi(eLatOut, phirot, rlarot, polphi, pollam, polgam)
USE DATAPOOL
IMPLICIT NONE
real(rkind), intent(out) :: eLatOut
real(rkind), intent(in) :: phirot, rlarot, polphi, pollam, polgam
real(rkind) :: zsinpol
real(rkind) :: zcospol
real(rkind) :: zphis
real(rkind) :: zrlas
real(rkind) :: zarg
real(rkind) :: zgam
zsinpol = sin(DEGRAD * polphi)
zcospol = cos(DEGRAD * polphi)
zphis = DEGRAD * phirot
if (rlarot .gt. 180.) THEN
zrlas = rlarot - 360.
ELSE
zrlas = rlarot
END IF
zrlas = DEGRAD * zrlas;
if (ABS(polgam) .gt. 0) THEN
zgam = DEGRAD * polgam;
zarg = zsinpol*sin(zphis) + zcospol*cos(zphis) * ( cos(zrlas)*cos(zgam) - sin(zgam)*sin(zrlas))
ELSE
zarg = zcospol * cos(zphis) * cos(zrlas) + zsinpol * sin(zphis)
END IF
eLatOut = RADDEG * asin(zarg);
END SUBROUTINE
!****************************************************************************
!* Reading grid information from a GRIB file (case 2) *
!****************************************************************************
SUBROUTINE rlarot2rla(eLonOut, phirot, rlarot, polphi, pollam, polgam)
USE DATAPOOL
IMPLICIT NONE
real(rkind), intent(out) :: eLonOut
real(rkind), intent(in) :: phirot, rlarot, polphi, pollam, polgam
!
real(rkind) :: zsinpol
real(rkind) :: zcospol
real(rkind) :: zphis
real(rkind) :: zrlas
real(rkind) :: zgam, zlampol
real(rkind) :: zarg1, zarg2
zsinpol = sin(DEGRAD * polphi)
zcospol = cos(DEGRAD * polphi)
zphis = DEGRAD * phirot
if (rlarot .gt. 180.) THEN
zrlas = rlarot - 360.
ELSE
zrlas = rlarot
END IF
zrlas = DEGRAD * zrlas
zlampol = DEGRAD * pollam
IF (ABS(polgam) .gt. 0) THEN
zgam = DEGRAD * polgam;
zarg1 = sin (zlampol) * &
& (- zsinpol*cos(zphis) * (cos(zrlas)*cos(zgam) - sin(zrlas)*sin(zgam)) &
& + zcospol * sin(zphis)) &
& - cos (zlampol)*cos(zphis) * (sin(zrlas)*cos(zgam) + cos(zrlas)*sin(zgam))
zarg2 = cos (zlampol) * &
& (- zsinpol*cos(zphis) * (cos(zrlas)*cos(zgam) - sin(zrlas)*sin(zgam)) &
& + zcospol * sin(zphis)) &
& + sin (zlampol)*cos(zphis) * (sin(zrlas)*cos(zgam) + cos(zrlas)*sin(zgam))
ELSE
zarg1 = sin (zlampol) * (-zsinpol * cos(zrlas) * cos(zphis) + &
& zcospol * sin(zphis)) - &
& cos (zlampol) * sin(zrlas) * cos(zphis)
zarg2 = cos (zlampol) * (-zsinpol * cos(zrlas) * cos(zphis) + &
& zcospol * sin(zphis)) + &
& sin (zlampol) * sin(zrlas) * cos(zphis)
END IF
if (zarg2 .eq. 0) zarg2=1.0e-20
eLonOut = RADDEG * atan2(zarg1,zarg2);
END SUBROUTINE
!****************************************************************************
!* Reading grid information from a GRIB file (case 2) *
!****************************************************************************
SUBROUTINE READ_GRID_INFO_FROM_GRIB_TYPE2(TheInfo, eGrib)
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
type(FD_FORCING_GRID), intent(out) :: TheInfo
integer, intent(in) :: eGrib
real(rkind) :: latitudeOfSouthernPoleInDegrees, longitudeOfSouthernPoleInDegrees
real(rkind) :: angleOfRotationInDegrees
real(rkind) :: latitudeOfFirstGridPointInDegrees, longitudeOfFirstGridPointInDegrees
real(rkind) :: latitudeOfLastGridPointInDegrees, longitudeOfLastGridPointInDegrees
real(rkind) :: iDirectionIncrementInDegrees, jDirectionIncrementInDegrees
real(rkind) :: pollat_sp, pollon_sp, polgam, zstartlon_tot, zstartlat_tot
real(rkind) :: zendlon_tot, zendlat_tot, dlon, dlat
real(rkind) :: eLonR, eLatR, eLonOut, eLatOut
real(rkind) :: pollat, pollon
real(rkind) :: startlon_tot, startlat_tot
integer :: nx_dim, ny_dim, iX, iY
!
call grib_get(eGrib, 'latitudeOfSouthernPoleInDegrees',latitudeOfSouthernPoleInDegrees)
call grib_get(eGrib, 'longitudeOfSouthernPoleInDegrees',longitudeOfSouthernPoleInDegrees)
call grib_get(eGrib, 'angleOfRotationInDegrees',angleOfRotationInDegrees)
call grib_get(eGrib, 'latitudeOfFirstGridPointInDegrees',latitudeOfFirstGridPointInDegrees)
call grib_get(eGrib, 'longitudeOfFirstGridPointInDegrees',longitudeOfFirstGridPointInDegrees)
call grib_get(eGrib, 'latitudeOfLastGridPointInDegrees',latitudeOfLastGridPointInDegrees)
call grib_get(eGrib, 'longitudeOfLastGridPointInDegrees',longitudeOfLastGridPointInDegrees)
call grib_get(eGrib, 'iDirectionIncrementInDegrees',iDirectionIncrementInDegrees)
call grib_get(eGrib, 'jDirectionIncrementInDegrees',jDirectionIncrementInDegrees)
pollat_sp=latitudeOfSouthernPoleInDegrees
pollon_sp=longitudeOfSouthernPoleInDegrees
polgam=angleOfRotationInDegrees
zstartlon_tot=longitudeOfFirstGridPointInDegrees
zstartlat_tot=latitudeOfFirstGridPointInDegrees
zendlon_tot=longitudeOfLastGridPointInDegrees
zendlat_tot=latitudeOfLastGridPointInDegrees
dlon=iDirectionIncrementInDegrees
dlat=jDirectionIncrementInDegrees
!
! Now reading the mapped grid
!
call grib_get(eGrib,"Nx", nx_dim)
call grib_get(eGrib,"Ny", ny_dim)
IF (WRITESTATFLAG == 1) THEN
WRITE(STAT%FHNDL, *) 'nx_dim = ', nx_dim
WRITE(STAT%FHNDL, *) 'ny_dim = ', ny_dim
END IF
TheInfo % nx_dim = nx_dim
TheInfo % ny_dim = ny_dim
pollat= - pollat_sp
pollon= pollon_sp - 180.
startlon_tot=zstartlon_tot
startlat_tot=zstartlat_tot
allocate(TheInfo % LON(nx_dim, ny_dim), TheInfo % LAT(nx_dim, ny_dim), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 47')
DO iX=1,nx_dim
DO iY=1,ny_dim
eLonR = startlon_tot + MyREAL(iX-1)*dlon
eLatR = startlat_tot + MyREAL(iY-1)*dlat
CALL phirot2phi(eLatOut, eLatR, eLonR, pollat, pollon, polgam)
CALL rlarot2rla(eLonOut, eLatR, eLonR, pollat, pollon, polgam)
TheInfo % LON(iX,iY) = eLonOut
TheInfo % LAT(iX,iY) = eLatOut
END DO
END DO
END SUBROUTINE
!****************************************************************************
!* Reading grid information from a GRIB file (case 1) *
!****************************************************************************
SUBROUTINE READ_GRID_INFO_FROM_GRIB_TYPE3(TheInfo, eGrib)
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
type(FD_FORCING_GRID), intent(out) :: TheInfo
integer, intent(in) :: eGrib
REAL(rkind), allocatable :: LON_serial(:), LAT_serial(:), DATA_Serial(:)
integer nx_dim, ny_dim, idx, eProd
integer status, iX, iY
!
call grib_get(eGrib,"Nx", nx_dim)
call grib_get(eGrib,"Ny", ny_dim)
IF (WRITESTATFLAG == 1) THEN
WRITE(STAT%FHNDL, *) 'nx_dim = ', nx_dim
WRITE(STAT%FHNDL, *) 'ny_dim = ', ny_dim
END IF
TheInfo % nx_dim = nx_dim
TheInfo % ny_dim = ny_dim
allocate(TheInfo % LON(nx_dim, ny_dim), TheInfo % LAT(nx_dim, ny_dim), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 47')
eProd=nx_dim*ny_dim
allocate(LON_serial(eProd), LAT_serial(eProd), DATA_serial(eProd), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 48')
call grib_get_data(eGrib, LAT_serial, LON_serial, DATA_serial, status)
idx=0
DO iY=1,ny_dim
DO iX=1,nx_dim
idx=idx+1
TheInfo % LON(iX,iY)=LON_serial(idx)
TheInfo % LAT(iX,iY)=LAT_serial(idx)
END DO
END DO
DEALLOCATE(LON_serial, LAT_serial, DATA_serial)
END SUBROUTINE
!****************************************************************************
!* Reading grid information from a GRIB file *
!****************************************************************************
SUBROUTINE READ_GRID_INFO_FROM_GRIB(TheInfo, TheFile, shortName, GRIB_TYPE)
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
type(FD_FORCING_GRID), intent(out) :: TheInfo
character(len=140), intent(in) :: TheFile
character(len=20), intent(in) :: shortName
integer, intent(in) :: GRIB_TYPE
!
integer ifile, i, n
logical WeFound
integer, allocatable :: igrib(:)
character(len=100) eShortName
integer eProd
integer status, idx
integer iX, iY
integer nx_dim, ny_dim
!
CALL TEST_FILE_EXIST_DIE("Missing grib file: ", TRIM(TheFile))
CALL GRIB_OPEN_FILE(ifile, TRIM(TheFile), 'r')
call grib_count_in_file(ifile,n)
allocate(igrib(n))
!
WeFound=.FALSE.;
DO i=1,n
call grib_new_from_file(ifile, igrib(i))
call grib_get(igrib(i), 'shortName', eShortName)
IF ((TRIM(eShortName) .eq. shortName).and.(WeFound .eqv. .FALSE.)) THEN
IF (GRIB_FILE_TYPE .eq. 1) THEN
CALL READ_GRID_INFO_FROM_GRIB_TYPE1(TheInfo, igrib(i))
END IF
IF (GRIB_FILE_TYPE .eq. 2) THEN
CALL READ_GRID_INFO_FROM_GRIB_TYPE2(TheInfo, igrib(i))
END IF
IF (GRIB_FILE_TYPE .eq. 3) THEN
CALL READ_GRID_INFO_FROM_GRIB_TYPE3(TheInfo, igrib(i))
END IF
WeFound=.TRUE.
END IF
END DO
IF (WeFound .eqv. .FALSE.) THEN
Print *, 'Failed to find the wind variable in the grib file'
CALL WWM_ABORT("Wind has not been found in grib file")
END IF
IF (WRITESTATFLAG == 1) WRITE(STAT%FHNDL, *) 'WeFound=', WeFound
CALL GRIB_CLOSE_FILE(ifile)
deallocate(igrib)
END SUBROUTINE
!****************************************************************************
!* This is functionality for reading GRIB file wind input *
!* Specific supported cases (or wished): ECMWF (IFS), COSMO, DHMZ (ALADIN) *
!****************************************************************************
SUBROUTINE INIT_GRIB_WIND
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
INTEGER IT
INTEGER ifile, i, n
integer, allocatable :: igrib(:)
integer WeFound
REAL(rkind) :: eTimeMjd
integer IPROC, eInt(1)
integer iX, iY
integer nbtime_mjd
integer eProd
character(len=20) shortName
REAL(rkind), allocatable :: wind_time_mjd(:)
REAL cf_scale_factor, cf_add_offset
TYPE(FD_FORCING_GRID) :: TheInfo
integer GRIB_TYPE
! PRint *, 'Begin of INIT_GRIB_WIND'
IF (WRITEWINDBGFLAG == 1) THEN
WRITE(WINDBG%FHNDL, *) 'GRIB_FILE_TYPE=', GRIB_FILE_TYPE
WRITE(WINDBG%FHNDL, *) 'MULTIPLE_IN_WIND=', MULTIPLE_IN_WIND
END IF
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
OPEN(WIN%FHNDL,FILE=WIN%FNAME,STATUS='OLD',ACTION='READ', IOSTAT = ISTAT)
NUM_GRIB_FILES = 0
DO
READ( WIN%FHNDL, *, IOSTAT = ISTAT )
IF ( ISTAT /= 0 ) EXIT
NUM_GRIB_FILES = NUM_GRIB_FILES + 1
END DO
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'NUM_GRIB_FILES=', NUM_GRIB_FILES
REWIND (WIN%FHNDL)
!
ALLOCATE(GRIB_FILE_NAMES(NUM_GRIB_FILES), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 18')
DO IT = 1, NUM_GRIB_FILES
READ(WIN%FHNDL, *) GRIB_FILE_NAMES(IT)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) IT, TRIM(GRIB_FILE_NAMES(IT))
END DO
CLOSE (WIN%FHNDL)
IF (WRITEWINDBGFLAG == 1) FLUSH(WINDBG%FHNDL)
!
nbtime_mjd=NUM_GRIB_FILES
allocate(wind_time_mjd(nbtime_mjd), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 48')
DO IT=1, nbTime_mjd
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL, *) '---------------------------------------'
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL, *) 'IT=', IT, 'file = ', TRIM(GRIB_FILE_NAMES(IT))
CALL READ_TIME_OF_GRIB_FILE(eTimeMjd, GRIB_FILE_NAMES(IT), USE_STEPRANGE)
wind_time_mjd(IT)=eTimeMjd
END DO
IF (WRITEWINDBGFLAG == 1) FLUSH(WINDBG%FHNDL)
cf_scale_factor=ONE
cf_add_offset=ZERO
!
! Now the longitude/latitude to read.
!
shortName='10u'
GRIB_TYPE = GRIB_FILE_TYPE
IT=1
CALL READ_GRID_INFO_FROM_GRIB(TheInfo, GRIB_FILE_NAMES(IT), shortName, GRIB_TYPE)
NDX_WIND_FD = TheInfo % nx_dim
NDY_WIND_FD = TheInfo % ny_dim
allocate(UWIND_FD(NDX_WIND_FD, NDY_WIND_FD), VWIND_FD(NDX_WIND_FD, NDY_WIND_FD), stat=istat)
CALL COMPUTE_CF_COEFFICIENTS(TheInfo)
DEALLOCATE(TheInfo % LON, TheInfo % LAT)
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (.NOT. MULTIPLE_IN_WIND) THEN
IF (myrank .eq. 0) THEN
eInt(1)=nbtime_mjd
DO IPROC=2,nproc
CALL MPI_SEND(eInt,1,itype, iProc-1, 811, comm, ierr)
END DO
DO IPROC=2,nproc
CALL MPI_SEND(WIND_TIME_MJD,nbtime_mjd,rtype, iProc-1, 812, comm, ierr)
END DO
ELSE
CALL MPI_RECV(eInt,1,itype, 0, 811, comm, istatus, ierr)
nbtime_mjd=eInt(1)
ALLOCATE(WIND_TIME_MJD(nbtime_mjd), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_wind, allocate error 3')
CALL MPI_RECV(WIND_TIME_MJD,nbtime_mjd,rtype, 0, 812, comm, istatus, ierr)
END IF
END IF
# endif
eVAR_WIND % cf_scale_factor = cf_scale_factor
eVAR_WIND % cf_add_offset = cf_add_offset
eVAR_WIND % nbTime= nbtime_mjd
allocate(eVAR_WIND % ListTime(nbtime_mjd))
eVAR_WIND % ListTime = WIND_TIME_MJD
! PRint *, 'End of INIT_GRIB_WIND'
END SUBROUTINE
!****************************************************************************
!* The read subroutine *
!****************************************************************************
SUBROUTINE READ_GRIB_WIND(IT, outwind)
USE DATAPOOL
USE GRIB_API
IMPLICIT NONE
integer, intent(in) :: IT
REAL(rkind), INTENT(out) :: outwind(MNP,2)
REAL(rkind) :: outTotal(MNP_WIND,2)
REAL(rkind) :: Vtotal(np_total)
REAL(rkind) :: Vlocal(MNP)
INTEGER ifile, irec, n, iret
integer, allocatable :: igrib(:)
integer WeFoundU, WeFoundV
integer i, j, idx
character(len=100) eShortName
real(rkind) valueU(NDX_WIND_FD*NDY_WIND_FD)
real(rkind) valueV(NDX_WIND_FD*NDY_WIND_FD)
! Print *, 'NDX_WIND_FD=', NDX_WIND_FD
! Print *, 'NDY_WIND_FD=', NDY_WIND_FD
!
! Print *, 'Begin of READ_GRIB_WIND'
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND .or. (myrank .eq. 0)) THEN
# endif
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'IT=', IT, 'file = ', TRIM(GRIB_FILE_NAMES(IT))
! Print *, 'GRIB_FILE_NAMES(IT)=', TRIM(GRIB_FILE_NAMES(IT))
CALL TEST_FILE_EXIST_DIE("Missing grib file: ", TRIM(GRIB_FILE_NAMES(IT)))
CALL GRIB_OPEN_FILE(ifile, TRIM(GRIB_FILE_NAMES(IT)), 'r')
call grib_count_in_file(ifile,n)
IF (WRITEWINDBGFLAG == 1) WRITE(WINDBG%FHNDL,*) 'n=', n
allocate(igrib(n))
WeFoundU=0
WeFoundV=0
DO irec=1,n
call grib_new_from_file(ifile, igrib(irec), iret)
! Print *, 'Step 1'
call grib_get(igrib(irec), 'shortName', eShortName)
! Print *, 'Step 2'
! Print *, 'eShortName=', eShortName
IF ((TRIM(eShortName) .eq. '10u').and.(WeFoundU .eq. 0)) THEN
WeFoundU=1
! Print *, 'Step 3'
CALL grib_get(igrib(irec), 'values', valueU)
! Print *, 'Step 4'
END IF
IF ((TRIM(eShortName) .eq. '10v').and.(WeFoundV .eq. 0)) THEN
WeFoundV=1
! Print *, 'Step 5'
CALL grib_get(igrib(irec), 'values', valueV)
! Print *, 'Step 6'
END IF
END DO
idx=0
DO J=1,NDY_WIND_FD
DO I=1,NDX_WIND_FD
idx=idx+1
UWIND_FD(I,J)=valueU(idx)
VWIND_FD(I,J)=valueV(idx)
END DO
END DO
CALL GRIB_CLOSE_FILE(ifile)
CALL KERNEL_INTERP_UV_WINDFD(outTotal)
# ifdef MPI_PARALL_GRID
END IF
# endif
# ifdef MPI_PARALL_GRID
IF (MULTIPLE_IN_WIND) THEN
outwind=outTotal
ELSE
Vtotal=outTotal(:,1)
CALL SCATTER_ONED_ARRAY(Vtotal, Vlocal)
outwind(:,1)=Vlocal
!
Vtotal=outTotal(:,2)
CALL SCATTER_ONED_ARRAY(Vtotal, Vlocal)
outwind(:,2)=Vlocal
END IF
# else
outwind=outTotal
# endif
! Print *, 'End of READ_GRIB_WIND'
END SUBROUTINE
#endif