#include PROJECT_HEADER
#include GLOBAL_DEFS
MODULE CWSTWVCP
#if defined COAWST_COUPLING && defined MCT_LIB
!
!svn $Id: waves_coupler.F 756 2008-09-14 20:18:28Z jcwarner $
!==================================================== John C. Warner ===
! Copyright (c) 2002-2008 The ROMS/TOMS Group Hernan G. Arango !
! Licensed under a MIT/X style license !
! See License_ROMS.txt !
!=======================================================================
! !
! This module is used to communicate and exchange data between SWAN !
! other coupled model(s) using the Model Coupling Toolkit (MCT). !
! !
!=======================================================================
!
! Componenet model registry.
!
USE m_MCTWorld, ONLY : MCTWorld_init => init
USE m_MCTWorld, ONLY : MCTWorld_clean => clean
!
! Domain decompositin descriptor datatype and assocoiated methods.
!
USE m_GlobalSegMap, ONLY : GlobalSegMap
USE m_GlobalSegMap, ONLY : GlobalSegMap_init => init
USE m_GlobalSegMap, ONLY : GlobalSegMap_lsize => lsize
USE m_GlobalSegMap, ONLY : GlobalSegMap_clean => clean
USE m_GlobalSegMap, ONLY : GlobalSegMap_Ordpnts => OrderedPoints
!
! Field storage data types and associated methods.
!
USE m_AttrVect, ONLY : AttrVect
USE m_AttrVect, ONLY : AttrVect_init => init
USE m_AttrVect, ONLY : AttrVect_zero => zero
USE m_AttrVect, ONLY : AttrVect_clean => clean
USE m_AttrVect, ONLY : AttrVect_indxR => indexRA
USE m_AttrVect, ONLY : AttrVect_importRAttr => importRAttr
USE m_AttrVect, ONLY : AttrVect_exportRAttr => exportRAttr
!
! Intercomponent communitcations scheduler.
!
USE m_Router, ONLY : Router
USE m_Router, ONLY : Router_init => init
USE m_Router, ONLY : Router_clean => clean
!
! Intercomponent transfer.
!
USE m_Transfer, ONLY : MCT_isend => isend
USE m_Transfer, ONLY : MCT_irecv => irecv
USE m_Transfer, ONLY : MCT_waitr => waitrecv
USE m_Transfer, ONLY : MCT_waits => waitsend
!
# if defined MCT_INTERP_WV2AT
!
! Sparse Matrix DataType and associated methods.
!
USE m_SparseMatrix, ONLY : SparseMatrix
USE m_SparseMatrix, ONLY : SparseMatrix_init => init
USE m_SparseMatrix, ONLY : SparseMatrix_importGRowInd => &
& importGlobalRowIndices
USE m_SparseMatrix, ONLY : SparseMatrix_importGColInd => &
& importGlobalColumnIndices
USE m_SparseMatrix, ONLY : SparseMatrix_importMatrixElts => &
& importMatrixElements
USE m_SparseMatrix, only : SparseMatrix_lsize => lsize
USE m_SparseMatrix, only : SparseMatrix_clean => clean
USE m_SparseMatrixPlus, ONLY : SparseMatrixPlus
USE m_SparseMatrixPlus, ONLY : SparseMatrixPlus_init => init
USE m_SparseMatrixPlus, ONLY : SparseMatrixPlus_clean => clean
!
! Decompose matrix by row.
!
USE m_SparseMatrixPlus, ONLY : Xonly
!
! Matrix-Vector multiply methods.
!
USE m_MatAttrVectMul, ONLY : MCT_MatVecMul => sMatAvMult
# endif
implicit none
!
! PRIVATE
PUBLIC :: INIT_WVCP
PUBLIC :: initialize_wav_routers
PUBLIC :: COAWST_CPL
# ifdef ROMS_COUPLING
PUBLIC :: wav2ocn_coupling
PUBLIC :: wavfocn_coupling
# endif
# ifdef WRF_COUPLING
PUBLIC :: wav2atm_coupling
PUBLIC :: wavfatm_coupling
# endif
PUBLIC :: finalize_wav_coupling
!
! Declarations.
!
TYPE T_GlobalSegMap_G
TYPE(GlobalSegMap) :: GSMapSWAN ! GloabalSegMap variables
END TYPE T_GlobalSegMap_G
TYPE (T_GlobalSegMap_G), ALLOCATABLE :: GlobalSegMap_G(:)
TYPE T_AttrVect_G
# ifdef ROMS_COUPLING
TYPE(AttrVect) :: wav2ocn_AV ! AttrVect variables
TYPE(AttrVect) :: ocn2wav_AV
# endif
# ifdef WRF_COUPLING
TYPE(AttrVect) :: atm2wav_AV ! AttrVec variables
TYPE(AttrVect) :: wav2atm_AV ! AttrVec variables
# endif
END TYPE T_AttrVect_G
TYPE (T_AttrVect_G), ALLOCATABLE :: AttrVect_G(:)
# ifdef ROMS_COUPLING
TYPE T_Router_O
type(Router) :: SWANtoROMS ! Router variables
END TYPE T_Router_O
TYPE (T_Router_O), ALLOCATABLE :: Router_O(:,:)
# endif
# ifdef WRF_COUPLING
TYPE T_GSMapInterp_A
TYPE(GlobalSegMap) :: GSMapWRF ! GloabalSegMap variables
END TYPE T_GSMapInterp_A
TYPE (T_GSMapInterp_A), ALLOCATABLE :: GSMapInterp_A(:,:)
TYPE T_Router_A
type(Router) :: SWANtoWRF ! Router variables
END TYPE T_Router_A
TYPE (T_Router_A), ALLOCATABLE :: Router_A(:,:)
# endif
# ifdef MCT_INTERP_WV2AT
TYPE T_AV2_A
TYPE(AttrVect) :: wav2atm_AV2 ! AttrVect variables
TYPE(AttrVect) :: atm2wav_AV2
END TYPE T_AV2_A
TYPE (T_AV2_A), ALLOCATABLE :: AV2_A(:,:)
TYPE(SparseMatrix) :: sMatW ! Sparse matrix elements
TYPE(SparseMatrix) :: sMatA ! Sparse matrix elements
TYPE T_SMPlus_G
TYPE(SparseMatrixPlus) :: A2WMatPlus ! Sparse matrix plus elements
TYPE(SparseMatrixPlus) :: W2AMatPlus ! Sparse matrix plus elements
END TYPE T_SMPlus_G
TYPE (T_SMPlus_G), ALLOCATABLE :: SMPlus_G(:,:)
# endif
CONTAINS
# ifdef COAWST_COUPLING
!***********************************************************************
! *
SUBROUTINE COAWST_CPL (first)
! *
!***********************************************************************
!
USE MCTWW3PA
!
IMPLICIT NONE
!
INTEGER, intent(in) :: first
INTEGER io, iw, ia, offset, run_couple
# ifdef ROMS_COUPLING
!
! Send data to ocn model.
DO iw=1,Nwav_grids
DO io=1,Nocn_grids
run_couple=1
! IF ((first.eq.1).and.(iics(iw).eq.0)) run_couple=0
IF (MOD(first, nWAV2OCN(1,1)).ne.0) run_couple=0
write(*,*) 'in coawst cpl ', first, run_couple
IF (run_couple.eq.1) THEN
CALL WAV2OCN_COUPLING (iw, io)
ELSE
GOTO 40
END IF
END DO
END DO
40 CONTINUE
!
! Call to get data from ocn model.
DO iw=1,Nwav_grids
DO io=1,Nocn_grids
run_couple=1
! IF ((first.eq.1).and.(iics(iw).eq.0)) run_couple=0
IF (MOD(first, nWAVFOCN(1,1)).ne.0) run_couple=0
IF (run_couple.eq.1) THEN
CALL WAVFOCN_COUPLING (iw, io)
ELSE
GOTO 50
END IF
END DO
END DO
50 CONTINUE
# endif
# ifdef WRF_COUPLING
!
! Call to get data from atm model.
DO iw=1,Nwav_grids
! CALL INIT_POINTERS(iw)
! CALL INIT_COUPLING_POINTERS(iw)
DO ia=1,Natm_grids
run_couple=1
! IF ((first.eq.1).and.(iics(iw).eq.0)) run_couple=0
IF (MOD(first, nWAVFATM(1,1)).ne.0) run_couple=0
IF (run_couple.eq.1) THEN
CALL WAVFATM_COUPLING (iw, ia)
ELSE
GOTO 55
END IF
END DO
END DO
55 CONTINUE
!
! Send data to atm model.
DO iw=1,Nwav_grids
! CALL INIT_POINTERS(iw)
! CALL INIT_COUPLING_POINTERS(iw)
DO ia=1,Natm_grids
run_couple=1
! IF ((first.eq.1).and.(iics(iw).eq.0)) run_couple=0
IF (MOD(first, nWAV2ATM(1,1)).ne.0) run_couple=0
IF (run_couple.eq.1) THEN
CALL WAV2ATM_COUPLING (iw, ia)
ELSE
GOTO 45
END IF
END DO
END DO
45 CONTINUE
# endif
RETURN
END SUBROUTINE COAWST_CPL
# endif
SUBROUTINE INIT_WVCP (ng)
!
!=======================================================================
! !
! Initialize waves and ocean models coupling stream. This is the !
! training phase use to constuct MCT parallel interpolators and !
! stablish communication patterns. !
! !
!=======================================================================
!
USE MCTWW3PA
USE W3GDATMD, ONLY: NX, NY, NSEA
# ifdef MCT_INTERP_WV2AT
! USE mod_coupler_iounits
# endif
!
include 'mpif.h'
integer, intent(in) :: ng
!
! Local variable declarations.
!
integer :: MyError, MyRank
integer :: gsmsize, Nprocs
integer :: i, j, io, ia, Isize, Jsize, Asize
integer :: nRows, nCols, num_sparse_elems
integer :: cid, cad
character (len=70) :: nc_name
character (len=20) :: to_add
character (len=120) :: wostring
character (len=120) :: owstring
real :: cff
! integer, dimension(2) :: src_grid_dims, dst_grid_dims
integer, allocatable :: start(:), length(:)
!
!-----------------------------------------------------------------------
! Begin initialization phase.
!-----------------------------------------------------------------------
!
! Get communicator local rank and size.
!
CALL mpi_comm_rank (WAV_COMM_WORLD, MyRank, MyError)
CALL mpi_comm_size (WAV_COMM_WORLD, Nprocs, MyError)
!
! Initialize MCT coupled model registry.
!
IF (ng.eq.1) THEN
ALLOCATE(GlobalSegMap_G(Nwav_grids))
ALLOCATE(AttrVect_G(Nwav_grids))
# ifdef MCT_INTERP_WV2AT
ALLOCATE(SMPlus_G(Nwav_grids,Natm_grids))
ALLOCATE(AV2_A(Nwav_grids,Natm_grids))
ALLOCATE(GSMapInterp_A(Nwav_grids,Natm_grids))
# endif
END IF
!
WAVid=wavids(ng)
IF (Nwav_grids.gt.1) THEN
CALL MCTWorld_init (N_mctmodels, MPI_COMM_WORLD, &
& WAV_COMM_WORLD,myids=wavids)
ELSE
CALL MCTWorld_init (N_mctmodels, MPI_COMM_WORLD, &
& WAV_COMM_WORLD,WAVid)
END IF
!
! Initialize a Global Segment Map for non-haloed transfer of data for
! SWAN. Determine non-haloed start and length arrays for this
! processor. For now, this will set up a tiled exchange that is
! not identical to thw ww3 tiling.
!
Jsize=INT(NY/Nprocs)
IF (MyRank.eq.Nprocs-1) THEN
Jsize=NY-Jsize*(Nprocs-1)
ENDIF
IF (.not.allocated(start)) THEN
allocate ( start(1) )
END IF
IF (.not.allocated(length)) THEN
allocate ( length(1) )
END IF
start=(MyRank*INT(NY/Nprocs))*NX+1
length=Jsize*NX
gsmsize=length(1)
!
CALL GlobalSegMap_init (GlobalSegMap_G(ng)%GSMapSWAN, start, &
& length, 0, WAV_COMM_WORLD, WAVid)
deallocate (start)
deallocate (length)
# ifdef MCT_INTERP_WV2AT
!
! If wave grid and atm grids are different sizes, then
! develop sparse matrices for interpolation.
!
35 FORMAT(a3,i1,a7,i1,a11)
DO ia=1,Natm_grids
!!!!!!!!!!!!!!!!!!!!!!
! First work on atm to wave.
!!!!!!!!!!!!!!!!!!!!!!
!
IF (MyRank.eq.0) THEN
! IF (scrip_opt.eq.1) THEN
write(nc_name,35) 'atm',ia,'_to_wav',ng,'_weights.nc'
! ELSE
! nc_name=A2Wname(ia,ng)
! END IF
call get_sparse_matrix (ng, nc_name, num_sparse_elems, &
& src_grid_dims, dst_grid_dims)
!
! Init the sparse matrix.
!
nRows=dst_grid_dims(1)*dst_grid_dims(2)
nCols=src_grid_dims(1)*src_grid_dims(2)
!
! Create sparse matrix.
!
! Sparse rows is the dst address. Multiply the interp weights
! by the dst masking.
!
DO i=1,num_sparse_elems
j=sparse_rows(i)
cff=REAL(dst_grid_imask(j),m8)
sparse_weights(i)=sparse_weights(i)*cff
END DO
call SparseMatrix_init(sMatA,nRows,nCols,num_sparse_elems)
call SparseMatrix_importGRowInd(sMatA, sparse_rows, &
& size(sparse_rows))
call SparseMatrix_importGColInd(sMatA, sparse_cols, &
& size(sparse_cols))
call SparseMatrix_importMatrixElts(sMatA, sparse_weights, &
& size(sparse_weights))
!
! Deallocate arrays.
!
deallocate ( sparse_rows )
deallocate ( sparse_cols )
deallocate ( sparse_weights )
deallocate ( dst_grid_imask )
!!!!!!!!!!!!!!!!!!!!!!
! Second work on waves to atm.
!!!!!!!!!!!!!!!!!!!!!!
!
! IF (scrip_opt.eq.1) THEN
write(nc_name,35) 'wav',ng,'_to_atm',ia,'_weights.nc'
! ELSE
! nc_name=W2Aname(ng,ia)
! END IF
call get_sparse_matrix (ng, nc_name, num_sparse_elems, &
& src_grid_dims, dst_grid_dims)
!
! Init the sparse matrix.
!
nRows=dst_grid_dims(1)*dst_grid_dims(2)
nCols=src_grid_dims(1)*src_grid_dims(2)
!
! Create sparse matrix.
!
DO i=1,num_sparse_elems
j=sparse_rows(i)
cff=REAL(dst_grid_imask(j),m8)
sparse_weights(i)=sparse_weights(i)*cff
END DO
!
! Load the dst grid as a coupling mask.
!
allocate(W2A_CPLMASK(ng,ia)%dst_mask(nRows))
DO i=1,nRows
W2A_CPLMASK(ng,ia)%dst_mask(i)=dst_grid_imask(i)
END DO
!
call SparseMatrix_init(sMatW,nRows,nCols,num_sparse_elems)
call SparseMatrix_importGRowInd(sMatW, sparse_rows, &
& size(sparse_rows))
call SparseMatrix_importGColInd(sMatW, sparse_cols, &
& size(sparse_cols))
call SparseMatrix_importMatrixElts(sMatW, sparse_weights, &
& size(sparse_weights))
!
! Deallocate arrays.
!
deallocate ( sparse_rows )
deallocate ( sparse_cols )
deallocate ( sparse_weights )
deallocate ( dst_grid_imask )
END IF
CALL mpi_bcast(dst_grid_dims, 2, MPI_INTEGER, 0, &
& WAV_COMM_WORLD, MyError)
!
! scatter dst_grid_imask to be used as cpl_mask
!
IF (MyRank.ne.0) THEN
nRows=dst_grid_dims(1)*dst_grid_dims(2)
allocate(W2A_CPLMASK(ng,ia)%dst_mask(nRows))
END IF
CALL mpi_bcast(W2A_CPLMASK(ng,ia)%dst_mask, nRows, &
& MPI_INTEGER, 0, &
& WAV_COMM_WORLD, MyError)
!
! Initialize a Global Segment Map for non-haloed transfer of data
! for the atmosphere model.
!
Isize=INT(dst_grid_dims(1)/Nprocs)
IF (MyRank.eq.Nprocs-1) THEN
Isize=dst_grid_dims(1)-Isize*(Nprocs-1)
ENDIF
IF (.not.allocated(start)) THEN
allocate ( start(1) )
END IF
IF (.not.allocated(length)) THEN
allocate ( length(1) )
END IF
start=(MyRank*INT(dst_grid_dims(1)/Nprocs))*dst_grid_dims(2)+1
length=Isize*dst_grid_dims(2)
!
CALL GlobalSegMap_init (GSMapInterp_A(ng,ia)%GSMapWRF, &
& start, length, 0, WAV_COMM_WORLD, WAVid)
deallocate (start)
deallocate (length)
call mpi_barrier(WAV_COMM_WORLD, MyError)
!
! Create ATM sparse matrix plus for interpolation.
! Specify matrix decomposition to be by row.
!
call SparseMatrixPlus_init(SMPlus_G(ng,ia)%A2WMatPlus, sMatA, &
& GSMapInterp_A(ng,ia)%GSMapWRF, &
& GlobalSegMap_G(ng)%GSMapSWAN, &
& Xonly,0,WAV_COMM_WORLD, WAVid)
call SparseMatrix_clean(sMatA)
!
! Create Wave sparse matrix plus for interpolation.
! Specify matrix decomposition to be by row.
!
call SparseMatrixPlus_init(SMPlus_G(ng,ia)%W2AMatPlus, sMatW, &
& GlobalSegMap_G(ng)%GSMapSWAN, &
& GSMapInterp_A(ng,ia)%GSMapWRF, &
& Xonly,0,WAV_COMM_WORLD, WAVid)
call SparseMatrix_clean(sMatW)
END DO
# endif
!
! Initialize attribute vector holding the export data code strings of
! the wave model.
!
cad=LEN(wostring)
DO i=1,cad
wostring(i:i)=''
END DO
cid=1
!
# ifdef ROMS_COUPLING
to_add='DISBOT'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':DISSURF'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':DISWCAP'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':HSIGN'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':RTP'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':TMBOT'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':UBOT'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':DIRE'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
to_add=':DIRN'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':WLEN'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':WLENP'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':QB'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':WDSPR'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':WQP'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
#if defined VEGETATION && defined VEG_SWAN_COUPLING \
&& defined VEG_STREAMING
to_add=':DISVEG'
cad=LEN_TRIM(to_add)
write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
#endif
!
! Finalize and remove trailing spaces from the wostring
! for the rlist.
!
cad=LEN_TRIM(wostring)
wostring=wostring(1:cad)
!
CALL AttrVect_init(AttrVect_G(ng)%wav2ocn_AV, &
& rList=TRIM(wostring),lsize=gsmsize)
CALL AttrVect_zero(AttrVect_G(ng)%wav2ocn_AV)
!
! Initialize attribute vector holding the export data code string of
! the ocean model.
!
cad=LEN(owstring)
DO i=1,cad
owstring(i:i)=''
END DO
cid=1
!
to_add='DEPTH'
cad=LEN_TRIM(to_add)
write(owstring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':WLEV'
cad=LEN_TRIM(to_add)
write(owstring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':VELX'
cad=LEN_TRIM(to_add)
write(owstring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':VELY'
cad=LEN_TRIM(to_add)
write(owstring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
!
to_add=':ZO'
cad=LEN_TRIM(to_add)
write(owstring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
#if defined VEGETATION && defined VEG_SWAN_COUPLING
!
to_add=':VEGDENS'
cad=LEN_TRIM(to_add)
write(owstring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
#endif
# ifdef ICE_MODEL
!
to_add=':SEAICE'
cad=LEN_TRIM(to_add)
write(owstring(cid:cid+cad-1),'(a)') to_add(1:cad)
cid=cid+cad
# endif
!
! Finalize and remove trailing spaces from the owstring
! for the rlist.
!
cad=LEN_TRIM(owstring)
owstring=owstring(1:cad)
!
CALL AttrVect_init (AttrVect_G(ng)%ocn2wav_AV, &
& rList=TRIM(owstring),lsize=gsmsize)
CALL AttrVect_zero (AttrVect_G(ng)%ocn2wav_AV)
!
# endif
# if defined WRF_COUPLING
!
! Initialize attribute vector holding the export data code string of
! the atmosphere model.
!
Asize=GlobalSegMap_lsize(GlobalSegMap_G(ng)%GSMapSWAN, &
& WAV_COMM_WORLD)
CALL AttrVect_init (AttrVect_G(ng)%atm2wav_AV, rlist="U10:V10", &
& lsize=Asize)
CALL AttrVect_zero (AttrVect_G(ng)%atm2wav_AV)
!
! Initialize atribute vector holding wave data to atm.
!
CALL AttrVect_init(AttrVect_G(ng)%wav2atm_AV, &
& rList="HSIGN:WLENP:RTP", &
& lsize=Asize)
CALL AttrVect_zero(AttrVect_G(ng)%wav2atm_AV)
!
# if defined MCT_INTERP_WV2AT
DO ia=1,Natm_grids
! Initialize attribute vector holding the export data code strings of
! the atm model. The Asize is the number of grid point on this
! processor.
!
Asize=GlobalSegMap_lsize(GSMapInterp_A(ng,ia)%GSMapWRF, &
& WAV_COMM_WORLD)
CALL AttrVect_init (AV2_A(ng,ia)%atm2wav_AV2, rlist="U10:V10", &
& lsize=Asize)
CALL AttrVect_zero (AV2_A(ng,ia)%atm2wav_AV2)
!
! Initialize attribute vector holding the export data code string of
! the wave model.
!
CALL AttrVect_init(AV2_A(ng,ia)%wav2atm_AV2, &
& rList="HSIGN:WLENP:RTP:CPL_MASK", &
& lsize=Asize)
CALL AttrVect_zero(AV2_A(ng,ia)%wav2atm_AV2)
END DO
# endif
# endif
RETURN
END SUBROUTINE INIT_WVCP
SUBROUTINE INITIALIZE_WAV_ROUTERS
!
!=======================================================================
! !
! Initialize waves routers for wave model. !
! !
!=======================================================================
!
USE MCTWW3PA
! USE ww3_iounits
! USE M_PARALL
!
! include 'mpif.h'
!
! Local variable declarations.
!
integer :: MyError, MyRank
integer :: ng, iw, ia
!
! Initialize MCT Routers.
!
# ifdef ROMS_COUPLING
ALLOCATE(Router_O(Nwav_grids,Nocn_grids))
!
! Initialize a router to the ocean model component.
!
DO ng=1,Nocn_grids
DO iw=1,Nwav_grids
OCNid=ocnids(ng)
CALL Router_init (OCNid, GlobalSegMap_G(iw)%GSMapSWAN, &
& WAV_COMM_WORLD, Router_O(iw,ng)%SWANtoROMS)
END DO
END DO
# endif
# ifdef WRF_COUPLING
ALLOCATE(Router_A(Nwav_grids,Natm_grids))
!
! Initialize a router to the atmosphere model component.
!
DO ia=1,Natm_grids
DO iw=1,Nwav_grids
ATMid=atmids(ia)
# if defined MCT_INTERP_WV2AT
CALL Router_init (ATMid, GSMapInterp_A(iw,ia)%GSMapWRF, &
& WAV_COMM_WORLD, Router_A(iw,ia)%SWANtoWRF)
# else
CALL Router_init (ATMid, GlobalSegMap_G(iw)%GSMapSWAN, &
& WAV_COMM_WORLD, Router_A(iw,ia)%SWANtoWRF)
# endif
END DO
END DO
# endif
RETURN
END SUBROUTINE INITIALIZE_WAV_ROUTERS
# ifdef ROMS_COUPLING
SUBROUTINE WAV2OCN_COUPLING (iw, io)
!
!=======================================================================
! !
! This subroutine reads and writes the coupling data streams between !
! ocean and wave models. Currently, the following data streams are !
! processed: !
! !
! Fields exported to the OCEAN model: !
! !
! * Wave direction (degrees) !
! * Significant wave height (m) !
! * Average wave length (m) !
! * Surface wave relative peak period (s) !
! * Bottom wave period (s) !
! * Percent of breakig waves (nondimensional) !
! * Wave energy dissipation (W/m2) !
! * Wave bottom orbital velocity (m/s) !
! !
! Fields imported from the OCEAN Model: !
! !
! * Bathymetry, bottom elevation (m) !
! * Free-surface, water surface elevation (m) !
! * Depth integrated u-momentum (m/s) !
! * Depth integrated v-momentum (m/s) !
! !
!=======================================================================
!
USE CONSTANTS, ONLY: PI
USE W3GDATMD, ONLY: NX, NY, NSEA, NSEAL, MAPSF
USE MCTWW3PA
USE W3ADATMD, ONLY: HS, PHIBBL, PHIOC, FP0, T0M1, UBA
USE W3ADATMD, ONLY: THM, WLM
! USE W3ODATMD, ONLY: QB
!
implicit none
!/MPI INCLUDE "mpif.h"
!
! Imported variable declarations.
!
integer :: Numcouple, iw, io
integer :: IP, IX, IY
!
! Local variable declarations.
!
integer :: MyStatus, MyError, MySize, MyRank
integer :: i, id, j, gsmsize, ierr, indx, Tag
integer :: Istr, Iend, Jstr, Jend
integer :: start, length, grdsize
integer :: Isize, Jsize, Nprocs
integer, pointer :: points(:)
real, pointer :: SND_BUF(:), RCV_BUF(:)
real(m8) :: fac
real(m8), pointer :: avdata(:)
real(m8), pointer :: DIRE(:)
real(m8), pointer :: DIRN(:)
!
!-----------------------------------------------------------------------
! Send wave fields to ROMS.
!-----------------------------------------------------------------------
!
CALL MPI_COMM_RANK (WAV_COMM_WORLD, MyRank, MyError)
CALL MPI_COMM_SIZE (WAV_COMM_WORLD, Nprocs, MyError)
!
! Get the number of grid point on this processor.
!
gsmsize=GlobalSegMap_lsize(GlobalSegMap_G(iw)%GSMapSWAN, &
& WAV_COMM_WORLD)
!
! Allocate attribute vector array used to export/import data.
!
allocate ( avdata(gsmsize),stat=ierr )
allocate ( DIRE(gsmsize),stat=ierr )
allocate ( DIRN(gsmsize),stat=ierr )
avdata=0.0_m8
DIRE=0.0_m8
DIRN=0.0_m8
!
grdsize=NX*NY
allocate ( SND_BUF(grdsize),stat=ierr )
SND_BUF=0.0
allocate ( RCV_BUF(grdsize),stat=ierr )
RCV_BUF=0.0
!
! Ask for points in this tile.
!
CALL GlobalSegMap_Ordpnts (GlobalSegMap_G(iw)%GSMapSWAN, &
& MyRank, points)
!
! Determine grid tiling for exchanges.
!
Jsize=INT(NY/Nprocs)
IF (MyRank.eq.Nprocs-1) THEN
Jsize=NY-Jsize*(Nprocs-1)
ENDIF
start=(MyRank*INT(NY/Nprocs))*NX+1
length=Jsize*NX
!
! Load WW3 data into MCT storage buffers.
! The data is exported using ROMS definition for real kind m8=r8.
!
!-------------------------------------------------------------------
! DISBOT: Dissipation Bottom
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=PHIBBL(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "DISBOT",avdata)
!-------------------------------------------------------------------
! DISSURF: Dissipation surface
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=PHIOC(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "DISSURF",avdata)
!-------------------------------------------------------------------
! DISWCAP: Dissipation white capping
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=PHIOC(i)*0. ! jcw need this
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "DISWCAP",avdata)
!-------------------------------------------------------------------
! HS: Signfiicant wave height
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=HS(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "HSIGN",avdata)
!-------------------------------------------------------------------
! RTP: Peak surface period
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=FP0(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
fac=2.0_m8*PI
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=fac/MAX(REAL(RCV_BUF(i),m8),0.001_m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "RTP",avdata)
!-------------------------------------------------------------------
! TMBOT: Mean bottom period
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=T0M1(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "TMBOT",avdata)
!-------------------------------------------------------------------
! UBOT: bottom orbitral velocity
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=UBA(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "UBOT",avdata)
!-------------------------------------------------------------------
! DIRE and DIRN: mean direction, break into 2 components
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=THM(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
! fac=PI/180.0
fac=1.0_m8
DO IP=1,gsmsize
DIRE(IP)=1.0*SIN(avdata(IP)*fac)
DIRN(IP)=1.0*COS(avdata(IP)*fac)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "DIRE",DIRE)
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "DIRN",DIRN)
!-------------------------------------------------------------------
! WLEN: mean wave length
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=WLM(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "WLEN",avdata)
!-------------------------------------------------------------------
! WLENP: peak wave length
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=WLM(i) ! jcw need to make this peak
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "WLENP",avdata)
!-------------------------------------------------------------------
! QB: percent wave break
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
! SND_BUF(IP)=QB(i)
SND_BUF(IP)=0.
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "QB",avdata)
!-------------------------------------------------------------------
! WDSPR: wave dir spread
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
! SND_BUF(IP)=QB(i)*0. ! jcw need this
SND_BUF(IP)=0. ! jcw need this
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "WDSPR",avdata)
!-------------------------------------------------------------------
! WQP:
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
! SND_BUF(IP)=QB(i)*0. ! jcw need this
SND_BUF(IP)=0. ! jcw need this
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2ocn_AV, &
& "WQP",avdata)
!
!-----------------------------------------------------------------------
! Send wave fields bundle to ocean model, ROMS.
!-----------------------------------------------------------------------
!
# ifdef WAVES_OCEAN
Tag=io*100+0*10+iw
CALL MCT_isend (AttrVect_G(iw)%wav2ocn_AV, &
& Router_O(iw,io)%SWANtoROMS, Tag)
CALL MCT_waits (Router_O(iw,io)%SWANtoROMS)
IF (MyRank.EQ.0) THEN
WRITE (SCREEN,36)' == WW3 grid ',iw, &
& ' sent wave data to ROMS grid ', io
36 FORMAT (a14,i2,a29,i2)
END IF
IF (MyError.ne.0) THEN
WRITE (*,*)'coupling send fail swancplr, Error= ', MyError
CALL FINALIZE_WAV_COUPLING(iw)
END IF
# endif
deallocate (avdata, points, DIRE, DIRN)
deallocate (SND_BUF, RCV_BUF)
!
RETURN
END SUBROUTINE WAV2OCN_COUPLING
# endif
# ifdef WRF_COUPLING
SUBROUTINE WAV2ATM_COUPLING (iw, ia)
!
!=======================================================================
! !
! This subroutine reads and writes the coupling data streams between !
! atm and wave models. Currently, the following data streams are !
! processed: !
! !
! Fields exported to the ATM model: !
! !
! * Significant wave height (m) !
! * Surface wave relative peak period (s) !
! * Surface wave length (m) !
! !
! Fields imported from the ATM Model: !
! !
! * Wind Speed (m/s) !
! !
!=======================================================================
!
USE CONSTANTS, ONLY: PI
USE W3GDATMD, ONLY: NX, NY, NSEA, NSEAL, MAPSF
USE MCTWW3PA
USE W3ADATMD, ONLY: HS, PHIBBL, PHIOC, FP0, T0M1, UBA
USE W3ADATMD, ONLY: WLM
!
implicit none
!/MPI INCLUDE "mpif.h"
!
! Imported variable declarations.
!
integer :: Numcouple, iw, ia
integer :: IP, IX, IY
!
! Local variable declarations.
!
integer :: MyStatus, MyError, MySize, MyRank
integer :: i, id, j, gsmsize, ierr, indx, Tag
integer :: Istr, Iend, Jstr, Jend, Asize
integer :: Isize, Jsize, INDXG, Nprocs, OFFSET
integer :: start, length, grdsize
integer, pointer :: points(:)
real, pointer :: SND_BUF(:), RCV_BUF(:)
real(m8) :: fac
real(m8), pointer :: avdata(:)
#ifdef MCT_INTERP_WV2AT
integer, pointer :: indices(:)
real(m8), pointer :: Amask(:)
#endif
!
!-----------------------------------------------------------------------
! Send wave fields to WRF.
!-----------------------------------------------------------------------
!
CALL MPI_COMM_RANK (WAV_COMM_WORLD, MyRank, MyError)
CALL MPI_COMM_SIZE (WAV_COMM_WORLD, Nprocs, MyError)
!
! Get the number of grid point on this processor.
!
gsmsize=GlobalSegMap_lsize(GlobalSegMap_G(iw)%GSMapSWAN, &
& WAV_COMM_WORLD)
!
! Allocate attribute vector array used to export/import data.
!
allocate ( avdata(gsmsize),stat=ierr )
avdata=0.0_m8
!
grdsize=NX*NY
allocate ( SND_BUF(grdsize),stat=ierr )
SND_BUF=0.0
allocate ( RCV_BUF(grdsize),stat=ierr )
RCV_BUF=0.0
!
! Ask for points in this tile.
!
CALL GlobalSegMap_Ordpnts (GlobalSegMap_G(iw)%GSMapSWAN, &
& MyRank, points)
!
! Determine grid tiling for exchanges.
!
Jsize=INT(NY/Nprocs)
IF (MyRank.eq.Nprocs-1) THEN
Jsize=NY-Jsize*(Nprocs-1)
ENDIF
start=(MyRank*INT(NY/Nprocs))*NX+1
length=Jsize*NX
!
! Load WW3 data into MCT storage buffers.
! The data is exported using WRF definition for real kind m8=r8.
!-------------------------------------------------------------------
! HS: Signfiicant wave height
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=HS(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2atm_AV, &
& "HSIGN",avdata)
!-------------------------------------------------------------------
! RTP: Peak surface period
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=FP0(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
fac=2.0_m8*PI
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=fac/MAX(REAL(RCV_BUF(i),m8),0.001_m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2atm_AV, &
& "RTP",avdata)
!-------------------------------------------------------------------
! WLEN: mean wave length
!
! Fill wet parts of array SND_BUF that is NXxNY length.
! The local variable is only 1:NSEAL(M) long.
!
SND_BUF=0.0
DO i=1,NSEAL
IP=(MyRank+1)+(i-1)*Nprocs
IX = MAPSF(IP,1)
IY = MAPSF(IP,2)
IP=(IY-1)*NX+IX
SND_BUF(IP)=WLM(i)
END DO
!
! Gather up all the data.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Now extract the section of data from this tile
! and fill the mct array.
!
IP=0
DO i=start,start+length-1
IP=IP+1
avdata(IP)=REAL(RCV_BUF(i),m8)
END DO
CALL AttrVect_importRAttr (AttrVect_G(iw)%wav2atm_AV, &
& "WLENP",avdata)
!-------------------------------------------------------------------
! Send wave fields bundle to atm model, WRF.
!-----------------------------------------------------------------------
!
35 FORMAT (a14,i2,a23,i5)
!
! Send fields to atmosphere model.
!
Tag=ia*10+iw
# ifdef MCT_INTERP_WV2AT
CALL MCT_MatVecMul(AttrVect_G(iw)%wav2atm_AV, &
& SMPlus_G(iw,ia)%W2AMatPlus, &
& AV2_A(iw,ia)%wav2atm_AV2)
!
! Now add in the CPL_MASK before we send it over to wrf.
! Get the number of grid points on this processor.
!
Asize=GlobalSegMap_lsize (GSMapInterp_A(iw,ia)%GSMapWRF, &
& WAV_COMM_WORLD)
allocate (Amask(Asize))
Amask=0.0
!
! Ask for points in this tile.
!
CALL GlobalSegMap_Ordpnts (GSMapInterp_A(iw,ia)%GSMapWRF, &
& MyRank, points)
!
! Load the dst grid cpl mask into the attr vect.
!
DO i=1,Asize
Amask(i)=REAL(W2A_CPLMASK(iw,ia)%dst_mask(points(i)),m8)
END DO
deallocate (points)
CALL AttrVect_importRAttr (AV2_A(iw,ia)%wav2atm_AV2, "CPL_MASK",&
& Amask, Asize)
CALL MCT_isend (AV2_A(iw,ia)%wav2atm_AV2, &
& Router_A(iw,ia)%SWANtoWRF, Tag)
# else
CALL MCT_isend (AttrVect_G(iw)%wav2atm_AV, &
& Router_A(iw,ia)%SWANtoWRF, Tag)
# endif
CALL MCT_waits (Router_A(iw,ia)%SWANtoWRF)
IF (MyRank.EQ.0) THEN
WRITE (SCREEN,35) '== WW3 grid ',iw,' sent data to WRF grid ' &
& ,ia
END IF
IF (MyError.ne.0) THEN
WRITE (*,*) 'coupl fail swancplr, MyStatus= ', MyError
CALL FINALIZE_WAV_COUPLING(iw)
END IF
#ifdef MCT_INTERP_WV2AT
deallocate (Amask)
if (associated (indices)) then
deallocate (indices)
endif
#endif
deallocate (avdata)
deallocate (SND_BUF, RCV_BUF)
!
RETURN
END SUBROUTINE WAV2ATM_COUPLING
# endif
# ifdef ROMS_COUPLING
SUBROUTINE WAVFOCN_COUPLING (ng, io)
!
!=======================================================================
! !
! This subroutine reads and writes the coupling data streams between !
! ocean and wave models. Currently, the following data streams are !
! processed: !
! !
! Fields exported to the OCEAN model: !
! !
! * Wave direction (degrees) !
! * Significant wave height (m) !
! * Average wave length (m) !
! * Surface wave relative peak period (s) !
! * Bottom wave period (s) !
! * Percent of breakig waves (nondimensional) !
! * Wave energy dissipation (W/m2) !
! * Wave bottom orbital velocity (m/s) !
! !
! Fields imported from the OCEAN Model: !
! !
! * Bathymetry, bottom elevation (m) !
! * Free-surface, water surface elevation (m) !
! * Depth integrated u-momentum (m/s) !
! * Depth integrated v-momentum (m/s) !
! !
!=======================================================================
!
USE W3GDATMD, ONLY: NX, NY, NSEA, NSEAL, MAPSF
USE W3SERVMD
USE W3IDATMD
USE MCTWW3PA
USE W3WDATMD, ONLY: WLV
USE W3ADATMD, ONLY: CX, CY
!
implicit none
!
!/MPI INCLUDE "mpif.h"
!
! Imported variable declarations.
!
integer :: ng, io
!
! Local variable declarations.
!
integer :: IP, IX, IY, grdsize
integer :: MyStatus, MyError, MySize, MyRank, Nprocs
integer :: i, id, j, gsmsize, ierr, indx, Tag
integer :: Istr, Iend, Jstr, Jend, start, length
integer :: Isize, Jsize, INDXG, OFFSET, IAPROC
real :: cff, cffmin, cffmax
# if !defined BBL_MODEL
real :: cff2
# endif
real, parameter :: Large = 1.0E+20
real, dimension(2) :: range
real, pointer :: SND_BUF(:), RCV_BUF(:)
real(m8), pointer :: avdata(:)
!
!-----------------------------------------------------------------------
! Send wave fields to ROMS.
!-----------------------------------------------------------------------
!
CALL MPI_COMM_RANK (WAV_COMM_WORLD, MyRank, MyError)
CALL MPI_COMM_SIZE (WAV_COMM_WORLD, Nprocs, MyError)
IAPROC=MyRank+1
grdsize=NX*NY
!
! Get the number of grid point on this processor.
!
gsmsize=GlobalSegMap_lsize(GlobalSegMap_G(ng)%GSMapSWAN, &
& WAV_COMM_WORLD)
!
! Allocate attribute vector array used to export/import data.
!
allocate ( avdata(gsmsize),stat=ierr )
avdata=0.0_m8
allocate ( SND_BUF(grdsize),stat=ierr )
SND_BUF=0.
allocate ( RCV_BUF(grdsize),stat=ierr )
RCV_BUF=0.
!
!-----------------------------------------------------------------------
! Receive from ROMS: Depth, Water Level, VELX, and VELY.
!-----------------------------------------------------------------------
!
! Schedule receiving field from ocean model.
!
Tag=io*100+0*10+ng
CALL MCT_irecv (AttrVect_G(ng)%ocn2wav_AV, &
& Router_O(ng,io)%SWANtoROMS, Tag)
!
! Wait to make sure the OCN data has arrived.
!
CALL MCT_waitr (AttrVect_G(ng)%ocn2wav_AV, &
& Router_O(ng,io)%SWANtoROMS)
!
IF (MyRank.EQ.0) THEN
WRITE (SCREEN,35) ' == WW3 grid ',ng, &
& ' recv data from ROMS grid ', io
END IF
IF (MyError.ne.0) THEN
WRITE (*,*) 'coupling fail ww3cplr, MyStatus= ', MyError
CALL FINALIZE_WAV_COUPLING(ng)
END IF
35 FORMAT (a14,i2,a26,i2)
!
! Compute local non-halo data size.
!
Jsize=INT(NY/Nprocs)
IF (MyRank.eq.Nprocs-1) THEN
Jsize=NY-Jsize*(Nprocs-1)
ENDIF
start=(MyRank*INT(NY/Nprocs))*NX+1
length=Jsize*NX
!
40 FORMAT (a36,1x,2(1pe14.6))
# ifdef WAVES_OCEAN
!
! Bottom elevation....................................................
!
CALL AttrVect_exportRAttr (AttrVect_G(ng)%ocn2wav_AV, &
& "DEPTH",avdata,gsmsize)
range(1)= Large
range(2)=-Large
SND_BUF=0.0
IP=0
DO i=start,start+length-1
IP=IP+1
range(1)=MIN(range(1),REAL(avdata(IP)))
range(2)=MAX(range(2),REAL(avdata(IP)))
SND_BUF(i)=REAL(avdata(IP))
END DO
CALL MPI_ALLREDUCE(range(1), cffmin, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
CALL MPI_ALLREDUCE(range(2), cffmax, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
IF (MyRank.eq.0) THEN
write(SCREEN,40) 'ROMStoWW3 Min/Max DEPTH (m): ', &
& cffmin, cffmax
END IF
!
! now scatter data to all nodes.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Scatter to array Depth.
!
IF (MyRank.eq.0) THEN
write(SCREEN,*) 'Depth from ROMS to WW3 not coded yet.'
END IF
!
! Water surface elevation............................................
!
CALL AttrVect_exportRAttr (AttrVect_G(ng)%ocn2wav_AV, &
& "WLEV",avdata,gsmsize)
range(1)= Large
range(2)=-Large
SND_BUF=0.0
IP=0
DO i=start,start+length-1
IP=IP+1
range(1)=MIN(range(1),REAL(avdata(IP)))
range(2)=MAX(range(2),REAL(avdata(IP)))
SND_BUF(i)=REAL(avdata(IP))
END DO
CALL MPI_ALLREDUCE(range(1), cffmin, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
CALL MPI_ALLREDUCE(range(2), cffmax, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
IF (MyRank.eq.0) THEN
write(SCREEN,40) 'ROMStoWW3 Min/Max WLEV (m): ', &
& cffmin, cffmax
END IF
!
! now scatter data to all nodes.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Scatter to array WLV.
!
DO i=1,NSEA
IX = MAPSF(i,1)
IY = MAPSF(i,2)
IP=(IY-1)*NX+IX
WLV(i)=RCV_BUF(IP)
END DO
!
! Depth-integrated u-velocity........................................
!
CALL AttrVect_exportRAttr (AttrVect_G(ng)%ocn2wav_AV, &
& "VELX",avdata,gsmsize)
range(1)= Large
range(2)=-Large
SND_BUF=0.0
IP=0
DO i=start,start+length-1
IP=IP+1
range(1)=MIN(range(1),REAL(avdata(IP)))
range(2)=MAX(range(2),REAL(avdata(IP)))
SND_BUF(i)=REAL(avdata(IP))
END DO
CALL MPI_ALLREDUCE(range(1), cffmin, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
CALL MPI_ALLREDUCE(range(2), cffmax, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
IF (MyRank.eq.0) THEN
write(SCREEN,40) 'ROMStoWW3 Min/Max VELX (ms-1): ', &
& cffmin, cffmax
END IF
!
! now scatter data to all nodes.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Scatter to array VELX.
!
DO i=1,NSEA
IX = MAPSF(i,1)
IY = MAPSF(i,2)
IP=(IY-1)*NX+IX
CX(i)=RCV_BUF(IP)
END DO
!
! Depth-integrated v-velocity........................................
!
CALL AttrVect_exportRAttr (AttrVect_G(ng)%ocn2wav_AV, &
& "VELY",avdata,gsmsize)
range(1)= Large
range(2)=-Large
SND_BUF=0.0
IP=0
DO i=start,start+length-1
IP=IP+1
range(1)=MIN(range(1),REAL(avdata(IP)))
range(2)=MAX(range(2),REAL(avdata(IP)))
SND_BUF(i)=REAL(avdata(IP))
END DO
CALL MPI_ALLREDUCE(range(1), cffmin, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
CALL MPI_ALLREDUCE(range(2), cffmax, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
IF (MyRank.eq.0) THEN
write(SCREEN,40) 'ROMStoWW3 Min/Max VELY (ms-1): ', &
& cffmin, cffmax
END IF
!
! now scatter data to all nodes.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Scatter to array VELY.
!
DO i=1,NSEA
IX = MAPSF(i,1)
IY = MAPSF(i,2)
IP=(IY-1)*NX+IX
CY(i)=RCV_BUF(IP)
END DO
!
! Bottom roughness...................................................
!
CALL AttrVect_exportRAttr(AttrVect_G(ng)%ocn2wav_AV, &
& "ZO",avdata,gsmsize)
range(1)= Large
range(2)=-Large
SND_BUF=0.0
IP=0
DO i=start,start+length-1
IP=IP+1
range(1)=MIN(range(1),REAL(avdata(IP)))
range(2)=MAX(range(2),REAL(avdata(IP)))
SND_BUF(i)=REAL(avdata(IP))
END DO
CALL MPI_ALLREDUCE(range(1), cffmin, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
CALL MPI_ALLREDUCE(range(2), cffmax, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
IF (MyRank.eq.0) THEN
write(SCREEN,40) 'ROMStoWW3 Min/Max BottZ0 (m): ', &
& cffmin, cffmax
END IF
!
! now scatter data to all nodes.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Scatter to array Z0.
!
IF (MyRank.eq.0) THEN
write(SCREEN,*) 'Z0 from ROMS to WW3 not coded yet.'
END IF
# endif
!
deallocate (avdata, SND_BUF, RCV_BUF)
!
RETURN
END SUBROUTINE WAVFOCN_COUPLING
# endif
# ifdef WRF_COUPLING
SUBROUTINE WAVFATM_COUPLING (iw, ia)
!
!=======================================================================
! !
! This subroutine reads and writes the coupling data streams between !
! ocean and wave models. Currently, the following data streams are !
! processed: !
! !
! !
! Fields imported from the ATM Model: !
! !
! * Wind E and N (m/s) !
! !
!=======================================================================
!
USE W3GDATMD, ONLY: NX, NY, NSEA, NSEAL, MAPSF
USE W3SERVMD
USE W3IDATMD
USE MCTWW3PA
USE W3ADATMD, ONLY: UA, UD
!
implicit none
!/MPI INCLUDE "mpif.h"
!
! Imported variable declarations.
!
integer :: iw, ia
!
! Local variable declarations.
!
integer :: IP, IX, IY, grdsize
integer :: MyStatus, MyError, MySize, MyRank, Nprocs
integer :: i, id, j, gsmsize, ierr, indx, Tag
integer :: Istr, Iend, Jstr, Jend, start, length
integer :: Isize, Jsize, INDXG, OFFSET, IAPROC
# ifdef MCT_INTERP_WV2AT
integer, pointer :: indices(:)
# endif
real :: cff, cffmin, cffmax
real, parameter :: Large = 1.0E+20
real, dimension(2) :: range
real, pointer :: SND_BUF(:), RCV_BUF(:)
real, pointer :: WND_U10(:), WND_V10(:)
real(m8), pointer :: avdata(:)
!
!-----------------------------------------------------------------------
! Get wind data from atm.
!-----------------------------------------------------------------------
!
CALL MPI_COMM_RANK (WAV_COMM_WORLD, MyRank, MyError)
CALL MPI_COMM_SIZE (WAV_COMM_WORLD, Nprocs, MyError)
IAPROC=MyRank+1
grdsize=NX*NY
!
! Get the number of grid point on this processor.
!
gsmsize=GlobalSegMap_lsize(GlobalSegMap_G(iw)%GSMapSWAN, &
& WAV_COMM_WORLD)
!
! Allocate attribute vector array used to export/import data.
!
allocate ( avdata(gsmsize),stat=ierr )
avdata=0.0_m8
MyError=0
allocate ( SND_BUF(grdsize),stat=ierr )
SND_BUF=0.
allocate ( RCV_BUF(grdsize),stat=ierr )
RCV_BUF=0.
allocate ( WND_U10(grdsize),stat=ierr )
WND_U10=0.
allocate ( WND_V10(grdsize),stat=ierr )
WND_V10=0.
!
!-----------------------------------------------------------------------
! RCV the data from WRF.
!-----------------------------------------------------------------------
!
35 FORMAT (a14,i2,a24,i2)
!
! Receive fields from atmosphere model.
!
Tag=0*100+ia*10+iw
# ifdef MCT_INTERP_WV2AT
CALL MCT_irecv (AV2_A(iw,ia)%atm2wav_AV2, &
& Router_A(iw,ia)%SWANtoWRF, Tag)
! Wait to make sure the WRF data has arrived.
CALL MCT_waitr (AV2_A(iw,ia)%atm2wav_AV2, &
& Router_A(iw,ia)%SWANtoWRF)
CALL MCT_MatVecMul(AV2_A(iw,ia)%atm2wav_AV2, &
& SMPlus_G(iw,ia)%A2WMatPlus, &
& AttrVect_G(iw)%atm2wav_AV)
# else
CALL MCT_irecv (AttrVect_G(iw)%atm2wav_AV, &
& Router_A(iw,ia)%SWANtoWRF, Tag)
! Wait to make sure the WRF data has arrived.
CALL MCT_waitr (AttrVect_G(iw)%atm2wav_AV, &
& Router_A(iw,ia)%SWANtoWRF)
# endif
IF (MyRank.EQ.0) THEN
WRITE (SCREEN,35)'== WW3 grid ',iw,' recv data from WRF grid' &
& ,ia
END IF
IF (MyError.ne.0) THEN
WRITE (*,*) 'coupling fail swancplr, MyStatus= ', MyError
CALL FINALIZE_WAV_COUPLING(iw)
END IF
!
! Compute local non-halo data size.
!
Jsize=INT(NY/Nprocs)
IF (MyRank.eq.Nprocs-1) THEN
Jsize=NY-Jsize*(Nprocs-1)
ENDIF
start=(MyRank*INT(NY/Nprocs))*NX+1
length=Jsize*NX
!
# ifdef AIR_WAVES
!
! U-wind.
!
CALL AttrVect_exportRAttr(AttrVect_G(iw)%atm2wav_AV,"U10", &
& avdata,gsmsize)
range(1)= Large
range(2)=-Large
SND_BUF=0.0
IP=0
DO i=start,start+length-1
IP=IP+1
range(1)=MIN(range(1),REAL(avdata(IP)))
range(2)=MAX(range(2),REAL(avdata(IP)))
SND_BUF(i)=REAL(avdata(IP))
END DO
CALL MPI_ALLREDUCE(range(1), cffmin, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
CALL MPI_ALLREDUCE(range(2), cffmax, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
IF (MyRank.eq.0) THEN
write(SCREEN,40) 'WRFtoWW3 Min/Max U10 (ms-1): ', &
& cffmin, cffmax
END IF
!
! now scatter data to all nodes.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Scatter to array WND_U10 as temporary for now.
!
DO i=1,grdsize
WND_U10(i)=RCV_BUF(i)
END DO
! DO i=1,NSEA
! IX = MAPSF(i,1)
! IY = MAPSF(i,2)
! IP=(IY-1)*NX+IX
! WND_SPD(i)=RCV_BUF(IP)
! END DO
!
! V-wind.
!
CALL AttrVect_exportRAttr(AttrVect_G(iw)%atm2wav_AV,"V10", &
& avdata,gsmsize)
range(1)= Large
range(2)=-Large
SND_BUF=0.0
IP=0
DO i=start,start+length-1
IP=IP+1
range(1)=MIN(range(1),REAL(avdata(IP)))
range(2)=MAX(range(2),REAL(avdata(IP)))
SND_BUF(i)=REAL(avdata(IP))
END DO
CALL MPI_ALLREDUCE(range(1), cffmin, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
CALL MPI_ALLREDUCE(range(2), cffmax, 1, MPI_REAL, &
MPI_SUM, WAV_COMM_WORLD, MyError)
IF (MyRank.eq.0) THEN
write(SCREEN,40) 'WRFtoWW3 Min/Max V10 (ms-1): ', &
& cffmin, cffmax
END IF
!
! now scatter data to all nodes.
!
CALL MPI_ALLREDUCE(SND_BUF, RCV_BUF, grdsize, &
& MPI_REAL, MPI_SUM, WAV_COMM_WORLD, MyError)
!
! Scatter to array WND_V10 as temporary for now.
!
DO i=1,grdsize
WND_V10(i)=RCV_BUF(i)
END DO
!
! Now we need to combine wnd speed and dir and scatter.
!
DO i=1,NSEA
IX = MAPSF(i,1)
IY = MAPSF(i,2)
IP=(IY-1)*NX+IX
UA(i)=SQRT(WND_U10(IP)**2+WND_V10(IP)**2+0.000001)
UD(i)=ATAN2(WND_V10(IP),WND_U10(IP))
END DO
40 FORMAT (a36,1x,2(1pe14.6))
# endif
!
deallocate (avdata)
deallocate (SND_BUF, RCV_BUF)
deallocate (WND_U10, WND_V10)
# ifdef MCT_INTERP_WV2AT
if (associated (indices)) then
deallocate (indices)
endif
# endif
!
RETURN
END SUBROUTINE WAVFATM_COUPLING
# endif
SUBROUTINE FINALIZE_WAV_COUPLING(iw)
!
!=======================================================================
! ===
! This routines terminates execution during coupling error. ===
! ===
!=======================================================================
USE MCTWW3PA
!
! Local variable declarations.
!
integer :: iw, io, ia, MyError
!
!-----------------------------------------------------------------------
! Deallocate MCT environment.
!-----------------------------------------------------------------------
!
# ifdef ROMS_COUPLING
DO io=1,Nocn_grids
CALL Router_clean (Router_O(iw,io)%SWANtoROMS, MyError)
END DO
CALL AttrVect_clean (AttrVect_G(iw)%wav2ocn_AV, MyError)
# endif
# ifdef WRF_COUPLING
DO ia=1,Natm_grids
CALL Router_clean (Router_A(iw,ia)%SWANtoWRF, MyError)
END DO
CALL AttrVect_clean (AttrVect_G(iw)%atm2wav_AV, MyError)
# endif
CALL GlobalSegMap_clean (GlobalSegMap_G(iw)%GSMapSWAN, MyError)
END SUBROUTINE FINALIZE_WAV_COUPLING
#endif
END MODULE CWSTWVCP