Changeset 2865


Ignore:
Timestamp:
2011-09-27T14:33:01+02:00 (9 years ago)
Author:
davestorkey
Message:
  1. Updates for dynspg_exp option.
  2. Implement time_offset functionality in obc_dta.
  3. Add option to specify boundaries in the namelist.
  4. Re-activate obc_vol option.
  5. Update to namelist control of tidal harmonics.
Location:
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC
Files:
13 edited

Legend:

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

    r2818 r2865  
    3131   USE obc_oce         ! ocean open boundary conditions 
    3232   USE obcdta          ! ocean open boundary conditions 
    33    USE obcdyn3d        ! ocean open boundary conditions 
     33   USE obcdyn          ! ocean open boundary conditions 
    3434   USE obcvol          ! ocean open boundary condition (obc_vol routines) 
    3535   USE in_out_manager  ! I/O manager 
     
    153153# if defined key_obc 
    154154      !                                !* OBC open boundaries 
    155       IF( .NOT. lk_dynspg_flt ) THEN 
    156  
    157          CALL obc_dyn3d( kt ) 
     155      IF( lk_dynspg_exp ) CALL obc_dyn( kt ) 
     156      IF( lk_dynspg_ts )  CALL obc_dyn( kt, dyn3d_only=.true. ) 
    158157 
    159158!!$!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called 
     
    163162!!$         ! 
    164163!!$         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask ) 
    165  
    166       ENDIF 
    167164      ! 
    168165# endif 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r2814 r2865  
    355355         !                                                !* Update the forcing (OBC and tides) 
    356356         !                                                !  ------------------ 
    357          IF( lk_obc )   CALL obc_dta ( kt, jit=jn  ) 
     357         IF( lk_obc )   CALL obc_dta ( kt, jit=jn, time_offset=+1 ) 
    358358 
    359359         !                                                !* after ssh_e 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90

    r2831 r2865  
    5757   LOGICAL                    ::   ln_mask_file             !: =T read obcmask from file 
    5858   LOGICAL                    ::   ln_vol                   !: =T volume correction              
    59    LOGICAL, DIMENSION(jp_obc) ::   ln_clim                  !: =T obc data files contain climatological data (time-cyclic) 
    6059   ! 
    6160   INTEGER                    ::   nb_obc                   !: number of open boundary sets 
    62    INTEGER, DIMENSION(jp_obc) ::   nn_rimwidth              !: boundary rim width 
    63    INTEGER, DIMENSION(jp_obc) ::   nn_dtactl           !: = 0 use the initial state as obc dta ;  
    64                                                             !: = 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 
     61   INTEGER, DIMENSION(jp_obc) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme 
    6862   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
    6963   !                                                        !  = 1 the volume will be constant during all the integration. 
     
    7165   INTEGER, DIMENSION(jp_obc) ::   nn_dyn2d_dta           !: = 0 use the initial state as obc dta ;  
    7266                                                            !: = 1 read it in a NetCDF file 
     67                                                            !: = 2 read tidal harmonic forcing from a NetCDF file 
     68                                                            !: = 3 read external data AND tidal harmonic forcing from NetCDF files 
    7369   INTEGER, DIMENSION(jp_obc) ::   nn_dyn3d                 ! Choice of boundary condition for baroclinic velocities  
    7470   INTEGER, DIMENSION(jp_obc) ::   nn_dyn3d_dta           !: = 0 use the initial state as obc dta ;  
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90

    r2831 r2865  
    5555CONTAINS 
    5656 
    57       SUBROUTINE obc_dta( kt, jit ) 
     57      SUBROUTINE obc_dta( kt, jit, time_offset ) 
    5858      !!---------------------------------------------------------------------- 
    5959      !!                   ***  SUBROUTINE obc_dta  *** 
     
    6969      INTEGER, INTENT( in )           ::   kt    ! ocean time-step index  
    7070      INTEGER, INTENT( in ), OPTIONAL ::   jit   ! subcycle time-step index (for timesplitting option) 
     71      INTEGER, INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
     72                                                        ! is present then units = subcycle timesteps. 
     73                                                        ! time_offset = 0 => get data at "now" time level 
     74                                                        ! time_offset = -1 => get data at "before" time level 
     75                                                        ! time_offset = +1 => get data at "after" time level 
     76                                                        ! etc. 
    7177      !! 
    7278      INTEGER     ::  ib_obc, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices 
     
    202208            IF( PRESENT(jit) ) THEN 
    203209               ! Update barotropic boundary conditions only 
    204                ! jit is optional argument for fld_read 
    205                IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) THEN 
    206                   IF( nn_tides(ib_obc) .eq. 1 ) THEN 
     210               ! jit is optional argument for fld_read and tide_update 
     211               IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN 
     212                  IF( nn_dyn2d_dta(ib_obc) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    207213                     dta_obc(ib_obc)%ssh(:) = 0.0 
    208214                     dta_obc(ib_obc)%u2d(:) = 0.0 
    209215                     dta_obc(ib_obc)%v2d(:) = 0.0 
    210                   ELSE 
     216                  ENDIF 
     217                  IF( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) THEN ! update external data 
    211218                     jend = jstart + 2 
    212                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit ) 
     219                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset ) 
    213220                  ENDIF 
    214                ENDIF 
    215                IF( nn_tides(ib_obc) .gt. 0 ) THEN 
    216                   CALL tide_update( kt=kt, jit=jit, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) ) 
     221                  IF( nn_dyn2d_dta(ib_obc) .ge. 2 ) THEN ! update tidal harmonic forcing 
     222                     CALL tide_update( kt=kt, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc), jit=jit, time_offset=time_offset ) 
     223                  ENDIF 
    217224               ENDIF 
    218225            ELSE 
    219                IF( nb_obc_fld(ib_obc) .gt. 0 ) THEN  
    220                   jend = jstart + nb_obc_fld(ib_obc) - 1 
    221                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), timeshift=1 ) 
    222                ENDIF 
    223                IF( nn_tides(ib_obc) .eq. 1 ) THEN 
     226               IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    224227                  dta_obc(ib_obc)%ssh(:) = 0.0 
    225228                  dta_obc(ib_obc)%u2d(:) = 0.0 
    226229                  dta_obc(ib_obc)%v2d(:) = 0.0 
    227230               ENDIF 
    228                IF( nn_tides(ib_obc) .gt. 0 ) THEN 
    229                   !!! THINK ABOUT kt, jit VALUES !!! 
    230                   CALL tide_update( kt=kt, jit=0, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) ) 
     231               IF( nb_obc_fld(ib_obc) .gt. 0 ) THEN ! update external data 
     232                  jend = jstart + nb_obc_fld(ib_obc) - 1 
     233                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 
     234               ENDIF 
     235               IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .ge. 2 ) THEN ! update tidal harmonic forcing 
     236                  CALL tide_update( kt=kt, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc), time_offset=time_offset ) 
    231237               ENDIF 
    232238            ENDIF 
     
    238244 
    239245            IF( ln_full_vel_array(ib_obc) .and.                                             &  
    240            &    ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn3d_dta(ib_obc) .eq. 1 ) ) THEN  
     246           &    ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 .or. nn_dyn3d_dta(ib_obc) .eq. 1 ) ) THEN  
    241247 
    242248               igrd = 2                      ! zonal velocity 
     
    326332#endif 
    327333                              ) 
     334         IF(nn_dta(ib_obc) .gt. 1) nn_dta(ib_obc) = 1 
    328335      END DO 
    329336 
     
    333340      nb_obc_fld(:) = 0 
    334341      DO ib_obc = 1, nb_obc          
    335          IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) THEN 
     342         IF( nn_dyn2d(ib_obc) .gt. 0 .and. ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) THEN 
    336343            nb_obc_fld(ib_obc) = nb_obc_fld(ib_obc) + 3 
    337344         ENDIF 
     
    408415            ! Only read in necessary fields for this set. 
    409416            ! Important that barotropic variables come first. 
    410             IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) THEN  
    411  
    412                IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1) THEN 
     417            IF( nn_dyn2d(ib_obc) .gt. 0 .and. ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) THEN  
     418 
     419               IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN 
    413420                  jfld = jfld + 1 
    414421                  blf_i(jfld) = bn_ssh 
     
    419426               ENDIF 
    420427 
    421                IF( .not. ln_full_vel_array(ib_obc) .and. nn_tides(ib_obc) .ne. 1 ) THEN 
     428               IF( .not. ln_full_vel_array(ib_obc) ) THEN 
    422429 
    423430                  jfld = jfld + 1 
     
    449456            ! baroclinic velocities 
    450457            IF( ( nn_dyn3d(ib_obc) .gt. 0 .and. nn_dyn3d_dta(ib_obc) .eq. 1 ) .or. & 
    451                   ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) ) THEN 
     458           &      ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 .and.  & 
     459           &        ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) ) THEN 
    452460 
    453461               jfld = jfld + 1 
     
    583591 
    584592         IF (nn_dyn2d(ib_obc) .gt. 0) THEN 
    585             IF( nn_dyn2d_dta(ib_obc) .eq. 0 .or. ln_full_vel_array(ib_obc) .or. nn_tides(ib_obc) .eq. 1 ) THEN 
     593            IF( nn_dyn2d_dta(ib_obc) .eq. 0 .or. nn_dyn2d_dta(ib_obc) .eq. 2 .or. ln_full_vel_array(ib_obc) ) THEN 
    586594               IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN 
    587595                  ilen0(1:3) = nblen(1:3) 
     
    614622         ENDIF 
    615623         IF ( ( nn_dyn3d(ib_obc) .gt. 0 .and. nn_dyn3d_dta(ib_obc) .eq. 1 ).or. & 
    616            &  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) ) THEN 
     624           &  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 .and.   & 
     625           &    ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) ) THEN 
    617626            jfld = jfld + 1 
    618627            dta_obc(ib_obc)%u3d => bf(jfld)%fnow(:,1,:) 
     
    646655                  ilen0(1:3) = nblenrim(1:3) 
    647656               ENDIF 
    648                ALLOCATE( dta_obc(ib_obc)%ssh(ilen0(1)) ) 
    649                ALLOCATE( dta_obc(ib_obc)%u2d(ilen0(1)) ) 
    650                ALLOCATE( dta_obc(ib_obc)%v2d(ilen0(1)) ) 
     657               ALLOCATE( dta_obc(ib_obc)%frld(ilen0(1)) ) 
     658               ALLOCATE( dta_obc(ib_obc)%hicif(ilen0(1)) ) 
     659               ALLOCATE( dta_obc(ib_obc)%hsnif(ilen0(1)) ) 
    651660            ELSE 
    652661               jfld = jfld + 1 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90

    r2800 r2865  
    4141CONTAINS 
    4242 
    43    SUBROUTINE obc_dyn( kt ) 
     43   SUBROUTINE obc_dyn( kt, dyn3d_only ) 
    4444      !!---------------------------------------------------------------------- 
    4545      !!                  ***  SUBROUTINE obc_dyn  *** 
     
    5151      USE wrk_nemo, ONLY: wrk_2d_7, wrk_2d_8      ! 2D workspace 
    5252      !! 
    53       INTEGER, INTENT( in ) :: kt     ! Main time step counter 
     53      INTEGER, INTENT( in )           :: kt               ! Main time step counter 
     54      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities 
    5455      !! 
    5556      INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter 
     57      LOGICAL               :: ll_dyn2d, ll_dyn3d   
    5658      !! 
    5759 
     
    6062      END IF 
    6163 
     64      ll_dyn2d = .true. 
     65      ll_dyn3d = .true. 
     66 
     67      IF( PRESENT(dyn3d_only) ) THEN 
     68         IF( dyn3d_only ) ll_dyn2d = .false. 
     69      ENDIF 
     70 
    6271      !------------------------------------------------------- 
    6372      ! Set pointers 
    6473      !------------------------------------------------------- 
    6574 
    66       pssh => sshn_b 
     75      pssh => sshn 
    6776      phur => hur 
    6877      phvr => hvr 
     
    92101      !------------------------------------------------------- 
    93102 
    94       CALL obc_dyn2d( kt ) 
     103      IF( ll_dyn2d ) CALL obc_dyn2d( kt ) 
    95104 
    96       CALL obc_dyn3d( kt ) 
     105      IF( ll_dyn3d ) CALL obc_dyn3d( kt ) 
    97106 
    98107      !------------------------------------------------------- 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90

    r2818 r2865  
    5555            CYCLE 
    5656         CASE(jp_frs) 
    57             CALL obc_dyn2d_frs( idx_obc(ib_obc), dta_obc(ib_obc), kt ) 
     57            CALL obc_dyn2d_frs( idx_obc(ib_obc), dta_obc(ib_obc) ) 
    5858         CASE(jp_flather) 
    5959            CALL obc_dyn2d_fla( idx_obc(ib_obc), dta_obc(ib_obc) ) 
     
    6565   END SUBROUTINE obc_dyn2d 
    6666 
    67    SUBROUTINE obc_dyn2d_frs( idx, dta, kt ) 
     67   SUBROUTINE obc_dyn2d_frs( idx, dta ) 
    6868      !!---------------------------------------------------------------------- 
    6969      !!                  ***  SUBROUTINE obc_dyn2d_frs  *** 
     
    7676      !!               topography. Tellus, 365-382. 
    7777      !!---------------------------------------------------------------------- 
    78       INTEGER,         INTENT(in) ::   kt 
    7978      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    8079      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90

    r2831 r2865  
    5151      !! ** Input   :  obc_init.nc, input file for unstructured open boundaries 
    5252      !!----------------------------------------------------------------------       
    53       INTEGER  ::   ib_obc, ii, ij, ik, igrd, ib, ir       ! dummy loop indices 
     53      ! namelist variables 
     54      !------------------- 
     55      INTEGER, PARAMETER          :: jp_nseg = 100 
     56      INTEGER                     :: nobcsege, nobcsegw, nobcsegn, nobcsegs  
     57      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft 
     58      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft 
     59      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft 
     60      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft 
     61 
     62      ! local variables 
     63      !------------------- 
     64      INTEGER  ::   ib_obc, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices 
    5465      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers 
    5566      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       - 
     
    7283         &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
    7384#endif 
    74          &             nn_tides, ln_vol, ln_clim, nn_volctl,               & 
     85         &             ln_vol, nn_volctl,                                  & 
    7586         &             nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out,             & 
    7687         &             nn_dmp3d_in, nn_dmp3d_out 
     88      !! 
     89      NAMELIST/namobc_index/ nobcsege, jpieob, jpjedt, jpjeft,             & 
     90                             nobcsegw, jpiwob, jpjwdt, jpjwft,             & 
     91                             nobcsegn, jpjnob, jpindt, jpinft,             & 
     92                             nobcsegs, jpjsob, jpisdt, jpisft 
     93 
    7794      !!---------------------------------------------------------------------- 
    7895 
     
    109126#endif 
    110127      ln_vol            = .false. 
    111       ln_clim(:)        = .false. 
    112       nn_tides(:)       =  0  ! default to no tidal forcing 
    113128      nn_volctl         = -1  ! uninitialised flag 
    114129      nn_rimwidth(:)    = -1  ! uninitialised flag 
     
    137152        IF(lwp) WRITE(numout,*) ' '  
    138153        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_obc,'------'  
     154 
     155        IF( ln_coords_file(ib_obc) ) THEN 
     156           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_obc)) 
     157        ELSE 
     158           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.' 
     159        ENDIF 
     160        IF(lwp) WRITE(numout,*) 
    139161 
    140162        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
     
    149171              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for obc data'         
    150172              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    151               CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be 0 or 1' ) 
     173              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file' 
     174              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files' 
     175              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
    152176           END SELECT 
    153177        ENDIF 
     
    201225#endif 
    202226 
    203         IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth 
     227        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_obc) 
    204228        IF(lwp) WRITE(numout,*) 
    205  
    206         SELECT CASE( nn_tides(ib_obc) )  
    207         CASE(0) 
    208           IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing' 
    209           IF(lwp) WRITE(numout,*) 
    210         CASE(1) 
    211           IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution' 
    212           IF(lwp) WRITE(numout,*) 
    213         CASE(2) 
    214           IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions' 
    215           IF(lwp) WRITE(numout,*) 
    216         CASE DEFAULT 
    217           CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' ) 
    218         END SELECT 
    219229 
    220230      ENDDO 
     
    240250      ! Work out global dimensions of boundary data 
    241251      ! --------------------------------------------- 
     252      REWIND( numnam )                     
    242253      DO ib_obc = 1, nb_obc 
    243254 
     
    245256         IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Work out size of global arrays from namelist parameters 
    246257  
    247  
    248            !! 1. Read parameters from namelist 
    249            !! 2. Work out global size of boundary data arrays nblendta and jpbdta 
     258            ! No REWIND here because may need to read more than one namobc_index namelist. 
     259            READ  ( numnam, namobc_index ) 
     260 
     261            ! Automatic boundary definition: if nobcsegX = -1 
     262            ! set boundary to whole side of model domain. 
     263            IF( nobcsege == -1 ) THEN 
     264               nobcsege = 1 
     265               jpieob(1) = jpiglo - 1 
     266               jpjedt(1) = 2 
     267               jpjeft(1) = jpjglo - 1 
     268            ENDIF 
     269            IF( nobcsegw == -1 ) THEN 
     270               nobcsegw = 1 
     271               jpiwob(1) = 2 
     272               jpjwdt(1) = 2 
     273               jpjwft(1) = jpjglo - 1 
     274            ENDIF 
     275            IF( nobcsegn == -1 ) THEN 
     276               nobcsegn = 1 
     277               jpjnob(1) = jpjglo - 1 
     278               jpindt(1) = 2 
     279               jpinft(1) = jpiglo - 1 
     280            ENDIF 
     281            IF( nobcsegs == -1 ) THEN 
     282               nobcsegs = 1 
     283               jpjsob(1) = 2 
     284               jpisdt(1) = 2 
     285               jpisft(1) = jpiglo - 1 
     286            ENDIF 
     287 
     288            nblendta(:,ib_obc) = 0 
     289            DO iseg = 1, nobcsege 
     290               igrd = 1 
     291               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjeft(iseg) - jpjedt(iseg) + 1                
     292               igrd = 2 
     293               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjeft(iseg) - jpjedt(iseg) + 1                
     294               igrd = 3 
     295               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjeft(iseg) - jpjedt(iseg)                
     296            ENDDO 
     297            DO iseg = 1, nobcsegw 
     298               igrd = 1 
     299               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     300               igrd = 2 
     301               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     302               igrd = 3 
     303               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjwft(iseg) - jpjwdt(iseg)                
     304            ENDDO 
     305            DO iseg = 1, nobcsegn 
     306               igrd = 1 
     307               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpinft(iseg) - jpindt(iseg) + 1                
     308               igrd = 2 
     309               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpinft(iseg) - jpindt(iseg)                
     310               igrd = 3 
     311               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpinft(iseg) - jpindt(iseg) + 1 
     312            ENDDO 
     313            DO iseg = 1, nobcsegs 
     314               igrd = 1 
     315               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpisft(iseg) - jpisdt(iseg) + 1                
     316               igrd = 2 
     317               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpisft(iseg) - jpisdt(iseg) 
     318               igrd = 3 
     319               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpisft(iseg) - jpisdt(iseg) + 1                
     320            ENDDO 
     321 
     322            nblendta(:,ib_obc) = nblendta(:,ib_obc) * nn_rimwidth(ib_obc) 
     323            jpbdta = MAXVAL(nblendta(:,ib_obc))                
    250324 
    251325 
     
    263337         ENDIF  
    264338 
    265       ENDDO 
     339      ENDDO ! ib_obc 
    266340 
    267341      ! Allocate arrays 
     
    274348      ! Calculate global boundary index arrays or read in from file 
    275349      !------------------------------------------------------------ 
     350      REWIND( numnam )                     
    276351      DO ib_obc = 1, nb_obc 
    277352 
    278353         IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Calculate global index arrays from namelist parameters 
    279354 
    280             !! Calculate global index arrays from namelist parameters 
     355            ! No REWIND here because may need to read more than one namobc_index namelist. 
     356            READ  ( numnam, namobc_index ) 
     357 
     358            ! Automatic boundary definition: if nobcsegX = -1 
     359            ! set boundary to whole side of model domain. 
     360            IF( nobcsege == -1 ) THEN 
     361               nobcsege = 1 
     362               jpieob(1) = jpiglo - 1 
     363               jpjedt(1) = 2 
     364               jpjeft(1) = jpjglo - 1 
     365            ENDIF 
     366            IF( nobcsegw == -1 ) THEN 
     367               nobcsegw = 1 
     368               jpiwob(1) = 2 
     369               jpjwdt(1) = 2 
     370               jpjwft(1) = jpjglo - 1 
     371            ENDIF 
     372            IF( nobcsegn == -1 ) THEN 
     373               nobcsegn = 1 
     374               jpjnob(1) = jpjglo - 1 
     375               jpindt(1) = 2 
     376               jpinft(1) = jpiglo - 1 
     377            ENDIF 
     378            IF( nobcsegs == -1 ) THEN 
     379               nobcsegs = 1 
     380               jpjsob(1) = 2 
     381               jpisdt(1) = 2 
     382               jpisft(1) = jpiglo - 1 
     383            ENDIF 
     384 
     385            ! ------------ T points ------------- 
     386            igrd = 1   
     387            icount = 0 
     388            DO ir = 1, nn_rimwidth(ib_obc) 
     389               ! east 
     390               DO iseg = 1, nobcsege 
     391                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     392                     icount = icount + 1 
     393                     nbidta(icount, igrd, ib_obc) = jpieob(iseg) - ir + 1 
     394                     nbjdta(icount, igrd, ib_obc) = ij 
     395                     nbrdta(icount, igrd, ib_obc) = ir 
     396                  ENDDO 
     397               ENDDO 
     398               ! west 
     399               DO iseg = 1, nobcsegw 
     400                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     401                     icount = icount + 1 
     402                     nbidta(icount, igrd, ib_obc) = jpiwob(iseg) + ir - 1 
     403                     nbjdta(icount, igrd, ib_obc) = ij 
     404                     nbrdta(icount, igrd, ib_obc) = ir 
     405                  ENDDO 
     406               ENDDO 
     407               ! north 
     408               DO iseg = 1, nobcsegn 
     409                  DO ii = jpindt(iseg), jpinft(iseg) 
     410                     icount = icount + 1 
     411                     nbidta(icount, igrd, ib_obc) = ii 
     412                     nbjdta(icount, igrd, ib_obc) = jpjnob(iseg) - ir + 1 
     413                     nbrdta(icount, igrd, ib_obc) = ir 
     414                  ENDDO 
     415               ENDDO 
     416               ! south 
     417               DO iseg = 1, nobcsegs 
     418                  DO ii = jpisdt(iseg), jpisft(iseg) 
     419                     icount = icount + 1 
     420                     nbidta(icount, igrd, ib_obc) = ii 
     421                     nbjdta(icount, igrd, ib_obc) = jpjsob(iseg) + ir + 1 
     422                     nbrdta(icount, igrd, ib_obc) = ir 
     423                  ENDDO 
     424               ENDDO 
     425            ENDDO 
     426 
     427            ! ------------ U points ------------- 
     428            igrd = 2   
     429            icount = 0 
     430            DO ir = 1, nn_rimwidth(ib_obc) 
     431               ! east 
     432               DO iseg = 1, nobcsege 
     433                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     434                     icount = icount + 1 
     435                     nbidta(icount, igrd, ib_obc) = jpieob(iseg) - ir 
     436                     nbjdta(icount, igrd, ib_obc) = ij 
     437                     nbrdta(icount, igrd, ib_obc) = ir 
     438                  ENDDO 
     439               ENDDO 
     440               ! west 
     441               DO iseg = 1, nobcsegw 
     442                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     443                     icount = icount + 1 
     444                     nbidta(icount, igrd, ib_obc) = jpiwob(iseg) + ir - 1 
     445                     nbjdta(icount, igrd, ib_obc) = ij 
     446                     nbrdta(icount, igrd, ib_obc) = ir 
     447                  ENDDO 
     448               ENDDO 
     449               ! north 
     450               DO iseg = 1, nobcsegn 
     451                  DO ii = jpindt(iseg), jpinft(iseg) - 1 
     452                     icount = icount + 1 
     453                     nbidta(icount, igrd, ib_obc) = ii 
     454                     nbjdta(icount, igrd, ib_obc) = jpjnob(iseg) - ir + 1 
     455                     nbrdta(icount, igrd, ib_obc) = ir 
     456                  ENDDO 
     457               ENDDO 
     458               ! south 
     459               DO iseg = 1, nobcsegs 
     460                  DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     461                     icount = icount + 1 
     462                     nbidta(icount, igrd, ib_obc) = ii 
     463                     nbjdta(icount, igrd, ib_obc) = jpjsob(iseg) + ir + 1 
     464                     nbrdta(icount, igrd, ib_obc) = ir 
     465                  ENDDO 
     466               ENDDO 
     467            ENDDO 
     468 
     469            ! ------------ V points ------------- 
     470            igrd = 3   
     471            icount = 0 
     472            DO ir = 1, nn_rimwidth(ib_obc) 
     473               ! east 
     474               DO iseg = 1, nobcsege 
     475                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     476                     icount = icount + 1 
     477                     nbidta(icount, igrd, ib_obc) = jpieob(iseg) - ir + 1 
     478                     nbjdta(icount, igrd, ib_obc) = ij 
     479                     nbrdta(icount, igrd, ib_obc) = ir 
     480                  ENDDO 
     481               ENDDO 
     482               ! west 
     483               DO iseg = 1, nobcsegw 
     484                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     485                     icount = icount + 1 
     486                     nbidta(icount, igrd, ib_obc) = jpiwob(iseg) + ir - 1 
     487                     nbjdta(icount, igrd, ib_obc) = ij 
     488                     nbrdta(icount, igrd, ib_obc) = ir 
     489                  ENDDO 
     490               ENDDO 
     491               ! north 
     492               DO iseg = 1, nobcsegn 
     493                  DO ii = jpindt(iseg), jpinft(iseg) 
     494                     icount = icount + 1 
     495                     nbidta(icount, igrd, ib_obc) = ii 
     496                     nbjdta(icount, igrd, ib_obc) = jpjnob(iseg) - ir 
     497                     nbrdta(icount, igrd, ib_obc) = ir 
     498                  ENDDO 
     499               ENDDO 
     500               ! south 
     501               DO iseg = 1, nobcsegs 
     502                  DO ii = jpisdt(iseg), jpisft(iseg) 
     503                     icount = icount + 1 
     504                     nbidta(icount, igrd, ib_obc) = ii 
     505                     nbjdta(icount, igrd, ib_obc) = jpjsob(iseg) + ir + 1 
     506                     nbrdta(icount, igrd, ib_obc) = ir 
     507                  ENDDO 
     508               ENDDO 
     509            ENDDO 
    281510 
    282511         ELSE            ! Read global index arrays from boundary coordinates file. 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90

    r2814 r2865  
    9494      REWIND(numnam) 
    9595      DO ib_obc = 1, nb_obc 
    96          IF( nn_tides(ib_obc) .gt. 0 ) THEN 
     96         IF( nn_dyn2d_dta(ib_obc) .ge. 2 ) THEN 
    9797 
    9898            td => tides(ib_obc) 
     
    260260            ENDIF ! date correction 
    261261            ! 
    262          ENDIF ! nn_tides .gt. 0 
     262         ENDIF ! nn_dyn2d_dta(ib_obc) .ge. 2 
    263263         ! 
    264264      END DO ! loop on ib_obc 
     
    267267 
    268268 
    269    SUBROUTINE tide_update ( kt, jit, idx, dta, td ) 
     269   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) 
    270270      !!---------------------------------------------------------------------- 
    271271      !!                 ***  SUBROUTINE tide_update  *** 
     
    276276      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter 
    277277!!gm doctor jit ==> kit 
    278       INTEGER, INTENT( in )          ::   jit     ! Barotropic timestep counter (for timesplitting option) 
    279278      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices 
    280279      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data 
    281280      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. 
    282288      !! 
    283289      INTEGER                          :: itide, igrd, ib      ! dummy loop indices 
     290      INTEGER                          :: time_add             ! time offset in units of timesteps 
    284291      REAL(wp)                         :: z_arg, z_sarg       
    285292      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 
    286293      !!---------------------------------------------------------------------- 
    287294 
     295      time_add = 0 
     296      IF( PRESENT(time_offset) ) THEN 
     297         time_add = time_offset 
     298      ENDIF 
     299          
    288300      ! Note tide phase speeds are in deg/hour, so we need to convert the 
    289301      ! 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 
     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 
    295306      ENDIF 
    296307 
     
    301312      END DO 
    302313 
    303       ! 
    304314      DO itide = 1, td%ncpt 
    305315         igrd=1                              ! SSH on tracer grid. 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctra.F90

    r2797 r2865  
    5252            CALL obc_tra_frs( idx_obc(ib_obc), dta_obc(ib_obc), kt ) 
    5353         CASE DEFAULT 
    54             CALL ctl_stop( 'obc_tra : unrecognised option for open boundaries for T an S' ) 
     54            CALL ctl_stop( 'obc_tra : unrecognised option for open boundaries for T and S' ) 
    5555         END SELECT 
    5656      ENDDO 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90

    r2797 r2865  
    1010   !!---------------------------------------------------------------------- 
    1111#if   defined key_obc   &&   defined key_dynspg_flt 
    12 !!$   !!---------------------------------------------------------------------- 
    13 !!$   !!   'key_obc'            AND      unstructured open boundary conditions 
    14 !!$   !!   'key_dynspg_flt'                              filtered free surface 
    15 !!$   !!---------------------------------------------------------------------- 
    16 !!$   USE oce             ! ocean dynamics and tracers  
    17 !!$   USE dom_oce         ! ocean space and time domain  
    18 !!$   USE phycst          ! physical constants 
    19 !!$   USE obc_oce         ! ocean open boundary conditions 
    20 !!$   USE lib_mpp         ! for mppsum 
    21 !!$   USE in_out_manager  ! I/O manager 
    22 !!$   USE sbc_oce         ! ocean surface boundary conditions 
    23 !!$ 
    24 !!$   IMPLICIT NONE 
    25 !!$   PRIVATE 
    26 !!$ 
    27 !!$   PUBLIC obc_vol        ! routine called by dynspg_flt.h90 
    28 !!$ 
    29 !!$   !! * Substitutions 
    30 !!$#  include "domzgr_substitute.h90" 
    31 !!$   !!---------------------------------------------------------------------- 
    32 !!$   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    33 !!$   !! $Id$  
    34 !!$   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    35 !!$   !!---------------------------------------------------------------------- 
    36 !!$CONTAINS 
    37 !!$ 
    38 !!$   SUBROUTINE obc_vol( kt ) 
    39 !!$      !!---------------------------------------------------------------------- 
    40 !!$      !!                      ***  ROUTINE obcvol  *** 
    41 !!$      !! 
    42 !!$      !! ** Purpose :   This routine is called in dynspg_flt to control  
    43 !!$      !!      the volume of the system. A correction velocity is calculated 
    44 !!$      !!      to correct the total transport through the unstructured OBC.  
    45 !!$      !!      The total depth used is constant (H0) to be consistent with the  
    46 !!$      !!      linear free surface coded in OPA 8.2 
    47 !!$      !! 
    48 !!$      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating 
    49 !!$      !!      the total transport through all open boundaries (trans_obc) minus 
    50 !!$      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral  
    51 !!$      !!      surface (obcsurftot) of the unstructured boundary.  
    52 !!$      !!         zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot) 
    53 !!$      !!      with z_cflxemp => sum of (Evaporation minus Precipitation) 
    54 !!$      !!                       over all the domain in m3/s at each time step. 
    55 !!$      !!      z_cflxemp < 0 when precipitation dominate 
    56 !!$      !!      z_cflxemp > 0 when evaporation dominate 
    57 !!$      !! 
    58 !!$      !!      There are 2 options (user's desiderata):  
    59 !!$      !!         1/ The volume changes according to E-P, this is the default 
    60 !!$      !!            option. In this case the cumulate E-P flux are setting to 
    61 !!$      !!            zero (z_cflxemp=0) to calculate the correction velocity. So 
    62 !!$      !!            it will only balance the flux through open boundaries. 
    63 !!$      !!            (set nn_volctl to 0 in tne namelist for this option) 
    64 !!$      !!         2/ The volume is constant even with E-P flux. In this case 
    65 !!$      !!            the correction velocity must balance both the flux  
    66 !!$      !!            through open boundaries and the ones through the free 
    67 !!$      !!            surface.  
    68 !!$      !!            (set nn_volctl to 1 in tne namelist for this option) 
    69 !!$      !!---------------------------------------------------------------------- 
    70 !!$      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    71 !!$      !! 
    72 !!$      INTEGER  ::   ji, jj, jk, jb, jgrd 
    73 !!$      INTEGER  ::   ii, ij 
    74 !!$      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst 
    75 !!$      !!----------------------------------------------------------------------------- 
    76 !!$ 
    77 !!$      IF( ln_vol ) THEN 
    78 !!$ 
    79 !!$      IF( kt == nit000 ) THEN  
    80 !!$         IF(lwp) WRITE(numout,*) 
    81 !!$         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC' 
    82 !!$         IF(lwp) WRITE(numout,*)'~~~~~~~' 
    83 !!$      END IF  
    84 !!$ 
    85 !!$      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 
    86 !!$      ! ----------------------------------------------------------------------- 
    87 !!$      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0 
    88 !!$      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain 
    89 !!$ 
    90 !!$      ! Transport through the unstructured open boundary 
    91 !!$      ! ------------------------------------------------ 
    92 !!$      zubtpecor = 0.e0 
    93 !!$      jgrd = 2                               ! cumulate u component contribution first  
    94 !!$      DO jb = 1, nblenrim(jgrd) 
    95 !!$         DO jk = 1, jpkm1 
    96 !!$            ii = nbi(jb,jgrd) 
    97 !!$            ij = nbj(jb,jgrd) 
    98 !!$            zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
    99 !!$         END DO 
    100 !!$      END DO 
    101 !!$      jgrd = 3                               ! then add v component contribution 
    102 !!$      DO jb = 1, nblenrim(jgrd) 
    103 !!$         DO jk = 1, jpkm1 
    104 !!$            ii = nbi(jb,jgrd) 
    105 !!$            ij = nbj(jb,jgrd) 
    106 !!$            zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
    107 !!$         END DO 
    108 !!$      END DO 
    109 !!$      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain 
    110 !!$ 
    111 !!$      ! The normal velocity correction 
    112 !!$      ! ------------------------------ 
    113 !!$      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot  
    114 !!$      ELSE                   ;   zubtpecor =   zubtpecor             / obcsurftot 
    115 !!$      END IF 
    116 !!$ 
    117 !!$      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 
    118 !!$      ! ------------------------------------------------------------- 
    119 !!$      ztranst = 0.e0 
    120 !!$      jgrd = 2                               ! correct u component 
    121 !!$      DO jb = 1, nblenrim(jgrd) 
    122 !!$         DO jk = 1, jpkm1 
    123 !!$            ii = nbi(jb,jgrd) 
    124 !!$            ij = nbj(jb,jgrd) 
    125 !!$            ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk) 
    126 !!$            ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
    127 !!$         END DO 
    128 !!$      END DO 
    129 !!$      jgrd = 3                              ! correct v component 
    130 !!$      DO jb = 1, nblenrim(jgrd) 
    131 !!$         DO jk = 1, jpkm1 
    132 !!$            ii = nbi(jb,jgrd) 
    133 !!$            ij = nbj(jb,jgrd) 
    134 !!$            va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk) 
    135 !!$            ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
    136 !!$         END DO 
    137 !!$      END DO 
    138 !!$      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain 
    139 !!$  
    140 !!$      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 
    141 !!$      ! ------------------------------------------------------ 
    142 !!$      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
    143 !!$         IF(lwp) WRITE(numout,*) 
    144 !!$         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt 
    145 !!$         IF(lwp) WRITE(numout,*)'~~~~~~~ ' 
    146 !!$         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)' 
    147 !!$         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', obcsurftot, '(m2)' 
    148 !!$         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)' 
    149 !!$         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)' 
    150 !!$      END IF  
    151 !!$      ! 
    152 !!$      END IF ! ln_vol 
    153 !!$ 
    154 !!$   END SUBROUTINE obc_vol 
    155 !!$ 
    156 !!$#else 
     12   !!---------------------------------------------------------------------- 
     13   !!   'key_obc'            AND      unstructured open boundary conditions 
     14   !!   'key_dynspg_flt'                              filtered free surface 
     15   !!---------------------------------------------------------------------- 
     16   USE oce             ! ocean dynamics and tracers  
     17   USE dom_oce         ! ocean space and time domain  
     18   USE phycst          ! physical constants 
     19   USE obc_oce         ! ocean open boundary conditions 
     20   USE lib_mpp         ! for mppsum 
     21   USE in_out_manager  ! I/O manager 
     22   USE sbc_oce         ! ocean surface boundary conditions 
     23 
     24   IMPLICIT NONE 
     25   PRIVATE 
     26 
     27   PUBLIC obc_vol        ! routine called by dynspg_flt.h90 
     28 
     29   !! * Substitutions 
     30#  include "domzgr_substitute.h90" 
     31   !!---------------------------------------------------------------------- 
     32   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     33   !! $Id$  
     34   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     35   !!---------------------------------------------------------------------- 
     36CONTAINS 
     37 
     38   SUBROUTINE obc_vol( kt ) 
     39      !!---------------------------------------------------------------------- 
     40      !!                      ***  ROUTINE obcvol  *** 
     41      !! 
     42      !! ** Purpose :   This routine is called in dynspg_flt to control  
     43      !!      the volume of the system. A correction velocity is calculated 
     44      !!      to correct the total transport through the unstructured OBC.  
     45      !!      The total depth used is constant (H0) to be consistent with the  
     46      !!      linear free surface coded in OPA 8.2 
     47      !! 
     48      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating 
     49      !!      the total transport through all open boundaries (trans_obc) minus 
     50      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral  
     51      !!      surface (obcsurftot) of the unstructured boundary.  
     52      !!         zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot) 
     53      !!      with z_cflxemp => sum of (Evaporation minus Precipitation) 
     54      !!                       over all the domain in m3/s at each time step. 
     55      !!      z_cflxemp < 0 when precipitation dominate 
     56      !!      z_cflxemp > 0 when evaporation dominate 
     57      !! 
     58      !!      There are 2 options (user's desiderata):  
     59      !!         1/ The volume changes according to E-P, this is the default 
     60      !!            option. In this case the cumulate E-P flux are setting to 
     61      !!            zero (z_cflxemp=0) to calculate the correction velocity. So 
     62      !!            it will only balance the flux through open boundaries. 
     63      !!            (set nn_volctl to 0 in tne namelist for this option) 
     64      !!         2/ The volume is constant even with E-P flux. In this case 
     65      !!            the correction velocity must balance both the flux  
     66      !!            through open boundaries and the ones through the free 
     67      !!            surface.  
     68      !!            (set nn_volctl to 1 in tne namelist for this option) 
     69      !!---------------------------------------------------------------------- 
     70      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     71      !! 
     72      INTEGER  ::   ji, jj, jk, jb, jgrd 
     73      INTEGER  ::   ib_obc, ii, ij 
     74      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst 
     75      TYPE(OBC_INDEX), POINTER :: idx 
     76      !!----------------------------------------------------------------------------- 
     77 
     78      IF( ln_vol ) THEN 
     79 
     80      IF( kt == nit000 ) THEN  
     81         IF(lwp) WRITE(numout,*) 
     82         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC' 
     83         IF(lwp) WRITE(numout,*)'~~~~~~~' 
     84      END IF  
     85 
     86      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 
     87      ! ----------------------------------------------------------------------- 
     88      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0 
     89      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain 
     90 
     91      ! Transport through the unstructured open boundary 
     92      ! ------------------------------------------------ 
     93      zubtpecor = 0.e0 
     94      DO ib_obc = 1, nb_obc 
     95         idx => idx_obc(ib_obc) 
     96 
     97         jgrd = 2                               ! cumulate u component contribution first  
     98         DO jb = 1, idx%nblenrim(jgrd) 
     99            DO jk = 1, jpkm1 
     100               ii = idx%nbi(jb,jgrd) 
     101               ij = idx%nbj(jb,jgrd) 
     102               zubtpecor = zubtpecor + idx%flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     103            END DO 
     104         END DO 
     105         jgrd = 3                               ! then add v component contribution 
     106         DO jb = 1, idx%nblenrim(jgrd) 
     107            DO jk = 1, jpkm1 
     108               ii = idx%nbi(jb,jgrd) 
     109               ij = idx%nbj(jb,jgrd) 
     110               zubtpecor = zubtpecor + idx%flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
     111            END DO 
     112         END DO 
     113 
     114      END DO 
     115      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain 
     116 
     117      ! The normal velocity correction 
     118      ! ------------------------------ 
     119      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot  
     120      ELSE                   ;   zubtpecor =   zubtpecor             / obcsurftot 
     121      END IF 
     122 
     123      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 
     124      ! ------------------------------------------------------------- 
     125      ztranst = 0.e0 
     126      DO ib_obc = 1, nb_obc 
     127         idx => idx_obc(ib_obc) 
     128 
     129         jgrd = 2                               ! correct u component 
     130         DO jb = 1, idx%nblenrim(jgrd) 
     131            DO jk = 1, jpkm1 
     132               ii = idx%nbi(jb,jgrd) 
     133               ij = idx%nbj(jb,jgrd) 
     134               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb) * zubtpecor * umask(ii,ij,jk) 
     135               ztranst = ztranst + idx%flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     136            END DO 
     137         END DO 
     138         jgrd = 3                              ! correct v component 
     139         DO jb = 1, idx%nblenrim(jgrd) 
     140            DO jk = 1, jpkm1 
     141               ii = idx%nbi(jb,jgrd) 
     142               ij = idx%nbj(jb,jgrd) 
     143               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb) * zubtpecor * vmask(ii,ij,jk) 
     144               ztranst = ztranst + idx%flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
     145            END DO 
     146         END DO 
     147 
     148      END DO 
     149      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain 
     150  
     151      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 
     152      ! ------------------------------------------------------ 
     153      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 
     154         IF(lwp) WRITE(numout,*) 
     155         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt 
     156         IF(lwp) WRITE(numout,*)'~~~~~~~ ' 
     157         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)' 
     158         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', obcsurftot, '(m2)' 
     159         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)' 
     160         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)' 
     161      END IF  
     162      ! 
     163      END IF ! ln_vol 
     164 
     165   END SUBROUTINE obc_vol 
     166 
     167#else 
    157168   !!---------------------------------------------------------------------- 
    158169   !!   Dummy module                   NO Unstruct Open Boundary Conditions 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r2818 r2865  
    104104CONTAINS 
    105105 
    106    SUBROUTINE fld_read( kt, kn_fsbc, sd, map, jit, timeshift ) 
     106   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, jit, time_offset ) 
    107107      !!--------------------------------------------------------------------- 
    108108      !!                    ***  ROUTINE fld_read  *** 
     
    121121      TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping index 
    122122      INTEGER  , INTENT(in   ), OPTIONAL     ::   jit       ! subcycle timestep for timesplitting option 
    123       INTEGER  , INTENT(in   ), OPTIONAL     ::   timeshift ! provide fields at time other than "now" 
     123      INTEGER  , INTENT(in   ), OPTIONAL     ::   time_offset ! provide fields at time other than "now" 
     124                                                              ! time_offset = -1 => fields at "before" time level 
     125                                                              ! time_offset = +1 => fields at "after" time levels 
     126                                                              ! etc. 
    124127      !! 
    125128      INTEGER  ::   imf        ! size of the structure sd 
     
    128131      INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend 
    129132      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 
     133      INTEGER  ::   time_add   ! local time_offset variable 
    130134      LOGICAL  ::   llnxtyr    ! open next year  file? 
    131135      LOGICAL  ::   llnxtmth   ! open next month file? 
     
    143147      ENDIF 
    144148 
     149      time_add = 0 
     150      IF( PRESENT(time_offset) ) THEN 
     151         time_add = time_offset 
     152      ENDIF 
     153          
    145154      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar  
    146155      IF( present(jit) ) THEN  
    147156         ! ignore kn_fsbc in this case 
    148          isecsbc = nsec_year + nsec1jan000 + jit*rdt/REAL(nn_baro,wp) 
    149       ELSE IF( present(timeshift) ) THEN 
    150          isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + timeshift * rdttra(1)  ! middle of sbc time step 
     157         isecsbc = nsec_year + nsec1jan000 + (jit+time_add)*rdt/REAL(nn_baro,wp)  
    151158      ELSE 
    152          isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) ! middle of sbc time step 
     159         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + time_add * rdttra(1) ! middle of sbc time step 
    153160      ENDIF 
    154161      imf = SIZE( sd ) 
     
    263270                  clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   & 
    264271                     &    "', records b/a: ', i4.4, '/', i4.4, ' (days ', f7.2,'/', f7.2, ')')" 
    265                   WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   & 
     272                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &             
    266273                     & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday 
     274                  WRITE(numout, *) 'time_add is : ',time_add 
    267275               ENDIF 
    268276               ! temporal interpolation weights 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r2715 r2865  
    3535   !! free surface                                      !  before  ! now    ! after  ! 
    3636   !! ------------                                      !  fields  ! fields ! trends ! 
    37    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m] 
    38    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m] 
    39    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m] 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::            sshf_n          !: sea surface height at f-point [m] 
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m] 
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m] 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m] 
     40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::            sshf_n          !: sea surface height at f-point [m] 
    4141   ! 
    4242   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   spgu, spgv               !: horizontal surface pressure gradient 
  • branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/step.F90

    r2797 r2865  
    9797      IF( lk_dtasal  )   CALL dta_sal( kstp )         ! update 3D salinity data 
    9898                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
    99       IF( lk_obc     )   CALL obc_dta( kstp )        ! update dynamic and tracer data at open boundaries 
     99      IF( lk_obc     )   CALL obc_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries 
    100100 
    101101      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.