source: branches/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90 @ 6370

Last change on this file since 6370 was 6370, checked in by deazer, 6 years ago

Added restart switch for harmonic analysis reading

File size: 23.5 KB
Line 
1MODULE harmonic_analysis
2   !!======================================================================
3   !!                       ***  MODULE  example  ***
4   !! Ocean physics:  On line harmonic analyser
5   !!                 
6   !!=====================================================================
7
8#if defined key_harm_ana
9   !!----------------------------------------------------------------------
10   !!   'key_harm_ana'  :                Calculate harmonic analysis
11   !!----------------------------------------------------------------------
12   !!   harm_ana        :
13   !!   harm_ana_init   :
14   !!----------------------------------------------------------------------
15
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE iom
19   USE in_out_manager  ! I/O units
20   USE phycst          ! physical constants
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE bdy_par         ! Unstructured boundary parameters
23   USE bdy_oce         ! ocean open boundary conditions
24   USE bdytides        ! tidal bdy forcing
25   USE daymod          ! calendar
26   USE tideini
27   USE restart
28   USE ioipsl, ONLY : ju2ymds    ! for calendar
29
30   IMPLICIT NONE
31   PRIVATE
32
33   !! *  Routine accessibility
34   PUBLIC harm_ana    ! routine called in step.F90 module
35
36   !! * Module variables
37   INTEGER, PARAMETER ::  nharm_max  =   jptides_max ! max number of harmonics to be analysed
38   INTEGER, PARAMETER ::  nhm_max    =   2*jptides_max+1 
39   INTEGER, PARAMETER ::  nvab       = 2 ! number of 3D variables
40   INTEGER            ::  nharm
41   INTEGER            ::  nhm 
42   INTEGER ::                 & !!! ** toto namelist (namtoto) **
43      nflag  =  1                ! default value of nflag
44   REAL(wp), DIMENSION(nharm_max) ::                & 
45      om_tide                     ! tidal frequencies ( rads/sec)
46   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:)   ::                & 
47      bzz,c,x    ! work arrays
48   REAL(wp) :: cca,ssa,zm,bt
49   REAL(wp) :: ccau,ccav,ssau,ssav
50   REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf   ! nodel/phase corrections used by diaharmana
51!
52   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::   &
53      bssh,bubar,bvbar        ! work array for ssh anaylsis
54   REAL(wp), ALLOCATABLE,SAVE, DIMENSION(:,:,:) ::   &
55      cosampz, cosampu,cosampv,      &
56      sinampz, sinampu,sinampv
57   REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a
58!
59   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
60      gout,hout        ! arrays for output
61   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
62      huout,hvout,guout,gvout        ! arrays for output
63   REAL(wp), PUBLIC ::   fjulday_startharm       !: Julian Day since start of harmonic analysis
64
65   !! * Substitutions
66
67   !!----------------------------------------------------------------------
68   !!  OPA 9.0 , LOCEAN-IPSL (2005)
69   !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
70   !! or  TOP 1.0 , LOCEAN-IPSL (2005)
71   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
72   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
73   !!----------------------------------------------------------------------
74
75CONTAINS
76
77   SUBROUTINE harm_ana( kt )
78      !!----------------------------------------------------------------------
79      !!                    ***  ROUTINE harm_ana  ***
80      !!
81      !! ** Purpose :   Harmonic analyser
82      !!
83      !! ** Method  :   
84      !!
85      !! ** Action  : - first action (share memory array/varible modified
86      !!                in this routine
87      !!              - second action .....
88      !!              - .....
89      !!
90      !! References :
91      !!   Give references if exist otherwise suppress these lines
92      !!
93      !! History :
94      !!   9.0  !  03-08  (Autor Names)  Original code
95      !!        !  02-08  (Author names)  brief description of modifications
96      !!----------------------------------------------------------------------
97      !! * Modules used
98     
99      !! * arguments
100      INTEGER, INTENT( in  ) ::   & 
101         kt                          ! describe it!!!
102
103      !! * local declarations
104      INTEGER ::   ji, jj, jk,  &        ! dummy loop arguments
105                   ih,i1,i2,iv,jgrid
106
107
108      !!--------------------------------------------------------------------
109
110
111
112      IF( kt == nit000  )   CALL harm_ana_init    ! Initialization (first time-step only)
113
114! this bit done every time step
115        !nharm=ncpt_anal
116        nharm=nb_harmo
117        nhm=2*nb_harmo+1
118       
119        c(1)=1.0
120
121        do ih=1,nb_harmo
122             c(2*ih)=    cos((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
123             c(2*ih+1)=  sin((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
124        enddo
125
126
127         do ji=1,jpi
128          do jj=1,jpj
129            do ih=1,nhm
130            bssh (ih,ji,jj)=bssh (ih,ji,jj)+c(ih)*sshn(ji,jj)
131            bubar(ih,ji,jj)=bubar(ih,ji,jj)+c(ih)*un_b(ji,jj)
132            bvbar(ih,ji,jj)=bvbar(ih,ji,jj)+c(ih)*vn_b(ji,jj)
133            enddo
134          enddo
135        enddo
136
137!        IF(lwp) WRITE(666,'(4E15.5)') adatrj*86400._wp, sshn(10,10),bssh(2,10,10),bssh(3,10,10)
138!
139         do i1=1,nhm
140           do i2=1,nhm
141             cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
142            enddo
143         enddo
144!
145! At End of run
146
147      IF (kt ==  nitend  ) then
148
149      IF(lwp) WRITE(numout,*)
150      IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
151      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
152      CALL harm_rst_write(kt)     ! Dump out data for a restarted run
153
154
155
156!
157        cosampu=0.0_wp
158        sinampu=0.0_wp
159        cosampv=0.0_wp
160        sinampv=0.0_wp
161        cosampz=0.0_wp
162        sinampz=0.0_wp
163!
164     do jgrid=1,3 ! elevation, Ubar,Vbar
165        do ji=1,jpi
166            do jj=1,jpj
167               bt=1.0
168               do ih=1,nhm
169                   if(jgrid .eq. 1) then
170                     bzz(ih)=bssh(ih,ji,jj)
171                   endif
172                   if(jgrid .eq. 2) then
173                     bzz(ih)=bubar(ih,ji,jj)
174                   endif
175                   if(jgrid .eq. 3) then
176                     bzz(ih)=bvbar(ih,ji,jj)
177                   endif
178                   bt=bt*bzz(ih)
179               enddo
180
181               do i1=1,nhm
182                   do i2=1,nhm
183                       a(i1,i2)=cc(i1,i2)
184                   enddo
185               enddo
186!  now do gaussian elimination of the system
187!  a * x = b
188!  the matrix x is (a0,a1,b1,a2,b2 ...)
189!  the matrix a and rhs b solved here for x
190              x=0.0d0
191              if(bt.ne.0.) then
192                  if(lwp .and. ji == 10 .and. ji == 10) then
193              endif
194              call gelim(a,bzz,x,nhm)
195
196              do ih=1,nb_harmo
197                if(jgrid .eq. 1) then
198                   cosampz(ih,ji,jj)=x(ih*2)
199                   sinampz(ih,ji,jj)=x(ih*2+1)
200                endif
201                if(jgrid .eq. 2) then
202                   cosampu(ih,ji,jj)=x(ih*2)
203                   sinampu(ih,ji,jj)=x(ih*2+1)
204                endif
205                if(jgrid .eq. 3) then
206                   cosampv(ih,ji,jj)=x(ih*2)
207                   sinampv(ih,ji,jj)=x(ih*2+1)
208                endif
209              enddo
210                if(jgrid .eq. 1) then
211                  cosampz(0,ji,jj)=x(1)
212                  sinampz(0,ji,jj)=0.0_wp
213                endif
214                if(jgrid .eq. 2) then
215                  cosampu(0,ji,jj)=x(1)
216                  sinampu(0,ji,jj)=0.0_wp
217                endif
218                if(jgrid .eq. 3) then
219                  cosampv(0,ji,jj)=x(1)
220                  sinampv(0,ji,jj)=0.0_wp
221                endif
222              endif
223           enddo
224        enddo
225      enddo ! jgrid
226!
227!
228        CALL harm_ana_out     ! output analysis (last time step)
229       
230    ENDIF
231   END SUBROUTINE harm_ana
232
233   SUBROUTINE harm_ana_init
234      !!----------------------------------------------------------------------
235      !!                  ***  ROUTINE harm_ana_init  ***
236      !!                   
237      !! ** Purpose :   initialization of ....
238      !!
239      !! ** Method  :   blah blah blah ...
240      !!
241      !! ** input   :   Namlist namexa
242      !!
243      !! ** Action  :   ... 
244      !!
245      !! history :
246      !!   9.0  !  03-08  (Autor Names)  Original code
247      !!----------------------------------------------------------------------
248      !! * local declarations
249      INTEGER ::   ji, jj, jk, jit,ih   ! dummy loop indices
250
251      !!----------------------------------------------------------------------
252      IF(lwp) WRITE(numout,*)
253      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
254      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
255     ! DO ALLOCATIONS
256
257      ALLOCATE( bubar(nb_harmo*2+1,jpi,jpj) )
258      ALLOCATE( bvbar(nb_harmo*2+1,jpi,jpj) )
259      ALLOCATE( bssh (nb_harmo*2+1,jpi,jpj) )
260
261      ALLOCATE( cosampz(0:nb_harmo*2+1,jpi,jpj))
262      ALLOCATE( cosampu(0:nb_harmo*2+1,jpi,jpj))
263      ALLOCATE( cosampv(0:nb_harmo*2+1,jpi,jpj))
264
265      ALLOCATE( sinampz(0:nb_harmo*2+1,jpi,jpj))
266      ALLOCATE( sinampu(0:nb_harmo*2+1,jpi,jpj))
267      ALLOCATE( sinampv(0:nb_harmo*2+1,jpi,jpj))
268
269      ALLOCATE( cc(nb_harmo*2+1,nb_harmo*2+1) )
270
271      ALLOCATE( a(nb_harmo*2+1,nb_harmo*2+1) )
272
273      ALLOCATE( anau(nb_harmo) )
274      ALLOCATE( anav(nb_harmo) )
275      ALLOCATE( anaf(nb_harmo) )
276
277      ALLOCATE( bzz(nb_harmo*2+1) )
278      ALLOCATE( x(nb_harmo*2+1) )
279      ALLOCATE( c(nb_harmo*2+1) )
280
281      ALLOCATE( gout(jpi,jpj))
282      ALLOCATE( hout(jpi,jpj))
283
284      ALLOCATE( guout(jpi,jpj))
285      ALLOCATE( huout(jpi,jpj))
286
287      ALLOCATE( gvout(jpi,jpj))
288      ALLOCATE( hvout(jpi,jpj))
289
290        ! END ALLOCATE
291
292                 do ih=1,nb_harmo
293                  om_tide(ih)= omega_tide(ih) 
294                 enddo
295       if(ln_harmana_read) then
296               WRITE(numout,*) "Reading previous harmonic data from previous run"
297               ! Need to read in  bssh bz, cc anau anav and anaf
298               call harm_rst_read  ! This reads in from the previous day
299                                   ! Currrently the data in in assci format
300       else
301               WRITE(numout,*) "Starting harmonic analysis from Fresh "
302         bssh(:,:,:)  = 0.0_wp
303         bubar(:,:,:) = 0.0_wp
304         bvbar(:,:,:) = 0.0_wp
305
306         cc=0.0_wp
307
308         anau(:) = 0.0_wp
309         anav(:) = 0.0_wp
310         anaf(:) = 0.0_wp
311         bzz(:) = 0.0_wp
312         x(:) = 0.0_wp
313         c(:) = 0.0_wp
314
315         DO ih = 1, nb_harmo
316            anau(ih) = utide(ih)
317            anav(ih) = v0tide(ih)
318            anaf(ih) = ftide(ih)
319
320         END DO
321         fjulday_startharm=fjulday !Set this at very start and store
322      endif
323
324   END SUBROUTINE harm_ana_init
325!
326   SUBROUTINE gelim (a,b,x,n)
327      !!----------------------------------------------------------------------
328      !!                    ***  ROUTINE harm_ana  ***
329      !!
330      !! ** Purpose :   Guassian elimination
331      !!
332      !! ** Method  :   
333      !!
334      !! ** Action  : - first action (share memory array/varible modified
335      !!                in this routine
336      !!              - second action .....
337      !!              - .....
338      !!
339      !! References :
340      !!   Give references if exist otherwise suppress these lines
341      !!
342      !! History :
343        implicit none
344!
345        integer  :: n
346   REAL(WP) :: b(nb_harmo*2+1),a(nb_harmo*2+1,nb_harmo*2+1)
347   REAL(WP) :: x(nb_harmo*2+1)
348        INTEGER  :: row,col,prow,pivrow,rrow,ntemp
349   REAL(WP) ::  atemp
350   REAL(WP) ::  pivot
351   REAL(WP) ::  m
352   REAL(WP) ::  tempsol(nb_harmo*2+1)
353
354
355        do row=1,n-1
356           pivrow=row
357           pivot=a(row,n-row+1)
358           do prow=row+1,n
359              if (abs(a(prow,n-row+1)).gt.abs(pivot)  ) then
360                 pivot=a(prow,n-row+1)
361                 pivrow=prow
362              endif
363           enddo
364!  swap row and prow
365           if ( pivrow .ne. row ) then
366              atemp=b(pivrow)
367              b(pivrow)=b(row)
368              b(row)=atemp
369              do col=1,n
370                 atemp=a(pivrow,col)
371                 a(pivrow,col)=a(row,col)
372                 a(row,col)=atemp
373              enddo
374           endif
375
376           do rrow=row+1,n
377              if (a(row,row).ne.0) then
378   
379                 m=-a(rrow,n-row+1)/a(row,n-row+1)
380                 do col=1,n
381                    a(rrow,col)=m*a(row,col)+a(rrow,col)
382                 enddo
383                 b(rrow)=m*b(row)+b(rrow)
384              endif
385           enddo
386        enddo
387!  back substitution now
388
389        x(1)=b(n)/a(n,1)
390        do row=n-1,1,-1
391           x(n-row+1)=b(row)
392           do col=1,(n-row)
393              x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) 
394           enddo
395
396           x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
397        enddo
398
399        return
400        END SUBROUTINE gelim
401
402   SUBROUTINE harm_ana_out
403      !!----------------------------------------------------------------------
404      !!                  ***  ROUTINE harm_ana_init  ***
405      !!                   
406      !! ** Purpose :   initialization of ....
407      !!
408      !! ** Method  :   blah blah blah ...
409      !!
410      !! ** input   :   Namlist namexa
411      !!
412      !! ** Action  :   ... 
413      !!
414      !! history :
415      !!   9.0  !  03-08  (Autor Names)  Original code
416      !!----------------------------------------------------------------------
417        USE dianam          ! build name of file (routine)
418 
419      !! * local declarations
420      INTEGER ::   ji, jj, jk, jgrid,jit, ih, jjk   ! dummy loop indices
421      INTEGER :: nh_T
422      INTEGER :: nid_harm
423      CHARACTER (len=40) ::           &
424         clhstnamt, clop1, clop2          ! temporary names
425      CHARACTER (len=40) ::           &
426         clhstnamu, clhstnamv   ! temporary names
427      REAL(wp) ::   &
428         zsto1, zsto2, zout, zmax,            &  ! temporary scalars
429         zjulian, zdt, zmdi 
430         
431        do jgrid=1,3
432          do ih=1,nb_harmo
433            hout = 0.0
434            gout = 0.0
435            do jj=1,nlcj
436                do ji=1,nlci
437                    if(jgrid .eq. 1) then ! SSH
438                       cca=cosampz(ih,ji,jj)
439                       ssa=sinampz(ih,ji,jj)
440                    endif
441                    if(jgrid .eq. 2) then ! UBAR
442                       cca=cosampu(ih,ji,jj)
443                       ssa=sinampu(ih,ji,jj)
444                    endif
445                    if(jgrid .eq. 3) then ! VBAR
446                       cca=cosampv(ih,ji,jj)
447                       ssa=sinampv(ih,ji,jj)
448                    endif
449
450                    hout(ji,jj)=sqrt(cca**2+ssa**2)
451       
452                    if (cca.ne.0.0) then
453                         gout(ji,jj)=(180.0/rpi)*atan(ssa/cca)
454                    else
455                        if (ssa.ne.0.0) then
456                            if (ssa.gt. 0.0) then
457                                 gout(ji,jj)=90.0
458                            else
459                                 gout(ji,jj)=270.0
460                            endif
461                        else
462                            gout(ji,jj)=0.0
463                        endif
464                    endif
465
466                    if (gout(ji,jj).lt.0) then
467                        gout(ji,jj)=gout(ji,jj)+180.0
468                    endif
469                    if (ssa.lt.0) then
470                        gout(ji,jj)=gout(ji,jj)+180.0
471                    endif
472                   
473                    if (hout(ji,jj).ne.0) then
474                        hout(ji,jj)=hout(ji,jj)/anaf(ih)
475                    endif
476                    if (gout(ji,jj).ne.0) then
477                                                !Convert Rad to degree and take
478                                                !modulus
479                         gout(ji,jj)=gout(ji,jj)+MOD( (anau(ih)+anav(ih))/rad , 360.0)
480                        if (gout(ji,jj).gt.360.0) then
481                            gout(ji,jj)=gout(ji,jj)-360.0
482                        else if (gout(ji,jj).lt.0) then
483                            gout(ji,jj)=gout(ji,jj)+360.0
484                        endif
485                    endif
486                enddo
487            enddo
488!NETCDF OUTPUT
489            if(jgrid==1) then
490             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'amp'
491             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'phase'
492             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp', hout(:,:) )
493             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase', gout(:,:) )
494            endif
495            if(jgrid==2) then
496             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_u', hout(:,:) )
497             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_u', gout(:,:) )
498            endif
499            if(jgrid==3) then
500             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_v', hout(:,:) )
501             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_v', gout(:,:) )
502            endif
503          enddo
504        enddo
505 
506!
507   END SUBROUTINE harm_ana_out
508
509   SUBROUTINE harm_rst_write(kt)
510      !!----------------------------------------------------------------------
511      !!                  ***  ROUTINE harm_ana_init  ***
512      !!                   
513      !! ** Purpose :  To write out cummulated Tidal Harmomnic data to file for
514      !!               restarting
515      !!
516      !! ** Method  :   restart files will be dated by default
517      !!
518      !! ** input   :   
519      !!
520      !! ** Action  :   ... 
521      !!
522      !! history :
523      !!   0.0  !  01-16  (Enda O'Dea)  Original code
524      !! ASSUMES  dated file for rose  , can change later to be more generic
525      !!----------------------------------------------------------------------
526      INTEGER, INTENT(in) ::   kt     ! ocean time-step
527      INTEGER             ::   iyear, imonth, iday,ih
528      REAL (wp)           ::   zsec
529      !!
530      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
531      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
532      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
533      CHARACTER(LEN=250)  ::   clfinal   ! full name
534
535      CALL ju2ymds( fjulday , iyear, imonth, iday, zsec)
536
537      WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
538      ! create the file
539      clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana"
540      clpath = TRIM(cn_ocerst_outdir)
541      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
542      WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
543
544
545
546
547!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
548      ! Open a file to write the data into
549         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
550         open(66,file=TRIM(clfinal))
551         write(66,'(1e15.5)') cc
552         !These Three which contain the most data can be moved to a regular
553         !restart file
554
555            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bssh'     , bssh(1,:,:)       )
556            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bubar'    , bubar(1,:,:)       )
557            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bvbar'    , bvbar(1,:,:)       )
558         do ih=1,nb_harmo
559            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
560            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
561
562            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
563            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
564
565            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
566            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
567         enddo
568
569         write(66,'(1e15.5)') anau
570         write(66,'(1e15.5)') anav
571         write(66,'(1e15.5)') anaf
572         write(66,'(1e25.20)') fjulday_startharm
573         close(66)
574 
575   END SUBROUTINE harm_rst_write
576
577   SUBROUTINE harm_rst_read
578      !!----------------------------------------------------------------------
579      !!                  ***  ROUTINE harm_ana_init  ***
580      !!                   
581      !! ** Purpose :  To read in  cummulated Tidal Harmomnic data to file for
582      !!               restarting
583      !!
584      !! ** Method  :   
585      !!
586      !! ** input   :   
587      !!
588      !! ** Action  :   ... 
589      !!
590      !! history :
591      !!   0.0  !  01-16  (Enda O'Dea)  Original code
592      !! ASSUMES  dated file for rose  , can change later to be more generic
593      !!----------------------------------------------------------------------
594      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
595      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
596      CHARACTER(LEN=250)  ::   clfinal   ! full name
597      INTEGER             ::   ih
598
599      ! create the file
600      clname = "restart_harm_ana"
601      clpath = TRIM(cn_ocerst_indir)
602      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
603      WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
604
605
606
607
608!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
609      ! Open a file to read the data ifrom
610         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
611         open(66,file=TRIM(clfinal),status='old')
612         read(66,'(1e15.5)') cc
613!2D regular fields can be read from normal restart this saves space and handy to
614!view in netcdf format also.
615         CALL iom_get( numror,jpdom_autoglo, 'Mean_bssh'     , bssh(1,:,:)       )
616         CALL iom_get( numror,jpdom_autoglo, 'Mean_bubar'    , bubar(1,:,:)       )
617         CALL iom_get( numror,jpdom_autoglo, 'Mean_bvbar'    , bvbar(1,:,:)       )
618         
619         do ih=1,nb_harmo
620
621            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
622            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
623
624            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
625            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
626
627            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
628            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
629         enddo
630         read(66,'(1e15.5)') anau
631         read(66,'(1e15.5)') anav
632         read(66,'(1e15.5)') anaf
633         read(66,'(1e25.20)') fjulday_startharm
634      WRITE(numout,*) 'Checking anaf is correct'
635      WRITE(numout,*) anaf
636      WRITE(numout,*) '---end of anaf checking----'
637
638         close(66)
639
640   END SUBROUTINE harm_rst_read
641
642   !!======================================================================
643#else
644!!---------------------------------------------------------------------------------
645!!   Dummy module                                   NO harmonic Analysis
646!!---------------------------------------------------------------------------------
647        CONTAINS
648           SUBROUTINE harm_rst_write(kt)     ! Dummy routine
649           END SUBROUTINE harm_rst_write
650           SUBROUTINE harm_rst_read    ! Dummy routine
651           END SUBROUTINE harm_rst_read
652           SUBROUTINE harm_ana_out      ! Dummy routine
653           END SUBROUTINE harm_ana_out
654           SUBROUTINE harm_ana_init
655           END SUBROUTINE harm_ana_init
656           SUBROUTINE harm_ana( kt )
657           END SUBROUTINE harm_ana_init
658           SUBROUTINE gelim (a,b,x,n)
659           END SUBROUTINE gelim (a,b,x,n)
660           
661#endif
662
663END MODULE harmonic_analysis
Note: See TracBrowser for help on using the repository browser.