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 @ 15575

Last change on this file since 15575 was 15575, checked in by hadjt, 2 years ago

Improved version of Diaharm_fast, for testing. using sec2start based on adatrj (* 86400._wp)

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