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

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

DIA/diaharm_fast.F90

Nemo compile key dependence removed

Commented code added to help start tidal ellipse paramter calculation (how to copy U and V onto the T grid.

File size: 44.8 KB
Line 
1MODULE diaharm_fast 
2   !!======================================================================
3   !!                       ***  MODULE  example  ***
4   !! Ocean physics:  On line harmonic analyser
5   !!                 
6   !!=====================================================================
7
8   !!----------------------------------------------------------------------
9   !!   'key_harm_ana'  :                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   !JT
86   LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
87   INTEGER ::   nb_ana        ! Number of harmonics to analyse
88   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
89   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_all ! INDEX within the full set of constituents (tide.h90)
90   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_sub ! INDEX within the subset of constituents pass in input
91
92   !! * Substitutions
93
94   !!----------------------------------------------------------------------
95   !!  OPA 9.0 , LOCEAN-IPSL (2005)
96   !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
97   !! or  TOP 1.0 , LOCEAN-IPSL (2005)
98   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
99   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
100   !!----------------------------------------------------------------------
101
102CONTAINS
103
104   SUBROUTINE dia_harm_fast( kt )
105      !!----------------------------------------------------------------------
106      !!                    ***  ROUTINE harm_ana  ***
107      !!
108      !! ** Purpose :   Harmonic analyser
109      !!
110      !! ** Method  :   
111      !!
112      !! ** Action  : - first action (share memory array/varible modified
113      !!                in this routine
114      !!              - second action .....
115      !!              - .....
116      !!
117      !! References :
118      !!   Give references if exist otherwise suppress these lines
119      !!
120      !! History :
121      !!   9.0  !  03-08  (Autor Names)  Original code
122      !!        !  02-08  (Author names)  brief description of modifications
123      !!----------------------------------------------------------------------
124      !! * Modules used
125     
126      !! * arguments
127      INTEGER, INTENT( in  ) ::   & 
128         kt                          ! describe it!!!
129
130      !! * local declarations
131      INTEGER  :: ji, jk, jj          ! dummy loop arguments
132      INTEGER  :: jh, i1, i2, jgrid
133      INTEGER  :: j2d, j3d
134      REAL(WP) :: sec2start,sec2start_old
135      CHARACTER (len=40) :: tmp_name
136      !!--------------------------------------------------------------------
137
138      !JT IF( nn_timing == 1 )   CALL timing_start( 'dia_harm_fast' )
139      IF( ln_timing )   CALL timing_start( 'dia_harm_fast' )
140      IF( kt == nit000   )   CALL harm_ana_init(kt)    ! Initialization (first time-step only)
141
142
143      IF ( ln_diaharm_update_nodal_daily ) THEN
144         !IF (MOD(kt,nint(86400./rdt)) == 0) THEN
145
146
147         IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN     
148            DO jh = 1, nb_ana
149               !JT anau(jh) = 3.141579*utide ( ntide_sub(jh) )/180.
150               !JT anav(jh) = 3.141579*v0tide( ntide_sub(jh) )/180.
151               anau(jh) = utide ( ntide_sub(jh) )
152               anav(jh) = v0tide( ntide_sub(jh) )
153               anaf(jh) = ftide ( ntide_sub(jh) )
154            END DO
155
156            IF(lwp) WRITE(numout,*)
157            IF(lwp) WRITE(numout,*) 'harm_ana : update nodes?',ln_diaharm_update_nodal_daily
158            IF(lwp) WRITE(numout,*) 'harm_ana : date, time',ndastp, nhour, nminute
159
160         ENDIF
161      ENDIF
162
163       
164    DO jh = 1, nb_ana       
165
166        tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anau_utide'
167        IF( iom_use(TRIM(tmp_name)) )  THEN
168        !    IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anau(jh) )
169            CALL iom_put( TRIM(tmp_name), anau(jh) )
170        !ELSE
171        !    IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
172        ENDIF   
173
174        tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anav_v0tide'
175        IF( iom_use(TRIM(tmp_name)) )  THEN
176        !    IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anav(jh) )
177            CALL iom_put( TRIM(tmp_name), anav(jh) )
178        !ELSE
179        !    IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
180        ENDIF
181
182        tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anaf_ftide'
183        IF( iom_use(TRIM(tmp_name)) )  THEN
184        !    IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anaf(jh) )
185            CALL iom_put( TRIM(tmp_name), anaf(jh) )
186        !ELSE
187        !    IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
188        ENDIF
189
190     END DO
191       
192
193     IF ( ln_diaharm_fast .and. ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
194
195          ! this bit done every time step
196          nhm=2*nb_ana+1
197          c(1) = 1.0
198
199          sec2start_old = nint( (fjulday-fjulday_startharm)*86400._wp ) 
200          sec2start = nsec_day - NINT(0.5_wp * rdt)
201          !IF(lwp) WRITE(numout,*) "ztime NEW", kt, sec2start, fjulday_startharm
202
203          DO jh=1,nb_ana
204             c(2*jh  ) = anaf(jh)*cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
205             c(2*jh+1) = anaf(jh)*sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
206             
207             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)
208          ENDDO 
209
210          !IF(lwp) WRITE(numout,*) "c init", c, "c end", sec2start, om_tide(1), anau(1), anav(1),"end nodal"
211
212
213          ! CUMULATE
214          DO ji=1,jpi         ! loop lon
215             DO jj=1,jpj      ! loop lat
216                DO jh=1,nhm   ! loop harmonic
217
218                   DO j2d=1,nvar_2d
219                      IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj)             ! analysis elevation
220                      IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj)             ! analysis depth average velocities
221                      IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj)
222                      !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
223                      !JT IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj)
224                      g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul
225                   ENDDO
226
227                   DO j3d=1,nvar_3d
228                      DO jk=1,jpkm1
229                         IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) *  rhd(ji,jj,jk)               * tmask(ji,jj,jk)   
230                         IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk) 
231                         IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk)
232                         IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) *   wn(ji,jj,jk)               * wmask(ji,jj,jk)
233                         g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul
234                      ENDDO
235                   ENDDO
236
237                ENDDO     ! end loop harmonic
238             ENDDO        ! end loop lat
239          ENDDO           ! end loop lon
240
241          ! Compute nodal factor cumulative cross-product
242          DO i1=1,nhm
243             DO i2=1,nhm
244                cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
245             ENDDO
246          ENDDO
247
248          ! Output RESTART
249          IF( kt == nitrst ) THEN
250             CALL harm_rst_write(kt) ! Dump out data for a restarted run
251          ENDIF
252
253          ! At End of run
254          IF ( kt ==  nitend ) THEN
255
256             IF(lwp) WRITE(numout,*)
257             IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
258             IF(lwp) WRITE(numout,*) '~~~~~~~~~'
259
260             IF( ln_diaharm_compute ) THEN
261
262                 ! INITIALISE TABLE TO 0
263                 IF ( nvar_2d .gt. 0 ) THEN
264                    g_cosamp2D = 0.0_wp
265                    g_sinamp2D = 0.0_wp
266                 ENDIF
267                 IF ( nvar_3d .gt. 0 ) THEN
268                    g_cosamp3D = 0.0_wp
269                    g_sinamp3D = 0.0_wp
270                 ENDIF
271
272                 ! FIRST OUTPUT 2D VARIABLES
273                 DO jgrid=1,nvar_2d    ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically
274                    DO ji=1,jpi        ! loop lon
275                       DO jj=1,jpj     ! loop lat
276                          bt = 1.0_wp; bzz(:) = 0.0_wp
277                          DO jh=1,nhm  ! loop harmonic
278                             bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid)
279                             bt = bt*bzz(jh)
280                          ENDDO
281                          ! Copy back original cumulated nodal factor
282                          a(:,:) = cc(:,:)
283    !                     now do gaussian elimination of the system
284    !                     a * x = b
285    !                     the matrix x is (a0,a1,b1,a2,b2 ...)
286    !                     the matrix a and rhs b solved here for x
287                          x=0.0_wp
288                          IF(bt.ne.0.) THEN
289                            CALL gelim( a, bzz, x, nhm )
290    !                       Backup output in variables
291                            DO jh=1,nb_ana
292                               g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2  )
293                               g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1)
294                            ENDDO
295                            g_cosamp2D   ( 0,ji,jj,jgrid) = x(1)
296                            g_sinamp2D   ( 0,ji,jj,jgrid) = 0.0_wp
297                          ENDIF     ! bt.ne.0.
298                       ENDDO        ! jj
299                    ENDDO           ! ji
300                 ENDDO              ! jgrid
301
302                 ! SECOND OUTPUT 3D VARIABLES
303                 DO jgrid=1,nvar_3d     ! loop number of 3d variables rho, U, V, W
304                    DO jk=1,jpkm1       ! loop over vertical level
305                       DO ji=1,jpi      ! loop over lon
306                          DO jj=1,jpj   ! loop over lat
307                             bt = 1.0_wp; bzz(:) = 0.0_wp
308                             DO jh=1,nhm
309                                bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid)
310                                bt = bt*bzz(jh)
311                             ENDDO
312                             ! Copy back original cumulated nodal factor
313                             a(:,:) = cc(:,:)                     
314    !                        now do gaussian elimination of the system
315    !                        a * x = b
316    !                        the matrix x is (a0,a1,b1,a2,b2 ...)
317    !                        the matrix a and rhs b solved here for x
318                             x=0.0_wp
319                             IF(bt.ne.0.) THEN
320                               CALL gelim( a, bzz, x, nhm )
321    !                          Backup output in variables
322                               DO jh=1,nb_ana
323                                  g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2  )
324                                  g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1)
325                               ENDDO
326                               g_cosamp3D   ( 0,ji,jj,jk,jgrid) = x(1)
327                               g_sinamp3D   ( 0,ji,jj,jk,jgrid) = 0.0_wp
328                            ENDIF     ! bt.ne.0.
329                          ENDDO       ! jj
330                       ENDDO          ! ji
331                    ENDDO             ! jk
332                 ENDDO                ! jgrid
333
334                 CALL harm_ana_out     ! output analysis (last time step)
335
336             ELSE    ! ln_harmana_compute = False
337                 IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step"
338
339             ENDIF   ! ln_harmana_compute
340          ENDIF      ! kt ==  nitend
341
342     ENDIF
343
344      !JT IF( nn_timing == 1 )   CALL timing_stop( 'dia_harm_fast' )
345      IF( ln_timing  )   CALL timing_stop( 'dia_harm_fast' )
346
347   END SUBROUTINE dia_harm_fast 
348
349   SUBROUTINE harm_ana_init( kt )   !JT
350      !!----------------------------------------------------------------------
351      !!                  ***  ROUTINE harm_ana_init  ***
352      !!                   
353      !! ** Purpose :   initialization of ....
354      !!
355      !! ** Method  :   blah blah blah ...
356      !!
357      !! ** input   :   Namlist namexa
358      !!
359      !! ** Action  :   ... 
360      !!
361      !! history :
362      !!   9.0  !  03-08  (Autor Names)  Original code
363      !!----------------------------------------------------------------------
364      !! * local declarations
365
366      INTEGER, INTENT(in) ::   kt     ! ocean time-step  !JT
367      !!
368      INTEGER ::   ji, jk, jh  ! dummy loop indices
369      INTEGER ::   ios                  ! Local integer output status for namelist read
370      INTEGER ::   k2d, k3d             ! dummy number of analysis
371
372
373     
374
375      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, &
376               & tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily
377      !!----------------------------------------------------------------------
378      !JT
379      ln_diaharm_fast = .FALSE.
380      ln_diaharm_multiyear = .FALSE.
381      nn_diaharm_multiyear = 20
382      !JT
383      lk_diaharm_2D    = .TRUE.   ! to run 2d
384      lk_diaharm_3D    = .TRUE.   ! to run 3d
385
386      ln_diaharm_store = .TRUE.
387
388      IF(lwp) WRITE(numout,*)
389      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
390      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
391
392      ! GET NAMELIST DETAILS
393      REWIND( numnam_ref )              ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
394      READ  ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
395901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist' )
396
397      REWIND( numnam_cfg )              ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
398      READ  ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
399902   IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist' )
400      IF(lwm) WRITE ( numond, nam_diaharm_fast )
401
402
403      !
404      IF(lwp) THEN
405         WRITE(numout,*) 'Tidal diagnostics_fast'
406         WRITE(numout,*) '   Fast Harmonic analysis ?:         ln_diaharm_fast= ', ln_diaharm_fast
407         WRITE(numout,*) '   Store output in restart?:         ln_diaharm_store= ', ln_diaharm_store
408         WRITE(numout,*) '   Compute analysis?:         ln_diaharm_compute= ', ln_diaharm_compute
409         WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
410         WRITE(numout,*) '   SSH harmonic analysis: ln_ana_ssh = ', ln_ana_ssh
411         WRITE(numout,*) '   Barotropic Velocities harmonic analysis: ln_ana_uvbar = ', ln_ana_uvbar
412         WRITE(numout,*) '   Bed Friction for harmonic analysis (not implemented): ln_ana_bfric = ', ln_ana_bfric
413         WRITE(numout,*) '   Density harmonic analysis: ln_ana_rho = ', ln_ana_rho
414         WRITE(numout,*) '   3D velocities harmonic analysis: ln_ana_uv3d = ', ln_ana_uv3d
415         WRITE(numout,*) '   Vertical Velocities harmonic analysis: ln_ana_w3d = ', ln_ana_w3d
416         WRITE(numout,*) '   Names of harmonics: tname = ', tname
417         WRITE(numout,*) '   Max number of harmonics: jpmax_harmo = ', jpmax_harmo
418         WRITE(numout,*) '   Number of Harmonics: nb_harmo = ', nb_harmo
419         WRITE(numout,*) '   Multi-year harmonic analysis: ln_diaharm_multiyear = ', ln_diaharm_multiyear
420         WRITE(numout,*) '   Multi-year harmonic analysis - number of years: nn_diaharm_multiyear = ', nn_diaharm_multiyear
421         WRITE(numout,*) '   Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily
422         WRITE(numout,*) '   Number of Harmonics: nyear, nmonth = ', nyear, nmonth
423
424      ENDIF
425      ! JT
426
427
428      IF ( ln_diaharm_multiyear ) THEN
429          ln_diaharm_store = .True.
430          IF(lwp) WRITE(numout,*) '   Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
431          IF ((mod(nyear,nn_diaharm_multiyear) == 0) .AND. ( nmonth == 1 )) THEN ! Jan, year = 1980,2000,2020,2040, restart tidal calculation
432              ln_diaharm_read_restart = .FALSE.
433              IF(lwp) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
434          ELSE
435              ln_diaharm_read_restart = .TRUE.
436              IF(lwp) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
437          ENDIF
438
439
440
441          IF(lwp) WRITE(numout,*) '   Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
442          IF ((mod(nyear,nn_diaharm_multiyear) == (nn_diaharm_multiyear - 1)) .AND. ( nmonth == 12 )) THEN ! Dec year = 1999,2019,2039,2040, restart tidal calculation
443              ln_diaharm_compute = .TRUE.
444              IF(lwp) WRITE(numout,*) '   Compute analysis?:         ln_diaharm_compute= ', ln_diaharm_compute
445          ELSE
446              ln_diaharm_compute = .FALSE.
447              IF(lwp) WRITE(numout,*) '   Compute analysis?:         ln_diaharm_compute= ', ln_diaharm_compute
448          ENDIF
449
450      ENDIF
451      IF ( kt < 10 ) THEN
452        ln_diaharm_read_restart = .FALSE.
453        IF(lwp) WRITE(numout,*) '   kt = ',kt
454        IF(lwp) WRITE(numout,*) '   kt < 10, so setting ln_diaharm_read_restart to .FALSE.'
455        IF(lwp) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart       
456      ENDIF
457
458      ! JT
459
460
461      ! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
462      nb_ana = 0
463      DO jk=1,jpmax_harmo
464         DO ji=1,nb_harmo
465            IF(TRIM(tname(jk)) == TRIM(Wave( ntide(ji) )%cname_tide) ) THEN
466               nb_ana=nb_ana+1
467            ENDIF
468         END DO
469      END DO
470      !
471      IF(lwp) THEN
472         WRITE(numout,*) '        Namelist nam_diaharm_fast'
473         WRITE(numout,*) '        nb_ana    = ', nb_ana
474         CALL flush(numout)
475      ENDIF
476      !
477      IF (nb_ana > nharm_max) THEN
478        IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
479        IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
480        nstop = nstop + 1
481      ENDIF
482
483
484      ALLOCATE(ntide_all(nb_ana))
485      ALLOCATE(ntide_sub(nb_ana))
486
487      DO jk=1,nb_ana
488       DO ji=1,nb_harmo
489          IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
490             ntide_sub(jk) = ji
491             ntide_all(jk) = ntide(ji)
492             EXIT
493          END IF
494       END DO
495      END DO
496
497     ALLOCATE( anau(nb_ana) )
498     ALLOCATE( anav(nb_ana) )
499     ALLOCATE( anaf(nb_ana) )
500
501
502    IF( ln_diaharm_fast ) THEN
503
504          ! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
505          nvar_2d = 0; nvar_3d = 0
506          IF ( ln_ana_ssh   ) nvar_2d = nvar_2d + 1       ! analysis elevation
507          IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2       ! analysis depth-averaged velocity
508          IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2       ! analysis bottom friction
509               
510          IF ( ln_ana_rho   ) nvar_3d = nvar_3d + 1       ! analysis density
511          IF ( ln_ana_uv3d  ) nvar_3d = nvar_3d + 2       ! analysis 3D horizontal velocities
512          IF ( ln_ana_w3d   ) nvar_3d = nvar_3d + 1       ! analysis 3D vertical velocity
513
514          ! CHECK IF SOMETHING TO RUN
515          IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE.   ! no 2d to run
516          IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE.   ! no 3d to run
517    !      IF ( nvar_2d .gt. 0 .and. nvar_3d .gt. 0 ) lk_diaharm_fast = .FALSE.
518    !      IF ( .NOT. ln_diaharm_store ) lk_diaharm_fast = .FALSE.
519
520          IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
521
522             ! DO ALLOCATIONS
523             IF ( lk_diaharm_2D ) THEN
524                ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj,    nvar_2d) )
525                ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
526                ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
527                ALLOCATE( g_out2D (jpi,jpj) )
528                ALLOCATE( h_out2D (jpi,jpj) )
529                ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
530             ENDIF
531     
532             IF ( lk_diaharm_3D ) THEN
533                ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
534                ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
535                ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
536                ALLOCATE( g_out3D (jpi,jpj,jpk) )
537                ALLOCATE( h_out3D (jpi,jpj,jpk) )
538                ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
539             ENDIF
540
541             ALLOCATE( cc(nb_ana*2+1,nb_ana*2+1) )
542             ALLOCATE( a (nb_ana*2+1,nb_ana*2+1) )
543             ALLOCATE( bzz(nb_ana*2+1) )
544             ALLOCATE( x  (nb_ana*2+1) )
545             ALLOCATE( c  (nb_ana*2+1) )
546             ! END ALLOCATE
547
548             ! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
549             ! MAKES THINGS EASIER AND FASTER LATER
550             ! !!! UGLY !!!
551             jh = 1; k2d = 0; 
552             IF ( ln_ana_ssh   ) THEN
553                k2d = k2d + 1; m_posi_2d(k2d) = jh
554                IF(lwp) WRITE(numout,*) "   - ssh harmonic analysis activated (ln_ana_ssh)"
555             ENDIF
556             jh = jh + 1
557             IF ( ln_ana_uvbar ) THEN
558                k2d = k2d + 1; m_posi_2d(k2d) = jh
559                jh  = jh  + 1 
560                k2d = k2d + 1; m_posi_2d(k2d) = jh
561                IF(lwp) WRITE(numout,*) "   - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
562             ELSE
563                jh  = jh  + 1
564             ENDIF
565             jh = jh + 1
566             IF ( ln_ana_bfric ) THEN
567                k2d = k2d + 1; m_posi_2d(k2d) = jh
568                jh  = jh  + 1; 
569                k2d = k2d + 1; m_posi_2d(k2d) = jh
570                IF(lwp) WRITE(numout,*) "   - bottom friction harmonic analysis activated (ln_ana_vbfr)"
571             ELSE
572                jh  = jh  + 1
573             ENDIF
574
575             ! and for 3D
576             jh = 1; k3d = 0; 
577             IF ( ln_ana_rho  ) THEN
578                k3d = k3d + 1; m_posi_3d(k3d) = jh
579                IF(lwp) WRITE(numout,*) "   - 3D density harmonic analysis activated (ln_ana_rho)"
580             ENDIF
581             jh = jh + 1
582             IF ( ln_ana_uv3d )  THEN
583                k3d = k3d + 1; m_posi_3d(k3d) = jh
584                jh  = jh  + 1 
585                k3d = k3d + 1; m_posi_3d(k3d) = jh
586                IF(lwp) WRITE(numout,*) "   - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
587             ELSE
588                jh  = jh  + 1
589             ENDIF
590             jh = jh + 1
591             IF ( ln_ana_w3d ) THEN
592                k3d = k3d + 1; m_posi_3d(k3d) = jh
593                IF(lwp) WRITE(numout,*) "   - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
594             ENDIF
595
596             ! SELECT AND STORE FREQUENCIES
597             IF(lwp)    WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
598             DO jh=1,nb_ana
599                om_tide(jh) = omega_tide( ntide_sub(jh) ) 
600                IF(lwp) WRITE(numout,*) '        - ',tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr"
601             ENDDO
602
603             ! READ RESTART IF
604             IF ( ln_diaharm_read_restart ) THEN
605                IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run. kt = ",kt
606                ! Need to read in  bssh bz, cc anau anav and anaf
607                call harm_rst_read  ! This reads in from the previous day
608                                    ! Currrently the data in in assci format
609             ELSE
610
611                IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh. kt = ",kt 
612     
613                IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,:  ) = 0.0_wp
614                IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
615                cc           = 0.0_wp
616                a    (:,:)   = 0.0_wp ! NB
617                bzz  (:)     = 0.0_wp
618                x    (:)     = 0.0_wp
619                c    (:)     = 0.0_wp
620                anau (:)     = 0.0_wp
621                anav (:)     = 0.0_wp
622                anaf (:)     = 0.0_wp
623
624                DO jh = 1, nb_ana
625                   anau(jh) = utide ( ntide_sub(jh) )
626                   anav(jh) = v0tide( ntide_sub(jh) )
627                   anaf(jh) = ftide ( ntide_sub(jh) )
628                END DO
629
630                fjulday_startharm=fjulday !Set this at very start and store
631                !JT this is a mistake - only works on daily cycles, should use fjulnsec_dayday
632
633                IF (lwp) THEN
634                   WRITE(numout,*) '--------------------------'
635                   WRITE(numout,*) '   - Output anaf for check'
636                   WRITE(numout,*) 'ANA F', anaf
637                   WRITE(numout,*) 'ANA U', anau
638                   WRITE(numout,*) 'ANA V', anav
639                   WRITE(numout,*) 'fjulday',fjulday
640                   WRITE(numout,*) 'fjulday_startharm',fjulday_startharm
641                   WRITE(numout,*) 'nsec_day',nsec_day
642                   WRITE(numout,*) 'kt',kt
643                   WRITE(numout,*) '--------------------------'
644                ENDIF
645
646             ENDIF
647
648          ELSE
649
650             IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
651
652          ENDIF
653      ENDIF
654
655   END SUBROUTINE harm_ana_init
656!
657   SUBROUTINE gelim (a,b,x,n)
658      !!----------------------------------------------------------------------
659      !!                    ***  ROUTINE harm_ana  ***
660      !!
661      !! ** Purpose :   Guassian elimination
662      !!
663      !!
664      !! ** Action  : - first action (share memory array/varible modified
665      !!                in this routine
666      !!              - second action .....
667      !!              - .....
668      !!
669      !! References :
670      !!   Give references if exist otherwise suppress these lines
671      !!
672      !! History :
673        implicit none
674!
675        integer  :: n
676        REAL(WP) :: b(nb_ana*2+1), a(nb_ana*2+1,nb_ana*2+1)
677        REAL(WP) :: x(nb_ana*2+1)
678        INTEGER  :: row,col,prow,pivrow,rrow
679        REAL(WP) :: atemp
680        REAL(WP) :: pivot
681        REAL(WP) :: m
682
683        do row=1,n-1
684           pivrow=row
685           pivot=a(row,n-row+1)
686           do prow=row+1,n
687              if (abs(a(prow,n-row+1)).gt.abs(pivot)  ) then
688                 pivot=a(prow,n-row+1)
689                 pivrow=prow
690              endif
691           enddo
692!  swap row and prow
693           if ( pivrow .ne. row ) then
694              atemp=b(pivrow)
695              b(pivrow)=b(row)
696              b(row)=atemp
697              do col=1,n
698                 atemp=a(pivrow,col)
699                 a(pivrow,col)=a(row,col)
700                 a(row,col)=atemp
701              enddo
702           endif
703
704           do rrow=row+1,n
705              if (a(row,row).ne.0) then
706   
707                 m=-a(rrow,n-row+1)/a(row,n-row+1)
708                 do col=1,n
709                    a(rrow,col)=m*a(row,col)+a(rrow,col)
710                 enddo
711                 b(rrow)=m*b(row)+b(rrow)
712              endif
713           enddo
714        enddo
715!  back substitution now
716
717        x(1)=b(n)/a(n,1)
718        do row=n-1,1,-1
719           x(n-row+1)=b(row)
720           do col=1,(n-row)
721              x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) 
722           enddo
723
724           x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
725        enddo
726
727        return
728   END SUBROUTINE gelim
729
730   SUBROUTINE harm_ana_out
731      !!----------------------------------------------------------------------
732      !!                  ***  ROUTINE harm_ana_init  ***
733      !!                   
734      !! ** Purpose :   initialization of ....
735      !!
736      !! ** Method  :   blah blah blah ...
737      !!
738      !! ** input   :   Namlist namexa
739      !!
740      !! ** Action  :   ... 
741      !!
742      !! history :
743      !!   9.0  !  03-08  (Autor Names)  Original code
744      !!----------------------------------------------------------------------
745        USE dianam          ! build name of file (routine)
746 
747      !! * local declarations
748      INTEGER :: ji, jj, jk, jgrid, jh    ! dummy loop indices
749!      INTEGER :: nh_T
750!      INTEGER :: nid_harm
751!      CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names
752!      CHARACTER (len=40) :: clhstnamu, clhstnamv    ! temporary names
753      CHARACTER (len=40) :: suffix
754      CHARACTER (len=40) :: tmp_name
755!      REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi  ! temporary scalars
756
757      do jgrid=1,nvar_2d
758          do jh=1,nb_ana
759             h_out2D = 0.0
760             g_out2D = 0.0
761             do jj=1,nlcj
762                do ji=1,nlci
763                   cca=g_cosamp2D(jh,ji,jj,jgrid)
764                   ssa=g_sinamp2D(jh,ji,jj,jgrid)
765                   h_out2D(ji,jj)=sqrt(cca**2+ssa**2)
766                   IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
767                      g_out2D(ji,jj)= 0.0_wp
768                   ELSE
769                      g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)       
770                   ENDIF
771                   IF (h_out2D(ji,jj).ne.0) THEN
772                       h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
773                   ENDIF
774                   IF (g_out2D(ji,jj).ne.0) THEN  !Correct and take modulus
775                       g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
776                       if (g_out2D(ji,jj).gt.360.0) then
777                           g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
778                       else if (g_out2D(ji,jj).lt.0.0) then
779                           g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
780                       endif
781                   ENDIF
782                enddo
783             enddo
784             !
785             ! NETCDF OUTPUT
786             suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
787             IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
788
789             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//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(h_out2D)
792                IF(lwp) WRITE(numout,*) "harm_ana_out_names", tmp_name,tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr"
793                CALL iom_put( TRIM(tmp_name), h_out2D(:,:) )
794             ELSE
795                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
796             ENDIF
797
798             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
799             IF( iom_use(TRIM(tmp_name)) )  THEN
800                IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)
801                CALL iom_put( TRIM(tmp_name), g_out2D(:,:) )
802             ELSE
803                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
804             ENDIF
805
806          enddo
807
808         suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
809         tmp_name='TA_'//TRIM(suffix)//'_off'
810         IF( iom_use(TRIM(tmp_name)) )  THEN
811            IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name)
812            CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid))
813         ELSE
814            IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
815         ENDIF
816
817      enddo
818!
819! DO THE SAME FOR 3D VARIABLES
820!
821      do jgrid=1,nvar_3d
822          do jh=1,nb_ana
823             h_out3D = 0.0
824             g_out3D = 0.0
825             DO jk=1,jpkm1
826                do jj=1,nlcj
827                   do ji=1,nlci
828                      cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
829                      ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
830                      h_out3D(ji,jj,jk)=sqrt(cca**2+ssa**2)
831                      IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
832                         g_out3D(ji,jj,jk) = 0.0_wp
833                      ELSE
834                         g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
835                      ENDIF
836                      IF (h_out3D(ji,jj,jk).ne.0) THEN
837                          h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
838                      ENDIF
839                      IF (g_out3D(ji,jj,jk).ne.0) THEN  !Correct and take modulus
840                          g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
841                          if      (g_out3D(ji,jj,jk).gt.360.0) then
842                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
843                          else if (g_out3D(ji,jj,jk).lt.0.0) then
844                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
845                          endif
846                      ENDIF
847                   enddo    ! ji
848                enddo       ! jj
849             ENDDO          ! jk
850             !
851             ! NETCDF OUTPUT
852             suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
853             IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
854
855             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
856             IF( iom_use(TRIM(tmp_name)) )  THEN
857                IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)
858                CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) )
859             ELSE
860                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
861             ENDIF
862
863             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
864             IF( iom_use(TRIM(tmp_name)) )  THEN
865                IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)
866                CALL iom_put(tmp_name, g_out3D(:,:,:) )
867             ELSE
868                IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
869             ENDIF
870
871          enddo             ! jh
872
873         suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
874         tmp_name='TA_'//TRIM(suffix)//'_off'
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
885
886! to output tidal parameters, u and v on t grid
887!
888!                                  !==  standard Cd  ==!
889!         DO jj = 2, jpjm1
890!            DO ji = 2, jpim1
891!               imk = k_mk(ji,jj)    ! ocean bottom level at t-points
892!               zut = un(ji,jj,imk) + un(ji-1,jj,imk)     ! 2 x velocity at t-point
893!               zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk)
894!               !                                                           ! here pCd0 = mask*boost * drag
895!               pCdU(ji,jj) = - pCd0(ji,jj) * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  )
896!            END DO
897!         END DO
898
899
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.