Changeset 2814


Ignore:
Timestamp:
2011-07-27T14:41:28+02:00 (10 years ago)
Author:
davestorkey
Message:
  1. Implement tidal harmonics forcing (UKMO version) in new structure.
  2. Other bug fixes and updates.
Location:
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r2800 r2814  
    2626   USE zdfbfr          ! bottom friction 
    2727   USE dynvor          ! vorticity term 
     28   USE obc_par         ! for lk_obc 
    2829   USE obc_oce         ! Lateral open boundary condition 
    2930   USE obcdta          ! open boundary condition data      
     
    528529      !                                                    ! ==================== ! 
    529530 
    530 #if defined key_obc 
    531       IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ????? 
    532       IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:) 
    533       IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:) 
    534       IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:) 
    535 #endif 
    536  
    537531      ! ----------------------------------------------------------------------------- 
    538532      ! Phase 3. update the general trend with the barotropic trend 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/FLO/floblk.F90

    r2715 r2814  
    345345      ! more time.       
    346346# if defined key_obc 
    347       DO jfl = 1, jpnfl 
    348          IF( lp_obc_east ) THEN 
    349             IF( jped <=  zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <=  zgifl(jfl) ) THEN 
    350                zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
    351                zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
    352                zagefl(jfl) = rdt 
    353             END IF 
    354          END IF 
    355          IF( lp_obc_west ) THEN 
    356             IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >=  zgifl(jfl) ) THEN 
    357                zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
    358                zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
    359                zagefl(jfl) = rdt 
    360             END IF 
    361          END IF 
    362          IF( lp_obc_north ) THEN 
    363             IF( jpnd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >=  zgjfl(jfl) ) THEN 
    364                zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
    365                zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
    366                zagefl(jfl) = rdt 
    367             END IF 
    368          END IF 
    369          IF( lp_obc_south ) THEN 
    370             IF( jpsd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND.  njsob >= zgjfl(jfl) ) THEN 
    371                zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
    372                zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
    373                zagefl(jfl) = rdt 
    374             END IF 
    375          END IF 
    376       END DO 
     347!!!!!!!! NEED TO SORT THIS OUT !!!!!!!! 
     348!!$      DO jfl = 1, jpnfl 
     349!!$         IF( lp_obc_east ) THEN 
     350!!$            IF( jped <=  zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <=  zgifl(jfl) ) THEN 
     351!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
     352!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
     353!!$               zagefl(jfl) = rdt 
     354!!$            END IF 
     355!!$         END IF 
     356!!$         IF( lp_obc_west ) THEN 
     357!!$            IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >=  zgifl(jfl) ) THEN 
     358!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
     359!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
     360!!$               zagefl(jfl) = rdt 
     361!!$            END IF 
     362!!$         END IF 
     363!!$         IF( lp_obc_north ) THEN 
     364!!$            IF( jpnd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >=  zgjfl(jfl) ) THEN 
     365!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
     366!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
     367!!$               zagefl(jfl) = rdt 
     368!!$            END IF 
     369!!$         END IF 
     370!!$         IF( lp_obc_south ) THEN 
     371!!$            IF( jpsd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND.  njsob >= zgjfl(jfl) ) THEN 
     372!!$               zgifl (jfl) = INT(zgifl(jfl)) + 0.5 
     373!!$               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5 
     374!!$               zagefl(jfl) = rdt 
     375!!$            END IF 
     376!!$         END IF 
     377!!$      END DO 
    377378#endif 
    378379 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90

    r2800 r2814  
    5656   !                                                        !: =F read obc coordinates from namelist 
    5757   LOGICAL                    ::   ln_mask_file             !: =T read obcmask from file 
    58    LOGICAL, DIMENSION(jp_obc) ::   ln_tides                 !: =T apply tidal harmonic forcing along open boundaries 
    5958   LOGICAL                    ::   ln_vol                   !: =T volume correction              
    6059   LOGICAL, DIMENSION(jp_obc) ::   ln_clim                  !: =T obc data files contain climatological data (time-cyclic) 
     
    6463   INTEGER, DIMENSION(jp_obc) ::   nn_dtactl           !: = 0 use the initial state as obc dta ;  
    6564                                                            !: = 1 read it in a NetCDF file 
     65   INTEGER, DIMENSION(jp_obc) ::   nn_tides                 !: = 0 no tidal harmonic forcing 
     66                                                            !: = 1 apply ONLY tidal harmonic forcing for barotropic solution 
     67                                                            !: = 2 ADD tidal harmonic forcing to other barotropic boundary data 
    6668   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
    6769   !                                                        !  = 1 the volume will be constant during all the integration. 
     
    98100   !!---------------------------------------------------------------------- 
    99101 
    100    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)      ::   dta_global        !: workspace for reading in global data arrays 
    101    TYPE(OBC_INDEX), DIMENSION(jp_obc), TARGET    ::   idx_obc           !: obc indices (local process) 
    102    TYPE(OBC_DATA) , DIMENSION(jp_obc)            ::   dta_obc           !: obc external data (local process) 
     102   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays 
     103   TYPE(OBC_INDEX), DIMENSION(jp_obc), TARGET      ::   idx_obc           !: obc indices (local process) 
     104   TYPE(OBC_DATA) , DIMENSION(jp_obc)              ::   dta_obc           !: obc external data (local process) 
    103105 
    104106   !!---------------------------------------------------------------------- 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90

    r2800 r2814  
    6565      !!---------------------------------------------------------------------- 
    6666      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 
    67       USE wrk_nemo, ONLY: wrk_2d_9, wrk_2d_10      ! 2D workspace 
     67      USE wrk_nemo, ONLY: wrk_2d_22, wrk_2d_23   ! 2D workspace 
    6868      !! 
    6969      INTEGER, INTENT( in )           ::   kt    ! ocean time-step index  
     
    7676      !!--------------------------------------------------------------------------- 
    7777 
    78       IF(wrk_in_use(2, 9,10) ) THEN 
     78      IF(wrk_in_use(2, 22,23) ) THEN 
    7979         CALL ctl_stop('obc_dta: ERROR: requested workspace arrays are unavailable.')   ;   RETURN 
    8080      END IF 
     
    8787         ! Calculate depth-mean currents 
    8888         !----------------------------- 
    89          pu2d => wrk_2d_9 
    90          pu2d => wrk_2d_10 
     89         pu2d => wrk_2d_22 
     90         pu2d => wrk_2d_23 
    9191 
    9292         pu2d(:,:) = 0.e0 
     
    206206               ! Update barotropic boundary conditions only 
    207207               ! jit is optional argument for fld_read 
    208                IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN 
     208               IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_tides(ib_obc) .ne. 1 ) THEN 
    209209                  jend = jstart + 2 
    210210                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit ) 
    211211               ENDIF 
     212               IF( nn_tides(ib_obc) .gt. 0 ) THEN 
     213                  CALL tide_update( kt=kt, jit=jit, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) ) 
     214               ENDIF 
    212215            ELSE 
    213                jend = jstart + nb_obc_fld(ib_obc) - 1 
    214                CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend ), map=nbmap_ptr(jstart:jend), timeshift=1 ) 
     216               IF( nb_obc_fld(ib_obc) .gt. 0 ) THEN  
     217                  jend = jstart + nb_obc_fld(ib_obc) - 1 
     218                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), timeshift=1 ) 
     219               ENDIF 
     220               IF( nn_tides(ib_obc) .gt. 0 ) THEN 
     221                  !!! THINK ABOUT kt, jit VALUES !!! 
     222                  CALL tide_update( kt=kt, jit=0, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) ) 
     223               ENDIF 
    215224            ENDIF 
    216225            jstart = jend+1 
     
    257266      END DO  ! ib_obc 
    258267 
    259       IF(wrk_not_released(2, 9,10) )    CALL ctl_stop('obc_dta: ERROR: failed to release workspace arrays.') 
     268      IF(wrk_not_released(2, 22,23) )    CALL ctl_stop('obc_dta: ERROR: failed to release workspace arrays.') 
    260269 
    261270      END SUBROUTINE obc_dta 
     
    383392            IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN  
    384393 
    385                IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN 
     394               IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1) THEN 
    386395                  jfld = jfld + 1 
    387396                  blf_i(jfld) = bn_ssh 
     
    392401               ENDIF 
    393402 
    394                IF( .not. ln_full_vel_array(ib_obc) ) THEN 
     403               IF( .not. ln_full_vel_array(ib_obc) .and. nn_tides(ib_obc) .ne. 1 ) THEN 
    395404 
    396405                  jfld = jfld + 1 
     
    604613 
    605614            ! nn_dtactl = 1 
    606             ! Set boundary data arrays to point to relevant "fnow" arrays 
    607             !----------------------------------------------------------- 
     615            ! Set boundary data arrays to point to "fnow" arrays 
     616            !--------------------------------------------------- 
    608617            IF (nn_dyn2d(ib_obc) .gt. 0) THEN 
    609                IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN 
     618               IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1 ) THEN 
    610619                  jfld = jfld + 1 
    611620                  dta_obc(ib_obc)%ssh => bf(jfld)%fnow(:,1,1) 
    612621               ENDIF 
    613                IF( ln_full_vel_array(ib_obc) ) THEN 
     622               IF( ln_full_vel_array(ib_obc) .or. nn_tides(ib_obc) .eq. 1 ) THEN 
    614623                  ! In this case we need space but we aren't reading it  
    615624                  ! directly from the external file.  
    616625                  IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN 
    617                      ilen0(2) = nblen(2) 
    618                      ilen0(3) = nblen(3) 
     626                     ilen0(:) = nblen(:) 
    619627                  ELSE 
    620                      ilen0(2) = nblenrim(2) 
    621                      ilen0(3) = nblenrim(3) 
     628                     ilen0(:) = nblenrim(:) 
    622629                  ENDIF 
     630                  IF( nn_tides(ib_obc) .eq. 1 ) ALLOCATE( dta_obc(ib_obc)%ssh(ilen0(1)) ) 
    623631                  ALLOCATE( dta_obc(ib_obc)%u2d(ilen0(2)) ) 
    624632                  ALLOCATE( dta_obc(ib_obc)%v2d(ilen0(3)) ) 
     
    631639            ENDIF 
    632640            IF (nn_dyn3d(ib_obc) .gt. 0 .or. & 
    633               &   ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) THEN 
     641              &  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) ) THEN 
    634642               jfld = jfld + 1 
    635643               dta_obc(ib_obc)%u3d => bf(jfld)%fnow(:,1,:) 
     
    666674CONTAINS 
    667675   SUBROUTINE obc_dta( kt, jit )              ! Empty routine 
     676      INTEGER, INTENT( in )           ::   kt     
     677      INTEGER, INTENT( in ), OPTIONAL ::   jit    
    668678      WRITE(*,*) 'obc_dta: You should not have seen this print! error?', kt 
    669679   END SUBROUTINE obc_dta 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90

    r2800 r2814  
    2424   USE phycst          ! physical constants 
    2525   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    26    USE obctides        ! for tidal harmonic forcing at boundary 
    2726   USE in_out_manager  ! 
    2827 
     
    3029   PRIVATE 
    3130 
    32    PUBLIC   obc_dyn2d     ! routine called in dynspg_ts (time splitting case ONLY) 
     31   PUBLIC   obc_dyn2d     ! routine called in dynspg_ts and dyn_nxt and dynspg_flt 
    3332 
    3433   !!---------------------------------------------------------------------- 
     
    7776      !!               topography. Tellus, 365-382. 
    7877      !!---------------------------------------------------------------------- 
    79       INTEGER,         INTENT(in) :: kt 
     78      INTEGER,         INTENT(in) ::   kt 
    8079      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    8180      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     
    114113      !!              
    115114      !!              - Apply Flather boundary conditions on normal barotropic velocities  
    116       !!                (ln_dyn_fla=.true. or ln_tides=.true.) 
    117115      !! 
    118116      !! ** WARNINGS about FLATHER implementation: 
     
    143141      ! ---------------------------------!  
    144142      
    145 !!!!!!!!!!!! SOME WORK TO DO ON THE TIDES HERE !!!!!!!!!!!!! 
    146  
    147143!!! REPLACE spgu with nemo_wrk work space 
    148144 
     
    154150         ij = idx%nbj(jb,igrd) 
    155151         spgu(ii, ij) = dta%ssh(jb) 
    156 !!$         IF( ln_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(jb) 
    157152      END DO 
    158153      ! 
     
    167162         ! 
    168163         zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
    169 !!$         zforc = dta%u2d(jb) + utide(jb) 
    170164         zforc = dta%u2d(jb) 
    171165         pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1)  
     
    181175         ! 
    182176         zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
    183 !!$         zforc = dta%v2d(jb) + vtide(jb) 
    184177         zforc = dta%v2d(jb) 
    185178         pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90

    r2797 r2814  
    2222   USE dom_oce         ! ocean space and time domain 
    2323   USE obc_oce         ! unstructured open boundary conditions 
    24    USE obctides        ! tides at open boundaries initialization (tide_init routine) 
    2524   USE in_out_manager  ! I/O units 
    2625   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    7372         &             nn_ice_lim2,                                        & 
    7473#endif 
    75          &             ln_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl,    & 
     74         &             nn_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl,    & 
    7675         &             nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out,             & 
    7776         &             nn_dmp3d_in, nn_dmp3d_out 
     
    105104      nn_ice_lim2(:)    = 0 
    106105#endif 
    107       ln_tides(:)       = .false. 
    108106      ln_vol            = .false. 
    109107      ln_clim(:)        = .false. 
    110108      nn_dtactl(:)      = -1  ! uninitialised flag 
     109      nn_tides(:)       =  0  ! default to no tidal forcing 
    111110      nn_volctl         = -1  ! uninitialised flag 
    112111      nn_rimwidth(:)    = -1  ! uninitialised flag 
     
    182181        IF(lwp) WRITE(numout,*) 
    183182 
    184         IF( ln_tides(ib_obc) ) THEN 
    185           IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 
     183        SELECT CASE( nn_tides(ib_obc) )  
     184        CASE(0) 
     185          IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing' 
    186186          IF(lwp) WRITE(numout,*) 
    187         ENDIF 
    188  
    189 !!$        ! Presumably will need to read in a separate namelist for each boundary that includes tides??? 
    190 !!$        IF( ln_tides )   CALL tide_init( ib_obc )      ! Read tides namelist  
    191  
     187        CASE(1) 
     188          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution' 
     189          IF(lwp) WRITE(numout,*) 
     190        CASE(2) 
     191          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions' 
     192          IF(lwp) WRITE(numout,*) 
     193        CASE DEFAULT 
     194          CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' ) 
     195        END SELECT 
    192196 
    193197      ENDDO 
     
    345349         ENDDO  
    346350 
    347          ! Compute rim weights 
    348          ! ------------------- 
     351         ! Compute rim weights for FRS scheme 
     352         ! ---------------------------------- 
    349353         DO igrd = 1, jpbgrd 
    350354            DO ib = 1, idx_obc(ib_obc)%nblen(igrd) 
     
    537541         IF( lk_mpp )   CALL mpp_sum( obcsurftot )      ! sum over the global domain 
    538542      END IF    
    539  
    540       ! Read in tidal constituents and adjust for model start time 
    541       ! ---------------------------------------------------------- 
    542 !!$      IF( ln_tides )   CALL tide_data 
    543543      ! 
    544544      ! Tidy up 
     
    550550#else 
    551551   !!--------------------------------------------------------------------------------- 
    552    !!   Dummy module                                   NO unstructured open boundaries 
     552   !!   Dummy module                                   NO open boundaries 
    553553   !!--------------------------------------------------------------------------------- 
    554554CONTAINS 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90

    r2797 r2814  
    1010   !!---------------------------------------------------------------------- 
    1111#if defined key_obc 
    12 !!$   !!---------------------------------------------------------------------- 
    13 !!$   !!   'key_obc'     Unstructured Open Boundary Condition 
    14 !!$   !!---------------------------------------------------------------------- 
    15 !!$   !!   PUBLIC 
    16 !!$   !!      tide_init     : read of namelist 
    17 !!$   !!      tide_data     : read in and initialisation of tidal constituents at boundary 
    18 !!$   !!      tide_update   : calculation of tidal forcing at each timestep 
    19 !!$   !!   PRIVATE 
    20 !!$   !!      uvset         :\ 
    21 !!$   !!      vday          : | Routines to correct tidal harmonics forcing for 
    22 !!$   !!      shpen         : | start time of integration 
    23 !!$   !!      ufset         : | 
    24 !!$   !!      vset          :/ 
    25 !!$   !!---------------------------------------------------------------------- 
    26 !!$   USE oce             ! ocean dynamics and tracers  
    27 !!$   USE dom_oce         ! ocean space and time domain 
    28 !!$   USE iom 
    29 !!$   USE in_out_manager  ! I/O units 
    30 !!$   USE phycst          ! physical constants 
    31 !!$   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    32 !!$   USE obc_par         ! Unstructured boundary parameters 
    33 !!$   USE obc_oce         ! ocean open boundary conditions 
    34 !!$   USE daymod          ! calendar 
    35 !!$ 
    36 !!$   IMPLICIT NONE 
    37 !!$   PRIVATE 
    38 !!$ 
    39 !!$   PUBLIC   tide_init     ! routine called in obcini 
    40 !!$   PUBLIC   tide_data     ! routine called in obcini 
    41 !!$   PUBLIC   tide_update   ! routine called in obcdyn 
    42 !!$ 
    43 !!$   LOGICAL, PUBLIC            ::   ln_tide_date          !: =T correct tide phases and amplitude for model start date 
    44 !!$   INTEGER, PUBLIC, PARAMETER ::   jptides_max = 15      !: Max number of tidal contituents 
    45 !!$   INTEGER, PUBLIC            ::   ntide                 !: Actual number of tidal constituents 
    46 !!$ 
    47 !!$   CHARACTER(len=80), PUBLIC                         ::   filtide    !: Filename root for tidal input files 
    48 !!$   CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) ::   tide_cpt   !: Names of tidal components used. 
    49 !!$ 
    50 !!$   INTEGER , PUBLIC, DIMENSION(jptides_max) ::   nindx        !: ??? 
    51 !!$   REAL(wp), PUBLIC, DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr) 
    52 !!$    
    53 !!$   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   ! Tidal constituents : SSH 
    54 !!$   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     ! Tidal constituents : U 
    55 !!$   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     ! Tidal constituents : V 
    56 !!$    
    57 !!$   !!---------------------------------------------------------------------- 
    58 !!$   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    59 !!$   !! $Id: obctides.F90 2528 2010-12-27 17:33:53Z rblod $  
    60 !!$   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    61 !!$   !!---------------------------------------------------------------------- 
    62 !!$CONTAINS 
    63 !!$ 
    64 !!$   SUBROUTINE tide_init 
    65 !!$      !!---------------------------------------------------------------------- 
    66 !!$      !!                    ***  SUBROUTINE tide_init  *** 
    67 !!$      !!                      
    68 !!$      !! ** Purpose : - Read in namelist for tides 
    69 !!$      !! 
    70 !!$      !!---------------------------------------------------------------------- 
    71 !!$      INTEGER ::   itide                  ! dummy loop index  
    72 !!$      !! 
    73 !!$      NAMELIST/namobc_tide/ln_tide_date, filtide, tide_cpt, tide_speed 
    74 !!$      !!---------------------------------------------------------------------- 
    75 !!$ 
    76 !!$      IF(lwp) WRITE(numout,*) 
    77 !!$      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' 
    78 !!$      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
    79 !!$ 
    80 !!$      ! Namelist namobc_tide : tidal harmonic forcing at open boundaries 
    81 !!$      ln_tide_date = .false. 
    82 !!$      filtide(:) = '' 
    83 !!$      tide_speed(:) = 0.0 
    84 !!$      tide_cpt(:) = '' 
    85 !!$      REWIND( numnam )                                ! Read namelist parameters 
    86 !!$      READ  ( numnam, namobc_tide ) 
    87 !!$      !                                               ! Count number of components specified 
    88 !!$      ntide = jptides_max 
    89 !!$      DO itide = 1, jptides_max 
    90 !!$        IF( tide_cpt(itide) == '' ) THEN 
    91 !!$           ntide = itide-1 
    92 !!$           exit 
    93 !!$        ENDIF 
    94 !!$      END DO 
    95 !!$ 
    96 !!$      !                                               ! find constituents in standard list 
    97 !!$      DO itide = 1, ntide 
    98 !!$         nindx(itide) = 0 
    99 !!$         IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1 
    100 !!$         IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2 
    101 !!$         IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3 
    102 !!$         IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4 
    103 !!$         IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5 
    104 !!$         IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6 
    105 !!$         IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7 
    106 !!$         IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8 
    107 !!$         IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9 
    108 !!$         IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10 
    109 !!$         IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11 
    110 !!$         IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12 
    111 !!$         IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13 
    112 !!$         IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14 
    113 !!$         IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15 
    114 !!$         IF( nindx(itide) == 0  .AND. lwp ) THEN 
    115 !!$            WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 
    116 !!$            CALL ctl_warn( ctmp1 ) 
    117 !!$         ENDIF 
    118 !!$      END DO 
    119 !!$      !                                               ! Parameter control and print 
    120 !!$      IF( ntide < 1 ) THEN  
    121 !!$         CALL ctl_stop( '          Did not find any tidal components in namelist namobc_tide' ) 
    122 !!$      ELSE 
    123 !!$         IF(lwp) WRITE(numout,*) '          Namelist namobc_tide : tidal harmonic forcing at open boundaries' 
    124 !!$         IF(lwp) WRITE(numout,*) '             tidal components specified ', ntide 
    125 !!$         IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:ntide) 
    126 !!$         IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
    127 !!$         IF(lwp) WRITE(numout,*) '                ', tide_speed(1:ntide) 
    128 !!$      ENDIF 
    129 !!$ 
    130 !!$      ! Initialisation of tidal harmonics arrays 
    131 !!$      sshtide(:) = 0.e0 
    132 !!$      utide  (:) = 0.e0 
    133 !!$      vtide  (:) = 0.e0 
    134 !!$      ! 
    135 !!$   END SUBROUTINE tide_init 
    136 !!$ 
    137 !!$ 
    138 !!$   SUBROUTINE tide_data 
    139 !!$      !!---------------------------------------------------------------------- 
    140 !!$      !!                    ***  SUBROUTINE tide_data  *** 
    141 !!$      !!                     
    142 !!$      !! ** Purpose : - Read in tidal harmonics data and adjust for the start  
    143 !!$      !!                time of the model run.  
    144 !!$      !! 
    145 !!$      !!---------------------------------------------------------------------- 
    146 !!$      INTEGER ::   itide, igrd, ib        ! dummy loop indices 
    147 !!$      CHARACTER(len=80) :: clfile         ! full file name for tidal input file  
    148 !!$      INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read) 
    149 !!$      INTEGER, DIMENSION(6) :: lendta=0   ! length of data in the file (note may be different from nblendta!) 
    150 !!$      REAL(wp) ::  z_arg, z_atde, z_btde, z1t, z2t            
    151 !!$      REAL(wp), DIMENSION(jpbdta,1) ::   zdta   ! temporary array for data fields 
    152 !!$      REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc 
    153 !!$      !!------------------------------------------------------------------------------ 
    154 !!$ 
    155 !!$      ! Open files and read in tidal forcing data 
    156 !!$      ! ----------------------------------------- 
    157 !!$ 
    158 !!$      ipj = 1 
    159 !!$ 
    160 !!$      DO itide = 1, ntide 
    161 !!$         !                                                              ! SSH fields 
    162 !!$         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
    163 !!$         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
    164 !!$         CALL iom_open( clfile, inum ) 
    165 !!$         igrd = 4 
    166 !!$         IF( nblendta(igrd) <= 0 ) THEN  
    167 !!$            idvar = iom_varid( inum,'z1' ) 
    168 !!$            IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 
    169 !!$            nblendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
    170 !!$            WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd) 
    171 !!$         ENDIF 
    172 !!$         ipi = nblendta(igrd) 
    173 !!$         CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) ) 
    174 !!$         DO ib = 1, nblenrim(igrd) 
    175 !!$            ssh1(ib,itide) = zdta(nbmap(ib,igrd),1) 
    176 !!$         END DO 
    177 !!$         CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) ) 
    178 !!$         DO ib = 1, nblenrim(igrd) 
    179 !!$            ssh2(ib,itide) = zdta(nbmap(ib,igrd),1) 
    180 !!$         END DO 
    181 !!$         CALL iom_close( inum ) 
    182 !!$         ! 
    183 !!$         !                                                              ! U fields 
    184 !!$         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
    185 !!$         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
    186 !!$         CALL iom_open( clfile, inum ) 
    187 !!$         igrd = 5 
    188 !!$         IF( lendta(igrd) <= 0 ) THEN  
    189 !!$            idvar = iom_varid( inum,'u1' ) 
    190 !!$            lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
    191 !!$            WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd) 
    192 !!$         ENDIF 
    193 !!$         ipi = lendta(igrd) 
    194 !!$         CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) ) 
    195 !!$         DO ib = 1, nblenrim(igrd) 
    196 !!$            u1(ib,itide) = zdta(nbmap(ib,igrd),1) 
    197 !!$         END DO 
    198 !!$         CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) ) 
    199 !!$         DO ib = 1, nblenrim(igrd) 
    200 !!$            u2(ib,itide) = zdta(nbmap(ib,igrd),1) 
    201 !!$         END DO 
    202 !!$         CALL iom_close( inum ) 
    203 !!$         ! 
    204 !!$         !                                                              ! V fields 
    205 !!$         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
    206 !!$         if(lwp) write(numout,*) 'Reading data from file ', clfile 
    207 !!$         CALL iom_open( clfile, inum ) 
    208 !!$         igrd = 6 
    209 !!$         IF( lendta(igrd) <= 0 ) THEN  
    210 !!$            idvar = iom_varid( inum,'v1' ) 
    211 !!$            lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 
    212 !!$            WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd) 
    213 !!$         ENDIF 
    214 !!$         ipi = lendta(igrd) 
    215 !!$         CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) ) 
    216 !!$         DO ib = 1, nblenrim(igrd) 
    217 !!$            v1(ib,itide) = zdta(nbmap(ib,igrd),1) 
    218 !!$         END DO 
    219 !!$         CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) ) 
    220 !!$         DO ib=1, nblenrim(igrd) 
    221 !!$            v2(ib,itide) = zdta(nbmap(ib,igrd),1) 
    222 !!$         END DO 
    223 !!$         CALL iom_close( inum ) 
    224 !!$         ! 
    225 !!$      END DO ! end loop on tidal components 
    226 !!$ 
    227 !!$      IF( ln_tide_date ) THEN      ! correct for date factors 
    228 !!$ 
    229 !!$!! used nmonth, nyear and nday from daymod.... 
    230 !!$         ! Calculate date corrects for 15 standard consituents 
    231 !!$         ! This is the initialisation step, so nday, nmonth, nyear are the  
    232 !!$         ! initial date/time of the integration. 
    233 !!$           print *, nday,nmonth,nyear 
    234 !!$           nyear  = int(ndate0 / 10000  )                           ! initial year 
    235 !!$           nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month 
    236 !!$           nday   = int(ndate0 - nyear * 10000 - nmonth * 100) 
    237 !!$ 
    238 !!$         CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 
    239 !!$ 
    240 !!$         IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 
    241 !!$ 
    242 !!$         DO itide = 1, ntide       ! loop on tidal components 
    243 !!$            ! 
    244 !!$            IF( nindx(itide) /= 0 ) THEN 
    245 !!$!!gm use rpi  and rad global variable   
    246 !!$               z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 
    247 !!$               z_atde=z_ftc(nindx(itide))*cos(z_arg) 
    248 !!$               z_btde=z_ftc(nindx(itide))*sin(z_arg) 
    249 !!$               IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide),   & 
    250 !!$                  &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide)) 
    251 !!$            ELSE 
    252 !!$               z_atde = 1.0_wp 
    253 !!$               z_btde = 0.0_wp 
    254 !!$            ENDIF 
    255 !!$            !                                         !  elevation          
    256 !!$            igrd = 4 
    257 !!$            DO ib = 1, nblenrim(igrd)                 
    258 !!$               z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide) 
    259 !!$               z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide) 
    260 !!$               ssh1(ib,itide) = z1t 
    261 !!$               ssh2(ib,itide) = z2t 
    262 !!$            END DO 
    263 !!$            !                                         !  u        
    264 !!$            igrd = 5 
    265 !!$            DO ib = 1, nblenrim(igrd)                 
    266 !!$               z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide) 
    267 !!$               z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide) 
    268 !!$               u1(ib,itide) = z1t 
    269 !!$               u2(ib,itide) = z2t 
    270 !!$            END DO 
    271 !!$            !                                         !  v        
    272 !!$            igrd = 6 
    273 !!$            DO ib = 1, nblenrim(igrd)                 
    274 !!$               z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide) 
    275 !!$               z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide) 
    276 !!$               v1(ib,itide) = z1t 
    277 !!$               v2(ib,itide) = z2t 
    278 !!$            END DO 
    279 !!$            ! 
    280 !!$         END DO   ! end loop on tidal components 
    281 !!$         ! 
    282 !!$      ENDIF ! date correction 
    283 !!$      ! 
    284 !!$   END SUBROUTINE tide_data 
    285 !!$ 
    286 !!$ 
    287 !!$   SUBROUTINE tide_update ( kt, jit ) 
    288 !!$      !!---------------------------------------------------------------------- 
    289 !!$      !!                 ***  SUBROUTINE tide_update  *** 
    290 !!$      !!                 
    291 !!$      !! ** Purpose : - Add tidal forcing to ssh_obc, ubt_obc and vbt_obc arrays.  
    292 !!$      !!                 
    293 !!$      !!---------------------------------------------------------------------- 
    294 !!$      INTEGER, INTENT( in ) ::   kt      ! Main timestep counter 
    295 !!$!!gm doctor jit ==> kit 
    296 !!$      INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option) 
    297 !!$      !! 
    298 !!$      INTEGER  ::   itide, igrd, ib      ! dummy loop indices 
    299 !!$      REAL(wp) ::   z_arg, z_sarg            !             
    300 !!$      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
    301 !!$      !!---------------------------------------------------------------------- 
    302 !!$ 
    303 !!$      ! Note tide phase speeds are in deg/hour, so we need to convert the 
    304 !!$      ! elapsed time in seconds to hours by dividing by 3600.0 
    305 !!$      IF( jit == 0 ) THEN   
    306 !!$         z_arg = kt * rdt * rad / 3600.0 
    307 !!$      ELSE                              ! we are in a barotropic subcycle (for timesplitting option) 
    308 !!$!         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0 
    309 !!$         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 
    310 !!$      ENDIF 
    311 !!$ 
    312 !!$      DO itide = 1, ntide 
    313 !!$         z_sarg = z_arg * tide_speed(itide) 
    314 !!$         z_cost(itide) = COS( z_sarg ) 
    315 !!$         z_sist(itide) = SIN( z_sarg ) 
    316 !!$      END DO 
    317 !!$ 
    318 !!$      ! summing of tidal constituents into BDY arrays 
    319 !!$      sshtide(:) = 0.0 
    320 !!$      utide (:) = 0.0 
    321 !!$      vtide (:) = 0.0 
    322 !!$      ! 
    323 !!$      DO itide = 1, ntide 
    324 !!$         igrd=4                              ! SSH on tracer grid. 
    325 !!$         DO ib = 1, nblenrim(igrd) 
    326 !!$            sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide) 
    327 !!$            !    if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide) 
    328 !!$         END DO 
    329 !!$         igrd=5                              ! U grid 
    330 !!$         DO ib=1, nblenrim(igrd) 
    331 !!$            utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide) 
    332 !!$            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide) 
    333 !!$         END DO 
    334 !!$         igrd=6                              ! V grid 
    335 !!$         DO ib=1, nblenrim(igrd) 
    336 !!$            vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide) 
    337 !!$            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide) 
    338 !!$         END DO 
    339 !!$      END DO 
    340 !!$      ! 
    341 !!$   END SUBROUTINE tide_update 
    342 !!$ 
    343 !!$!!gm  doctor naming of dummy argument variables!!!   and all local variables 
    344 !!$ 
    345 !!$   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 
    346 !!$      !!---------------------------------------------------------------------- 
    347 !!$      !!                   ***  SUBROUTINE uvset  *** 
    348 !!$      !!                      
    349 !!$      !! ** Purpose : - adjust tidal forcing for date factors 
    350 !!$      !!                 
    351 !!$      !!---------------------------------------------------------------------- 
    352 !!$      implicit none 
    353 !!$      INTEGER, INTENT( in ) ::   ihs     ! Start time hours 
    354 !!$      INTEGER, INTENT( in ) ::   iday    ! start time days 
    355 !!$      INTEGER, INTENT( in ) ::   imnth   ! start time month 
    356 !!$      INTEGER, INTENT( in ) ::   iyr     ! start time year 
    357 !!$      !! 
    358 !!$!!gm  nc is jptides_max ???? 
    359 !!$      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents 
    360 !!$      CHARACTER(len=8), DIMENSION(nc) :: cname 
    361 !!$      INTEGER                 ::   year, vd, ivdy, ndc, i, k 
    362 !!$      REAL(wp)                ::   ss, h, p, en, p1, rtd 
    363 !!$      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction 
    364 !!$      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction 
    365 !!$      REAL(wp), DIMENSION(nc) ::   u, v, zig 
    366 !!$      !! 
    367 !!$      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   & 
    368 !!$         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   & 
    369 !!$         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      / 
    370 !!$      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   & 
    371 !!$         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   & 
    372 !!$         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   / 
    373 !!$      !!---------------------------------------------------------------------- 
    374 !!$! 
    375 !!$!     ihs             -  start time gmt on ... 
    376 !!$!     iday/imnth/iyr  -  date   e.g.   12/10/87 
    377 !!$! 
    378 !!$      CALL vday(iday,imnth,iyr,ivdy) 
    379 !!$      vd = ivdy 
    380 !!$! 
    381 !!$!rp   note change of year number for d. blackman shpen 
    382 !!$!rp   if(iyr.ge.1000) year=iyr-1900 
    383 !!$!rp   if(iyr.lt.1000) year=iyr 
    384 !!$      year = iyr 
    385 !!$! 
    386 !!$!.....year  =  year of required data 
    387 !!$!.....vd    =  day of required data..set up for 0000gmt day year 
    388 !!$      ndc = nc 
    389 !!$!.....ndc   =  number of constituents allowed 
    390 !!$!!gm use rpi ? 
    391 !!$      rtd = 360.0 / 6.2831852 
    392 !!$      DO i = 1, ndc 
    393 !!$         zig(i) = zig(i)*rtd 
    394 !!$         ! sigo(i)= zig(i) 
    395 !!$      END DO 
    396 !!$ 
    397 !!$!!gm try to avoid RETURN  in F90 
    398 !!$      IF( year == 0 )   RETURN 
    399 !!$       
    400 !!$         CALL shpen( year, vd, ss, h , p , en, p1 ) 
    401 !!$         CALL ufset( p   , en, u , f ) 
    402 !!$         CALL vset ( ss  , h , p , en, p1, v ) 
    403 !!$         ! 
    404 !!$         DO k = 1, nc 
    405 !!$            z_vplu(k) = v(k) + u(k) 
    406 !!$            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 
    407 !!$            DO WHILE( z_vplu(k) < 0    ) 
    408 !!$               z_vplu(k) = z_vplu(k) + 360.0 
    409 !!$            END DO 
    410 !!$            DO WHILE( z_vplu(k) > 360. ) 
    411 !!$               z_vplu(k) = z_vplu(k) - 360.0 
    412 !!$            END DO 
    413 !!$         END DO 
    414 !!$         ! 
    415 !!$      END SUBROUTINE uvset 
    416 !!$ 
    417 !!$ 
    418 !!$      SUBROUTINE vday( iday, imnth, iy, ivdy ) 
    419 !!$      !!---------------------------------------------------------------------- 
    420 !!$      !!                   *** SUBROUTINE vday  *** 
    421 !!$      !!                   
    422 !!$      !! ** Purpose : - adjust tidal forcing for date factors 
    423 !!$      !!                 
    424 !!$      !!---------------------------------------------------------------------- 
    425 !!$      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ???? 
    426 !!$      INTEGER, INTENT(  out) ::   ivdy              ! ??? 
    427 !!$      !!  
    428 !!$      INTEGER ::   iyr 
    429 !!$      !!------------------------------------------------------------------------------ 
    430 !!$ 
    431 !!$!!gm   nday_year in day mode is the variable compiuted here, no? 
    432 !!$!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year 
    433 !!$ 
    434 !!$      !calculate day number in year from day/month/year 
    435 !!$       if(imnth.eq.1) ivdy=iday 
    436 !!$       if(imnth.eq.2) ivdy=iday+31 
    437 !!$       if(imnth.eq.3) ivdy=iday+59 
    438 !!$       if(imnth.eq.4) ivdy=iday+90 
    439 !!$       if(imnth.eq.5) ivdy=iday+120 
    440 !!$       if(imnth.eq.6) ivdy=iday+151 
    441 !!$       if(imnth.eq.7) ivdy=iday+181 
    442 !!$       if(imnth.eq.8) ivdy=iday+212 
    443 !!$       if(imnth.eq.9) ivdy=iday+243 
    444 !!$       if(imnth.eq.10) ivdy=iday+273 
    445 !!$       if(imnth.eq.11) ivdy=iday+304 
    446 !!$       if(imnth.eq.12) ivdy=iday+334 
    447 !!$       iyr=iy 
    448 !!$       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
    449 !!$       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 
    450 !!$       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
    451 !!$      ! 
    452 !!$   END SUBROUTINE vday 
    453 !!$ 
    454 !!$!!doctor norme for dummy arguments 
    455 !!$ 
    456 !!$   SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 
    457 !!$      !!---------------------------------------------------------------------- 
    458 !!$      !!                    ***  SUBROUTINE shpen  *** 
    459 !!$      !!                    
    460 !!$      !! ** Purpose : - calculate astronomical arguments for tides 
    461 !!$      !!                this version from d. blackman 30 nove 1990 
    462 !!$      !! 
    463 !!$      !!---------------------------------------------------------------------- 
    464 !!$!!gm add INTENT in, out or inout....    DOCTOR name.... 
    465 !!$!!gm please do not use variable name with a single letter:  impossible to search in a code 
    466 !!$      INTEGER  ::   year,vd 
    467 !!$      REAL(wp) ::   s,h,p,en,p1 
    468 !!$      !! 
    469 !!$      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd 
    470 !!$      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat 
    471 !!$      !! 
    472 !!$      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           & 
    473 !!$         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    & 
    474 !!$         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    & 
    475 !!$         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    & 
    476 !!$         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    & 
    477 !!$         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    & 
    478 !!$         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    & 
    479 !!$         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    & 
    480 !!$         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    & 
    481 !!$         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    & 
    482 !!$         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    & 
    483 !!$         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    & 
    484 !!$         &       52.57 / 
    485 !!$      !!---------------------------------------------------------------------- 
    486 !!$ 
    487 !!$      cycle = 360.0 
    488 !!$      ilc = 0 
    489 !!$      icent = year / 100 
    490 !!$      yr = year - icent * 100 
    491 !!$      t = icent - 20 
    492 !!$! 
    493 !!$!     for the following equations 
    494 !!$!     time origin is fixed at 00 hr of jan 1st,2000. 
    495 !!$!     see notes by cartwright 
    496 !!$! 
    497 !!$!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90) 
    498 !!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    499 !!$      it = icent - 20 
    500 !!$      if (it) 1,2,2 
    501 !!$    1 iday = it/4 -it 
    502 !!$      go to 3 
    503 !!$    2 iday = (it+3)/4 - it 
    504 !!$! 
    505 !!$!     t is in julian century 
    506 !!$!     correction in gegorian calander where only century year divisible 
    507 !!$!     by 4 is leap year. 
    508 !!$! 
    509 !!$    3 continue 
    510 !!$! 
    511 !!$      td = 0.0 
    512 !!$! 
    513 !!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    514 !!$      if (yr) 4,5,4 
    515 !!$! 
    516 !!$    4 iyd = 365*yr 
    517 !!$      ild = (yr-1)/4 
    518 !!$      if((icent - (icent/4)*4) .eq. 0) ilc = 1 
    519 !!$      td = iyd + ild + ilc 
    520 !!$! 
    521 !!$    5 td = td + iday + vd -1.0 - 0.5 
    522 !!$      t = t + (td/36525.0) 
    523 !!$! 
    524 !!$      ipos=year-1899 
    525 !!$      if (ipos .lt. 0) go to 7 
    526 !!$      if (ipos .gt. 83) go to 6 
    527 !!$! 
    528 !!$      delta = (delt(ipos+1)+delt(ipos))/2.0 
    529 !!$      go to 7 
    530 !!$! 
    531 !!$    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5 
    532 !!$! 
    533 !!$    7 deltat = delta * 1.0e-6 
    534 !!$! 
    535 !!$!!gm   precision of the computation   :  example for s it should be replace by: 
    536 !!$!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results 
    537 !!$      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 
    538 !!$      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat 
    539 !!$      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat 
    540 !!$      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat 
    541 !!$      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t 
    542 !!$      ! 
    543 !!$      nn = s / cycle 
    544 !!$      s = s - nn * cycle 
    545 !!$      IF( s < 0.e0 )   s = s + cycle 
    546 !!$      ! 
    547 !!$      nn = h / cycle 
    548 !!$      h = h - cycle * nn 
    549 !!$      IF( h < 0.e0 )   h = h + cycle 
    550 !!$      ! 
    551 !!$      nn = p / cycle 
    552 !!$      p = p - cycle * nn 
    553 !!$      IF( p < 0.e0)   p = p + cycle 
    554 !!$      ! 
    555 !!$      nn = en / cycle 
    556 !!$      en = en - cycle * nn 
    557 !!$      IF( en < 0.e0 )   en = en + cycle 
    558 !!$      en = cycle - en 
    559 !!$      ! 
    560 !!$      nn = p1 / cycle 
    561 !!$      p1 = p1 - nn * cycle 
    562 !!$      ! 
    563 !!$   END SUBROUTINE shpen 
    564 !!$ 
    565 !!$ 
    566 !!$   SUBROUTINE ufset( p, cn, b, a ) 
    567 !!$      !!---------------------------------------------------------------------- 
    568 !!$      !!                    ***  SUBROUTINE ufset  *** 
    569 !!$      !!                     
    570 !!$      !! ** Purpose : - calculate nodal parameters for the tides 
    571 !!$      !!                 
    572 !!$      !!---------------------------------------------------------------------- 
    573 !!$!!gm  doctor naming of dummy argument variables!!!   and all local variables 
    574 !!$!!gm  nc is jptides_max ???? 
    575 !!$      integer nc 
    576 !!$      parameter (nc=15) 
    577 !!$      REAL(wp) p,cn 
    578 !!$      !! 
    579 !!$       
    580 !!$!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z" 
    581 !!$      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 
    582 !!$      REAL(wp) ::   a(nc), b(nc) 
    583 !!$      INTEGER  ::   ndc, k 
    584 !!$      !!---------------------------------------------------------------------- 
    585 !!$ 
    586 !!$      ndc = nc 
    587 !!$ 
    588 !!$!    a=f       ,  b =u 
    589 !!$!    t is  zero as compared to tifa. 
    590 !!$!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module 
    591 !!$      rad = 6.2831852d0/360.0 
    592 !!$      pw = p  * rad 
    593 !!$      nw = cn * rad 
    594 !!$      w1 = cos(   nw ) 
    595 !!$      w2 = cos( 2*nw ) 
    596 !!$      w3 = cos( 3*nw ) 
    597 !!$      w4 = sin(   nw ) 
    598 !!$      w5 = sin( 2*nw ) 
    599 !!$      w6 = sin( 3*nw ) 
    600 !!$      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   & 
    601 !!$         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw ) 
    602 !!$      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   & 
    603 !!$         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw ) 
    604 !!$      ! 
    605 !!$      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 
    606 !!$      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 
    607 !!$      !   q1 
    608 !!$      a(2) = a(1) 
    609 !!$      b(2) = b(1) 
    610 !!$      !   o1 
    611 !!$      a(3) = 1.0 
    612 !!$      b(3) = 0.0 
    613 !!$      !   p1 
    614 !!$      a(4) = 1.0 
    615 !!$      b(4) = 0.0 
    616 !!$      !   s1 
    617 !!$      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 
    618 !!$      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 
    619 !!$      !   k1 
    620 !!$      a(6) =1.0004 -0.0373*w1+ 0.0002*w2 
    621 !!$      b(6) = -0.0374*w4 
    622 !!$      !  2n2 
    623 !!$      a(7) = a(6) 
    624 !!$      b(7) = b(6) 
    625 !!$      !  mu2 
    626 !!$      a(8) = a(6) 
    627 !!$      b(8) = b(6) 
    628 !!$      !   n2 
    629 !!$      a(9) = a(6) 
    630 !!$      b(9) = b(6) 
    631 !!$      !  nu2 
    632 !!$      a(10) = a(6) 
    633 !!$      b(10) = b(6) 
    634 !!$      !   m2 
    635 !!$      a(11) = SQRT( w7 * w7 + w8 * w8 ) 
    636 !!$      b(11) = ATAN( w8 / w7 ) 
    637 !!$!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ???? 
    638 !!$      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992 
    639 !!$      !   l2 
    640 !!$      a(12) = 1.0 
    641 !!$      b(12) = 0.0 
    642 !!$      !   t2 
    643 !!$      a(13)= a(12) 
    644 !!$      b(13)= b(12) 
    645 !!$      !   s2 
    646 !!$      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 
    647 !!$      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 
    648 !!$      !   k2 
    649 !!$      a(15) = a(6)*a(6) 
    650 !!$      b(15) = 2*b(6) 
    651 !!$      !   m4 
    652 !!$!!gm  old coding,  remove GOTO and label of lines 
    653 !!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    654 !!$      DO 40 k = 1,ndc 
    655 !!$         b(k) = b(k)/rad 
    656 !!$32       if (b(k)) 34,35,35 
    657 !!$34       b(k) = b(k) + 360.0 
    658 !!$         go to 32 
    659 !!$35       if (b(k)-360.0) 40,37,37 
    660 !!$37       b(k) = b(k)-360.0 
    661 !!$         go to 35 
    662 !!$40    continue 
    663 !!$      ! 
    664 !!$   END SUBROUTINE ufset 
    665 !!$    
    666 !!$ 
    667 !!$   SUBROUTINE vset( s,h,p,en,p1,v) 
    668 !!$      !!---------------------------------------------------------------------- 
    669 !!$      !!                    ***  SUBROUTINE vset  *** 
    670 !!$      !!                    
    671 !!$      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 
    672 !!$      !!                 
    673 !!$      !!---------------------------------------------------------------------- 
    674 !!$!!gm  doctor naming of dummy argument variables!!!   and all local variables 
    675 !!$!!gm  nc is jptides_max ???? 
    676 !!$!!gm  en argument is not used: suppress it ? 
    677 !!$      integer nc 
    678 !!$      parameter (nc=15) 
    679 !!$      real(wp) s,h,p,en,p1 
    680 !!$      real(wp)   v(nc) 
    681 !!$      !! 
    682 !!$      integer ndc, k 
    683 !!$      !!---------------------------------------------------------------------- 
    684 !!$ 
    685 !!$      ndc = nc 
    686 !!$      !   v s  are computed here. 
    687 !!$      v(1) =-3*s +h +p +270      ! Q1 
    688 !!$      v(2) =-2*s +h +270.0       ! O1 
    689 !!$      v(3) =-h +270              ! P1 
    690 !!$      v(4) =180                  ! S1 
    691 !!$      v(5) =h +90.0              ! K1 
    692 !!$      v(6) =-4*s +2*h +2*p       ! 2N2 
    693 !!$      v(7) =-4*(s-h)             ! MU2 
    694 !!$      v(8) =-3*s +2*h +p         ! N2 
    695 !!$      v(9) =-3*s +4*h -p         ! MU2 
    696 !!$      v(10) =-2*s +2*h           ! M2 
    697 !!$      v(11) =-s +2*h -p +180     ! L2  
    698 !!$      v(12) =-h +p1              ! T2 
    699 !!$      v(13) =0                   ! S2 
    700 !!$      v(14) =h+h                 ! K2 
    701 !!$      v(15) =2*v(10)             ! M4 
    702 !!$! 
    703 !!$!!gm  old coding,  remove GOTO and label of lines 
    704 !!$!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
    705 !!$      do 72 k = 1, ndc 
    706 !!$69       if( v(k) )   70,71,71 
    707 !!$70       v(k) = v(k)+360.0 
    708 !!$         go to 69 
    709 !!$71       if( v(k) - 360.0 )   72,73,73 
    710 !!$73       v(k) = v(k)-360.0 
    711 !!$         go to 71 
    712 !!$72    continue 
    713 !!$      ! 
    714 !!$   END SUBROUTINE vset 
    715 !!$ 
    716 !!$#else 
     12   !!---------------------------------------------------------------------- 
     13   !!   'key_obc'     Open Boundary Condition 
     14   !!---------------------------------------------------------------------- 
     15   !!   PUBLIC 
     16   !!      tide_init     : read of namelist and initialisation of tidal harmonics data 
     17   !!      tide_update   : calculation of tidal forcing at each timestep 
     18   !!   PRIVATE 
     19   !!      uvset         :\ 
     20   !!      vday          : | Routines to correct tidal harmonics forcing for 
     21   !!      shpen         : | start time of integration 
     22   !!      ufset         : | 
     23   !!      vset          :/ 
     24   !!---------------------------------------------------------------------- 
     25   USE oce             ! ocean dynamics and tracers  
     26   USE dom_oce         ! ocean space and time domain 
     27   USE iom 
     28   USE in_out_manager  ! I/O units 
     29   USE phycst          ! physical constants 
     30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     31   USE obc_par         ! Unstructured boundary parameters 
     32   USE obc_oce         ! ocean open boundary conditions 
     33   USE daymod          ! calendar 
     34   USE fldread, ONLY: fld_map 
     35 
     36   IMPLICIT NONE 
     37   PRIVATE 
     38 
     39   PUBLIC   tide_init     ! routine called in nemo_init 
     40   PUBLIC   tide_update   ! routine called in obcdyn 
     41 
     42   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
     43      INTEGER                                ::   ncpt       !: Actual number of tidal components 
     44      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr) 
     45      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH 
     46      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U 
     47      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V 
     48   END TYPE TIDES_DATA 
     49 
     50   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents 
     51 
     52   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_obc), TARGET ::   tides                 !: External tidal harmonics data 
     53    
     54   !!---------------------------------------------------------------------- 
     55   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     56   !! $Id: obctides.F90 2528 2010-12-27 17:33:53Z rblod $  
     57   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     58   !!---------------------------------------------------------------------- 
     59CONTAINS 
     60 
     61   SUBROUTINE tide_init 
     62      !!---------------------------------------------------------------------- 
     63      !!                    ***  SUBROUTINE tide_init  *** 
     64      !!                      
     65      !! ** Purpose : - Read in namelist for tides and initialise external 
     66      !!                tidal harmonics data 
     67      !! 
     68      !!---------------------------------------------------------------------- 
     69      !! namelist variables 
     70      !!------------------- 
     71      LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date 
     72      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files 
     73      CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used. 
     74      REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr) 
     75      !! 
     76      INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components 
     77      INTEGER                                   ::   ib_obc, itide, ib  !: dummy loop indices 
     78      INTEGER                                   ::   inum, igrd 
     79      INTEGER, DIMENSION(3)                     ::   ilen0              !: length of boundary data (from OBC arrays) 
     80      CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file  
     81      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t            
     82      REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc 
     83      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data 
     84      !! 
     85      TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut    
     86      !! 
     87      NAMELIST/namobc_tide/ln_tide_date, filtide, tide_cpt, tide_speed 
     88      !!---------------------------------------------------------------------- 
     89 
     90      IF(lwp) WRITE(numout,*) 
     91      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' 
     92      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
     93 
     94      REWIND(numnam) 
     95      DO ib_obc = 1, nb_obc 
     96         IF( nn_tides(ib_obc) .gt. 0 ) THEN 
     97 
     98            td => tides(ib_obc) 
     99 
     100            ! Namelist namobc_tide : tidal harmonic forcing at open boundaries 
     101            ln_tide_date = .false. 
     102            filtide(:) = '' 
     103            tide_speed(:) = 0.0 
     104            tide_cpt(:) = '' 
     105 
     106            ! Don't REWIND here - may need to read more than one of these namelists. 
     107            READ  ( numnam, namobc_tide ) 
     108            !                                               ! Count number of components specified 
     109            td%ncpt = 0 
     110            DO itide = 1, jptides_max 
     111              IF( tide_cpt(itide) /= '' ) THEN 
     112                 td%ncpt = td%ncpt + 1 
     113              ENDIF 
     114            END DO 
     115 
     116            ! Fill in phase speeds from namelist 
     117            ALLOCATE( td%speed(td%ncpt) ) 
     118            td%speed = tide_speed(1:td%ncpt) 
     119 
     120            ! Find constituents in standard list 
     121            DO itide = 1, td%ncpt 
     122               nindx(itide) = 0 
     123               IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1 
     124               IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2 
     125               IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3 
     126               IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4 
     127               IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5 
     128               IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6 
     129               IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7 
     130               IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8 
     131               IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9 
     132               IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10 
     133               IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11 
     134               IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12 
     135               IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13 
     136               IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14 
     137               IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15 
     138               IF( nindx(itide) == 0  .AND. lwp ) THEN 
     139                  WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 
     140                  CALL ctl_warn( ctmp1 ) 
     141               ENDIF 
     142            END DO 
     143            !                                               ! Parameter control and print 
     144            IF( td%ncpt < 1 ) THEN  
     145               CALL ctl_stop( '          Did not find any tidal components in namelist namobc_tide' ) 
     146            ELSE 
     147               IF(lwp) WRITE(numout,*) '          Namelist namobc_tide : tidal harmonic forcing at open boundaries' 
     148               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt 
     149               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt) 
     150               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
     151               IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt) 
     152            ENDIF 
     153 
     154 
     155            ! Allocate space for tidal harmonics data -  
     156            ! get size from OBC data arrays 
     157            ! --------------------------------------- 
     158 
     159            ilen0(1) = SIZE( dta_obc(ib_obc)%ssh )  
     160            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 
     161 
     162            ilen0(2) = SIZE( dta_obc(ib_obc)%u2d )  
     163            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 
     164 
     165            ilen0(3) = SIZE( dta_obc(ib_obc)%v2d )  
     166            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 
     167 
     168            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 
     169 
     170 
     171            ! Open files and read in tidal forcing data 
     172            ! ----------------------------------------- 
     173 
     174            DO itide = 1, td%ncpt 
     175               !                                                              ! SSH fields 
     176               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
     177               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     178               CALL iom_open( clfile, inum ) 
     179               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,1) ) 
     180               td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1) 
     181               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,1) ) 
     182               td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     183               CALL iom_close( inum ) 
     184               !                                                              ! U fields 
     185               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
     186               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     187               CALL iom_open( clfile, inum ) 
     188               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,2) ) 
     189               td%u(:,itide,1) = dta_read(1:ilen0(2),1,1) 
     190               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,2) ) 
     191               td%u(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     192               CALL iom_close( inum ) 
     193               !                                                              ! V fields 
     194               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
     195               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     196               CALL iom_open( clfile, inum ) 
     197               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,3) ) 
     198               td%v(:,itide,1) = dta_read(1:ilen0(3),1,1) 
     199               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_obc(ib_obc)%nbmap(:,3) ) 
     200               td%v(:,itide,2) = dta_read(1:ilen0(3),1,1) 
     201               CALL iom_close( inum ) 
     202               ! 
     203            END DO ! end loop on tidal components 
     204 
     205            IF( ln_tide_date ) THEN      ! correct for date factors 
     206 
     207!! used nmonth, nyear and nday from daymod.... 
     208               ! Calculate date corrects for 15 standard consituents 
     209               ! This is the initialisation step, so nday, nmonth, nyear are the  
     210               ! initial date/time of the integration. 
     211                 print *, nday,nmonth,nyear 
     212                 nyear  = int(ndate0 / 10000  )                          ! initial year 
     213                 nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month 
     214                 nday   = int(ndate0 - nyear * 10000 - nmonth * 100) 
     215 
     216               CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 
     217 
     218               IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 
     219 
     220               DO itide = 1, td%ncpt       ! loop on tidal components 
     221                  ! 
     222                  IF( nindx(itide) /= 0 ) THEN 
     223!!gm use rpi  and rad global variable   
     224                     z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 
     225                     z_atde=z_ftc(nindx(itide))*cos(z_arg) 
     226                     z_btde=z_ftc(nindx(itide))*sin(z_arg) 
     227                     IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   & 
     228                        &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide)) 
     229                  ELSE 
     230                     z_atde = 1.0_wp 
     231                     z_btde = 0.0_wp 
     232                  ENDIF 
     233                  !                                         !  elevation          
     234                  igrd = 1 
     235                  DO ib = 1, ilen0(igrd)                 
     236                     z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2) 
     237                     z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1) 
     238                     td%ssh(ib,itide,1) = z1t 
     239                     td%ssh(ib,itide,2) = z2t 
     240                  END DO 
     241                  !                                         !  u        
     242                  igrd = 2 
     243                  DO ib = 1, ilen0(igrd)                 
     244                     z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2) 
     245                     z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1) 
     246                     td%u(ib,itide,1) = z1t 
     247                     td%u(ib,itide,2) = z2t 
     248                  END DO 
     249                  !                                         !  v        
     250                  igrd = 3 
     251                  DO ib = 1, ilen0(igrd)                 
     252                     z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2) 
     253                     z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1) 
     254                     td%v(ib,itide,1) = z1t 
     255                     td%v(ib,itide,2) = z2t 
     256                  END DO 
     257                  ! 
     258               END DO   ! end loop on tidal components 
     259               ! 
     260            ENDIF ! date correction 
     261            ! 
     262         ENDIF ! nn_tides .gt. 0 
     263         ! 
     264      END DO ! loop on ib_obc 
     265 
     266   END SUBROUTINE tide_init 
     267 
     268 
     269   SUBROUTINE tide_update ( kt, jit, idx, dta, td ) 
     270      !!---------------------------------------------------------------------- 
     271      !!                 ***  SUBROUTINE tide_update  *** 
     272      !!                 
     273      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
     274      !!                 
     275      !!---------------------------------------------------------------------- 
     276      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter 
     277!!gm doctor jit ==> kit 
     278      INTEGER, INTENT( in )          ::   jit     ! Barotropic timestep counter (for timesplitting option) 
     279      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices 
     280      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data 
     281      TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data 
     282      !! 
     283      INTEGER                          :: itide, igrd, ib      ! dummy loop indices 
     284      REAL(wp)                         :: z_arg, z_sarg       
     285      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
     286      !!---------------------------------------------------------------------- 
     287 
     288      ! Note tide phase speeds are in deg/hour, so we need to convert the 
     289      ! elapsed time in seconds to hours by dividing by 3600.0 
     290      IF( jit == 0 ) THEN   
     291         z_arg = kt * rdt * rad / 3600.0 
     292      ELSE                              ! we are in a barotropic subcycle (for timesplitting option) 
     293!         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0 
     294         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 
     295      ENDIF 
     296 
     297      DO itide = 1, td%ncpt 
     298         z_sarg = z_arg * td%speed(itide) 
     299         z_cost(itide) = COS( z_sarg ) 
     300         z_sist(itide) = SIN( z_sarg ) 
     301      END DO 
     302 
     303      ! 
     304      DO itide = 1, td%ncpt 
     305         igrd=1                              ! SSH on tracer grid. 
     306         DO ib = 1, idx%nblenrim(igrd) 
     307            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 
     308            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2) 
     309         END DO 
     310         igrd=2                              ! U grid 
     311         DO ib=1, idx%nblenrim(igrd) 
     312            dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) 
     313            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) 
     314         END DO 
     315         igrd=3                              ! V grid 
     316         DO ib=1, idx%nblenrim(igrd) 
     317            dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) 
     318            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) 
     319         END DO 
     320      END DO 
     321      ! 
     322   END SUBROUTINE tide_update 
     323 
     324!!gm  doctor naming of dummy argument variables!!!   and all local variables 
     325 
     326   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 
     327      !!---------------------------------------------------------------------- 
     328      !!                   ***  SUBROUTINE uvset  *** 
     329      !!                      
     330      !! ** Purpose : - adjust tidal forcing for date factors 
     331      !!                 
     332      !!---------------------------------------------------------------------- 
     333      implicit none 
     334      INTEGER, INTENT( in ) ::   ihs     ! Start time hours 
     335      INTEGER, INTENT( in ) ::   iday    ! start time days 
     336      INTEGER, INTENT( in ) ::   imnth   ! start time month 
     337      INTEGER, INTENT( in ) ::   iyr     ! start time year 
     338      !! 
     339!!gm  nc is jptides_max ???? 
     340      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents 
     341      CHARACTER(len=8), DIMENSION(nc) :: cname 
     342      INTEGER                 ::   year, vd, ivdy, ndc, i, k 
     343      REAL(wp)                ::   ss, h, p, en, p1, rtd 
     344      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction 
     345      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction 
     346      REAL(wp), DIMENSION(nc) ::   u, v, zig 
     347      !! 
     348      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   & 
     349         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   & 
     350         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      / 
     351      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   & 
     352         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   & 
     353         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   / 
     354      !!---------------------------------------------------------------------- 
     355! 
     356!     ihs             -  start time gmt on ... 
     357!     iday/imnth/iyr  -  date   e.g.   12/10/87 
     358! 
     359      CALL vday(iday,imnth,iyr,ivdy) 
     360      vd = ivdy 
     361! 
     362!rp   note change of year number for d. blackman shpen 
     363!rp   if(iyr.ge.1000) year=iyr-1900 
     364!rp   if(iyr.lt.1000) year=iyr 
     365      year = iyr 
     366! 
     367!.....year  =  year of required data 
     368!.....vd    =  day of required data..set up for 0000gmt day year 
     369      ndc = nc 
     370!.....ndc   =  number of constituents allowed 
     371!!gm use rpi ? 
     372      rtd = 360.0 / 6.2831852 
     373      DO i = 1, ndc 
     374         zig(i) = zig(i)*rtd 
     375         ! sigo(i)= zig(i) 
     376      END DO 
     377 
     378!!gm try to avoid RETURN  in F90 
     379      IF( year == 0 )   RETURN 
     380       
     381         CALL shpen( year, vd, ss, h , p , en, p1 ) 
     382         CALL ufset( p   , en, u , f ) 
     383         CALL vset ( ss  , h , p , en, p1, v ) 
     384         ! 
     385         DO k = 1, nc 
     386            z_vplu(k) = v(k) + u(k) 
     387            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 
     388            DO WHILE( z_vplu(k) < 0    ) 
     389               z_vplu(k) = z_vplu(k) + 360.0 
     390            END DO 
     391            DO WHILE( z_vplu(k) > 360. ) 
     392               z_vplu(k) = z_vplu(k) - 360.0 
     393            END DO 
     394         END DO 
     395         ! 
     396      END SUBROUTINE uvset 
     397 
     398 
     399      SUBROUTINE vday( iday, imnth, iy, ivdy ) 
     400      !!---------------------------------------------------------------------- 
     401      !!                   *** SUBROUTINE vday  *** 
     402      !!                   
     403      !! ** Purpose : - adjust tidal forcing for date factors 
     404      !!                 
     405      !!---------------------------------------------------------------------- 
     406      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ???? 
     407      INTEGER, INTENT(  out) ::   ivdy              ! ??? 
     408      !!  
     409      INTEGER ::   iyr 
     410      !!------------------------------------------------------------------------------ 
     411 
     412!!gm   nday_year in day mode is the variable compiuted here, no? 
     413!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year 
     414 
     415      !calculate day number in year from day/month/year 
     416       if(imnth.eq.1) ivdy=iday 
     417       if(imnth.eq.2) ivdy=iday+31 
     418       if(imnth.eq.3) ivdy=iday+59 
     419       if(imnth.eq.4) ivdy=iday+90 
     420       if(imnth.eq.5) ivdy=iday+120 
     421       if(imnth.eq.6) ivdy=iday+151 
     422       if(imnth.eq.7) ivdy=iday+181 
     423       if(imnth.eq.8) ivdy=iday+212 
     424       if(imnth.eq.9) ivdy=iday+243 
     425       if(imnth.eq.10) ivdy=iday+273 
     426       if(imnth.eq.11) ivdy=iday+304 
     427       if(imnth.eq.12) ivdy=iday+334 
     428       iyr=iy 
     429       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
     430       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 
     431       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 
     432      ! 
     433   END SUBROUTINE vday 
     434 
     435!!doctor norme for dummy arguments 
     436 
     437   SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 
     438      !!---------------------------------------------------------------------- 
     439      !!                    ***  SUBROUTINE shpen  *** 
     440      !!                    
     441      !! ** Purpose : - calculate astronomical arguments for tides 
     442      !!                this version from d. blackman 30 nove 1990 
     443      !! 
     444      !!---------------------------------------------------------------------- 
     445!!gm add INTENT in, out or inout....    DOCTOR name.... 
     446!!gm please do not use variable name with a single letter:  impossible to search in a code 
     447      INTEGER  ::   year,vd 
     448      REAL(wp) ::   s,h,p,en,p1 
     449      !! 
     450      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd 
     451      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat 
     452      !! 
     453      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           & 
     454         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    & 
     455         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    & 
     456         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    & 
     457         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    & 
     458         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    & 
     459         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    & 
     460         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    & 
     461         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    & 
     462         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    & 
     463         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    & 
     464         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    & 
     465         &       52.57 / 
     466      !!---------------------------------------------------------------------- 
     467 
     468      cycle = 360.0 
     469      ilc = 0 
     470      icent = year / 100 
     471      yr = year - icent * 100 
     472      t = icent - 20 
     473! 
     474!     for the following equations 
     475!     time origin is fixed at 00 hr of jan 1st,2000. 
     476!     see notes by cartwright 
     477! 
     478!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90) 
     479!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
     480      it = icent - 20 
     481      if (it) 1,2,2 
     482    1 iday = it/4 -it 
     483      go to 3 
     484    2 iday = (it+3)/4 - it 
     485! 
     486!     t is in julian century 
     487!     correction in gegorian calander where only century year divisible 
     488!     by 4 is leap year. 
     489! 
     490    3 continue 
     491! 
     492      td = 0.0 
     493! 
     494!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
     495      if (yr) 4,5,4 
     496! 
     497    4 iyd = 365*yr 
     498      ild = (yr-1)/4 
     499      if((icent - (icent/4)*4) .eq. 0) ilc = 1 
     500      td = iyd + ild + ilc 
     501! 
     502    5 td = td + iday + vd -1.0 - 0.5 
     503      t = t + (td/36525.0) 
     504! 
     505      ipos=year-1899 
     506      if (ipos .lt. 0) go to 7 
     507      if (ipos .gt. 83) go to 6 
     508! 
     509      delta = (delt(ipos+1)+delt(ipos))/2.0 
     510      go to 7 
     511! 
     512    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5 
     513! 
     514    7 deltat = delta * 1.0e-6 
     515! 
     516!!gm   precision of the computation   :  example for s it should be replace by: 
     517!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results 
     518      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 
     519      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat 
     520      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat 
     521      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat 
     522      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t 
     523      ! 
     524      nn = s / cycle 
     525      s = s - nn * cycle 
     526      IF( s < 0.e0 )   s = s + cycle 
     527      ! 
     528      nn = h / cycle 
     529      h = h - cycle * nn 
     530      IF( h < 0.e0 )   h = h + cycle 
     531      ! 
     532      nn = p / cycle 
     533      p = p - cycle * nn 
     534      IF( p < 0.e0)   p = p + cycle 
     535      ! 
     536      nn = en / cycle 
     537      en = en - cycle * nn 
     538      IF( en < 0.e0 )   en = en + cycle 
     539      en = cycle - en 
     540      ! 
     541      nn = p1 / cycle 
     542      p1 = p1 - nn * cycle 
     543      ! 
     544   END SUBROUTINE shpen 
     545 
     546 
     547   SUBROUTINE ufset( p, cn, b, a ) 
     548      !!---------------------------------------------------------------------- 
     549      !!                    ***  SUBROUTINE ufset  *** 
     550      !!                     
     551      !! ** Purpose : - calculate nodal parameters for the tides 
     552      !!                 
     553      !!---------------------------------------------------------------------- 
     554!!gm  doctor naming of dummy argument variables!!!   and all local variables 
     555!!gm  nc is jptides_max ???? 
     556      integer nc 
     557      parameter (nc=15) 
     558      REAL(wp) p,cn 
     559      !! 
     560       
     561!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z" 
     562      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 
     563      REAL(wp) ::   a(nc), b(nc) 
     564      INTEGER  ::   ndc, k 
     565      !!---------------------------------------------------------------------- 
     566 
     567      ndc = nc 
     568 
     569!    a=f       ,  b =u 
     570!    t is  zero as compared to tifa. 
     571!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module 
     572      rad = 6.2831852d0/360.0 
     573      pw = p  * rad 
     574      nw = cn * rad 
     575      w1 = cos(   nw ) 
     576      w2 = cos( 2*nw ) 
     577      w3 = cos( 3*nw ) 
     578      w4 = sin(   nw ) 
     579      w5 = sin( 2*nw ) 
     580      w6 = sin( 3*nw ) 
     581      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   & 
     582         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw ) 
     583      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   & 
     584         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw ) 
     585      ! 
     586      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 
     587      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 
     588      !   q1 
     589      a(2) = a(1) 
     590      b(2) = b(1) 
     591      !   o1 
     592      a(3) = 1.0 
     593      b(3) = 0.0 
     594      !   p1 
     595      a(4) = 1.0 
     596      b(4) = 0.0 
     597      !   s1 
     598      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 
     599      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 
     600      !   k1 
     601      a(6) =1.0004 -0.0373*w1+ 0.0002*w2 
     602      b(6) = -0.0374*w4 
     603      !  2n2 
     604      a(7) = a(6) 
     605      b(7) = b(6) 
     606      !  mu2 
     607      a(8) = a(6) 
     608      b(8) = b(6) 
     609      !   n2 
     610      a(9) = a(6) 
     611      b(9) = b(6) 
     612      !  nu2 
     613      a(10) = a(6) 
     614      b(10) = b(6) 
     615      !   m2 
     616      a(11) = SQRT( w7 * w7 + w8 * w8 ) 
     617      b(11) = ATAN( w8 / w7 ) 
     618!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ???? 
     619      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992 
     620      !   l2 
     621      a(12) = 1.0 
     622      b(12) = 0.0 
     623      !   t2 
     624      a(13)= a(12) 
     625      b(13)= b(12) 
     626      !   s2 
     627      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 
     628      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 
     629      !   k2 
     630      a(15) = a(6)*a(6) 
     631      b(15) = 2*b(6) 
     632      !   m4 
     633!!gm  old coding,  remove GOTO and label of lines 
     634!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
     635      DO 40 k = 1,ndc 
     636         b(k) = b(k)/rad 
     63732       if (b(k)) 34,35,35 
     63834       b(k) = b(k) + 360.0 
     639         go to 32 
     64035       if (b(k)-360.0) 40,37,37 
     64137       b(k) = b(k)-360.0 
     642         go to 35 
     64340    continue 
     644      ! 
     645   END SUBROUTINE ufset 
     646    
     647 
     648   SUBROUTINE vset( s,h,p,en,p1,v) 
     649      !!---------------------------------------------------------------------- 
     650      !!                    ***  SUBROUTINE vset  *** 
     651      !!                    
     652      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 
     653      !!                 
     654      !!---------------------------------------------------------------------- 
     655!!gm  doctor naming of dummy argument variables!!!   and all local variables 
     656!!gm  nc is jptides_max ???? 
     657!!gm  en argument is not used: suppress it ? 
     658      integer nc 
     659      parameter (nc=15) 
     660      real(wp) s,h,p,en,p1 
     661      real(wp)   v(nc) 
     662      !! 
     663      integer ndc, k 
     664      !!---------------------------------------------------------------------- 
     665 
     666      ndc = nc 
     667      !   v s  are computed here. 
     668      v(1) =-3*s +h +p +270      ! Q1 
     669      v(2) =-2*s +h +270.0       ! O1 
     670      v(3) =-h +270              ! P1 
     671      v(4) =180                  ! S1 
     672      v(5) =h +90.0              ! K1 
     673      v(6) =-4*s +2*h +2*p       ! 2N2 
     674      v(7) =-4*(s-h)             ! MU2 
     675      v(8) =-3*s +2*h +p         ! N2 
     676      v(9) =-3*s +4*h -p         ! MU2 
     677      v(10) =-2*s +2*h           ! M2 
     678      v(11) =-s +2*h -p +180     ! L2  
     679      v(12) =-h +p1              ! T2 
     680      v(13) =0                   ! S2 
     681      v(14) =h+h                 ! K2 
     682      v(15) =2*v(10)             ! M4 
     683! 
     684!!gm  old coding,  remove GOTO and label of lines 
     685!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90 
     686      do 72 k = 1, ndc 
     68769       if( v(k) )   70,71,71 
     68870       v(k) = v(k)+360.0 
     689         go to 69 
     69071       if( v(k) - 360.0 )   72,73,73 
     69173       v(k) = v(k)-360.0 
     692         go to 71 
     69372    continue 
     694      ! 
     695   END SUBROUTINE vset 
     696 
     697#else 
    717698   !!---------------------------------------------------------------------- 
    718699   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r2800 r2814  
    2424   IMPLICIT NONE 
    2525   PRIVATE    
     26  
     27   PUBLIC   fld_map    ! routine called by tides_init 
    2628 
    2729   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations 
     
    629631      !!                using a general mapping (for open boundaries) 
    630632      !!---------------------------------------------------------------------- 
     633#if defined key_obc 
    631634      USE obc_oce, ONLY:  dta_global         ! workspace to read in global data arrays 
     635#endif  
    632636 
    633637      INTEGER                   , INTENT(in ) ::   num     ! stream number 
    634638      CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name 
    635       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta     ! output field on model grid 
     639      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta   ! output field on model grid (2 dimensional) 
    636640      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice) 
    637641      INTEGER,  DIMENSION(:)    , INTENT(in ) ::   map     ! global-to-local mapping indices 
    638642      !! 
    639       INTEGER                  ::   ipi      ! length of boundary data on local process 
    640       INTEGER                  ::   ipj      ! length of dummy dimension ( = 1 ) 
    641       INTEGER                  ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    642       INTEGER                  ::   ilendta  ! length of data in file 
    643       INTEGER                  ::   idvar    ! variable ID 
    644       INTEGER                  ::   ib, ik   ! loop counters 
    645       INTEGER                  ::   ierr 
    646       !! 
    647       CHARACTER(len=80)                   :: zfile 
     643      INTEGER                                 ::   ipi      ! length of boundary data on local process 
     644      INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 ) 
     645      INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     646      INTEGER                                 ::   ilendta  ! length of data in file 
     647      INTEGER                                 ::   idvar    ! variable ID 
     648      INTEGER                                 ::   ib, ik   ! loop counters 
     649      INTEGER                                 ::   ierr 
     650      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read ! work space for global data 
    648651      !!--------------------------------------------------------------------- 
    649652             
     653#if defined key_obc 
     654      dta_read => dta_global 
     655#endif 
     656 
    650657      ipi = SIZE( dta, 1 ) 
    651658      ipj = 1 
     
    655662      ilendta = iom_file(num)%dimsz(1,idvar) 
    656663      IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta 
    657  
    658       CALL iom_get ( num, jpdom_unknown, clvar, dta_global(1:ilendta,1:ipj,1:ipk), nrec ) 
     664      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 
     665 
     666      SELECT CASE( ipk ) 
     667      CASE(1)    
     668         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec ) 
     669      CASE DEFAULT 
     670         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 
     671      END SELECT 
    659672      ! 
    660673      DO ib = 1, ipi 
    661674         DO ik = 1, ipk 
    662             dta(ib,1,ik) =  dta_global(map(ib),1,ik) 
     675            dta(ib,1,ik) =  dta_read(map(ib),1,ik) 
    663676         END DO 
    664677      END DO 
    665678 
    666679   END SUBROUTINE fld_map 
     680 
    667681 
    668682   SUBROUTINE fld_rot( kt, sd ) 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r2797 r2814  
    4747   USE obcini          ! open boundary cond. initialization (obc_init routine) 
    4848   USE obcdta          ! open boundary cond. initialization (obc_dta_init routine) 
     49   USE obctides        ! open boundary cond. initialization (tide_init routine) 
    4950   USE istate          ! initial state setting          (istate_init routine) 
    5051   USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine) 
     
    296297      IF( lk_obc        )   CALL     obc_init       ! Open boundaries initialisation 
    297298      IF( lk_obc        )   CALL     obc_dta_init   ! Open boundaries initialisation of external data arrays 
     299      IF( lk_obc        )   CALL     tide_init      ! Open boundaries initialisation of tidal harmonic forcing 
    298300 
    299301                            CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
Note: See TracChangeset for help on using the changeset viewer.