MODULE diaharm #if defined key_diaharm && defined key_tide !!================================================================================= !! *** MODULE diaharm *** !! Harmonic analysis of tidal constituents !!================================================================================= !! * Modules used USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain USE in_out_manager ! I/O units USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE ioipsl ! NetCDF IPSL library USE diadimg ! To write dimg USE phycst USE dynspg_oce USE dynspg_ts USE surdetermine USE daymod USE tide_mod USE iom IMPLICIT NONE PRIVATE INTEGER, PARAMETER :: nb_harmo_max=9 LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .TRUE. INTEGER :: & !! namelist variables nit000_han=1, & ! First time step used for harmonic analysis nitend_han=1, & ! Last time step used for harmonic analysis nstep_han=1, & ! Time step frequency for harmonic analysis nb_ana ! Number of harmonics to analyse CHARACTER (LEN=4), DIMENSION(nb_harmo_max) :: & tname ! Names of tidal constituents ('M2', 'K1',...) REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, vt, ut, ft REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta, & out_u, & out_v INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: name !! * Routine accessibility PUBLIC dia_harm ! routine called by step.F90 !!--------------------------------------------------------------------------------- !! !!--------------------------------------------------------------------------------- CONTAINS SUBROUTINE dia_harm_init !!---------------------------------------------------------------------- !! *** ROUTINE dia_harm_init *** !!---------------------------------------------------------------------- !! !! ** Purpose : Initialization of tidal harmonic analysis !! !! ** Method : Initialize frequency array and nodal factor for nit000_han !! !! History : !! 9.0 O. Le Galloudec and J. Chanut (Original) !!-------------------------------------------------------------------- !! * Local declarations INTEGER :: jh, nhan, jk, ji NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, nb_ana, tname !!---------------------------------------------------------------------- ! Read namelist parameters: ! ------------------------- REWIND ( numnam ) READ ( numnam, nam_diaharm ) IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' IF(lwp) WRITE(numout,*) 'First time step used for analysis: nit000_han= ', nit000_han IF(lwp) WRITE(numout,*) 'Last time step used for analysis: nitend_han= ', nitend_han IF(lwp) WRITE(numout,*) 'Time step frequency for harmonic analysis: nstep_han= ', nstep_han IF (nb_ana > nb_harmo_max) THEN IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. & & nb_ana must be lower than nb_harmo_max, stop' IF(lwp) WRITE(numout,*) 'nb_harmo_max= ', nb_harmo_max nstop = nstop + 1 ENDIF ! Basic checks on harmonic analysis time window: ! ---------------------------------------------- IF (nit000 > nit000_han) THEN IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. & & nit000_han must be greater than nit000, stop' IF(lwp) WRITE(numout,*) 'restart capability not implemented' nstop = nstop + 1 ENDIF IF (nitend < nitend_han) THEN IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. & & nitend_han must be lower than nitend, stop' IF(lwp) WRITE(numout,*) 'restart capability not implemented' nstop = nstop + 1 ENDIF IF (MOD(nitend_han-nit000_han+1,nstep_han).NE.0) THEN IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. & & analysis time span must be a multiple of nstep_han, stop' nstop = nstop + 1 END IF CALL tide_init_Wave ALLOCATE(name (nb_ana)) DO jk=1,nb_ana DO ji=1,jpmax_harmo IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN name(jk) = ji EXIT END IF END DO END DO ! Initialize frequency array: ! --------------------------- ALLOCATE(ana_freq(nb_ana)) ALLOCATE(vt (nb_ana)) ALLOCATE(ut (nb_ana)) ALLOCATE(ft (nb_ana)) CALL tide_harmo(ana_freq, vt, ut , ft, name ,nb_ana) IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency ' DO jh = 1, nb_ana IF(lwp) WRITE(numout,*) ' : ',tname(jh),' ',ana_freq(jh) END DO ! Initialize temporary arrays: ! ---------------------------- ALLOCATE( ana_temp(jpi,jpj,nb_ana*2,3)) ana_temp(:,:,:,:) = 0.e0 END SUBROUTINE dia_harm_init SUBROUTINE dia_harm ( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dia_harm *** !!---------------------------------------------------------------------- !! !! ** Purpose : Tidal harmonic analysis main routine !! !! ** Action : Sums ssh/u/v over time analysis [nit000_han,nitend_han] !! !! History : !! 9.0 O. Le Galloudec and J. Chanut (Original) !!-------------------------------------------------------------------- !! * Argument: INTEGER, INTENT( IN ) :: kt !! * Local declarations INTEGER :: ji, jj, jh, jc, nhc REAL(wp) :: ztime, ztemp IF ( kt .EQ. nit000 ) CALL dia_harm_init IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. & (MOD(kt,nstep_han).EQ.0) ) THEN ztime = kt*rdt nhc = 0 DO jh = 1,nb_ana DO jc = 1,2 nhc = nhc+1 ztemp =( MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) & +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) DO jj = 1,jpj DO ji = 1,jpi ! Elevation ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) & + ztemp*sshn(ji,jj)*tmask(ji,jj,1) #if defined key_dynspg_ts ! ubar ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) & + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) ! vbar ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) & + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) #endif END DO END DO END DO END DO END IF IF ( kt .EQ. nitend_han ) CALL dia_harm_end END SUBROUTINE dia_harm SUBROUTINE dia_harm_end !!---------------------------------------------------------------------- !! *** ROUTINE diaharm_end *** !!---------------------------------------------------------------------- !! !! ** Purpose : Compute the Real and Imaginary part of tidal constituents !! !! ** Action : Decompose the signal on the harmonic constituents !! !! History : !! 9.0 O. Le Galloudec and J. Chanut (Original) !!-------------------------------------------------------------------- !! * Local declarations INTEGER :: ji, jj, jh, jc, jn, nhan, jl INTEGER :: ksp, kun, keq REAL(wp) :: ztime, ztime_ini, ztime_end REAL(wp) :: X1,X2 REAL(wp), DIMENSION(jpi,jpj,nb_harmo_max,2) :: ana_amp IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis ninco = 2*nb_ana ksp = 0 keq = 0 DO jn = 1, nhan ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1) keq = keq + 1 kun = 0 DO jh = 1,nb_ana DO jc = 1,2 kun = kun + 1 ksp = ksp + 1 nisparse(ksp) = keq njsparse(ksp) = kun valuesparse(ksp)= & +( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) & +(1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) END DO END DO END DO nsparse=ksp ! Elevation: DO jj = 1, jpj DO ji = 1, jpi ! Fill input array kun=0 DO jh = 1,nb_ana DO jc = 1,2 kun = kun + 1 tmp4(kun)=ana_temp(ji,jj,kun,1) ENDDO ENDDO CALL SUR_DETERMINE(jj) ! Fill output array DO jh = 1, nb_ana ana_amp(ji,jj,jh,1)=tmp7((jh-1)*2+1) ana_amp(ji,jj,jh,2)=tmp7((jh-1)*2+2) END DO END DO END DO ALLOCATE(out_eta(jpi,jpj,2*nb_ana)) ALLOCATE(out_u (jpi,jpj,2*nb_ana)) ALLOCATE(out_v (jpi,jpj,2*nb_ana)) DO jj = 1, jpj DO ji = 1, jpi DO jh = 1, nb_ana X1=ana_amp(ji,jj,jh,1) X2=-ana_amp(ji,jj,jh,2) out_eta(ji,jj,jh)=X1 * tmask(ji,jj,1) out_eta(ji,jj,nb_ana+jh)=X2 * tmask(ji,jj,1) ENDDO ENDDO ENDDO ! ubar: DO jj = 1, jpj DO ji = 1, jpi ! Fill input array kun=0 DO jh = 1,nb_ana DO jc = 1,2 kun = kun + 1 tmp4(kun)=ana_temp(ji,jj,kun,2) ENDDO ENDDO CALL SUR_DETERMINE(jj+1) ! Fill output array DO jh = 1, nb_ana ana_amp(ji,jj,jh,1)=tmp7((jh-1)*2+1) ana_amp(ji,jj,jh,2)=tmp7((jh-1)*2+2) END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi DO jh = 1, nb_ana X1=ana_amp(ji,jj,jh,1) X2=-ana_amp(ji,jj,jh,2) out_u(ji,jj,jh) = X1 * umask(ji,jj,1) out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1) ENDDO ENDDO ENDDO ! vbar: DO jj = 1, jpj DO ji = 1, jpi ! Fill input array kun=0 DO jh = 1,nb_ana DO jc = 1,2 kun = kun + 1 tmp4(kun)=ana_temp(ji,jj,kun,3) ENDDO ENDDO CALL SUR_DETERMINE(jj+1) ! Fill output array DO jh = 1, nb_ana ana_amp(ji,jj,jh,1)=tmp7((jh-1)*2+1) ana_amp(ji,jj,jh,2)=tmp7((jh-1)*2+2) END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi DO jh = 1, nb_ana X1=ana_amp(ji,jj,jh,1) X2=-ana_amp(ji,jj,jh,2) out_v(ji,jj,jh)=X1 * vmask(ji,jj,1) out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1) ENDDO ENDDO ENDDO CALL dia_wri_harm ! Write results in files END SUBROUTINE dia_harm_end SUBROUTINE dia_wri_harm !!-------------------------------------------------------------------- !! *** ROUTINE dia_wri_harm *** !!-------------------------------------------------------------------- !! !! ** Purpose : Write tidal harmonic analysis results in a netcdf file !! !! !! History : !! 9.0 O. Le Galloudec and J. Chanut (Original) !!-------------------------------------------------------------------- !! * Local declarations CHARACTER(LEN=lc) :: cltext CHARACTER(LEN=lc) :: & cdfile_name_T , & ! name of the file created (T-points) cdfile_name_U , & ! name of the file created (U-points) cdfile_name_V ! name of the file created (V-points) INTEGER :: jh !!---------------------------------------------------------------------- #if defined key_dimgout cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc' cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc' cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc' #endif IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' #if defined key_dimgout IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T) IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U) IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V) #endif IF(lwp) WRITE(numout,*) ' ' ! A) Elevation !///////////// ! #if defined key_dimgout cltext='Elevation amplitude and phase' CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2') #else DO jh = 1, nb_ana CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) ) CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) ) END DO #endif ! B) ubar !///////// ! #if defined key_dimgout cltext='ubar amplitude and phase' CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2') #else DO jh = 1, nb_ana CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) ) CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) ) END DO #endif ! C) vbar !///////// ! #if defined key_dimgout cltext='vbar amplitude and phase' CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2') #else DO jh = 1, nb_ana CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh) ) CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,nb_ana+jh) ) END DO #endif END SUBROUTINE dia_wri_harm #else !!---------------------------------------------------------------------- !! Default case : Empty module !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE. CONTAINS SUBROUTINE dia_harm ( kt ) ! Empty routine INTEGER, INTENT( IN ) :: kt WRITE(*,*) 'dia_harm: you should not have seen this print' END SUBROUTINE dia_harm #endif !!====================================================================== END MODULE diaharm