New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diaharm_fast.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaharm_fast.F90 @ 15500

Last change on this file since 15500 was 15500, checked in by hadjt, 8 months ago

diaharm_fast.F90

Found an important bug. it wasn't set up for monthly cycles, so was only reseting the "seconds of the day" count once a month. Therefore, the v0tide was introducing a jump every day, so was fitting poorly.

It is still not working properly, as the amplitudes output are too high, and the phase maps are not correct.

File size: 45.2 KB
Line 
1MODULE diaharm_fast 
2   !!======================================================================
3   !!                       ***  MODULE  example  ***
4   !! Ocean physics:  On line harmonic analyser
5   !!                 
6   !!=====================================================================
7
8#if defined key_diaharm_fast 
9
10   !!----------------------------------------------------------------------
11   !!   'key_harm_ana'  :                Calculate harmonic analysis
12   !!----------------------------------------------------------------------
13   !!   harm_ana        :
14   !!   harm_ana_init   :
15   !!   NB: 2017-12 : add 3D harmonic analysis of velocities
16   !!                 integration of Maria Luneva's development
17   !!   'key_3Ddiaharm'
18   !!----------------------------------------------------------------------
19
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE iom
23   USE in_out_manager  ! I/O units
24   USE phycst          ! physical constants
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE bdy_oce         ! ocean open boundary conditions
27   USE bdytides        ! tidal bdy forcing
28   USE daymod          ! calendar
29   USE tideini
30   USE restart
31   USE ioipsl, ONLY : ju2ymds    ! for calendar
32   !
33   !
34   USE timing          ! preformance summary
35   USE zdf_oce
36
37   IMPLICIT NONE
38   PRIVATE
39
40   !! *  Routine accessibility
41   PUBLIC dia_harm_fast                                      ! routine called in step.F90 module
42   LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast  = .TRUE.   ! to be run or not
43   LOGICAL, PUBLIC :: lk_diaharm_2D   ! = .TRUE.   ! to run 2d
44   LOGICAL, PUBLIC :: lk_diaharm_3D   ! = .TRUE.   ! to run 3d
45
46   !! * Module variables
47   INTEGER, PARAMETER ::  nharm_max  = jpmax_harmo  ! max number of harmonics to be analysed
48   INTEGER, PARAMETER ::  nhm_max    = 2*nharm_max+1 
49   INTEGER, PARAMETER ::  nvab       = 2 ! number of 3D variables
50   INTEGER            ::  nharm
51   INTEGER            ::  nhm 
52   INTEGER ::                 & !!! ** toto namelist (namtoto) **
53      nflag  =  1                ! default value of nflag
54   REAL(wp), DIMENSION(nharm_max) ::                & 
55      om_tide                     ! tidal frequencies ( rads/sec)
56   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:)   ::                & 
57      bzz,c,x    ! work arrays
58   REAL(wp) :: cca,ssa,zm,bt,dd_cumul
59!
60   REAL(wp), PUBLIC ::   fjulday_startharm       !: Julian Day since start of harmonic analysis
61   REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf   ! nodel/phase corrections used by diaharmana
62   REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:)   :: cc,a
63!
64   INTEGER ::  nvar_2d, nvar_3d    !: number of 2d and 3d variables to analyse
65   INTEGER, ALLOCATABLE,DIMENSION(:) :: m_posi_2d, m_posi_3d
66
67!  Name of variables used in the restart
68   CHARACTER( LEN = 10 ), DIMENSION(5), PARAMETER :: m_varName2d = (/'ssh','u2d','v2d','ubfr','vbfr'/)
69   CHARACTER( LEN = 10 ), DIMENSION(4), PARAMETER :: m_varName3d = (/'rho','u3d','v3d','w3d'/)
70!
71   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:  ) :: g_cosamp2D, g_sinamp2D, g_cumul_var2D
72   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:) :: g_cosamp3D, g_sinamp3D, g_cumul_var3D 
73!
74   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:)       :: g_out2D,h_out2D  ! arrays for output
75   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:)     :: g_out3D,h_out3D  ! arrays for 3D output
76!
77!  NAMELIST
78   LOGICAL, PUBLIC :: ln_diaharm_store           !: =T  Stores data for harmonic Analysis
79   LOGICAL, PUBLIC :: ln_diaharm_compute         !: =T  Compute harmonic Analysis
80   LOGICAL, PUBLIC :: ln_diaharm_read_restart   !: =T  Read restart from a previous run
81   !JT
82   LOGICAL, PUBLIC :: ln_diaharm_multiyear   !: =T  Read restart from a previous run
83   INTEGER, PUBLIC :: nn_diaharm_multiyear   !: =T  Read restart from a previous run
84   LOGICAL, PUBLIC :: ln_diaharm_update_nodal_daily   !: =T  update the nodes every day
85   LOGICAL, PUBLIC :: ln_diaharm_fast
86
87   !JT
88   LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
89   INTEGER ::   nb_ana        ! Number of harmonics to analyse
90   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
91   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_all ! INDEX within the full set of constituents (tide.h90)
92   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_sub ! INDEX within the subset of constituents pass in input
93
94   !! * Substitutions
95
96   !!----------------------------------------------------------------------
97   !!  OPA 9.0 , LOCEAN-IPSL (2005)
98   !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
99   !! or  TOP 1.0 , LOCEAN-IPSL (2005)
100   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
101   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
102   !!----------------------------------------------------------------------
103
104CONTAINS
105
106   SUBROUTINE dia_harm_fast( kt )
107      !!----------------------------------------------------------------------
108      !!                    ***  ROUTINE harm_ana  ***
109      !!
110      !! ** Purpose :   Harmonic analyser
111      !!
112      !! ** Method  :   
113      !!
114      !! ** Action  : - first action (share memory array/varible modified
115      !!                in this routine
116      !!              - second action .....
117      !!              - .....
118      !!
119      !! References :
120      !!   Give references if exist otherwise suppress these lines
121      !!
122      !! History :
123      !!   9.0  !  03-08  (Autor Names)  Original code
124      !!        !  02-08  (Author names)  brief description of modifications
125      !!----------------------------------------------------------------------
126      !! * Modules used
127     
128      !! * arguments
129      INTEGER, INTENT( in  ) ::   & 
130         kt                          ! describe it!!!
131
132      !! * local declarations
133      INTEGER  :: ji, jk, jj          ! dummy loop arguments
134      INTEGER  :: jh, i1, i2, jgrid
135      INTEGER  :: j2d, j3d
136      REAL(WP) :: sec2start,sec2start_old
137      CHARACTER (len=40) :: tmp_name
138      !!--------------------------------------------------------------------
139
140      !JT IF( nn_timing == 1 )   CALL timing_start( 'dia_harm_fast' )
141      IF( ln_timing )   CALL timing_start( 'dia_harm_fast' )
142      IF( kt == nit000   )   CALL harm_ana_init(kt)    ! Initialization (first time-step only)
143
144
145      IF ( ln_diaharm_update_nodal_daily ) THEN
146         !IF (MOD(kt,nint(86400./rdt)) == 0) THEN
147
148
149         IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN     
150            DO jh = 1, nb_ana
151               !JT anau(jh) = 3.141579*utide ( ntide_sub(jh) )/180.
152               !JT anav(jh) = 3.141579*v0tide( ntide_sub(jh) )/180.
153               anau(jh) = utide ( ntide_sub(jh) )
154               anav(jh) = v0tide( ntide_sub(jh) )
155               anaf(jh) = ftide ( ntide_sub(jh) )
156            END DO
157
158            IF(lwp) WRITE(numout,*)
159            IF(lwp) WRITE(numout,*) 'harm_ana : update nodes?',ln_diaharm_update_nodal_daily
160            IF(lwp) WRITE(numout,*) 'harm_ana : date, time',ndastp, nhour, nminute
161
162         ENDIF
163      ENDIF
164
165       
166    DO jh = 1, nb_ana       
167
168        tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anau_utide'
169        IF( iom_use(TRIM(tmp_name)) )  THEN
170        !    IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anau(jh) )
171            CALL iom_put( TRIM(tmp_name), anau(jh) )
172        !ELSE
173        !    IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
174        ENDIF   
175
176        tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anav_v0tide'
177        IF( iom_use(TRIM(tmp_name)) )  THEN
178        !    IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anav(jh) )
179            CALL iom_put( TRIM(tmp_name), anav(jh) )
180        !ELSE
181        !    IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
182        ENDIF
183
184        tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anaf_ftide'
185        IF( iom_use(TRIM(tmp_name)) )  THEN
186        !    IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anaf(jh) )
187            CALL iom_put( TRIM(tmp_name), anaf(jh) )
188        !ELSE
189        !    IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
190        ENDIF
191
192     END DO
193       
194
195     IF ( ln_diaharm_fast .and. ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
196
197          ! this bit done every time step
198          nhm=2*nb_ana+1
199          c(1) = 1.0
200
201          sec2start_old = nint( (fjulday-fjulday_startharm)*86400._wp ) 
202          sec2start = nsec_day - NINT(0.5_wp * rdt)
203          !IF(lwp) WRITE(numout,*) "ztime NEW", kt, sec2start, fjulday_startharm
204
205          DO jh=1,nb_ana
206             c(2*jh  ) = anaf(jh)*cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
207             c(2*jh+1) = anaf(jh)*sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
208             
209             IF(lwp) WRITE(numout,*) 'diaharm_fast: analwave,',kt,tname(jh),sec2start,sec2start/3600.,sec2start_old,sec2start_old/3600,c(2*jh),c(2*jh+1),om_tide(jh),anau(jh),anav(jh)
210          ENDDO 
211
212          !IF(lwp) WRITE(numout,*) "c init", c, "c end", sec2start, om_tide(1), anau(1), anav(1),"end nodal"
213
214
215          ! CUMULATE
216          DO ji=1,jpi         ! loop lon
217             DO jj=1,jpj      ! loop lat
218                DO jh=1,nhm   ! loop harmonic
219
220                   DO j2d=1,nvar_2d
221                      IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj)             ! analysis elevation
222                      IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj)             ! analysis depth average velocities
223                      IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj)
224                      !JT IF ( m_posi_2d(j2d) .eq. 4 ) dd_cumul = c(jh) * bfrua(ji,jj) * un(ji,jj,mbku(ji,jj)) * ssumask(ji,jj) ! analysis bottom friction
225                      !JT IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj)
226                      g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul
227                   ENDDO
228
229                   DO j3d=1,nvar_3d
230                      DO jk=1,jpkm1
231                         IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) *  rhd(ji,jj,jk)               * tmask(ji,jj,jk)   
232                         IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk) 
233                         IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk)
234                         IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) *   wn(ji,jj,jk)               * wmask(ji,jj,jk)
235                         g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul
236                      ENDDO
237                   ENDDO
238
239                ENDDO     ! end loop harmonic
240             ENDDO        ! end loop lat
241          ENDDO           ! end loop lon
242
243          ! Compute nodal factor cumulative cross-product
244          DO i1=1,nhm
245             DO i2=1,nhm
246                cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
247             ENDDO
248          ENDDO
249
250          ! Output RESTART
251          IF( kt == nitrst ) THEN
252             CALL harm_rst_write(kt) ! Dump out data for a restarted run
253          ENDIF
254
255          ! At End of run
256          IF ( kt ==  nitend ) THEN
257
258             IF(lwp) WRITE(numout,*)
259             IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
260             IF(lwp) WRITE(numout,*) '~~~~~~~~~'
261
262             IF( ln_diaharm_compute ) THEN
263
264                 ! INITIALISE TABLE TO 0
265                 IF ( nvar_2d .gt. 0 ) THEN
266                    g_cosamp2D = 0.0_wp
267                    g_sinamp2D = 0.0_wp
268                 ENDIF
269                 IF ( nvar_3d .gt. 0 ) THEN
270                    g_cosamp3D = 0.0_wp
271                    g_sinamp3D = 0.0_wp
272                 ENDIF
273
274                 ! FIRST OUTPUT 2D VARIABLES
275                 DO jgrid=1,nvar_2d    ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically
276                    DO ji=1,jpi        ! loop lon
277                       DO jj=1,jpj     ! loop lat
278                          bt = 1.0_wp; bzz(:) = 0.0_wp
279                          DO jh=1,nhm  ! loop harmonic
280                             bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid)
281                             bt = bt*bzz(jh)
282                          ENDDO
283                          ! Copy back original cumulated nodal factor
284                          a(:,:) = cc(:,:)
285    !                     now do gaussian elimination of the system
286    !                     a * x = b
287    !                     the matrix x is (a0,a1,b1,a2,b2 ...)
288    !                     the matrix a and rhs b solved here for x
289                          x=0.0_wp
290                          IF(bt.ne.0.) THEN
291                            CALL gelim( a, bzz, x, nhm )
292    !                       Backup output in variables
293                            DO jh=1,nb_ana
294                               g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2  )
295                               g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1)
296                            ENDDO
297                            g_cosamp2D   ( 0,ji,jj,jgrid) = x(1)
298                            g_sinamp2D   ( 0,ji,jj,jgrid) = 0.0_wp
299                          ENDIF     ! bt.ne.0.
300                       ENDDO        ! jj
301                    ENDDO           ! ji
302                 ENDDO              ! jgrid
303
304                 ! SECOND OUTPUT 3D VARIABLES
305                 DO jgrid=1,nvar_3d     ! loop number of 3d variables rho, U, V, W
306                    DO jk=1,jpkm1       ! loop over vertical level
307                       DO ji=1,jpi      ! loop over lon
308                          DO jj=1,jpj   ! loop over lat
309                             bt = 1.0_wp; bzz(:) = 0.0_wp
310                             DO jh=1,nhm
311                                bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid)
312                                bt = bt*bzz(jh)
313                             ENDDO
314                             ! Copy back original cumulated nodal factor
315                             a(:,:) = cc(:,:)                     
316    !                        now do gaussian elimination of the system
317    !                        a * x = b
318    !                        the matrix x is (a0,a1,b1,a2,b2 ...)
319    !                        the matrix a and rhs b solved here for x
320                             x=0.0_wp
321                             IF(bt.ne.0.) THEN
322                               CALL gelim( a, bzz, x, nhm )
323    !                          Backup output in variables
324                               DO jh=1,nb_ana
325                                  g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2  )
326                                  g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1)
327                               ENDDO
328                               g_cosamp3D   ( 0,ji,jj,jk,jgrid) = x(1)
329                               g_sinamp3D   ( 0,ji,jj,jk,jgrid) = 0.0_wp
330                            ENDIF     ! bt.ne.0.
331                          ENDDO       ! jj
332                       ENDDO          ! ji
333                    ENDDO             ! jk
334                 ENDDO                ! jgrid
335
336                 CALL harm_ana_out     ! output analysis (last time step)
337
338             ELSE    ! ln_harmana_compute = False
339                 IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step"
340
341             ENDIF   ! ln_harmana_compute
342          ENDIF      ! kt ==  nitend
343
344     ENDIF
345
346      !JT IF( nn_timing == 1 )   CALL timing_stop( 'dia_harm_fast' )
347      IF( ln_timing  )   CALL timing_stop( 'dia_harm_fast' )
348
349   END SUBROUTINE dia_harm_fast 
350
351   SUBROUTINE harm_ana_init( kt )   !JT
352      !!----------------------------------------------------------------------
353      !!                  ***  ROUTINE harm_ana_init  ***
354      !!                   
355      !! ** Purpose :   initialization of ....
356      !!
357      !! ** Method  :   blah blah blah ...
358      !!
359      !! ** input   :   Namlist namexa
360      !!
361      !! ** Action  :   ... 
362      !!
363      !! history :
364      !!   9.0  !  03-08  (Autor Names)  Original code
365      !!----------------------------------------------------------------------
366      !! * local declarations
367
368      INTEGER, INTENT(in) ::   kt     ! ocean time-step  !JT
369      !!
370      INTEGER ::   ji, jk, jh  ! dummy loop indices
371      INTEGER ::   ios                  ! Local integer output status for namelist read
372      INTEGER ::   k2d, k3d             ! dummy number of analysis
373
374
375     
376
377      NAMELIST/nam_diaharm_fast/ ln_diaharm_fast, ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, &
378               & tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily
379      !!----------------------------------------------------------------------
380      !JT
381      ln_diaharm_fast = .FALSE.
382      ln_diaharm_multiyear = .FALSE.
383      nn_diaharm_multiyear = 20
384      !JT
385      lk_diaharm_2D    = .TRUE.   ! to run 2d
386      lk_diaharm_3D    = .TRUE.   ! to run 3d
387
388      ln_diaharm_store = .TRUE.
389
390      IF(lwp) WRITE(numout,*)
391      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
392      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
393
394      ! GET NAMELIST DETAILS
395      REWIND( numnam_ref )              ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
396      READ  ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
397901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist' )
398
399      REWIND( numnam_cfg )              ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
400      READ  ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
401902   IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist' )
402      IF(lwm) WRITE ( numond, nam_diaharm_fast )
403
404
405      !
406      IF(lwp) THEN
407         WRITE(numout,*) 'Tidal diagnostics_fast'
408         WRITE(numout,*) '   Fast Harmonic analysis ?:         ln_diaharm_fast= ', ln_diaharm_fast
409         WRITE(numout,*) '   Store output in restart?:         ln_diaharm_store= ', ln_diaharm_store
410         WRITE(numout,*) '   Compute analysis?:         ln_diaharm_compute= ', ln_diaharm_compute
411         WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
412         WRITE(numout,*) '   SSH harmonic analysis: ln_ana_ssh = ', ln_ana_ssh
413         WRITE(numout,*) '   Barotropic Velocities harmonic analysis: ln_ana_uvbar = ', ln_ana_uvbar
414         WRITE(numout,*) '   Bed Friction for harmonic analysis (not implemented): ln_ana_bfric = ', ln_ana_bfric
415         WRITE(numout,*) '   Density harmonic analysis: ln_ana_rho = ', ln_ana_rho
416         WRITE(numout,*) '   3D velocities harmonic analysis: ln_ana_uv3d = ', ln_ana_uv3d
417         WRITE(numout,*) '   Vertical Velocities harmonic analysis: ln_ana_w3d = ', ln_ana_w3d
418         WRITE(numout,*) '   Names of harmonics: tname = ', tname
419         WRITE(numout,*) '   Max number of harmonics: jpmax_harmo = ', jpmax_harmo
420         WRITE(numout,*) '   Number of Harmonics: nb_harmo = ', nb_harmo
421         WRITE(numout,*) '   Multi-year harmonic analysis: ln_diaharm_multiyear = ', ln_diaharm_multiyear
422         WRITE(numout,*) '   Multi-year harmonic analysis - number of years: nn_diaharm_multiyear = ', nn_diaharm_multiyear
423         WRITE(numout,*) '   Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily
424         WRITE(numout,*) '   Number of Harmonics: nyear, nmonth = ', nyear, nmonth
425
426      ENDIF
427      ! JT
428
429
430      IF ( ln_diaharm_multiyear ) THEN
431          ln_diaharm_store = .True.
432          IF(lwp) WRITE(numout,*) '   Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
433          IF ((mod(nyear,nn_diaharm_multiyear) == 0) .AND. ( nmonth == 1 )) THEN ! Jan, year = 1980,2000,2020,2040, restart tidal calculation
434              ln_diaharm_read_restart = .FALSE.
435              IF(lwp) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
436          ELSE
437              ln_diaharm_read_restart = .TRUE.
438              IF(lwp) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
439          ENDIF
440
441
442
443          IF(lwp) WRITE(numout,*) '   Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
444          IF ((mod(nyear,nn_diaharm_multiyear) == (nn_diaharm_multiyear - 1)) .AND. ( nmonth == 12 )) THEN ! Dec year = 1999,2019,2039,2040, restart tidal calculation
445              ln_diaharm_compute = .TRUE.
446              IF(lwp) WRITE(numout,*) '   Compute analysis?:         ln_diaharm_compute= ', ln_diaharm_compute
447          ELSE
448              ln_diaharm_compute = .FALSE.
449              IF(lwp) WRITE(numout,*) '   Compute analysis?:         ln_diaharm_compute= ', ln_diaharm_compute
450          ENDIF
451
452      ENDIF
453      IF ( kt < 10 ) THEN
454        ln_diaharm_read_restart = .FALSE.
455        IF(lwp) WRITE(numout,*) '   kt = ',kt
456        IF(lwp) WRITE(numout,*) '   kt < 10, so setting ln_diaharm_read_restart to .FALSE.'
457        IF(lwp) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart       
458      ENDIF
459
460      ! JT
461
462
463      ! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
464      nb_ana = 0
465      DO jk=1,jpmax_harmo
466         DO ji=1,nb_harmo
467            IF(TRIM(tname(jk)) == TRIM(Wave( ntide(ji) )%cname_tide) ) THEN
468               nb_ana=nb_ana+1
469            ENDIF
470         END DO
471      END DO
472      !
473      IF(lwp) THEN
474         WRITE(numout,*) '        Namelist nam_diaharm_fast'
475         WRITE(numout,*) '        nb_ana    = ', nb_ana
476         CALL flush(numout)
477      ENDIF
478      !
479      IF (nb_ana > nharm_max) THEN
480        IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
481        IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
482        nstop = nstop + 1
483      ENDIF
484
485
486      ALLOCATE(ntide_all(nb_ana))
487      ALLOCATE(ntide_sub(nb_ana))
488
489      DO jk=1,nb_ana
490       DO ji=1,nb_harmo
491          IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
492             ntide_sub(jk) = ji
493             ntide_all(jk) = ntide(ji)
494             EXIT
495          END IF
496       END DO
497      END DO
498
499     ALLOCATE( anau(nb_ana) )
500     ALLOCATE( anav(nb_ana) )
501     ALLOCATE( anaf(nb_ana) )
502
503
504    IF( ln_diaharm_fast ) THEN
505
506          ! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
507          nvar_2d = 0; nvar_3d = 0
508          IF ( ln_ana_ssh   ) nvar_2d = nvar_2d + 1       ! analysis elevation
509          IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2       ! analysis depth-averaged velocity
510          IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2       ! analysis bottom friction
511               
512          IF ( ln_ana_rho   ) nvar_3d = nvar_3d + 1       ! analysis density
513          IF ( ln_ana_uv3d  ) nvar_3d = nvar_3d + 2       ! analysis 3D horizontal velocities
514          IF ( ln_ana_w3d   ) nvar_3d = nvar_3d + 1       ! analysis 3D vertical velocity
515
516          ! CHECK IF SOMETHING TO RUN
517          IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE.   ! no 2d to run
518          IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE.   ! no 3d to run
519    !      IF ( nvar_2d .gt. 0 .and. nvar_3d .gt. 0 ) lk_diaharm_fast = .FALSE.
520    !      IF ( .NOT. ln_diaharm_store ) lk_diaharm_fast = .FALSE.
521
522          IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
523
524             ! DO ALLOCATIONS
525             IF ( lk_diaharm_2D ) THEN
526                ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj,    nvar_2d) )
527                ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
528                ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
529                ALLOCATE( g_out2D (jpi,jpj) )
530                ALLOCATE( h_out2D (jpi,jpj) )
531                ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
532             ENDIF
533     
534             IF ( lk_diaharm_3D ) THEN
535                ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
536                ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
537                ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
538                ALLOCATE( g_out3D (jpi,jpj,jpk) )
539                ALLOCATE( h_out3D (jpi,jpj,jpk) )
540                ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
541             ENDIF
542
543             ALLOCATE( cc(nb_ana*2+1,nb_ana*2+1) )
544             ALLOCATE( a (nb_ana*2+1,nb_ana*2+1) )
545             ALLOCATE( bzz(nb_ana*2+1) )
546             ALLOCATE( x  (nb_ana*2+1) )
547             ALLOCATE( c  (nb_ana*2+1) )
548             ! END ALLOCATE
549
550             ! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
551             ! MAKES THINGS EASIER AND FASTER LATER
552             ! !!! UGLY !!!
553             jh = 1; k2d = 0; 
554             IF ( ln_ana_ssh   ) THEN
555                k2d = k2d + 1; m_posi_2d(k2d) = jh
556                IF(lwp) WRITE(numout,*) "   - ssh harmonic analysis activated (ln_ana_ssh)"
557             ENDIF
558             jh = jh + 1
559             IF ( ln_ana_uvbar ) THEN
560                k2d = k2d + 1; m_posi_2d(k2d) = jh
561                jh  = jh  + 1 
562                k2d = k2d + 1; m_posi_2d(k2d) = jh
563                IF(lwp) WRITE(numout,*) "   - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
564             ELSE
565                jh  = jh  + 1
566             ENDIF
567             jh = jh + 1
568             IF ( ln_ana_bfric ) THEN
569                k2d = k2d + 1; m_posi_2d(k2d) = jh
570                jh  = jh  + 1; 
571                k2d = k2d + 1; m_posi_2d(k2d) = jh
572                IF(lwp) WRITE(numout,*) "   - bottom friction harmonic analysis activated (ln_ana_vbfr)"
573             ELSE
574                jh  = jh  + 1
575             ENDIF
576
577             ! and for 3D
578             jh = 1; k3d = 0; 
579             IF ( ln_ana_rho  ) THEN
580                k3d = k3d + 1; m_posi_3d(k3d) = jh
581                IF(lwp) WRITE(numout,*) "   - 3D density harmonic analysis activated (ln_ana_rho)"
582             ENDIF
583             jh = jh + 1
584             IF ( ln_ana_uv3d )  THEN
585                k3d = k3d + 1; m_posi_3d(k3d) = jh
586                jh  = jh  + 1 
587                k3d = k3d + 1; m_posi_3d(k3d) = jh
588                IF(lwp) WRITE(numout,*) "   - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
589             ELSE
590                jh  = jh  + 1
591             ENDIF
592             jh = jh + 1
593             IF ( ln_ana_w3d ) THEN
594                k3d = k3d + 1; m_posi_3d(k3d) = jh
595                IF(lwp) WRITE(numout,*) "   - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
596             ENDIF
597
598             ! SELECT AND STORE FREQUENCIES
599             IF(lwp)    WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
600             DO jh=1,nb_ana
601                om_tide(jh) = omega_tide( ntide_sub(jh) ) 
602                IF(lwp) WRITE(numout,*) '        - ',tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr"
603             ENDDO
604
605             ! READ RESTART IF
606             IF ( ln_diaharm_read_restart ) THEN
607                IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run. kt = ",kt
608                ! Need to read in  bssh bz, cc anau anav and anaf
609                call harm_rst_read  ! This reads in from the previous day
610                                    ! Currrently the data in in assci format
611             ELSE
612
613                IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh. kt = ",kt 
614     
615                IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,:  ) = 0.0_wp
616                IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
617                cc           = 0.0_wp
618                a    (:,:)   = 0.0_wp ! NB
619                bzz  (:)     = 0.0_wp
620                x    (:)     = 0.0_wp
621                c    (:)     = 0.0_wp
622                anau (:)     = 0.0_wp
623                anav (:)     = 0.0_wp
624                anaf (:)     = 0.0_wp
625
626                DO jh = 1, nb_ana
627                   anau(jh) = utide ( ntide_sub(jh) )
628                   anav(jh) = v0tide( ntide_sub(jh) )
629                   anaf(jh) = ftide ( ntide_sub(jh) )
630                END DO
631
632                fjulday_startharm=fjulday !Set this at very start and store
633                !JT this is a mistake - only works on daily cycles, should use fjulnsec_dayday
634
635                IF (lwp) THEN
636                   WRITE(numout,*) '--------------------------'
637                   WRITE(numout,*) '   - Output anaf for check'
638                   WRITE(numout,*) 'ANA F', anaf
639                   WRITE(numout,*) 'ANA U', anau
640                   WRITE(numout,*) 'ANA V', anav
641                   WRITE(numout,*) 'fjulday',fjulday
642                   WRITE(numout,*) 'fjulday_startharm',fjulday_startharm
643                   WRITE(numout,*) 'nsec_day',nsec_day
644                   WRITE(numout,*) 'kt',kt
645                   WRITE(numout,*) '--------------------------'
646                ENDIF
647
648             ENDIF
649
650          ELSE
651
652             IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
653
654          ENDIF
655      ENDIF
656
657   END SUBROUTINE harm_ana_init
658!
659   SUBROUTINE gelim (a,b,x,n)
660      !!----------------------------------------------------------------------
661      !!                    ***  ROUTINE harm_ana  ***
662      !!
663      !! ** Purpose :   Guassian elimination
664      !!
665      !!
666      !! ** Action  : - first action (share memory array/varible modified
667      !!                in this routine
668      !!              - second action .....
669      !!              - .....
670      !!
671      !! References :
672      !!   Give references if exist otherwise suppress these lines
673      !!
674      !! History :
675        implicit none
676!
677        integer  :: n
678        REAL(WP) :: b(nb_ana*2+1), a(nb_ana*2+1,nb_ana*2+1)
679        REAL(WP) :: x(nb_ana*2+1)
680        INTEGER  :: row,col,prow,pivrow,rrow
681        REAL(WP) :: atemp
682        REAL(WP) :: pivot
683        REAL(WP) :: m
684
685        do row=1,n-1
686           pivrow=row
687           pivot=a(row,n-row+1)
688           do prow=row+1,n
689              if (abs(a(prow,n-row+1)).gt.abs(pivot)  ) then
690                 pivot=a(prow,n-row+1)
691                 pivrow=prow
692              endif
693           enddo
694!  swap row and prow
695           if ( pivrow .ne. row ) then
696              atemp=b(pivrow)
697              b(pivrow)=b(row)
698              b(row)=atemp
699              do col=1,n
700                 atemp=a(pivrow,col)
701                 a(pivrow,col)=a(row,col)
702                 a(row,col)=atemp
703              enddo
704           endif
705
706           do rrow=row+1,n
707              if (a(row,row).ne.0) then
708   
709                 m=-a(rrow,n-row+1)/a(row,n-row+1)
710                 do col=1,n
711                    a(rrow,col)=m*a(row,col)+a(rrow,col)
712                 enddo
713                 b(rrow)=m*b(row)+b(rrow)
714              endif
715           enddo
716        enddo
717!  back substitution now
718
719        x(1)=b(n)/a(n,1)
720        do row=n-1,1,-1
721           x(n-row+1)=b(row)
722           do col=1,(n-row)
723              x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) 
724           enddo
725
726           x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
727        enddo
728
729        return
730   END SUBROUTINE gelim
731
732   SUBROUTINE harm_ana_out
733      !!----------------------------------------------------------------------
734      !!                  ***  ROUTINE harm_ana_init  ***
735      !!                   
736      !! ** Purpose :   initialization of ....
737      !!
738      !! ** Method  :   blah blah blah ...
739      !!
740      !! ** input   :   Namlist namexa
741      !!
742      !! ** Action  :   ... 
743      !!
744      !! history :
745      !!   9.0  !  03-08  (Autor Names)  Original code
746      !!----------------------------------------------------------------------
747        USE dianam          ! build name of file (routine)
748 
749      !! * local declarations
750      INTEGER :: ji, jj, jk, jgrid, jh    ! dummy loop indices
751!      INTEGER :: nh_T
752!      INTEGER :: nid_harm
753!      CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names
754!      CHARACTER (len=40) :: clhstnamu, clhstnamv    ! temporary names
755      CHARACTER (len=40) :: suffix
756      CHARACTER (len=40) :: tmp_name
757!      REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi  ! temporary scalars
758
759      do jgrid=1,nvar_2d
760          do jh=1,nb_ana
761             h_out2D = 0.0
762             g_out2D = 0.0
763             do jj=1,nlcj
764                do ji=1,nlci
765                   cca=g_cosamp2D(jh,ji,jj,jgrid)
766                   ssa=g_sinamp2D(jh,ji,jj,jgrid)
767                   h_out2D(ji,jj)=sqrt(cca**2+ssa**2)
768                   IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
769                      g_out2D(ji,jj)= 0.0_wp
770                   ELSE
771                      g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)       
772                   ENDIF
773                   IF (h_out2D(ji,jj).ne.0) THEN
774                       h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
775                   ENDIF
776                   IF (g_out2D(ji,jj).ne.0) THEN  !Correct and take modulus
777                       g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
778                       if (g_out2D(ji,jj).gt.360.0) then
779                           g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
780                       else if (g_out2D(ji,jj).lt.0.0) then
781                           g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
782                       endif
783                   ENDIF
784                enddo
785             enddo
786             !
787             ! NETCDF OUTPUT
788             suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
789             IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
790
791             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
792             IF( iom_use(TRIM(tmp_name)) )  THEN
793                IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D)
794                IF(lwp) WRITE(numout,*) "harm_ana_out_names", tmp_name,tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr"
795                CALL iom_put( TRIM(tmp_name), h_out2D(:,:) )
796             ELSE
797                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
798             ENDIF
799
800             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
801             IF( iom_use(TRIM(tmp_name)) )  THEN
802                IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)
803                CALL iom_put( TRIM(tmp_name), g_out2D(:,:) )
804             ELSE
805                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
806             ENDIF
807
808          enddo
809
810         !tmp_name=TRIM(suffix)//'_tidal_anal_offset'
811         !IF( iom_use(TRIM(tmp_name)) )  THEN
812         !   IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name)
813         !   CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid))
814         !ELSE
815         !   IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
816         !ENDIF
817
818      enddo
819!
820! DO THE SAME FOR 3D VARIABLES
821!
822      do jgrid=1,nvar_3d
823          do jh=1,nb_ana
824             h_out3D = 0.0
825             g_out3D = 0.0
826             DO jk=1,jpkm1
827                do jj=1,nlcj
828                   do ji=1,nlci
829                      cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
830                      ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
831                      h_out3D(ji,jj,jk)=sqrt(cca**2+ssa**2)
832                      IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
833                         g_out3D(ji,jj,jk) = 0.0_wp
834                      ELSE
835                         g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
836                      ENDIF
837                      IF (h_out3D(ji,jj,jk).ne.0) THEN
838                          h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
839                      ENDIF
840                      IF (g_out3D(ji,jj,jk).ne.0) THEN  !Correct and take modulus
841                          g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
842                          if      (g_out3D(ji,jj,jk).gt.360.0) then
843                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
844                          else if (g_out3D(ji,jj,jk).lt.0.0) then
845                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
846                          endif
847                      ENDIF
848                   enddo    ! ji
849                enddo       ! jj
850             ENDDO          ! jk
851             !
852             ! NETCDF OUTPUT
853             suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
854             IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
855
856             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
857             IF( iom_use(TRIM(tmp_name)) )  THEN
858                IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)
859                CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) )
860             ELSE
861                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
862             ENDIF
863
864             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
865             IF( iom_use(TRIM(tmp_name)) )  THEN
866                IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)
867                CALL iom_put(tmp_name, g_out3D(:,:,:) )
868             ELSE
869                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
870             ENDIF
871
872          enddo             ! jh
873
874         !tmp_name=TRIM(suffix)//'_tidal_anal_offset'
875         !IF( iom_use(TRIM(tmp_name)) )  THEN
876         !   IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name)
877         !   CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid))
878         !ELSE
879         !   IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
880         !ENDIF
881
882      enddo                 ! jgrid
883!
884   END SUBROUTINE harm_ana_out
885!
886   SUBROUTINE harm_rst_write(kt)
887      !!----------------------------------------------------------------------
888      !!                  ***  ROUTINE harm_ana_init  ***
889      !!                   
890      !! ** Purpose :  To write out cummulated Tidal Harmomnic data to file for
891      !!               restarting
892      !!
893      !! ** Method  :   restart files will be dated by default
894      !!
895      !! ** input   :   
896      !!
897      !! ** Action  :   ... 
898      !!
899      !! history :
900      !!   0.0  !  01-16  (Enda O'Dea)  Original code
901      !! ASSUMES  dated file for rose  , can change later to be more generic
902      !!----------------------------------------------------------------------
903      INTEGER, INTENT(in) ::   kt     ! ocean time-step
904      !!
905      INTEGER             ::   jh, j2d, j3d
906      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
907      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
908      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
909      CHARACTER(LEN=250)  ::   clfinal   ! full name
910
911      !restart file
912      DO j2d=1,nvar_2d
913         CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
914         DO jh=1,nb_ana
915            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2  , :, :, j2d ) )
916            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
917         ENDDO
918      ENDDO
919
920      DO j3d=1,nvar_3d
921         !JT CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
922         CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
923         DO jh=1,nb_ana
924            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2  , :, :, :, j3d ) )
925            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
926         ENDDO
927      ENDDO
928
929      IF(lwp) THEN
930        IF( kt > 999999999 ) THEN ; WRITE(clkt, *       ) kt
931        ELSE                      ; WRITE(clkt, '(i8.8)') kt
932        ENDIF
933        clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
934        clpath = TRIM(cn_ocerst_outdir)
935        IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
936        IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
937
938        WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
939        OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
940        WRITE(66) cc
941        WRITE(66) anau
942        WRITE(66) anav
943        WRITE(66) anaf
944        WRITE(66) fjulday_startharm
945        CLOSE(66)
946        WRITE(numout,*) '----------------------------'
947        WRITE(numout,*) '   harm_rst_write: DONE '
948        WRITE(numout,*) cc
949        WRITE(numout,*) anaf
950        WRITE(numout,*) anau
951        WRITE(numout,*) anav
952        WRITE(numout,*) fjulday_startharm
953        WRITE(numout,*) '----------------------------'
954      ENDIF
955 
956   END SUBROUTINE harm_rst_write
957
958   SUBROUTINE harm_rst_read
959      !!----------------------------------------------------------------------
960      !!                  ***  ROUTINE harm_ana_init  ***
961      !!                   
962      !! ** Purpose :  To read in  cummulated Tidal Harmomnic data to file for
963      !!               restarting
964      !!
965      !! ** Method  :   
966      !!
967      !! ** input   :   
968      !!
969      !! ** Action  :   ... 
970      !!
971      !! history :
972      !!   0.0  !  01-16  (Enda O'Dea)  Original code
973      !! ASSUMES  dated file for rose  , can change later to be more generic
974      !!----------------------------------------------------------------------
975      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
976      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
977      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
978      CHARACTER(LEN=250)  ::   clfinal   ! full name
979      INTEGER             ::   jh, j2d, j3d
980
981      IF( nit000 > 999999999 ) THEN ; WRITE(clkt, *       ) nit000-1
982      ELSE                      ; WRITE(clkt, '(i8.8)') nit000-1
983      ENDIF
984      clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
985      clpath = TRIM(cn_ocerst_outdir)
986      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
987
988      IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
989
990      DO j2d=1,nvar_2d
991         CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
992         IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) )
993         DO jh=1,nb_ana
994            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2  , :, :, j2d ) )
995            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
996         ENDDO
997      ENDDO
998
999      DO j3d=1,nvar_3d
1000         !JT  CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1001         CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1002         IF(lwp) WRITE(numout,*) "3D", j3d,  m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) )
1003
1004         DO jh=1,nb_ana
1005            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2  , :, :, :, j3d ) )
1006            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
1007         ENDDO
1008      ENDDO
1009
1010      WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
1011      OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
1012      READ(66) cc
1013      READ(66) anau
1014      READ(66) anav
1015      READ(66) anaf
1016      READ(66) fjulday_startharm
1017      CLOSE(66)
1018
1019      IF(lwp) THEN
1020        WRITE(numout,*) '----------------------------'
1021        WRITE(numout,*) '   Checking anaf is correct'
1022        WRITE(numout,*) cc
1023        WRITE(numout,*) anaf
1024        WRITE(numout,*) fjulday_startharm
1025        WRITE(numout,*) '----------------------------'
1026      ENDIF
1027 
1028   END SUBROUTINE harm_rst_read
1029
1030   !!======================================================================
1031#else
1032!!---------------------------------------------------------------------------------
1033!!   Dummy module                                   NO harmonic Analysis
1034!!---------------------------------------------------------------------------------
1035        LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast  = .FALSE.   ! to be run or not
1036
1037        CONTAINS
1038           SUBROUTINE harm_rst_write(kt)     ! Dummy routine
1039           END SUBROUTINE harm_rst_write
1040           SUBROUTINE harm_rst_read    ! Dummy routine
1041           END SUBROUTINE harm_rst_read
1042           SUBROUTINE harm_ana_out      ! Dummy routine
1043           END SUBROUTINE harm_ana_out
1044           SUBROUTINE harm_ana_init
1045           END SUBROUTINE harm_ana_init
1046           SUBROUTINE harm_ana( kt )
1047!--- NB : end call not properly written
1048           END SUBROUTINE harm_ana
1049!           END SUBROUTINE harm_ana_init
1050!--- END NB
1051           SUBROUTINE gelim (a,b,x,n)
1052!--- NB : end call not properly written
1053           END SUBROUTINE gelim
1054!           END SUBROUTINE gelim (a,b,x,n)
1055!--- END NB           
1056#endif
1057
1058END MODULE diaharm_fast 
Note: See TracBrowser for help on using the repository browser.