New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3116 for branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 – NEMO

Ignore:
Timestamp:
2011-11-15T21:55:40+01:00 (13 years ago)
Author:
cetlod
Message:

dev_NEMO_MERGE_2011: add in changes dev_NOC_UKMO_MERGE developments

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r2528 r3116  
    1111#if defined key_bdy 
    1212   !!---------------------------------------------------------------------- 
    13    !!   'key_bdy'     Unstructured Open Boundary Condition 
     13   !!   'key_bdy'     Open Boundary Condition 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   PUBLIC 
    16    !!      tide_init     : read of namelist 
    17    !!      tide_data     : read in and initialisation of tidal constituents at boundary 
     16   !!      tide_init     : read of namelist and initialisation of tidal harmonics data 
    1817   !!      tide_update   : calculation of tidal forcing at each timestep 
    1918   !!   PRIVATE 
     
    3332   USE bdy_oce         ! ocean open boundary conditions 
    3433   USE daymod          ! calendar 
     34   USE fldread, ONLY: fld_map 
    3535 
    3636   IMPLICIT NONE 
    3737   PRIVATE 
    3838 
    39    PUBLIC   tide_init     ! routine called in bdyini 
    40    PUBLIC   tide_data     ! routine called in bdyini 
     39   PUBLIC   tide_init     ! routine called in nemo_init 
    4140   PUBLIC   tide_update   ! routine called in bdydyn 
    4241 
    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 
     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_bdy), TARGET ::   tides                 !: External tidal harmonics data 
    5653    
    5754   !!---------------------------------------------------------------------- 
     
    6663      !!                    ***  SUBROUTINE tide_init  *** 
    6764      !!                      
    68       !! ** Purpose : - Read in namelist for tides 
    69       !! 
    70       !!---------------------------------------------------------------------- 
    71       INTEGER ::   itide                  ! dummy loop index  
     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_bdy, 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    
    7286      !! 
    7387      NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed 
     
    7892      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
    7993 
    80       ! Namelist nambdy_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, nambdy_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 nambdy_tide' ) 
    122       ELSE 
    123          IF(lwp) WRITE(numout,*) '          Namelist nambdy_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       ! 
     94      REWIND(numnam) 
     95      DO ib_bdy = 1, nb_bdy 
     96         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 
     97 
     98            td => tides(ib_bdy) 
     99 
     100            ! Namelist nambdy_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, nambdy_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 nambdy_tide' ) 
     146            ELSE 
     147               IF(lwp) WRITE(numout,*) '          Namelist nambdy_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_bdy(ib_bdy)%ssh )  
     160            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 
     161 
     162            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d )  
     163            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 
     164 
     165            ilen0(3) = SIZE( dta_bdy(ib_bdy)%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_bdy(ib_bdy)%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_bdy(ib_bdy)%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_bdy(ib_bdy)%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_bdy(ib_bdy)%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_bdy(ib_bdy)%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_bdy(ib_bdy)%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_dyn2d_dta(ib_bdy) .ge. 2 
     263         ! 
     264      END DO ! loop on ib_bdy 
     265 
    135266   END SUBROUTINE tide_init 
    136267 
    137268 
    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 ) 
     269   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) 
    288270      !!---------------------------------------------------------------------- 
    289271      !!                 ***  SUBROUTINE tide_update  *** 
    290272      !!                 
    291       !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.  
     273      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
    292274      !!                 
    293275      !!---------------------------------------------------------------------- 
    294       INTEGER, INTENT( in ) ::   kt      ! Main timestep counter 
     276      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter 
    295277!!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            !             
     278      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices 
     279      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data 
     280      TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data 
     281      INTEGER,INTENT(in),OPTIONAL    ::   jit     ! Barotropic timestep counter (for timesplitting option) 
     282      INTEGER,INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
     283                                                       ! is present then units = subcycle timesteps. 
     284                                                       ! time_offset = 0 => get data at "now" time level 
     285                                                       ! time_offset = -1 => get data at "before" time level 
     286                                                       ! time_offset = +1 => get data at "after" time level 
     287                                                       ! etc. 
     288      !! 
     289      INTEGER                          :: itide, igrd, ib      ! dummy loop indices 
     290      INTEGER                          :: time_add             ! time offset in units of timesteps 
     291      REAL(wp)                         :: z_arg, z_sarg       
    300292      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
    301293      !!---------------------------------------------------------------------- 
    302294 
     295      time_add = 0 
     296      IF( PRESENT(time_offset) ) THEN 
     297         time_add = time_offset 
     298      ENDIF 
     299          
    303300      ! Note tide phase speeds are in deg/hour, so we need to convert the 
    304301      ! 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 
     302      IF( PRESENT(jit) ) THEN   
     303         z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0 
     304      ELSE                               
     305         z_arg = (kt+time_add) * rdt * rad / 3600.0 
    310306      ENDIF 
    311307 
    312       DO itide = 1, ntide 
    313          z_sarg = z_arg * tide_speed(itide) 
     308      DO itide = 1, td%ncpt 
     309         z_sarg = z_arg * td%speed(itide) 
    314310         z_cost(itide) = COS( z_sarg ) 
    315311         z_sist(itide) = SIN( z_sarg ) 
    316312      END DO 
    317313 
    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) 
     314      DO itide = 1, td%ncpt 
     315         igrd=1                              ! SSH on tracer grid. 
     316         DO ib = 1, idx%nblenrim(igrd) 
     317            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 
     318            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2) 
    328319         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) 
     320         igrd=2                              ! U grid 
     321         DO ib=1, idx%nblenrim(igrd) 
     322            dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) 
     323            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) 
    333324         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) 
     325         igrd=3                              ! V grid 
     326         DO ib=1, idx%nblenrim(igrd) 
     327            dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) 
     328            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) 
    338329         END DO 
    339330      END DO 
Note: See TracChangeset for help on using the changeset viewer.