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.
Changeset 3953 – NEMO

Changeset 3953


Ignore:
Timestamp:
2013-07-03T13:41:32+02:00 (11 years ago)
Author:
gm
Message:

dev_r3858_NOC_ZTC, #863 : activate tide potential in filtered ssh case + style in tide modules

Location:
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r3294 r3953  
    11MODULE diaharm  
    2  
    3 #if defined key_diaharm && defined key_tide 
    4    !!================================================================================= 
     2   !!====================================================================== 
    53   !!                       ***  MODULE  diaharm  *** 
    64   !! Harmonic analysis of tidal constituents  
    7    !!================================================================================= 
    8    !! * Modules used 
     5   !!====================================================================== 
     6   !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code 
     7   !!---------------------------------------------------------------------- 
     8#if defined key_diaharm && defined key_tide 
     9   !!---------------------------------------------------------------------- 
     10   !!   'key_diaharm' 
     11   !!   'key_tide' 
     12   !!---------------------------------------------------------------------- 
    913   USE oce             ! ocean dynamics and tracers variables 
    1014   USE dom_oce         ! ocean space and time domain 
    11    USE in_out_manager  ! I/O units 
    12    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    13    USE ioipsl          ! NetCDF IPSL library 
    14    USE diadimg         ! To write dimg 
    1515   USE phycst 
    1616   USE dynspg_oce 
     
    1818   USE daymod 
    1919   USE tide_mod 
    20    USE iom  
     20   USE in_out_manager  ! I/O units 
     21   USE iom             ! I/0 library 
     22   USE ioipsl          ! NetCDF IPSL library 
     23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     24   USE diadimg         ! To write dimg 
    2125   USE timing          ! preformance summary 
    2226   USE wrk_nemo        ! working arrays 
     
    3034   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24 
    3135 
    32    INTEGER ::                            & !! namelist variables 
    33                          nit000_han = 1, & ! First time step used for harmonic analysis 
    34                          nitend_han = 1, & ! Last time step used for harmonic analysis 
    35                          nstep_han  = 1, & ! Time step frequency for harmonic analysis 
    36                          nb_ana            ! Number of harmonics to analyse 
    37  
    38    INTEGER , ALLOCATABLE, DIMENSION(:)       :: name 
    39    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp 
    40    REAL(wp), ALLOCATABLE, DIMENSION(:)       :: ana_freq, vt, ut, ft 
    41    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: out_eta, & 
    42                                                 out_u  , & 
    43                                                 out_v 
    44  
    45    INTEGER :: ninco, nsparse 
    46    INTEGER ,       DIMENSION(jpdimsparse)         :: njsparse, nisparse 
    47    INTEGER , SAVE, DIMENSION(jpincomax)           :: ipos1 
    48    REAL(wp),       DIMENSION(jpdimsparse)         :: valuesparse 
    49    REAL(wp),       DIMENSION(jpincomax)           :: ztmp4 , ztmp7 
    50    REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier 
    51    REAL(wp), SAVE, DIMENSION(jpincomax)           :: zpivot 
    52  
    53    CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   & 
    54        tname         ! Names of tidal constituents ('M2', 'K1',...) 
    55  
    56  
    57 !! * Routine accessibility 
    58    PUBLIC  dia_harm    ! routine called by step.F90 
    59  
    60    !!--------------------------------------------------------------------------------- 
    61    !!   
    62    !!--------------------------------------------------------------------------------- 
    63  
     36   !                            !!!namelist variables 
     37   INTEGER ::   nit000_han = 1   ! First time step used for harmonic analysis 
     38   INTEGER ::   nitend_han = 1   ! Last time step used for harmonic analysis 
     39   INTEGER ::   nstep_han  = 1   ! Time step frequency for harmonic analysis 
     40   INTEGER ::   nb_ana           ! Number of harmonics to analyse 
     41 
     42   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
     43   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ana_temp 
     44   REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut   , vt   , ft 
     45   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta , out_u, out_v 
     46 
     47   INTEGER ::   ninco, nsparse 
     48   INTEGER ,       DIMENSION(jpdimsparse)         ::   njsparse, nisparse 
     49   INTEGER , SAVE, DIMENSION(jpincomax)           ::   ipos1 
     50   REAL(wp),       DIMENSION(jpdimsparse)         ::   valuesparse 
     51   REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7 
     52   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) ::   ztmp3 , zpilier 
     53   REAL(wp), SAVE, DIMENSION(jpincomax)           ::   zpivot 
     54 
     55   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...) 
     56 
     57   PUBLIC   dia_harm   ! routine called by step.F90 
     58 
     59   !!---------------------------------------------------------------------- 
     60   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
     61   !! $Id:$ 
     62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     63   !!---------------------------------------------------------------------- 
    6464CONTAINS 
    6565 
     
    6767      !!---------------------------------------------------------------------- 
    6868      !!                 ***  ROUTINE dia_harm_init  *** 
    69       !!---------------------------------------------------------------------- 
    7069      !!          
    7170      !! ** Purpose :   Initialization of tidal harmonic analysis 
     
    7372      !! ** Method  :   Initialize frequency array and  nodal factor for nit000_han 
    7473      !! 
    75       !! History : 
    76       !!   9.0  O. Le Galloudec and J. Chanut (Original) 
    77       !!-------------------------------------------------------------------- 
    78       !! * Local declarations  
     74      !!-------------------------------------------------------------------- 
    7975      INTEGER :: jh, nhan, jk, ji 
    80  
     76      ! 
    8177      NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname 
    8278      !!---------------------------------------------------------------------- 
     
    9288      tname(:)='' 
    9389      ! 
    94       ! Read Namelist nam_diaharm 
    95       REWIND ( numnam ) 
    96       READ   ( numnam, nam_diaharm ) 
     90      REWIND( numnam )                   ! Read Namelist nam_diaharm 
     91      READ  ( numnam, nam_diaharm ) 
    9792      ! 
    9893      IF(lwp) THEN 
     
    10499      ! Basic checks on harmonic analysis time window: 
    105100      ! ---------------------------------------------- 
    106       IF (nit000 > nit000_han) THEN 
    107         IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nit000_han must be greater than nit000, stop' 
    108         IF(lwp) WRITE(numout,*) ' restart capability not implemented' 
    109         nstop = nstop + 1 
    110       ENDIF 
    111       IF (nitend < nitend_han) THEN 
    112         IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nitend_han must be lower than nitend, stop' 
    113         IF(lwp) WRITE(numout,*) ' restart capability not implemented' 
    114         nstop = nstop + 1 
    115       ENDIF 
    116  
    117       IF (MOD(nitend_han-nit000_han+1,nstep_han).NE.0) THEN 
    118         IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : analysis time span must be a multiple of nstep_han, stop' 
    119         nstop = nstop + 1 
    120       END IF 
    121  
    122       nb_ana=0 
     101      IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   & 
     102         &                                       ' restart capability not implemented' ) 
     103      IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   & 
     104         &                                       'restart capability not implemented' ) 
     105 
     106      IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   & 
     107         &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' ) 
     108 
     109      nb_ana = 0 
    123110      DO jk=1,jpmax_harmo 
    124111         DO ji=1,jpmax_harmo 
     
    153140      ! Initialize frequency array: 
    154141      ! --------------------------- 
    155       ALLOCATE(ana_freq(nb_ana)) 
    156       ALLOCATE(vt      (nb_ana)) 
    157       ALLOCATE(ut      (nb_ana)) 
    158       ALLOCATE(ft      (nb_ana)) 
    159  
    160       CALL tide_harmo(ana_freq, vt, ut , ft, name ,nb_ana) 
     142      ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) ) 
     143 
     144      CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana ) 
    161145 
    162146      IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency ' 
     
    168152      ! Initialize temporary arrays: 
    169153      ! ---------------------------- 
    170       ALLOCATE( ana_temp(jpi,jpj,nb_ana*2,3)) 
     154      ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    171155      ana_temp(:,:,:,:) = 0.e0 
    172156 
    173157   END SUBROUTINE dia_harm_init 
    174     
     158 
     159 
    175160   SUBROUTINE dia_harm ( kt ) 
    176161      !!---------------------------------------------------------------------- 
    177162      !!                 ***  ROUTINE dia_harm  *** 
    178       !!---------------------------------------------------------------------- 
    179163      !!          
    180164      !! ** Purpose :   Tidal harmonic analysis main routine 
     
    182166      !! ** Action  :   Sums ssh/u/v over time analysis [nit000_han,nitend_han] 
    183167      !! 
    184       !! History : 
    185       !!   9.0  O. Le Galloudec and J. Chanut (Original) 
    186       !!-------------------------------------------------------------------- 
    187       !! * Argument: 
     168      !!-------------------------------------------------------------------- 
    188169      INTEGER, INTENT( IN ) :: kt 
    189  
    190       !! * Local declarations 
     170      ! 
    191171      INTEGER  :: ji, jj, jh, jc, nhc 
    192172      REAL(wp) :: ztime, ztemp 
     
    194174      IF( nn_timing == 1 )   CALL timing_start('dia_harm') 
    195175 
    196       IF ( kt .EQ. nit000 ) CALL dia_harm_init 
     176      IF ( kt == nit000 ) CALL dia_harm_init 
    197177 
    198178      IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. & 
     
    211191              DO ji = 1,jpi 
    212192                ! Elevation 
    213                 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1)                & 
    214                                         + ztemp*sshn(ji,jj)*tmask(ji,jj,1)         
     193                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)         
    215194#if defined key_dynspg_ts 
    216                 ! ubar 
    217                 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2)                & 
    218                                         + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) 
    219                 ! vbar 
    220                 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3)                & 
    221                                         + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) 
     195                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) 
     196                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) 
    222197#endif 
    223198              END DO 
     
    229204      END IF 
    230205 
    231       IF ( kt .EQ. nitend_han ) CALL dia_harm_end 
     206      IF ( kt == nitend_han )  CALL dia_harm_end 
    232207 
    233208      IF( nn_timing == 1 )   CALL timing_stop('dia_harm') 
     
    235210   END SUBROUTINE dia_harm 
    236211 
     212 
    237213   SUBROUTINE dia_harm_end 
    238214      !!---------------------------------------------------------------------- 
    239215      !!                 ***  ROUTINE diaharm_end  *** 
    240       !!---------------------------------------------------------------------- 
    241216      !!          
    242217      !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents 
     
    244219      !! ** Action  :  Decompose the signal on the harmonic constituents  
    245220      !! 
    246       !! History : 
    247       !!   9.0  O. Le Galloudec and J. Chanut (Original) 
    248       !!-------------------------------------------------------------------- 
    249  
    250       !! * Local declarations 
     221      !!-------------------------------------------------------------------- 
    251222      INTEGER :: ji, jj, jh, jc, jn, nhan, jl 
    252223      INTEGER :: ksp, kun, keq 
     
    279250               nisparse(ksp) = keq 
    280251               njsparse(ksp) = kun 
    281                valuesparse(ksp)= & 
    282                    +(   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) & 
    283                    +(1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
     252               valuesparse(ksp) = (   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))   & 
     253                  &             + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) ) 
    284254            END DO 
    285255         END DO 
    286256      END DO 
    287257 
    288       nsparse=ksp 
     258      nsparse = ksp 
    289259 
    290260      ! Elevation: 
     
    292262         DO ji = 1, jpi 
    293263            ! Fill input array 
    294             kun=0 
    295             DO jh = 1,nb_ana 
    296                DO jc = 1,2 
     264            kun = 0 
     265            DO jh = 1, nb_ana 
     266               DO jc = 1, 2 
    297267                  kun = kun + 1 
    298268                  ztmp4(kun)=ana_temp(ji,jj,kun,1) 
    299                ENDDO 
    300             ENDDO 
     269               END DO 
     270            END DO 
    301271 
    302272            CALL SUR_DETERMINE(jj) 
     
    310280      END DO 
    311281 
    312       ALLOCATE(out_eta(jpi,jpj,2*nb_ana)) 
    313       ALLOCATE(out_u  (jpi,jpj,2*nb_ana)) 
    314       ALLOCATE(out_v  (jpi,jpj,2*nb_ana)) 
     282      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   &  
     283         &      out_u  (jpi,jpj,2*nb_ana),   & 
     284         &      out_v  (jpi,jpj,2*nb_ana)  ) 
    315285 
    316286      DO jj = 1, jpj 
    317287         DO ji = 1, jpi 
    318288            DO jh = 1, nb_ana  
    319                X1=ana_amp(ji,jj,jh,1) 
    320                X2=-ana_amp(ji,jj,jh,2) 
    321                out_eta(ji,jj,jh)=X1 * tmask(ji,jj,1) 
    322                out_eta(ji,jj,nb_ana+jh)=X2 * tmask(ji,jj,1) 
     289               X1 = ana_amp(ji,jj,jh,1) 
     290               X2 =-ana_amp(ji,jj,jh,2) 
     291               out_eta(ji,jj,jh       ) = X1 * tmask(ji,jj,1) 
     292               out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1) 
    323293            ENDDO 
    324294         ENDDO 
     
    398368   END SUBROUTINE dia_harm_end 
    399369 
     370 
    400371   SUBROUTINE dia_wri_harm 
    401372      !!-------------------------------------------------------------------- 
    402373      !!                 ***  ROUTINE dia_wri_harm  *** 
    403       !!-------------------------------------------------------------------- 
    404374      !!          
    405375      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file 
    406       !! 
    407       !! 
    408       !! History : 
    409       !!   9.0  O. Le Galloudec and J. Chanut (Original) 
    410       !!-------------------------------------------------------------------- 
    411  
    412       !! * Local declarations 
     376      !!-------------------------------------------------------------------- 
    413377      CHARACTER(LEN=lc) :: cltext 
    414378      CHARACTER(LEN=lc) ::   & 
     
    468432#else 
    469433      DO jh = 1, nb_ana 
    470       CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh) ) 
    471       CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,nb_ana+jh) ) 
     434         CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh       ) ) 
     435         CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,jh+nb_ana) ) 
    472436      END DO 
    473437#endif 
    474438 
    475439   END SUBROUTINE dia_wri_harm 
     440 
    476441 
    477442   SUBROUTINE SUR_DETERMINE(init) 
     
    482447   !!        
    483448   !!--------------------------------------------------------------------------------- 
    484    INTEGER, INTENT(in) :: init  
    485                 
     449   INTEGER, INTENT(in) ::   init  
     450   ! 
    486451   INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
    487452   REAL(wp)                        :: zval1, zval2, zx1 
     
    492457   CALL wrk_alloc( jpincomax , ipos2 , ipivot        ) 
    493458             
    494    IF( init==1 )THEN 
    495  
    496       IF( nsparse .GT. jpdimsparse ) & 
    497          CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
    498  
    499       IF( ninco .GT. jpincomax ) & 
    500          CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
    501  
    502       ztmp3(:,:)=0.e0 
    503  
     459   IF( init == 1 ) THEN 
     460      IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
     461      IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
     462      ! 
     463      ztmp3(:,:) = 0._wp 
     464      ! 
    504465      DO jk1_sd = 1, nsparse 
    505466         DO jk2_sd = 1, nsparse 
    506  
    507             nisparse(jk2_sd)=nisparse(jk2_sd) 
    508             njsparse(jk2_sd)=njsparse(jk2_sd) 
    509  
     467            nisparse(jk2_sd) = nisparse(jk2_sd) 
     468            njsparse(jk2_sd) = njsparse(jk2_sd) 
    510469            IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
    511470               ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
    512471                                                        + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
    513472            ENDIF 
    514  
    515          ENDDO 
    516       ENDDO 
     473         END DO 
     474      END DO 
    517475 
    518476      DO jj_sd = 1 ,ninco 
     
    584542   ENDDO 
    585543 
    586  
    587544   CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
    588545   CALL wrk_dealloc( jpincomax , ipos2 , ipivot        ) 
     
    590547  END SUBROUTINE SUR_DETERMINE 
    591548 
    592  
    593549#else 
    594550   !!---------------------------------------------------------------------- 
     
    597553   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE. 
    598554CONTAINS 
    599  
    600555   SUBROUTINE dia_harm ( kt )     ! Empty routine 
    601556      INTEGER, INTENT( IN ) :: kt   
    602557      WRITE(*,*) 'dia_harm: you should not have seen this print' 
    603558   END SUBROUTINE dia_harm 
    604  
    605  
    606 #endif 
     559#endif 
     560 
    607561   !!====================================================================== 
    608562END MODULE diaharm 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r3625 r3953  
    2222   USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine) 
    2323   USE dynadv         ! dynamics: vector invariant versus flux form 
     24   USE sbctide 
     25   USE updtide 
    2426   USE trdmod         ! ocean dynamics trends 
    2527   USE trdmod_oce     ! ocean variables trends 
     
    100102      ENDIF 
    101103 
    102       IF( ln_apr_dyn ) THEN                   !==  Atmospheric pressure gradient  ==! 
    103          zg_2 = grav * 0.5 
    104          DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh 
     104      IF(      ln_apr_dyn                                                &   ! atmos. pressure 
     105         .OR.  ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting) 
     106         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice 
     107         ! 
     108         DO jj = 2, jpjm1 
    105109            DO ji = fs_2, fs_jpim1   ! vector opt. 
    106                spgu(ji,jj) =  zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    & 
    107                   &                   + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
    108                spgv(ji,jj) =  zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    & 
    109                   &                   + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
    110             END DO 
    111          END DO 
    112          DO jk = 1, jpkm1                          ! Add the apg to the general trend 
     110               spgu(ji,jj) = 0._wp 
     111               spgv(ji,jj) = 0._wp 
     112            END DO 
     113         END DO          
     114         ! 
     115         IF( ln_apr_dyn ) THEN                !==  Atmospheric pressure gradient  ==! 
     116            zg_2 = grav * 0.5 
     117            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh 
     118               DO ji = fs_2, fs_jpim1   ! vector opt. 
     119                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    & 
     120                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
     121                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    & 
     122                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
     123               END DO 
     124            END DO 
     125         ENDIF 
     126         ! 
     127         !                                    !==  tide potential forcing term  ==! 
     128         IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case 
     129            ! 
     130            CALL upd_tide( kt )                      ! update tide potential 
     131            ! 
     132            DO jj = 2, jpjm1                         ! add tide potential forcing 
     133               DO ji = fs_2, fs_jpim1   ! vector opt. 
     134                  spgv(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     135                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     136               END DO  
     137            END DO 
     138         ENDIF 
     139         ! 
     140         IF( nn_ice_embd == 2 ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==! 
     141            CALL wrk_alloc( jpi, jpj, zpice ) 
     142            !                                             
     143            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc ) 
     144            zgrau0r     = - grav * r1_rau0 
     145            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r 
     146            DO jj = 2, jpjm1 
     147               DO ji = fs_2, fs_jpim1   ! vector opt. 
     148                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
     149                  spgv(ji,jj) = spgu(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
     150               END DO 
     151            END DO 
     152            ! 
     153            CALL wrk_dealloc( jpi, jpj, zpice )          
     154         ENDIF 
     155         ! 
     156         DO jk = 1, jpkm1                     !== Add all terms to the general trend 
    113157            DO jj = 2, jpjm1 
    114158               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    117161               END DO 
    118162            END DO 
    119          END DO 
    120       ENDIF 
    121  
    122       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: Pressure gradient due to snow-ice mass ==! 
    123          CALL wrk_alloc( jpi, jpj, zpice ) 
    124          !                                             
    125          zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc ) 
    126          zgrau0r     = - grav * r1_rau0 
    127          zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r 
    128          DO jj = 2, jpjm1 
    129             DO ji = fs_2, fs_jpim1   ! vector opt. 
    130                spgu(ji,jj) = ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
    131                spgv(ji,jj) = ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
    132             END DO 
    133          END DO 
    134          DO jk = 1, jpkm1                             ! Add the surface pressure trend to the general trend 
    135             DO jj = 2, jpjm1 
    136                DO ji = fs_2, fs_jpim1   ! vector opt. 
    137                   ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
    138                   va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
    139                END DO 
    140             END DO 
    141          END DO 
    142          ! 
    143          CALL wrk_dealloc( jpi, jpj, zpice ) 
    144       ENDIF 
    145  
     163         END DO          
     164      ENDIF 
    146165 
    147166      SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r3680 r3953  
    9191               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 
    9292            END DO  
    93          END DO  
     93         END DO 
     94         ! 
    9495         DO jk = 1, jpkm1                    ! Add it to the general trend 
    9596            DO jj = 2, jpjm1 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r3951 r3953  
    396396         !                                                !* Update the forcing (BDY and tides) 
    397397         !                                                !  ------------------ 
    398          IF( lk_obc )   CALL obc_dta_bt ( kt, jn   ) 
    399          IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
    400          IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn ) 
     398         IF( lk_obc                    )   CALL obc_dta_bt( kt, jn   ) 
     399         IF( lk_bdy                    )   CALL bdy_dta   ( kt, jit=jn, time_offset=1 ) 
     400         IF( ln_tide_pot .AND. lk_tide )   CALL upd_tide  ( kt, kit=jn, kbaro=nn_baro ) 
    401401 
    402402         !                                                !* after ssh_e 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90

    r3651 r3953  
    11MODULE sbctide 
    2   !!================================================================================= 
    3   !!                       ***  MODULE  sbctide  *** 
    4   !! Initialization of tidal forcing 
    5   !! History :  9.0  !  07  (O. Le Galloudec)  Original code 
    6   !!================================================================================= 
    7   !! * Modules used 
    8   USE oce             ! ocean dynamics and tracers variables 
    9   USE dom_oce         ! ocean space and time domain 
    10   USE in_out_manager  ! I/O units 
    11   USE ioipsl          ! NetCDF IPSL library 
    12   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    13   USE phycst 
    14   USE daymod 
    15   USE dynspg_oce 
    16   USE tideini 
    17   USE iom 
     2   !!====================================================================== 
     3   !!                       ***  MODULE  sbctide  *** 
     4   !! Initialization of tidal forcing 
     5   !!====================================================================== 
     6   !! History :  9.0  !  2007  (O. Le Galloudec)  Original code 
     7   !!---------------------------------------------------------------------- 
     8   USE oce             ! ocean dynamics and tracers variables 
     9   USE dom_oce         ! ocean space and time domain 
     10   USE phycst 
     11   USE daymod 
     12   USE dynspg_oce 
     13   USE tideini 
     14   ! 
     15   USE iom 
     16   USE in_out_manager  ! I/O units 
     17   USE ioipsl          ! NetCDF IPSL library 
     18   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1819 
    19   IMPLICIT NONE 
    20   PUBLIC 
     20   IMPLICIT NONE 
     21   PUBLIC 
    2122 
    22   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro 
     23   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   pot_astro   ! 
    2324 
    2425#if defined key_tide 
     26   !!---------------------------------------------------------------------- 
     27   !!   'key_tide' :                                        tidal potential 
     28   !!---------------------------------------------------------------------- 
     29   !!   sbc_tide            :  
     30   !!   tide_init_potential : 
     31   !!---------------------------------------------------------------------- 
    2532 
    26   LOGICAL, PUBLIC, PARAMETER ::   lk_tide  = .TRUE. 
    27   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot,phi_pot 
    28   !!--------------------------------------------------------------------------------- 
    29   !!   OPA 9.0 , LODYC-IPSL  (2003) 
    30   !!--------------------------------------------------------------------------------- 
     33   LOGICAL, PUBLIC, PARAMETER ::   lk_tide  = .TRUE. 
     34   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   amp_pot, phi_pot 
    3135 
     36   !!---------------------------------------------------------------------- 
     37   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
     38   !! $Id: $ 
     39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     40   !!---------------------------------------------------------------------- 
    3241CONTAINS 
    3342 
    34   SUBROUTINE sbc_tide ( kt ) 
    35     !!---------------------------------------------------------------------- 
    36     !!                 ***  ROUTINE sbc_tide  *** 
    37     !!----------------------------------------------------------------------       
    38     !! * Arguments 
    39     INTEGER, INTENT( in ) ::   kt     ! ocean time-step 
    40     !!---------------------------------------------------------------------- 
     43   SUBROUTINE sbc_tide( kt ) 
     44      !!---------------------------------------------------------------------- 
     45      !!                 ***  ROUTINE sbc_tide  *** 
     46      !!----------------------------------------------------------------------       
     47      INTEGER, INTENT( in ) ::   kt     ! ocean time-step 
     48      !!---------------------------------------------------------------------- 
    4149 
    42     IF ( kt == nit000 .AND. .NOT. lk_dynspg_ts )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' ) 
     50      IF( kt /= nit000 )   CALL tide_init( kt ) 
    4351 
    44     IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN 
     52      IF( nsec_day == NINT(0.5 * rdttra(1)) ) THEN      ! start a new day 
     53         ! 
     54         IF( kt == nit000 ) THEN 
     55            ALLOCATE( amp_pot(jpi,jpj,nb_harmo),                      & 
     56               &      phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj)   ) 
     57            ! 
     58            amp_pot(:,:,:) = 0._wp 
     59            phi_pot(:,:,:) = 0._wp 
     60            pot_astro(:,:) = 0._wp 
     61         ENDIF 
     62         ! 
     63         CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo ) 
     64         ! 
     65         kt_tide = kt 
     66         ! 
     67         IF(lwp) THEN 
     68            WRITE(numout,*) 
     69            WRITE(numout,*) 'sbc_tide : Update of the components and (re)Init. the potential at kt=', kt 
     70            WRITE(numout,*) '~~~~~~~~ ' 
     71            DO jk = 1, nb_harmo 
     72               WRITE(numout,*) Wave(ntide(jk))%cname_tide, utide(jk), ftide(jk), v0tide(jk), omega_tide(jk) 
     73            END DO 
     74         ENDIF 
     75         ! 
     76         IF( ln_tide_pot )   CALL tide_init_potential 
     77         ! 
     78      ENDIF 
    4579      ! 
    46       kt_tide = kt 
    47  
    48       IF(lwp) THEN 
    49          WRITE(numout,*) 
    50          WRITE(numout,*) 'sbc_tide : (re)Initialization of the tidal potential at kt=',kt 
    51          WRITE(numout,*) '~~~~~~~ ' 
    52       ENDIF 
    53  
    54       IF(lwp) THEN 
    55          IF ( kt == nit000 ) WRITE(numout,*) 'Apply astronomical potential : ln_tide_pot =', ln_tide_pot 
    56          CALL flush(numout) 
    57       ENDIF 
    58  
    59       IF ( kt == nit000 ) ALLOCATE(amp_pot(jpi,jpj,nb_harmo)) 
    60       IF ( kt == nit000 ) ALLOCATE(phi_pot(jpi,jpj,nb_harmo)) 
    61       IF ( kt == nit000 ) ALLOCATE(pot_astro(jpi,jpj)) 
    62  
    63       amp_pot(:,:,:) = 0.e0 
    64       phi_pot(:,:,:) = 0.e0 
    65       pot_astro(:,:) = 0.e0 
    66  
    67       IF ( ln_tide_pot ) CALL tide_init_potential 
    68       ! 
    69     ENDIF 
    70  
    71   END SUBROUTINE sbc_tide 
    72  
    73   SUBROUTINE tide_init_potential 
    74     !!---------------------------------------------------------------------- 
    75     !!                 ***  ROUTINE tide_init_potential  *** 
    76     !!---------------------------------------------------------------------- 
    77     !! * Local declarations 
    78     INTEGER  :: ji,jj,jk 
    79     REAL(wp) :: zcons,ztmp1,ztmp2,zlat,zlon 
     80   END SUBROUTINE sbc_tide 
    8081 
    8182 
    82     DO jk=1,nb_harmo 
    83        zcons=0.7*Wave(ntide(jk))%equitide*ftide(jk) 
    84        do ji=1,jpi 
    85           do jj=1,jpj 
    86              ztmp1 = amp_pot(ji,jj,jk)*COS(phi_pot(ji,jj,jk)) 
    87              ztmp2 = -amp_pot(ji,jj,jk)*SIN(phi_pot(ji,jj,jk)) 
    88              zlat = gphit(ji,jj)*rad !! latitude en radian 
    89              zlon = glamt(ji,jj)*rad !! longitude en radian 
    90              ! le potentiel est composé des effets des astres: 
    91              IF (Wave(ntide(jk))%nutide .EQ.1) THEN 
    92                 ztmp1= ztmp1 + zcons*(SIN(2.*zlat))*COS(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon) 
    93                 ztmp2= ztmp2 - zcons*(SIN(2.*zlat))*SIN(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon) 
    94              ENDIF 
    95              IF (Wave(ntide(jk))%nutide.EQ.2) THEN 
    96                 ztmp1= ztmp1 + zcons*(COS(zlat)**2)*COS(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon) 
    97                 ztmp2= ztmp2 - zcons*(COS(zlat)**2)*SIN(v0tide(jk)+utide(jk)+Wave(ntide(jk))%nutide*zlon) 
    98              ENDIF 
    99              amp_pot(ji,jj,jk)=SQRT(ztmp1**2+ztmp2**2) 
    100              phi_pot(ji,jj,jk)=ATAN2(-ztmp2/MAX(1.E-10,SQRT(ztmp1**2+ztmp2**2)),ztmp1/MAX(1.E-10,SQRT(ztmp1**2+ztmp2**2))) 
    101           enddo 
    102        enddo 
    103     END DO 
     83   SUBROUTINE tide_init_potential 
     84      !!---------------------------------------------------------------------- 
     85      !!                 ***  ROUTINE tide_init_potential  *** 
     86      !!---------------------------------------------------------------------- 
     87      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     88      REAL(wp) ::   zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zsc   ! local scalar 
     89      !!---------------------------------------------------------------------- 
    10490 
    105   END SUBROUTINE tide_init_potential 
     91      DO jk = 1, nb_harmo 
     92         zcons = 0.7 * Wave(ntide(jk))%equitide * ftide(jk) 
     93         DO ji = 1, jpi 
     94            DO jj = 1, jpj 
     95               ztmp1 =  amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) ) 
     96               ztmp2 = -amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) ) 
     97               zlat = gphit(ji,jj)*rad !! latitude en radian 
     98               zlon = glamt(ji,jj)*rad !! longitude en radian 
     99               ztmp = v0tide(jk) + utide(jk) + Wave(ntide(jk))%nutide * zlon 
     100               ! le potentiel est composé des effets des astres: 
     101               IF( Wave(ntide(jk))%nutide == 1 )   zcs = zcons * SIN( 2.*zlat ) 
     102               IF( Wave(ntide(jk))%nutide == 2 )   zcs = zcons * COS( zlat )**2 
     103               ztmp1 = ztmp1 + zcs * COS( ztmp ) 
     104               ztmp2 = ztmp2 - zcs * SIN( ztmp ) 
     105               zamp = SQRT( ztmp1*ztmp1 + ztmp2*ztmp2 ) 
     106               amp_pot(ji,jj,jk) = zamp 
     107               phi_pot(ji,jj,jk) = ATAN2( -ztmp2 / MAX( 1.e-10 , zamp ) ,   & 
     108                  &                        ztmp1 / MAX( 1.e-10,  zamp )   ) 
     109            END DO 
     110         END DO 
     111      END DO 
     112      ! 
     113   END SUBROUTINE tide_init_potential 
    106114 
    107115#else 
     
    116124  END SUBROUTINE sbc_tide 
    117125#endif 
     126 
    118127  !!====================================================================== 
    119  
    120128END MODULE sbctide 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90

    r3294 r3953  
    1 !! History :  9.0  !  07  (O. Le Galloudec)  Original code 
     1   !!---------------------------------------------------------------------- 
     2   !! History :  3.2  !  2007  (O. Le Galloudec)  Original code 
     3   !!---------------------------------------------------------------------- 
    24 
    3 ! Wave(1)= tide(name_tide,equitide,nutide,nt,ns,nh,np,np1,shift,nksi,nnu0,nnu1,nnu2,R,formula) 
    4   
    5    
    6   Wave(1)= tide('M2'     ,0.242297,2     ,2 ,-2,2 ,0 ,0  ,0    ,2   ,-2   ,0  ,0   ,0,78) 
    7   Wave(2)= tide('N2'     ,0.046313,2     ,2 ,-3,2 ,1 ,0  ,0    ,2   ,-2   ,0  ,0   ,0,78) 
    8   Wave(3)= tide('2N2'    ,0.006184,2     ,2 ,-4,2 ,2 ,0  ,0    ,2   ,-2   ,0  ,0   ,0,78) 
    9   Wave(4)= tide('S2'     ,0.113572,2     ,2 , 0,0 ,0 ,0  ,0    ,0   , 0   ,0  ,0   ,0,0) 
    10   Wave(5)= tide('K2'     ,0.030875,2     ,2 , 0,2 ,0 ,0  ,0    ,0   , 0   ,0  ,-2  ,0,235) 
    11  
    12   Wave(6)= tide('K1'     ,0.142408,1     ,1 , 0,1 ,0 ,0  ,-90  ,0   , 0   ,-1 ,0   ,0,227) 
    13   Wave(7)= tide('O1'     ,0.101266,1     ,1 ,-2,1 ,0 ,0  ,+90  ,2   ,-1   , 0 ,0   ,0,75) 
    14   Wave(8)= tide('Q1'     ,0.019387,1     ,1 ,-3,1 ,1 ,0  ,+90  ,2   ,-1   , 0 ,0   ,0,75) 
    15   Wave(9)= tide('P1'     ,0.047129,1     ,1 , 0,-1,0 ,0  ,+90  ,0   , 0   , 0 ,0   ,0,0)   
    16  
    17   Wave(10)= tide('M4'    ,0.000000,4     ,4 ,-4, 4,0 ,0  ,0    ,4   , -4  , 0 ,0   ,0,1) 
    18  
    19   Wave(11) = tide('Mf'   ,0.042017,0     ,0 , 2, 0,0 ,0  ,0    ,-2  , 0   , 0 ,0   ,0,74) 
    20   Wave(12) = tide('Mm'   ,0.022191,0     ,0 , 1,0 ,-1,0  ,0    ,0   , 0   , 0 ,0   ,0,73) 
    21   Wave(13) = tide('Msqm' ,0.000667,0     ,0 , 4,-2, 0,0  ,0    ,-2  , 0   , 0 ,0   ,0,74) 
    22   Wave(14) = tide('Mtm'  ,0.008049,0     ,0 , 3, 0,-1,0  ,0    ,-2  , 0   , 0 ,0   ,0,74) 
    23  
    24   Wave(15) = tide('S1'   ,0.000000,1     ,1,  0, 0, 0,0  ,0    , 0  , 0   , 0 ,0   ,0,0)    
    25   Wave(16) = tide('MU2'  ,0.005841,2     ,2, -4, 4, 0,0  ,0    ,2   ,-2   , 0, 0   ,0,78) 
    26   Wave(17) = tide('NU2'  ,0.009094,2     ,2, -3, 4,-1,0  ,0    ,2   ,-2   , 0, 0   ,0,78)  
    27   Wave(18) = tide('L2'   ,0.006694,2     ,2, -1, 2,-1,0  ,+180 ,2   ,-2   , 0, 0   ,0,215) 
    28   Wave(19) = tide('T2'   ,0.006614,2     ,2,  0,-1, 0,1  ,0    ,0   , 0   , 0, 0   ,0,0) 
     5   !             !! name_tide , equitide , nutide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !! 
     6   !             !!           !          !        !    !    !    !    !     !       !      !      !      !      !   !         !! 
     7   Wave( 1) = tide(  'M2'     , 0.242297 ,    2   ,  2 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   ) 
     8   Wave( 2) = tide(  'N2'     , 0.046313 ,    2   ,  2 , -3 ,  2 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   ) 
     9   Wave( 3) = tide( '2N2'     , 0.006184 ,    2   ,  2 , -4 ,  2 ,  2 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   ) 
     10   Wave( 4) = tide(  'S2'     , 0.113572 ,    2   ,  2 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   ) 
     11   Wave( 5) = tide(  'K2'     , 0.030875 ,    2   ,  2 ,  0 ,  2 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   , -2   , 0 ,   235   ) 
     12   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     13   Wave( 6) = tide(  'K1'     , 0.142408 ,    1   ,  1 ,  0 ,  1 ,  0 ,  0  ,  -90  ,  0   ,  0   , -1   ,  0   , 0 ,   227   ) 
     14   Wave( 7) = tide(  'O1'     , 0.101266 ,    1   ,  1 , -2 ,  1 ,  0 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   ) 
     15   Wave( 8) = tide(  'Q1'     , 0.019387 ,    1   ,  1 , -3 ,  1 ,  1 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   ) 
     16   Wave( 9) = tide(  'P1'     , 0.047129 ,    1   ,  1 ,  0 , -1 ,  0 ,  0  ,  +90  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    ) 
     17   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     18   Wave(10) = tide(  'M4'     , 0.000000 ,    4   ,  4 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    1    ) 
     19   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     20   Wave(11) = tide(  'Mf'     , 0.042017 ,    0   ,  0 ,  2 ,  0 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    ) 
     21   Wave(12) = tide(  'Mm'     , 0.022191 ,    0   ,  0 ,  1 ,  0 , -1 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,   73    ) 
     22   Wave(13) = tide(  'Msqm'   , 0.000667 ,    0   ,  0 ,  4 , -2 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    ) 
     23   Wave(14) = tide(  'Mtm'    , 0.008049 ,    0   ,  0 ,  3 ,  0 , -1 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    ) 
     24   !              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     25   Wave(15) = tide(  'S1'     , 0.000000 ,    1   ,  1 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )    
     26   Wave(16) = tide(  'MU2'    , 0.005841 ,    2   ,  2 , -4 ,  4 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    ) 
     27   Wave(17) = tide(  'NU2'    , 0.009094 ,    2   ,  2 , -3 ,  4 , -1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    )  
     28   Wave(18) = tide(  'L2'     , 0.006694 ,    2   ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,  215    ) 
     29   Wave(19) = tide(  'T2'     , 0.006614 ,    2   ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    ) 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90

    r3670 r3953  
    11MODULE tide_mod 
    2   !!================================================================================= 
    3   !!                       ***  MODULE  tide_mod  *** 
    4   !! Compute nodal modulations corrections and pulsations 
    5   !!================================================================================= 
    6   !!--------------------------------------------------------------------------------- 
    7   !!   OPA 9.0 , LODYC-IPSL  (2003) 
    8   !!--------------------------------------------------------------------------------- 
    9   USE dom_oce         ! ocean space and time domain 
    10   USE phycst 
    11   USE daymod 
    12  
    13   IMPLICIT NONE 
    14   PRIVATE 
    15  
    16   REAL(wp) :: sh_T, sh_s, sh_h, sh_p, sh_p1, & 
    17        sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R, & 
    18        sh_I, sh_x1ra, sh_N 
    19  
    20   INTEGER,PUBLIC, PARAMETER ::   & 
    21        jpmax_harmo = 19             ! maximum number of harmonic 
    22  
    23   TYPE, PUBLIC ::    tide 
    24      CHARACTER(LEN=4)  :: cname_tide 
    25      REAL(wp) :: equitide 
    26      INTEGER  :: nutide 
    27      INTEGER  ::  nt,ns,nh,np,np1,shift 
    28      INTEGER  ::  nksi,nnu0,nnu1,nnu2,R 
    29      INTEGER  :: nformula 
    30   END TYPE tide 
    31  
    32   TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) :: Wave 
    33  
    34   !! * Accessibility 
    35   PUBLIC tide_harmo 
    36   PUBLIC nodal_factort 
    37   PUBLIC tide_init_Wave 
    38  
     2   !!====================================================================== 
     3   !!                       ***  MODULE  tide_mod  *** 
     4   !! Compute nodal modulations corrections and pulsations 
     5   !!====================================================================== 
     6   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code 
     7   !!---------------------------------------------------------------------- 
     8   USE dom_oce        ! ocean space and time domain 
     9   USE phycst         ! physical constant 
     10   USE daymod         ! calendar 
     11 
     12   IMPLICIT NONE 
     13   PRIVATE 
     14 
     15   PUBLIC   tide_harmo       ! called by tideini and diaharm modules 
     16   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules 
     17 
     18   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic 
     19 
     20   TYPE, PUBLIC ::    tide 
     21      CHARACTER(LEN=4) ::   cname_tide 
     22      REAL(wp)         ::   equitide 
     23      INTEGER          ::   nutide 
     24      INTEGER          ::   nt, ns, nh, np, np1, shift 
     25      INTEGER          ::   nksi, nnu0, nnu1, nnu2, R 
     26      INTEGER          ::   nformula 
     27   END TYPE tide 
     28 
     29   TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) ::   Wave   !: 
     30 
     31   REAL(wp) ::   sh_T, sh_s, sh_h, sh_p, sh_p1             ! astronomic angles 
     32   REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   ! 
     33   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       ! 
     34 
     35   !!---------------------------------------------------------------------- 
     36   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
     37   !! $Id:$  
     38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     39   !!---------------------------------------------------------------------- 
    3940CONTAINS 
    4041 
    41   SUBROUTINE tide_init_Wave 
    42  
    43 #  include "tide.h90" 
    44  
    45   END SUBROUTINE tide_init_Wave 
    46  
    47   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc) 
    48  
    49     INTEGER, DIMENSION(kc), INTENT( in ) ::   & 
    50          ktide      ! Indice of tidal constituents 
    51  
    52     INTEGER, INTENT( in ) :: & 
    53          kc         ! Total number of tidal constituents 
    54  
    55     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    56          pomega      ! pulsation in radians/s 
    57  
    58     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    59          pvt, &      ! 
    60          put, &      ! 
    61          pcor         ! 
    62  
    63     CALL astronomic_angle 
    64     CALL tide_pulse(pomega, ktide ,kc) 
    65     CALL tide_vuf( pvt, put, pcor, ktide ,kc) 
    66  
    67   END SUBROUTINE tide_harmo 
    68  
    69   SUBROUTINE astronomic_angle 
    70  
    71     !!---------------------------------------------------------------------- 
    72     !! 
    73     !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian 
    74     !!  century (e.g. time in days divide by 36525) 
    75     !!---------------------------------------------------------------------- 
    76  
    77     REAL(wp) ::  cosI,p,q,t2,t4,sin2I,s2,tgI2,P1,sh_tgn2,at1,at2 
    78     REAL(wp) :: zqy,zsy,zday,zdj,zhfrac 
    79  
    80     zqy=AINT((nyear-1901.)/4.) 
    81     zsy=nyear-1900. 
    82  
    83     zdj=dayjul(nyear,nmonth,nday) 
    84     zday=zdj+zqy-1. 
    85  
    86     zhfrac=nsec_day/3600. 
    87  
    88     !---------------------------------------------------------------------- 
    89     !  Sh_n Longitude of ascending lunar node 
    90     !---------------------------------------------------------------------- 
    91  
    92     sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad 
    93     !---------------------------------------------------------------------- 
    94     ! T mean solar angle (Greenwhich time) 
    95     !---------------------------------------------------------------------- 
    96     sh_T=(180.+zhfrac*(360./24.))*rad 
    97     !---------------------------------------------------------------------- 
    98     ! h mean solar Longitude 
    99     !---------------------------------------------------------------------- 
    100  
    101     sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad 
    102     !---------------------------------------------------------------------- 
    103     ! s mean lunar Longitude 
    104     !---------------------------------------------------------------------- 
    105  
    106     sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad 
    107     !---------------------------------------------------------------------- 
    108     ! p1 Longitude of solar perigee 
    109     !---------------------------------------------------------------------- 
    110  
    111     sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad 
    112     !---------------------------------------------------------------------- 
    113     ! p Longitude of lunar perigee 
    114     !---------------------------------------------------------------------- 
    115  
    116     sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad 
    117  
    118     sh_N =mod(sh_N ,2*rpi) 
    119     sh_s =mod(sh_s ,2*rpi) 
    120     sh_h =mod(sh_h, 2*rpi) 
    121     sh_p =mod(sh_p, 2*rpi) 
    122     sh_p1=mod(sh_p1,2*rpi) 
    123  
    124     cosI=0.913694997 -0.035692561 *cos(sh_N) 
    125  
    126     sh_I=acos(cosI) 
    127  
    128     sin2I=sin(sh_I) 
    129     sh_tgn2=tan(sh_N/2.0) 
    130  
    131     at1=atan(1.01883*sh_tgn2) 
    132     at2=atan(0.64412*sh_tgn2) 
    133  
    134     sh_xi=-at1-at2+sh_N 
    135  
    136     if (sh_N > rpi) sh_xi=sh_xi-2.0*rpi 
    137  
    138     sh_nu=at1-at2 
    139  
    140     !---------------------------------------------------------------------- 
    141     ! For constituents l2 k1 k2 
    142     !---------------------------------------------------------------------- 
    143  
    144     tgI2=tan(sh_I/2.0) 
    145     P1=sh_p-sh_xi 
    146  
    147     t2=tgI2*tgI2 
    148     t4=t2*t2 
    149     sh_x1ra=sqrt(1.0-12.0*t2*cos(2.0*P1)+36.0*t4) 
    150  
    151     p=sin(2.0*P1) 
    152     q=1.0/(6.0*t2)-cos(2.0*P1) 
    153     sh_R=atan(p/q) 
    154  
    155     p=sin(2.0*sh_I)*sin(sh_nu) 
    156     q=sin(2.0*sh_I)*cos(sh_nu)+0.3347 
    157     sh_nuprim=atan(p/q) 
    158  
    159     s2=sin(sh_I)*sin(sh_I) 
    160     p=s2*sin(2.0*sh_nu) 
    161     q=s2*cos(2.0*sh_nu)+0.0727 
    162     sh_nusec=0.5*atan(p/q) 
    163  
    164   END SUBROUTINE astronomic_angle 
    165  
    166   SUBROUTINE tide_pulse( pomega, ktide ,kc) 
    167     !!---------------------------------------------------------------------- 
    168     !!                     ***  ROUTINE tide_pulse  *** 
    169     !!                       
    170     !! ** Purpose : Compute tidal frequencies 
    171     !! 
    172     !!---------------------------------------------------------------------- 
    173     !! * Arguments 
    174     INTEGER, DIMENSION(kc), INTENT( in ) ::   & 
    175          ktide      ! Indice of tidal constituents 
    176  
    177     INTEGER, INTENT( in ) :: & 
    178          kc         ! Total number of tidal constituents 
    179  
    180     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    181          pomega      ! pulsation in radians/s 
    182  
    183     !! * Local declarations 
    184     INTEGER :: jh 
    185     REAL(wp) :: zscale  =  36525*24.0 
    186     REAL(wp) :: zomega_T=  13149000.0 
    187     REAL(wp) :: zomega_s=    481267.892 
    188     REAL(wp) :: zomega_h=     36000.76892 
    189     REAL(wp) :: zomega_p=      4069.0322056 
    190     REAL(wp) :: zomega_n=      1934.1423972 
    191     REAL(wp) :: zomega_p1=        1.719175 
    192     !!---------------------------------------------------------------------- 
    193  
    194     DO jh=1,kc 
    195        pomega(jh) = zomega_T * Wave(ktide(jh))%nT & 
    196             + zomega_s * Wave(ktide(jh))%ns & 
    197             + zomega_h * Wave(ktide(jh))%nh & 
    198             + zomega_p * Wave(ktide(jh))%np & 
    199             + zomega_p1* Wave(ktide(jh))%np1 
    200        pomega(jh) = (pomega(jh)/zscale)*rad/3600. 
    201     END DO 
    202  
    203   END SUBROUTINE tide_pulse 
    204  
    205   SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc) 
    206     !!---------------------------------------------------------------------- 
    207     !!                     ***  ROUTINE tide_vuf  *** 
    208     !!                       
    209     !! ** Purpose : Compute nodal modulation corrections 
    210     !! 
    211     !! ** Outputs : 
    212     !!          vt: Pase of tidal potential relative to Greenwich (radians) 
    213     !!          ut: Phase correction u due to nodal motion (radians) 
    214     !!          ft: Nodal correction factor 
    215     !! 
    216     !! ** Inputs : 
    217     !!          tname: array of constituents names (dimension<=nc)  
    218     !!             nc: number of constituents 
    219     !!    
    220     !!---------------------------------------------------------------------- 
    221     !! * Arguments 
    222     INTEGER, DIMENSION(kc), INTENT( in ) ::   & 
    223          ktide      ! Indice of tidal constituents 
    224     INTEGER, INTENT( in ) :: & 
    225          kc         ! Total number of tidal constituents 
    226     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    227          pvt, &      ! 
    228          put, &      ! 
    229          pcor         ! 
    230     !! * Local declarations 
    231     INTEGER :: jh 
    232     !!---------------------------------------------------------------------- 
    233  
    234     DO jh =1,kc 
    235        !  Phase of the tidal potential relative to the Greenwhich  
    236        !  meridian (e.g. the position of the fictuous celestial body). Units are 
    237        !  radian: 
    238        pvt(jh) = sh_T *Wave(ktide(jh))%nT        & 
    239             +sh_s *Wave(ktide(jh))%ns        & 
    240             +sh_h *Wave(ktide(jh))%nh        & 
    241             +sh_p *Wave(ktide(jh))%np        & 
    242             +sh_p1*Wave(ktide(jh))%np1       & 
    243             +Wave(ktide(jh))%shift*rad 
     42   SUBROUTINE tide_init_Wave 
     43#     include "tide.h90" 
     44   END SUBROUTINE tide_init_Wave 
     45 
     46 
     47   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc) 
     48      !!---------------------------------------------------------------------- 
     49      !!---------------------------------------------------------------------- 
     50      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents 
     51      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents 
     52      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega           ! pulsation in radians/s 
     53      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   ! 
     54      !!---------------------------------------------------------------------- 
     55      ! 
     56      CALL astronomic_angle 
     57      CALL tide_pulse( pomega, ktide ,kc ) 
     58      CALL tide_vuf  ( pvt, put, pcor, ktide ,kc ) 
     59      ! 
     60   END SUBROUTINE tide_harmo 
     61 
     62 
     63   SUBROUTINE astronomic_angle 
     64      !!---------------------------------------------------------------------- 
     65      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian 
     66      !!  century (e.g. time in days divide by 36525) 
     67      !!---------------------------------------------------------------------- 
     68      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2 
     69      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac 
     70      !!---------------------------------------------------------------------- 
     71      ! 
     72      zqy = AINT( (nyear-1901.)/4. ) 
     73      zsy = nyear - 1900. 
     74      ! 
     75      zdj  = dayjul( nyear, nmonth, nday ) 
     76      zday = zdj + zqy - 1. 
     77      ! 
     78      zhfrac = nsec_day / 3600. 
     79      ! 
     80      !---------------------------------------------------------------------- 
     81      !  Sh_n Longitude of ascending lunar node 
     82      !---------------------------------------------------------------------- 
     83      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad 
     84      !---------------------------------------------------------------------- 
     85      ! T mean solar angle (Greenwhich time) 
     86      !---------------------------------------------------------------------- 
     87      sh_T=(180.+zhfrac*(360./24.))*rad 
     88      !---------------------------------------------------------------------- 
     89      ! h mean solar Longitude 
     90      !---------------------------------------------------------------------- 
     91      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad 
     92      !---------------------------------------------------------------------- 
     93      ! s mean lunar Longitude 
     94      !---------------------------------------------------------------------- 
     95      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad 
     96      !---------------------------------------------------------------------- 
     97      ! p1 Longitude of solar perigee 
     98      !---------------------------------------------------------------------- 
     99      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad 
     100      !---------------------------------------------------------------------- 
     101      ! p Longitude of lunar perigee 
     102      !---------------------------------------------------------------------- 
     103      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad 
     104 
     105      sh_N = MOD( sh_N ,2*rpi ) 
     106      sh_s = MOD( sh_s ,2*rpi ) 
     107      sh_h = MOD( sh_h, 2*rpi ) 
     108      sh_p = MOD( sh_p, 2*rpi ) 
     109      sh_p1= MOD( sh_p1,2*rpi ) 
     110 
     111      cosI = 0.913694997 -0.035692561 *cos(sh_N) 
     112 
     113      sh_I = ACOS( cosI ) 
     114 
     115      sin2I   = sin(sh_I) 
     116      sh_tgn2 = tan(sh_N/2.0) 
     117 
     118      at1=atan(1.01883*sh_tgn2) 
     119      at2=atan(0.64412*sh_tgn2) 
     120 
     121      sh_xi=-at1-at2+sh_N 
     122 
     123      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi 
     124 
     125      sh_nu = at1 - at2 
     126 
     127      !---------------------------------------------------------------------- 
     128      ! For constituents l2 k1 k2 
     129      !---------------------------------------------------------------------- 
     130 
     131      tgI2 = tan(sh_I/2.0) 
     132      P1   = sh_p-sh_xi 
     133 
     134      t2 = tgI2*tgI2 
     135      t4 = t2*t2 
     136      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 ) 
     137 
     138      p = sin(2.0*P1) 
     139      q = 1.0/(6.0*t2)-cos(2.0*P1) 
     140      sh_R = atan(p/q) 
     141 
     142      p = sin(2.0*sh_I)*sin(sh_nu) 
     143      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347 
     144      sh_nuprim = atan(p/q) 
     145 
     146      s2 = sin(sh_I)*sin(sh_I) 
     147      p  = s2*sin(2.0*sh_nu) 
     148      q  = s2*cos(2.0*sh_nu)+0.0727 
     149      sh_nusec = 0.5*atan(p/q) 
     150      ! 
     151   END SUBROUTINE astronomic_angle 
     152 
     153 
     154   SUBROUTINE tide_pulse( pomega, ktide ,kc ) 
     155      !!---------------------------------------------------------------------- 
     156      !!                     ***  ROUTINE tide_pulse  *** 
     157      !!                       
     158      !! ** Purpose : Compute tidal frequencies 
     159      !!---------------------------------------------------------------------- 
     160      INTEGER                , INTENT(in ) ::   kc       ! Total number of tidal constituents 
     161      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide    ! Indice of tidal constituents 
     162      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega   ! pulsation in radians/s 
     163      ! 
     164      INTEGER  ::   jh 
     165      REAL(wp) ::   zscale 
     166      REAL(wp) ::   zomega_T =  13149000.0_wp 
     167      REAL(wp) ::   zomega_s =    481267.892_wp 
     168      REAL(wp) ::   zomega_h =     36000.76892_wp 
     169      REAL(wp) ::   zomega_p =      4069.0322056_wp 
     170      REAL(wp) ::   zomega_n =      1934.1423972_wp 
     171      REAL(wp) ::   zomega_p1=         1.719175_wp 
     172      !!---------------------------------------------------------------------- 
     173      ! 
     174      zscale =  rad / ( 36525._wp * 86400._wp )  
     175      ! 
     176      DO jh = 1, kc 
     177         pomega(jh) = (  zomega_T * Wave( ktide(jh) )%nT   & 
     178            &          + zomega_s * Wave( ktide(jh) )%ns   & 
     179            &          + zomega_h * Wave( ktide(jh) )%nh   & 
     180            &          + zomega_p * Wave( ktide(jh) )%np   & 
     181            &          + zomega_p1* Wave( ktide(jh) )%np1  ) * zscale 
     182      END DO 
     183      ! 
     184   END SUBROUTINE tide_pulse 
     185 
     186 
     187   SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc ) 
     188      !!---------------------------------------------------------------------- 
     189      !!                     ***  ROUTINE tide_vuf  *** 
     190      !!                       
     191      !! ** Purpose : Compute nodal modulation corrections 
     192      !! 
     193      !! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians) 
     194      !!              ut: Phase correction u due to nodal motion (radians) 
     195      !!              ft: Nodal correction factor 
     196      !!---------------------------------------------------------------------- 
     197      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents 
     198      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents 
     199      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   ! 
     200      ! 
     201      INTEGER ::   jh   ! dummy loop index 
     202      !!---------------------------------------------------------------------- 
     203      ! 
     204      DO jh = 1, kc 
     205         !  Phase of the tidal potential relative to the Greenwhich  
     206         !  meridian (e.g. the position of the fictuous celestial body). Units are radian: 
     207         pvt(jh) = sh_T * Wave( ktide(jh) )%nT    & 
     208            &    + sh_s * Wave( ktide(jh) )%ns    & 
     209            &    + sh_h * Wave( ktide(jh) )%nh    & 
     210            &    + sh_p * Wave( ktide(jh) )%np    & 
     211            &    + sh_p1* Wave( ktide(jh) )%np1   & 
     212            &    +        Wave( ktide(jh) )%shift * rad 
     213         ! 
     214         !  Phase correction u due to nodal motion. Units are radian: 
     215         put(jh) = sh_xi     * Wave( ktide(jh) )%nksi   & 
     216            &    + sh_nu     * Wave( ktide(jh) )%nnu0   & 
     217            &    + sh_nuprim * Wave( ktide(jh) )%nnu1   & 
     218            &    + sh_nusec  * Wave( ktide(jh) )%nnu2   & 
     219            &    + sh_R      * Wave( ktide(jh) )%R 
     220 
     221         !  Nodal correction factor: 
     222         pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula ) 
     223      END DO 
     224      ! 
     225   END SUBROUTINE tide_vuf 
     226 
     227 
     228   RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf ) 
     229      !!---------------------------------------------------------------------- 
     230      !!---------------------------------------------------------------------- 
     231      INTEGER, INTENT(in) :: kformula 
     232      ! 
     233      REAL(wp) :: zf 
     234      REAL(wp) :: zs, zf1, zf2 
     235      !!---------------------------------------------------------------------- 
     236      ! 
     237      SELECT CASE( kformula ) 
     238      ! 
     239      CASE( 0 )                  !==  formule 0, solar waves 
     240         zf = 1.0 
     241         ! 
     242      CASE( 1 )                  !==  formule 1, compound waves (78 x 78) 
     243         zf=nodal_factort(78) 
     244         zf = zf * zf 
     245         ! 
     246      CASE ( 2 )                 !==  formule 2, compound waves (78 x 0)  ===  (78)  
     247       zf1= nodal_factort(78) 
     248       zf = nodal_factort( 0) 
     249       zf = zf1 * zf 
    244250       ! 
    245        !  Phase correction u due to nodal motion. Units are radian: 
    246        put(jh) = sh_xi    *Wave(ktide(jh))%nksi  & 
    247             +sh_nu    *Wave(ktide(jh))%nnu0  & 
    248             +sh_nuprim*Wave(ktide(jh))%nnu1  & 
    249             +sh_nusec *Wave(ktide(jh))%nnu2  & 
    250             +sh_R     *Wave(ktide(jh))%R 
    251  
    252        !  Nodal correction factor: 
    253        pcor(jh) = nodal_factort(Wave(ktide(jh))%nformula) 
    254     END DO 
    255  
    256   END SUBROUTINE tide_vuf 
    257  
    258   recursive function nodal_factort(kformula) result (zf) 
    259     !!---------------------------------------------------------------------- 
    260     INTEGER, INTENT(IN) :: kformula 
    261     REAL(wp) :: zf 
    262     REAL(wp) :: zs,zf1,zf2 
    263  
    264     SELECT CASE (kformula) 
    265  
    266        !!  formule 0, solar waves 
    267  
    268     case ( 0 ) 
    269        zf=1.0 
    270  
    271        !! formule 1, compound waves (78 x 78) 
    272  
    273     case ( 1 ) 
    274        zf=nodal_factort(78) 
    275        zf=zf*zf 
    276  
    277        !! formule 2, compound waves (78 x 0)  ===  (78)  
    278  
    279     case ( 2 ) 
    280        zf1=nodal_factort(78) 
    281        zf=nodal_factort(0) 
    282        zf=zf1*zf 
    283  
    284        !! formule 4,  compound waves (78 x 235)  
    285  
    286     case ( 4 ) 
    287        zf1=nodal_factort(78) 
    288        zf=nodal_factort(235) 
    289        zf=zf1*zf 
    290  
    291        !! formule 5,  compound waves (78 *78 x 235) 
    292  
    293     case ( 5 ) 
    294        zf1=nodal_factort(78) 
    295        zf=nodal_factort(235) 
    296        zf=zf*zf1*zf1 
    297  
    298        !! formule 6,  compound waves (78 *78 x 0) 
    299  
    300     case ( 6 ) 
    301        zf1=nodal_factort(78) 
    302        zf=nodal_factort(0) 
    303        zf=zf*zf1*zf1  
    304  
    305        !! formule 7, compound waves (75 x 75) 
    306  
    307     case ( 7 ) 
    308        zf=nodal_factort(75) 
    309        zf=zf*zf 
    310  
    311        !! formule 8,  compound waves (78 x 0 x 235) 
    312  
    313     case ( 8 ) 
    314        zf=nodal_factort(78) 
    315        zf1=nodal_factort(0) 
    316        zf2=nodal_factort(235) 
    317        zf=zf*zf1*zf2 
    318  
    319        !! formule 9,  compound waves (78 x 0 x 227) 
    320  
    321     case ( 9 ) 
    322        zf=nodal_factort(78) 
    323        zf1=nodal_factort(0) 
    324        zf2=nodal_factort(227) 
    325        zf=zf*zf1*zf2 
    326  
    327        !! formule 10,  compound waves (78 x 227) 
    328  
    329     case ( 10 ) 
    330        zf=nodal_factort(78) 
    331        zf1=nodal_factort(227) 
    332        zf=zf*zf1 
    333  
    334        !! formule 11,  compound waves (75 x 0) 
    335  
    336     case ( 11 ) 
    337        zf=nodal_factort(75) 
    338        zf=nodal_factort(0) 
    339        zf=zf*zf1 
    340  
    341        !! formule 12,  compound waves (78 x 78 x 78 x 0)  
    342  
    343     case ( 12 ) 
    344        zf1=nodal_factort(78) 
    345        zf=nodal_factort(0) 
    346        zf=zf*zf1*zf1*zf1 
    347  
    348        !! formule 13, compound waves (78 x 75) 
    349  
    350     case ( 13 ) 
    351        zf1=nodal_factort(78) 
    352        zf=nodal_factort(75) 
    353        zf=zf*zf1 
    354  
    355        !! formule 14, compound waves (235 x 0)  ===  (235) 
    356  
    357     case ( 14 ) 
    358        zf=nodal_factort(235) 
    359        zf1=nodal_factort(0) 
    360        zf=zf*zf1 
    361  
    362        !! formule 15, compound waves (235 x 75)  
    363  
    364     case ( 15 ) 
    365        zf=nodal_factort(235) 
    366        zf1=nodal_factort(75) 
    367        zf=zf*zf1 
    368  
    369        !! formule 16, compound waves (78 x 0 x 0)  ===  (78) 
    370  
    371     case ( 16 ) 
    372        zf=nodal_factort(78) 
    373        zf1=nodal_factort(0) 
    374        zf=zf*zf1*zf1 
    375  
    376        !! formule 17,  compound waves (227 x 0)  
    377  
    378     case ( 17 ) 
    379        zf1=nodal_factort(227) 
    380        zf=nodal_factort(0) 
    381        zf=zf*zf1 
    382  
    383        !! formule 18,  compound waves (78 x 78 x 78 ) 
    384  
    385     case ( 18 )  
    386        zf1=nodal_factort(78) 
    387        zf=zf1*zf1*zf1 
    388  
    389        !! formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78) 
    390  
    391     case ( 19 ) 
    392        zf=nodal_factort(78) 
    393        zf1=nodal_factort(0) 
    394        zf=zf*zf1*zf1 
    395  
    396        !! formule 73 
    397  
    398     case ( 73 ) 
    399        zs=sin(sh_I) 
    400        zf=(2./3.-zs*zs)/0.5021 
    401  
    402        !! formule 74 
    403  
    404     case ( 74 ) 
    405        zs=sin(sh_I) 
    406        zf=zs*zs/0.1578 
    407  
    408        !! formule 75 
    409  
    410     case ( 75 ) 
    411        zs=cos (sh_I/2) 
    412        zf=sin (sh_I)*zs*zs/0.3800 
    413  
    414        !! formule 76 
    415  
    416     case ( 76 ) 
    417        zf=sin (2*sh_I)/0.7214 
    418  
    419        !! formule 77 
    420  
    421     case ( 77 ) 
    422        zs=sin (sh_I/2) 
    423        zf=sin (sh_I)*zs*zs/0.0164 
    424  
    425        !! formule 78 
    426  
    427     case ( 78 ) 
    428        zs=cos (sh_I/2) 
    429        zf=zs*zs*zs*zs/0.9154 
    430  
    431        !! formule 79 
    432  
    433     case ( 79 ) 
    434        zs=sin(sh_I) 
    435        zf=zs*zs/0.1565 
    436  
    437        !! formule 144 
    438  
    439     case ( 144 ) 
    440        zs=sin (sh_I/2) 
    441        zf=(1-10*zs*zs+15*zs*zs*zs*zs)*cos(sh_I/2)/0.5873 
    442  
    443        !! formule 149 
    444  
    445     case ( 149 ) 
    446        zs=cos (sh_I/2) 
    447        zf=zs*zs*zs*zs*zs*zs/0.8758 
    448  
    449        !! formule 215 
    450  
    451     case ( 215 ) 
    452        zs=cos (sh_I/2) 
    453        zf=zs*zs*zs*zs/0.9154*sh_x1ra 
    454  
    455        !! formule 227  
    456  
    457     case ( 227 ) 
    458        zs=sin (2*sh_I) 
    459        zf=sqrt (0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006) 
    460  
    461        !! formule 235  
    462  
    463     case ( 235 ) 
    464        zs=sin (sh_I) 
    465        zf=sqrt (19.0444*zs*zs*zs*zs+2.7702*zs*zs*cos (2*sh_nu)+.0981) 
    466  
    467     END SELECT 
    468  
    469   end function nodal_factort 
    470  
    471   function dayjul(kyr,kmonth,kday) 
    472     ! 
    473     !*** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE) 
    474     ! 
    475     INTEGER,INTENT(IN) :: kyr,kmonth,kday 
    476     INTEGER,DIMENSION(12) ::  idayt,idays 
    477     INTEGER :: inc,ji 
    478     REAL(wp) :: dayjul,zyq 
    479  
    480     DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./ 
    481     idays(1)=0. 
    482     idays(2)=31. 
    483     inc=0. 
    484     zyq=MOD((kyr-1900.),4.) 
    485     IF(zyq .eq. 0.) inc=1. 
    486     DO ji=3,12 
    487        idays(ji)=idayt(ji)+inc 
    488     END DO 
    489     dayjul=idays(kmonth)+kday 
    490  
    491   END FUNCTION dayjul 
    492  
     251      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235)  
     252         zf1 = nodal_factort( 78) 
     253         zf  = nodal_factort(235) 
     254         zf  = zf1 * zf 
     255         ! 
     256      CASE ( 5 )                 !==  formule 5,  compound waves (78 *78 x 235) 
     257         zf1 = nodal_factort( 78) 
     258         zf  = nodal_factort(235) 
     259         zf  = zf * zf1 * zf1 
     260         ! 
     261      CASE ( 6 )                 !==  formule 6,  compound waves (78 *78 x 0) 
     262         zf1 = nodal_factort(78) 
     263         zf  = nodal_factort( 0) 
     264         zf  = zf * zf1 * zf1  
     265         ! 
     266      CASE( 7 )                  !==  formule 7, compound waves (75 x 75) 
     267         zf = nodal_factort(75) 
     268         zf = zf * zf 
     269         ! 
     270      CASE( 8 )                  !==  formule 8,  compound waves (78 x 0 x 235) 
     271         zf  = nodal_factort( 78) 
     272         zf1 = nodal_factort(  0) 
     273         zf2 = nodal_factort(235) 
     274         zf  = zf * zf1 * zf2 
     275         ! 
     276      CASE( 9 )                  !==  formule 9,  compound waves (78 x 0 x 227) 
     277         zf  = nodal_factort( 78) 
     278         zf1 = nodal_factort(  0) 
     279         zf2 = nodal_factort(227) 
     280         zf  = zf * zf1 * zf2 
     281         ! 
     282      CASE( 10 )                 !==  formule 10,  compound waves (78 x 227) 
     283         zf  = nodal_factort( 78) 
     284         zf1 = nodal_factort(227) 
     285         zf  = zf * zf1 
     286         ! 
     287      CASE( 11 )                 !==  formule 11,  compound waves (75 x 0) 
     288!!gm bug???? zf 2 fois ! 
     289         zf = nodal_factort(75) 
     290         zf = nodal_factort( 0) 
     291         zf = zf * zf1 
     292         ! 
     293      CASE( 12 )                 !==  formule 12,  compound waves (78 x 78 x 78 x 0)  
     294         zf1 = nodal_factort(78) 
     295         zf  = nodal_factort( 0) 
     296         zf  = zf * zf1 * zf1 * zf1 
     297         ! 
     298      CASE( 13 )                 !==  formule 13, compound waves (78 x 75) 
     299         zf1 = nodal_factort(78) 
     300         zf  = nodal_factort(75) 
     301         zf  = zf * zf1 
     302         ! 
     303      CASE( 14 )                 !==  formule 14, compound waves (235 x 0)  ===  (235) 
     304         zf  = nodal_factort(235) 
     305         zf1 = nodal_factort(  0) 
     306         zf  = zf * zf1 
     307         ! 
     308      CASE( 15 )                 !==  formule 15, compound waves (235 x 75)  
     309         zf  = nodal_factort(235) 
     310         zf1 = nodal_factort( 75) 
     311         zf  = zf * zf1 
     312         ! 
     313      CASE( 16 )                 !==  formule 16, compound waves (78 x 0 x 0)  ===  (78) 
     314         zf  = nodal_factort(78) 
     315         zf1 = nodal_factort( 0) 
     316         zf  = zf * zf1 * zf1 
     317         ! 
     318      CASE( 17 )                 !==  formule 17,  compound waves (227 x 0)  
     319         zf1 = nodal_factort(227) 
     320         zf  = nodal_factort(  0) 
     321         zf  = zf * zf1 
     322         ! 
     323      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 ) 
     324         zf1 = nodal_factort(78) 
     325         zf  = zf1 * zf1 * zf1 
     326         ! 
     327      CASE( 19 )                 !==  formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78) 
     328!!gm bug2 ==>>>   here identical to formule 16,  a third multiplication by zf1 is missing 
     329         zf  = nodal_factort(78) 
     330         zf1 = nodal_factort( 0) 
     331         zf = zf * zf1 * zf1 
     332         ! 
     333      CASE( 73 )                 !==  formule 73 
     334         zs = sin(sh_I) 
     335         zf = (2./3.-zs*zs)/0.5021 
     336         ! 
     337      CASE( 74 )                 !==  formule 74 
     338         zs = sin(sh_I) 
     339         zf = zs * zs / 0.1578 
     340         ! 
     341      CASE( 75 )                 !==  formule 75 
     342         zs = cos(sh_I/2) 
     343         zf = sin(sh_I) * zs * zs / 0.3800 
     344         ! 
     345      CASE( 76 )                 !==  formule 76 
     346         zf = sin(2*sh_I) / 0.7214 
     347         ! 
     348      CASE( 77 )                 !==  formule 77 
     349         zs = sin(sh_I/2) 
     350         zf = sin(sh_I) * zs * zs / 0.0164 
     351         ! 
     352      CASE( 78 )                 !==  formule 78 
     353         zs = cos(sh_I/2) 
     354         zf = zs * zs * zs * zs / 0.9154 
     355         ! 
     356      CASE( 79 )                 !==  formule 79 
     357         zs = sin(sh_I) 
     358         zf = zs * zs / 0.1565 
     359         ! 
     360      CASE( 144 )                !==  formule 144 
     361         zs = sin(sh_I/2) 
     362         zf = ( 1-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873 
     363         ! 
     364      CASE( 149 )                !==  formule 149 
     365         zs = cos(sh_I/2) 
     366         zf = zs*zs*zs*zs*zs*zs / 0.8758 
     367         ! 
     368      CASE( 215 )                !==  formule 215 
     369         zs = cos(sh_I/2) 
     370         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra 
     371         ! 
     372      CASE( 227 )                !==  formule 227  
     373         zs = sin(2*sh_I) 
     374         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 ) 
     375         ! 
     376      CASE ( 235 )               !==  formule 235  
     377         zs = sin(sh_I) 
     378         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 ) 
     379         ! 
     380      END SELECT 
     381      ! 
     382   END FUNCTION nodal_factort 
     383 
     384 
     385   FUNCTION dayjul( kyr, kmonth, kday ) 
     386      !!---------------------------------------------------------------------- 
     387      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE) 
     388      !!---------------------------------------------------------------------- 
     389      INTEGER,INTENT(in) ::   kyr, kmonth, kday 
     390      ! 
     391      INTEGER,DIMENSION(12) ::  idayt, idays 
     392      INTEGER  ::   inc, ji 
     393      REAL(wp) ::   dayjul, zyq 
     394      ! 
     395      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./ 
     396      !!---------------------------------------------------------------------- 
     397      ! 
     398      idays(1) = 0. 
     399      idays(2) = 31. 
     400      inc = 0. 
     401      zyq = MOD( kyr-1900. , 4. ) 
     402      IF( zyq == 0.)   inc = 1. 
     403      DO ji = 3, 12 
     404         idays(ji)=idayt(ji)+inc 
     405      END DO 
     406      dayjul = idays(kmonth) + kday 
     407      ! 
     408   END FUNCTION dayjul 
     409 
     410   !!====================================================================== 
    493411END MODULE tide_mod 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90

    r3651 r3953  
    11MODULE tideini 
    2   !!================================================================================= 
    3   !!                       ***  MODULE  tideini  *** 
    4   !! Initialization of tidal forcing 
    5   !! History :  9.0  !  07  (O. Le Galloudec)  Original code 
    6   !!================================================================================= 
    7   !! * Modules used 
    8   USE oce             ! ocean dynamics and tracers variables 
    9   USE dom_oce         ! ocean space and time domain 
    10   USE in_out_manager  ! I/O units 
    11   USE ioipsl          ! NetCDF IPSL library 
    12   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    13   USE phycst 
    14   USE daymod 
    15   USE dynspg_oce 
    16   USE tide_mod 
    17   USE iom 
     2   !!====================================================================== 
     3   !!                       ***  MODULE  tideini  *** 
     4   !! Initialization of tidal forcing 
     5   !!====================================================================== 
     6   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code 
     7   !!---------------------------------------------------------------------- 
     8   USE oce             ! ocean dynamics and tracers variables 
     9   USE dom_oce         ! ocean space and time domain 
     10   USE phycst 
     11   USE daymod 
     12   USE dynspg_oce 
     13   USE tide_mod 
     14   ! 
     15   USE iom 
     16   USE in_out_manager  ! I/O units 
     17   USE ioipsl          ! NetCDF IPSL library 
     18   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1819 
    19   IMPLICIT NONE 
    20   PUBLIC 
     20   IMPLICIT NONE 
     21   PUBLIC 
    2122 
    22   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::  & 
    23        omega_tide,  & 
    24        v0tide,      & 
    25        utide,       & 
    26        ftide 
     23   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   omega_tide   !: 
     24   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   v0tide       !: 
     25   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   utide        !: 
     26   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ftide        !: 
    2727 
    28   LOGICAL, PUBLIC :: ln_tide_pot = .false., ln_tide_ramp = .false. 
    29   REAL(wp), PUBLIC :: rdttideramp  
    30   INTEGER, PUBLIC :: nb_harmo 
    31   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide 
    32   INTEGER, PUBLIC :: kt_tide 
     28   LOGICAL , PUBLIC ::   ln_tide_pot  = .FALSE.   !: 
     29   LOGICAL , PUBLIC ::   ln_tide_ramp = .FALSE.   !: 
     30   INTEGER , PUBLIC ::   nb_harmo                 !: 
     31   INTEGER , PUBLIC ::   kt_tide                  !: 
     32   REAL(wp), PUBLIC ::   rdttideramp              !: 
     33    
     34   INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ntide   !: 
    3335 
    34   !!--------------------------------------------------------------------------------- 
    35   !!   OPA 9.0 , LODYC-IPSL  (2003) 
    36   !!--------------------------------------------------------------------------------- 
    37  
     36   !!---------------------------------------------------------------------- 
     37   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
     38   !! $Id: $ 
     39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     40   !!---------------------------------------------------------------------- 
    3841CONTAINS 
    3942    
    40   SUBROUTINE tide_init ( kt ) 
    41     !!---------------------------------------------------------------------- 
    42     !!                 ***  ROUTINE tide_init  *** 
    43     !!----------------------------------------------------------------------       
    44     !! * Local declarations 
    45     INTEGER  :: ji, jk 
    46     INTEGER, INTENT( in ) ::   kt     ! ocean time-step 
    47     CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 
    48     ! 
    49     NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname 
    50     !!---------------------------------------------------------------------- 
     43   SUBROUTINE tide_init ( kt ) 
     44      !!---------------------------------------------------------------------- 
     45      !!                 ***  ROUTINE tide_init  *** 
     46      !!----------------------------------------------------------------------       
     47      INTEGER  :: ji, jk 
     48      INTEGER, INTENT( in ) ::   kt     ! ocean time-step 
     49      CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 
     50      ! 
     51      NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname 
     52      !!---------------------------------------------------------------------- 
    5153 
    52     IF ( kt == nit000 ) THEN 
    53        ! 
    54        IF(lwp) THEN 
    55           WRITE(numout,*) 
    56           WRITE(numout,*) 'tide_init : Initialization of the tidal components' 
    57           WRITE(numout,*) '~~~~~~~~~ ' 
    58        ENDIF 
    59        ! 
    60        CALL tide_init_Wave 
    61        ! 
    62        clname(:)='' 
    63        ! 
    64        ! Read Namelist nam_tide 
    65        REWIND ( numnam ) 
    66        READ   ( numnam, nam_tide ) 
    67        ! 
    68        nb_harmo=0 
    69        DO jk=1,jpmax_harmo 
    70           DO ji=1,jpmax_harmo 
    71              IF(TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN 
    72                 nb_harmo=nb_harmo+1 
    73              ENDIF 
    74           END DO 
    75        ENDDO 
    76        ! 
    77        IF(lwp) THEN 
    78           WRITE(numout,*) '        Namelist nam_tide' 
    79           WRITE(numout,*) '        nb_harmo    = ', nb_harmo 
    80           WRITE(numout,*) '        ln_tide_ramp = ', ln_tide_ramp  
    81           WRITE(numout,*) '        rdttideramp = ', rdttideramp 
    82           IF (ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp)) & 
    83           & CALL ctl_stop('rdttideramp must be lower than run duration') 
    84           IF (ln_tide_ramp.AND.(rdttideramp<0.)) & 
    85           & CALL ctl_stop('rdttideramp must be positive') 
    86           CALL flush(numout) 
    87        ENDIF 
    88        ! 
    89        ALLOCATE(ntide(nb_harmo)) 
    90        DO jk=1,nb_harmo 
    91           DO ji=1,jpmax_harmo 
    92              IF (TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN 
    93                 ntide(jk) = ji 
    94                 EXIT 
    95              END IF 
    96           END DO 
    97        END DO 
    98        ! 
    99        ALLOCATE(omega_tide(nb_harmo)) 
    100        ALLOCATE(v0tide    (nb_harmo)) 
    101        ALLOCATE(utide     (nb_harmo)) 
    102        ALLOCATE(ftide     (nb_harmo)) 
    103        kt_tide = kt 
    104        ! 
    105     ENDIF 
    106  
    107     IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN 
    108        ! 
    109        IF(lwp) THEN 
    110           WRITE(numout,*) 
    111           WRITE(numout,*) 'tide_ini : Update of the tidal components at kt=',kt 
    112           WRITE(numout,*) '~~~~~~~~ ' 
    113        ENDIF 
    114        CALL tide_harmo(omega_tide, v0tide, utide, ftide, ntide, nb_harmo) 
    115        DO jk =1,nb_harmo 
    116          IF(lwp) WRITE(numout,*) Wave(ntide(jk))%cname_tide,utide(jk),ftide(jk),v0tide(jk),omega_tide(jk) 
    117          call flush(numout) 
    118        END DO 
    119        ! 
    120        kt_tide = kt 
    121        ! 
    122     ENDIF 
    123  
    124   END SUBROUTINE tide_init 
    125     
     54      IF( kt == nit000 ) THEN 
     55         ! 
     56         IF(lwp) THEN 
     57            WRITE(numout,*) 
     58            WRITE(numout,*) 'tide_init : Initialization of the tidal components' 
     59            WRITE(numout,*) '~~~~~~~~~ ' 
     60         ENDIF 
     61         ! 
     62         CALL tide_init_Wave 
     63         ! 
     64         clname(:)='' 
     65         ! 
     66         REWIND( numnam )            ! Read Namelist nam_tide 
     67         READ  ( numnam, nam_tide ) 
     68         ! 
     69         nb_harmo=0 
     70         DO jk = 1, jpmax_harmo 
     71            DO ji = 1,jpmax_harmo 
     72               IF( TRIM(clname(jk)) == Wave(ji)%cname_tide )   nb_harmo = nb_harmo + 1 
     73            END DO 
     74         END DO 
     75         ! 
     76         IF(lwp) THEN 
     77            WRITE(numout,*) '   Namelist nam_tide' 
     78            WRITE(numout,*) '      Apply astronomical potential : ln_tide_pot  =', ln_tide_pot 
     79            WRITE(numout,*) '                                     nb_harmo     = ', nb_harmo 
     80            WRITE(numout,*) '                                     ln_tide_ramp = ', ln_tide_ramp  
     81            WRITE(numout,*) '                                     rdttideramp  = ', rdttideramp 
     82         ENDIF 
     83         IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) )   & 
     84            &   CALL ctl_stop('rdttideramp must be lower than run duration') 
     85         IF( ln_tide_ramp.AND.(rdttideramp<0.) ) & 
     86            &   CALL ctl_stop('rdttideramp must be positive') 
     87         ! 
     88         IF( .NOT. lk_dynspg_ts )   CALL ctl_warn( 'sbc_tide : use of time splitting is recommended' ) 
     89         ! 
     90         ALLOCATE( ntide(nb_harmo) ) 
     91         DO jk = 1, nb_harmo 
     92            DO ji = 1, jpmax_harmo 
     93               IF( TRIM(clname(jk)) .eq. Wave(ji)%cname_tide ) THEN 
     94                  ntide(jk) = ji 
     95                  EXIT 
     96               END IF 
     97            END DO 
     98         END DO 
     99         ! 
     100         ALLOCATE( omega_tide(nb_harmo), v0tide    (nb_harmo),   & 
     101            &      utide     (nb_harmo), ftide     (nb_harmo)  ) 
     102         kt_tide = kt 
     103         ! 
     104      ENDIF 
     105      ! 
     106      IF( nsec_day == NINT(0.5 * rdttra(1)) ) THEN 
     107         ! 
     108         CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo ) 
     109         ! 
     110         kt_tide = kt 
     111         ! 
     112         IF(lwp) THEN 
     113            WRITE(numout,*) 
     114            WRITE(numout,*) 'tide_ini : Update of the tidal components at kt=', kt 
     115            WRITE(numout,*) '~~~~~~~~ ' 
     116            DO jk = 1, nb_harmo 
     117               WRITE(numout,*) Wave(ntide(jk))%cname_tide, utide(jk), ftide(jk), v0tide(jk), omega_tide(jk) 
     118            END DO 
     119         ENDIF 
     120         ! 
     121      ENDIF 
     122      ! 
     123   END SUBROUTINE tide_init 
     124      
     125   !!====================================================================== 
    126126END MODULE tideini 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90

    r3651 r3953  
    11MODULE updtide 
    2   !!================================================================================= 
    3   !!                       ***  MODULE  updtide  *** 
    4   !! Initialization of tidal forcing 
    5   !! History :  9.0  !  07  (O. Le Galloudec)  Original code 
    6   !!================================================================================= 
     2   !!====================================================================== 
     3   !!                       ***  MODULE  updtide  *** 
     4   !! Initialization of tidal forcing 
     5   !!====================================================================== 
     6   !! History :  9.0  !  07  (O. Le Galloudec)  Original code 
     7   !!---------------------------------------------------------------------- 
    78#if defined key_tide 
    8   !! * Modules used 
    9   USE oce             ! ocean dynamics and tracers variables 
    10   USE dom_oce         ! ocean space and time domain 
    11   USE in_out_manager  ! I/O units 
    12   USE phycst 
    13   USE sbctide 
    14   USE dynspg_oce 
    15   USE tideini, ONLY: ln_tide_ramp, rdttideramp 
     9   !!---------------------------------------------------------------------- 
     10   !!   'key_tide' :                                        tidal potential 
     11   !!---------------------------------------------------------------------- 
     12   !!   upd_tide       : update tidal potential 
     13   !!---------------------------------------------------------------------- 
     14   USE oce             ! ocean dynamics and tracers variables 
     15   USE dom_oce         ! ocean space and time domain 
     16   USE in_out_manager  ! I/O units 
     17   USE phycst          ! physical constant 
     18   USE sbctide         ! tide potential variable 
     19   USE tideini, ONLY: ln_tide_ramp, rdttideramp 
    1620 
    17   IMPLICIT NONE 
    18   PUBLIC 
     21   IMPLICIT NONE 
     22   PUBLIC 
    1923 
    20   !! * Routine accessibility 
    21   PUBLIC upd_tide 
    22   !!--------------------------------------------------------------------------------- 
    23   !!   OPA 9.0 , LODYC-IPSL  (2003) 
    24   !!--------------------------------------------------------------------------------- 
    25  
     24   PUBLIC   upd_tide   ! called in dynspg_... modules 
     25   
     26   !!---------------------------------------------------------------------- 
     27   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     28   !! $Id: sbcfwb.F90 3625 2012-11-21 13:19:18Z acc $ 
     29   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     30   !!---------------------------------------------------------------------- 
    2631CONTAINS 
    2732 
    28   SUBROUTINE upd_tide (kt,kit) 
    29     !!---------------------------------------------------------------------- 
    30     !!                 ***  ROUTINE upd_tide  *** 
    31     !!----------------------------------------------------------------------       
    32     !! * Local declarations 
     33   SUBROUTINE upd_tide( kt, kit, kbaro ) 
     34      !!---------------------------------------------------------------------- 
     35      !!                 ***  ROUTINE upd_tide  *** 
     36      !! 
     37      !! ** Purpose :   provide at each time step the astronomical potential 
     38      !! 
     39      !! ** Method  :   computed from pulsation and amplitude of all tide components 
     40      !! 
     41      !! ** Action  :   pot_astro   actronomical potential 
     42      !!----------------------------------------------------------------------       
     43      INTEGER, INTENT(in)           ::   kt      ! ocean time-step index 
     44      INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T only) 
     45      INTEGER, INTENT(in), OPTIONAL ::   kbaro   ! number of sub-time-step           (lk_dynspg_ts=T only) 
     46      ! 
     47      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     48      REAL(wp) ::   zt, zramp    ! local scalar 
     49      REAL(wp), DIMENSION(nb_harmo) ::   zwt  
     50      !!----------------------------------------------------------------------       
     51      ! 
     52      !                               ! tide pulsation at model time step (or sub-time-step) 
     53      zt = ( kt - kt_tide ) * rdt 
     54      IF( PRESENT( kit ) .AND. PRESENT( kbaro ) )   zt = zt + kit * rdt / REAL( kbaro, wp ) 
     55      zwt(:) = omega_tide(:) * zt 
    3356 
    34     INTEGER, INTENT( in ) ::   kt,kit      ! ocean time-step index 
    35     INTEGER  :: ji,jj,jk 
    36     REAL (wp) :: zramp 
    37     REAL (wp), DIMENSION(nb_harmo) :: zwt  
    38     !............................................................................... 
    39  
    40     pot_astro(:,:)=0.e0 
    41     zramp = 1.e0 
    42  
    43     IF (lk_dynspg_ts) THEN 
    44        zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp))) 
    45        IF (ln_tide_ramp) THEN 
    46           zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.) 
    47        ENDIF 
    48     ELSE 
    49        zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt 
    50        IF (ln_tide_ramp) THEN 
    51           zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.)  
    52        ENDIF   
    53     ENDIF 
    54  
    55     do jk=1,nb_harmo 
    56        do ji=1,jpi 
    57           do jj=1,jpj 
    58              pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))       
    59           enddo 
    60        enddo 
    61     enddo 
    62  
    63   END SUBROUTINE upd_tide 
     57      pot_astro(:,:) = 0._wp          ! update tidal potential (sum of all harmonics) 
     58      DO jk = 1, nb_harmo    
     59         pot_astro(:,:) = pot_astro(:,:) + amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )       
     60      END DO 
     61      ! 
     62      IF( ln_tide_ramp ) THEN         ! linear increase if asked 
     63         zt = ( kt - nit000 ) * rdt 
     64         IF( PRESENT( kit ) .AND. PRESENT( kbaro ) )   zt = zt + kit * rdt / REAL( kbaro, wp ) 
     65         zramp = MIN(  MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp  ) 
     66         pot_astro(:,:) = zramp * pot_astro(:,:) 
     67      ENDIF 
     68      ! 
     69   END SUBROUTINE upd_tide 
    6470 
    6571#else 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90

    r3865 r3953  
    8989      ! Update data, open boundaries, surface boundary condition (including sea-ice) 
    9090      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    91                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
    92       IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp ) 
     91                         CALL sbc     ( kstp )        ! Sea Boundary Condition (including sea-ice) 
    9392      IF( lk_tide    )   CALL sbc_tide( kstp ) 
    94       IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries 
    95       IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries 
    96       IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries 
     93      IF( lk_obc     )   CALL obc_dta ( kstp )        ! update dynamic and tracer data at open boundaries 
     94      IF( lk_obc     )   CALL obc_rad ( kstp )        ! compute phase velocities at open boundaries 
     95      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    9796 
    9897      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    117116      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz 
    118117      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz 
    119       IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value) 
     118      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value) 
    120119         avt (:,:,:) = rn_avt0 * tmask(:,:,:) 
    121120         avmu(:,:,:) = rn_avm0 * umask(:,:,:) 
Note: See TracChangeset for help on using the changeset viewer.