Changeset 3490 for branches/2012
- Timestamp:
- 2012-10-08T16:27:20+02:00 (11 years ago)
- Location:
- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r3367 r3490 28 28 INTEGER, POINTER, DIMENSION(:,:) :: nbmap 29 29 REAL , POINTER, DIMENSION(:,:) :: nbw 30 REAL , POINTER, DIMENSION(:,:) :: nbd 30 31 REAL , POINTER, DIMENSION(:) :: flagu 31 32 REAL , POINTER, DIMENSION(:) :: flagv … … 73 74 INTEGER, DIMENSION(jp_bdy) :: nn_tra_dta !: = 0 use the initial state as bdy dta ; 74 75 !: = 1 read it in a NetCDF file 75 LOGICAL, DIMENSION(jp_bdy) :: ln_tra_dmp 76 LOGICAL, DIMENSION(jp_bdy) :: ln_dyn3d_dmp 77 REAL, DIMENSION(jp_bdy) :: rn_time_dmp 78 LOGICAL :: ln_tra_dmp_tot 79 LOGICAL :: ln_dyn3d_dmp_tot 76 LOGICAL, DIMENSION(jp_bdy) :: ln_tra_dmp !: =T Tracer damping 77 LOGICAL, DIMENSION(jp_bdy) :: ln_dyn3d_dmp !: =T Baroclinic velocity damping 78 REAL, DIMENSION(jp_bdy) :: rn_time_dmp !: Damping time scale in days 80 79 81 80 #if defined key_lim2 … … 107 106 INTEGER, DIMENSION(jp_bdy) :: nn_dta !: =0 => *all* data is set to initial conditions 108 107 !: =1 => some data to be read in from data files 109 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global !: workspace for reading in global data arrays 108 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global !: workspace for reading in global data arrays (unstr. bdy) 109 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2 !: workspace for reading in global data arrays (struct. bdy) 110 110 TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET :: idx_bdy !: bdy indices (local process) 111 111 TYPE(OBC_DATA) , DIMENSION(jp_bdy) :: dta_bdy !: bdy external data (local process) -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r3367 r3490 192 192 IF( PRESENT(jit) ) THEN 193 193 ! Update barotropic boundary conditions only 194 ! jit is optional argument for fld_read and tide_update194 ! jit is optional argument for fld_read and bdytide_update 195 195 IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 196 196 IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays … … 200 200 ENDIF 201 201 IF (nn_tra(ib_bdy).ne.4) THEN 202 IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN 203 ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 204 jend = nb_bdy_fld(ib_bdy) 205 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 206 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset ) 207 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 208 ! If full velocities in boundary data then split into barotropic and baroclinic data 209 IF( ln_full_vel_array(ib_bdy) .and. & 210 & ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 211 igrd = 2 ! zonal velocity 212 dta_bdy(ib_bdy)%u2d(:) = 0.0 213 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 214 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 215 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 216 DO ik = 1, jpkm1 217 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 218 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 202 IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN 203 ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 204 jend = nb_bdy_fld(ib_bdy) 205 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 206 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset ) 207 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 208 ! If full velocities in boundary data then split into barotropic and baroclinic data 209 IF( ln_full_vel_array(ib_bdy) .and. & 210 & ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 211 igrd = 2 ! zonal velocity 212 dta_bdy(ib_bdy)%u2d(:) = 0.0 213 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 214 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 215 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 216 DO ik = 1, jpkm1 217 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 218 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 219 END DO 220 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 221 DO ik = 1, jpkm1 222 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 223 END DO 219 224 END DO 220 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 221 DO ik = 1, jpkm1 222 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 225 igrd = 3 ! meridional velocity 226 dta_bdy(ib_bdy)%v2d(:) = 0.0 227 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 228 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 229 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 230 DO ik = 1, jpkm1 231 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 232 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 233 END DO 234 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 235 DO ik = 1, jpkm1 236 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 237 END DO 223 238 END DO 224 END DO 225 igrd = 3 ! meridional velocity 226 dta_bdy(ib_bdy)%v2d(:) = 0.0 227 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 228 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 229 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 230 DO ik = 1, jpkm1 231 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 232 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 233 END DO 234 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 235 DO ik = 1, jpkm1 236 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 237 END DO 238 END DO 239 ENDIF 240 ENDIF 241 IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 242 CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), & 243 & jit=jit, time_offset=time_offset ) 244 ENDIF 239 ENDIF 240 ENDIF 241 IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 242 CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), & 243 & jit=jit, time_offset=time_offset ) 244 ENDIF 245 245 ENDIF 246 246 ENDIF … … 306 306 ENDIF 307 307 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 308 CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset )308 CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 309 309 ENDIF 310 310 ENDIF … … 378 378 ! Set nn_dta 379 379 DO ib_bdy = 1, nb_bdy 380 IF (nn_dyn2d_dta(ib_bdy).eq.1.OR.nn_dyn2d_dta(ib_bdy).eq.3) nn_dta(ib_bdy) = 1 381 IF (nn_dyn3d_dta(ib_bdy).ge.1) nn_dta(ib_bdy) = 1 382 IF (nn_tra_dta (ib_bdy).ge.1) nn_dta(ib_bdy) = 1 383 #if defined key_lim2 384 IF (nn_ice_lim2_dta(ib_bdy).eq.1) nn_dta(ib_bdy) = 1 385 #endif 386 IF (lwp) THEN 387 WRITE(numout,*) 'Bdy package number ',ib_bdy 388 IF (nn_dta(ib_bdy).eq.1) THEN 389 WRITE(numout,*) 'We use external data' 390 ELSE 391 WRITE(numout,*) 'We do not use external data' 392 ENDIF 393 ENDIF 380 nn_dta(ib_bdy) = MAX( nn_dyn2d_dta(ib_bdy) & 381 ,nn_dyn3d_dta(ib_bdy) & 382 ,nn_tra_dta(ib_bdy) & 383 #if defined key_lim2 384 ,nn_ice_lim2_dta(ib_bdy) & 385 #endif 386 ) 387 IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1 394 388 END DO 395 389 -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
r3367 r3490 14 14 !!---------------------------------------------------------------------- 15 15 USE timing ! Timing 16 USE wrk_nemo ! Memory Allocation 16 17 USE oce ! ocean dynamics and tracers 17 18 USE dom_oce ! ocean space and time domain … … 27 28 PUBLIC bdy_dyn3d_dmp ! routine called by step 28 29 30 !! * Substitutions 31 # include "domzgr_substitute.h90" 29 32 !!---------------------------------------------------------------------- 30 33 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 222 225 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp') 223 226 ! 227 !------------------------------------------------------- 228 ! Remove barotropic part from before velocity 229 !------------------------------------------------------- 230 CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 231 232 pu2d(:,:) = 0.e0 233 pv2d(:,:) = 0.e0 234 235 DO jk = 1, jpkm1 236 #if defined key_vvl 237 pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk) *umask(:,:,jk) 238 pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk) *vmask(:,:,jk) 239 #else 240 pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 241 pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 242 #endif 243 END DO 244 245 IF( lk_vvl ) THEN 246 pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 247 pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 248 ELSE 249 pu2d(:,:) = pv2d(:,:) * hur(:,:) 250 pv2d(:,:) = pu2d(:,:) * hvr(:,:) 251 ENDIF 252 224 253 DO ib_bdy=1, nb_bdy 225 254 IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN … … 228 257 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 229 258 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 230 zwgt = idx_bdy(ib_bdy)%nb w(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday)259 zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 231 260 DO jk = 1, jpkm1 232 ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) ) ) * umask(ii,ij,jk)261 ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk) 233 262 END DO 234 263 END DO … … 238 267 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 239 268 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 240 zwgt = idx_bdy(ib_bdy)%nb w(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday)269 zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 241 270 DO jk = 1, jpkm1 242 va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) ) ) * vmask(ii,ij,jk)271 va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk) 243 272 END DO 244 273 END DO 245 274 ENDIF 246 275 ENDDO 276 ! 277 CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 247 278 ! 248 279 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) ! Boundary points should be updated -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r3367 r3490 11 11 !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions 12 12 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 13 !! 3.4 ! 2012 (J. Chanut) straight open boundary case update 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_bdy … … 26 27 USE lib_mpp ! for mpp_sum 27 28 USE iom ! I/O 29 USE sbctide, ONLY: lk_tide ! Tidal forcing or not 30 USE phycst, ONLY: rday 28 31 29 32 IMPLICIT NONE … … 32 35 PUBLIC bdy_init ! routine called in nemo_init 33 36 37 INTEGER, PARAMETER :: jp_nseg = 100 38 INTEGER, PARAMETER :: nrimmax = 20 ! maximum rimwidth in structured 39 ! open boundary data files 40 ! Straight open boundary segment parameters: 41 INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs 42 INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge 43 INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw 44 INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn 45 INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs 34 46 !!---------------------------------------------------------------------- 35 47 !! NEMO/OPA 4.0 , NEMO Consortium (2011) … … 53 65 ! namelist variables 54 66 !------------------- 55 INTEGER, PARAMETER :: jp_nseg = 100 56 INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs 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 67 CHARACTER(LEN=80),DIMENSION(jpbgrd) :: clfile 68 CHARACTER(LEN=1) :: ctypebdy 69 INTEGER :: nbdyind, nbdybeg, nbdyend 61 70 62 71 ! local variables … … 66 75 INTEGER :: iw, ie, is, in, inum, id_dummy ! - - 67 76 INTEGER :: igrd_start, igrd_end, jpbdta ! - - 77 INTEGER :: jpbdtau, jpbdtas ! - - 78 INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - - 68 79 INTEGER, POINTER :: nbi, nbj, nbr ! short cuts 69 80 REAL , POINTER :: flagu, flagv ! - - 70 81 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 71 INTEGER, DIMENSION (2) :: kdimsz82 INTEGER, DIMENSION (2) :: kdimsz 72 83 INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays 73 84 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 74 85 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 75 REAL(wp), DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask 76 CHARACTER(LEN=80),DIMENSION(jpbgrd) :: clfile 77 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 86 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 78 87 !! 79 88 NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file, & … … 86 95 & ln_vol, nn_volctl, nn_rimwidth 87 96 !! 88 NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft, & 89 nbdysegw, jpiwob, jpjwdt, jpjwft, & 90 nbdysegn, jpjnob, jpindt, jpinft, & 91 nbdysegs, jpjsob, jpisdt, jpisft 97 NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 92 98 93 99 !!---------------------------------------------------------------------- … … 106 112 107 113 cgrid= (/'t','u','v'/) 108 114 109 115 ! ----------------------------------------- 110 116 ! Initialise and read namelist parameters … … 139 145 ! Check and write out namelist parameters 140 146 ! ----------------------------------------- 141 142 ln_tra_dmp_tot=.false.143 ln_dyn3d_dmp_tot=.false.144 147 ! ! control prints 145 148 IF(lwp) WRITE(numout,*) ' nambdy' … … 164 167 IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: ' 165 168 SELECT CASE( nn_dyn2d(ib_bdy) ) 166 CASE( 0 ); IF(lwp) WRITE(numout,*) ' no open boundary condition'167 CASE( 1 ); IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'168 CASE( 2) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition'169 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 170 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 171 CASE(jp_flather) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition' 169 172 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 170 173 END SELECT … … 177 180 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 178 181 END SELECT 182 IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN 183 CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' ) 184 ENDIF 179 185 ENDIF 180 186 IF(lwp) WRITE(numout,*) … … 182 188 IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: ' 183 189 SELECT CASE( nn_dyn3d(ib_bdy) ) 184 CASE( 0 ); IF(lwp) WRITE(numout,*) ' no open boundary condition'185 CASE( 1 ); IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'190 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 191 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 186 192 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 187 193 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' … … 195 201 END SELECT 196 202 ENDIF 197 IF(lwp) WRITE(numout,*)198 203 199 204 IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 200 205 IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 201 IF(lwp) WRITE(numout,*) 'No open boundary condition for the baroclinic velocitues : we force ln_dyn3d_dmp=.false.'206 IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 202 207 ln_dyn3d_dmp(ib_bdy)=.false. 208 ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN 209 CALL ctl_stop( 'Use FRS OR relaxation' ) 203 210 ELSE 204 IF(lwp) WRITE(numout,*) ' Damping of the baroclinic velocities along the boundaries'205 IF(lwp) WRITE(numout,*) ' Time damping :',rn_time_dmp(ib_bdy),' days'206 ln_dyn3d_dmp_tot=.true.211 IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone' 212 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 213 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 207 214 ENDIF 215 ELSE 216 IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities' 208 217 ENDIF 209 218 IF(lwp) WRITE(numout,*) … … 211 220 IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' 212 221 SELECT CASE( nn_tra(ib_bdy) ) 213 CASE( 0 ); IF(lwp) WRITE(numout,*) ' no open boundary condition'214 CASE( 1 ); IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'222 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 223 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 215 224 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 216 225 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Neumann conditions' … … 225 234 END SELECT 226 235 ENDIF 227 IF(lwp) WRITE(numout,*)228 236 229 237 IF ( ln_tra_dmp(ib_bdy) ) THEN 230 238 IF ( nn_tra(ib_bdy).EQ.0 ) THEN 231 IF(lwp) WRITE(numout,*) 'No open boundary condition for t he tracer : we force ln_tra_dmp=.false.'239 IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 232 240 ln_tra_dmp(ib_bdy)=.false. 241 ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN 242 CALL ctl_stop( 'Use FRS OR relaxation' ) 233 243 ELSE 234 IF(lwp) WRITE(numout,*) ' Damping of the scalars along the boundaries'235 IF(lwp) WRITE(numout,*) ' Time dumping :',rn_time_dmp(ib_bdy),' days'236 ln_tra_dmp_tot=.true.244 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone' 245 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 246 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 237 247 ENDIF 248 ELSE 249 IF(lwp) WRITE(numout,*) ' NO T/S relaxation' 238 250 ENDIF 239 251 IF(lwp) WRITE(numout,*) … … 256 268 #endif 257 269 258 IF(lwp) WRITE(numout,*) ' Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy)270 IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy) 259 271 IF(lwp) WRITE(numout,*) 260 272 261 273 ENDDO 262 274 263 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 264 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 265 IF(lwp) WRITE(numout,*) 266 SELECT CASE ( nn_volctl ) 267 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant' 268 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 269 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 270 END SELECT 271 IF(lwp) WRITE(numout,*) 272 ELSE 273 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 274 IF(lwp) WRITE(numout,*) 275 IF (nb_bdy .gt. 0) THEN 276 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 277 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 278 IF(lwp) WRITE(numout,*) 279 SELECT CASE ( nn_volctl ) 280 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant' 281 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 282 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 283 END SELECT 284 IF(lwp) WRITE(numout,*) 285 ELSE 286 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 287 IF(lwp) WRITE(numout,*) 288 ENDIF 275 289 ENDIF 276 290 … … 282 296 ! --------------------------------------------- 283 297 REWIND( numnam ) 284 jpbdta = 1298 285 299 nblendta(:,:) = 0 300 nbdysege = 0 301 nbdysegw = 0 302 nbdysegn = 0 303 nbdysegs = 0 304 icount = 0 ! count user defined segments 305 ! Dimensions below are used to allocate arrays to read external data 306 jpbdtas = 1 ! Maximum size of boundary data (structured case) 307 jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 308 286 309 DO ib_bdy = 1, nb_bdy 287 310 288 311 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 289 312 313 icount = icount + 1 290 314 ! No REWIND here because may need to read more than one nambdy_index namelist. 291 315 READ ( numnam, nambdy_index ) 292 316 293 ! Automatic boundary definition: if nbdysegX = -1294 ! set boundary to whole side of model domain.295 IF( nbdysege == -1 ) THEN296 nbdysege = 1297 jpieob(1) = jpiglo - 1298 jpjedt(1) = 2299 jpjeft(1) = jpjglo - 1300 ENDIF301 IF( nbdysegw == -1 ) THEN302 nbdysegw = 1303 jpiwob(1) = 2304 jpjwdt(1) = 2305 jpjwft(1) = jpjglo - 1306 ENDIF307 IF( nbdysegn == -1 ) THEN308 nbdysegn = 1309 jpjnob(1) = jpjglo - 1310 jpindt(1) = 2311 jpinft(1) = jpiglo - 1312 ENDIF313 IF( nbdysegs == -1 ) THEN314 nbdysegs = 1315 jpjsob(1) = 2316 jpisdt(1) = 2317 jpisft(1) = jpiglo - 1318 ENDIF319 320 DO iseg = 1, nbdysege321 igrd = 1322 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1323 igrd = 2324 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1325 igrd = 3326 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)327 ENDDO328 DO iseg = 1, nbdysegw329 igrd = 1330 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1331 igrd = 2332 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1333 igrd = 3334 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)335 ENDDO336 DO iseg = 1, nbdysegn337 igrd = 1338 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1339 igrd = 2340 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)341 igrd = 3342 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1343 END DO344 DO iseg = 1, nbdysegs 345 igrd = 1346 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1347 igrd = 2348 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg)349 igrd = 3 350 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1351 ENDDO352 353 nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)317 SELECT CASE ( TRIM(ctypebdy) ) 318 CASE( 'N' ) 319 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 320 nbdyind = jpjglo - 2 ! set boundary to whole side of model domain. 321 nbdybeg = 2 322 nbdyend = jpiglo - 1 323 ENDIF 324 nbdysegn = nbdysegn + 1 325 npckgn(nbdysegn) = ib_bdy ! Save bdy package number 326 jpjnob(nbdysegn) = nbdyind 327 jpindt(nbdysegn) = nbdybeg 328 jpinft(nbdysegn) = nbdyend 329 ! 330 CASE( 'S' ) 331 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 332 nbdyind = 2 ! set boundary to whole side of model domain. 333 nbdybeg = 2 334 nbdyend = jpiglo - 1 335 ENDIF 336 nbdysegs = nbdysegs + 1 337 npckgs(nbdysegs) = ib_bdy ! Save bdy package number 338 jpjsob(nbdysegs) = nbdyind 339 jpisdt(nbdysegs) = nbdybeg 340 jpisft(nbdysegs) = nbdyend 341 ! 342 CASE( 'E' ) 343 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 344 nbdyind = jpiglo - 2 ! set boundary to whole side of model domain. 345 nbdybeg = 2 346 nbdyend = jpjglo - 1 347 ENDIF 348 nbdysege = nbdysege + 1 349 npckge(nbdysege) = ib_bdy ! Save bdy package number 350 jpieob(nbdysege) = nbdyind 351 jpjedt(nbdysege) = nbdybeg 352 jpjeft(nbdysege) = nbdyend 353 ! 354 CASE( 'W' ) 355 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 356 nbdyind = 2 ! set boundary to whole side of model domain. 357 nbdybeg = 2 358 nbdyend = jpjglo - 1 359 ENDIF 360 nbdysegw = nbdysegw + 1 361 npckgw(nbdysegw) = ib_bdy ! Save bdy package number 362 jpiwob(nbdysegw) = nbdyind 363 jpjwdt(nbdysegw) = nbdybeg 364 jpjwft(nbdysegw) = nbdyend 365 ! 366 CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 367 END SELECT 368 369 ! For simplicity we assume that in case of straight bdy, arrays have the same length 370 ! (even if it is true that last tangential velocity points 371 ! are useless). This simplifies a little bit boundary data format (and agrees with format 372 ! used so far in obc package) 373 374 nblendta(1:jpbgrd,ib_bdy) = (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 375 jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 376 IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 377 & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 354 378 355 379 ELSE ! Read size of arrays in boundary coordinates file. 356 357 358 380 CALL iom_open( cn_coords_file(ib_bdy), inum ) 359 381 DO igrd = 1, jpbgrd 360 382 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 361 383 nblendta(igrd,ib_bdy) = kdimsz(1) 362 ENDDO 384 jpbdtau = MAX(jpbdtau, kdimsz(1)) 385 ENDDO 386 CALL iom_close( inum ) 363 387 364 388 ENDIF … … 366 390 ENDDO ! ib_bdy 367 391 368 jpbdta = MAXVAL(nblendta(:,:)) 369 370 ! Allocate arrays 371 !--------------- 372 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), & 373 & nbrdta(jpbdta, jpbgrd, nb_bdy) ) 374 375 ALLOCATE( dta_global(jpbdta, 1, jpk) ) 392 IF (nb_bdy>0) THEN 393 jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 394 395 ! Allocate arrays 396 !--------------- 397 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), & 398 & nbrdta(jpbdta, jpbgrd, nb_bdy) ) 399 400 ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 401 IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 402 ! 403 ENDIF 404 405 ! Now look for crossings in user (namelist) defined open boundary segments: 406 !-------------------------------------------------------------------------- 407 IF ( icount>0 ) CALL bdy_ctl_seg 376 408 377 409 ! Calculate global boundary index arrays or read in from file 378 !------------------------------------------------------------ 379 REWIND( numnam )410 !------------------------------------------------------------ 411 ! 1. Read global index arrays from boundary coordinates file. 380 412 DO ib_bdy = 1, nb_bdy 381 413 382 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 383 384 ! No REWIND here because may need to read more than one nambdy_index namelist. 385 READ ( numnam, nambdy_index ) 386 387 ! Automatic boundary definition: if nbdysegX = -1 388 ! set boundary to whole side of model domain. 389 IF( nbdysege == -1 ) THEN 390 nbdysege = 1 391 jpieob(1) = jpiglo - 1 392 jpjedt(1) = 2 393 jpjeft(1) = jpjglo - 1 394 ENDIF 395 IF( nbdysegw == -1 ) THEN 396 nbdysegw = 1 397 jpiwob(1) = 2 398 jpjwdt(1) = 2 399 jpjwft(1) = jpjglo - 1 400 ENDIF 401 IF( nbdysegn == -1 ) THEN 402 nbdysegn = 1 403 jpjnob(1) = jpjglo - 1 404 jpindt(1) = 2 405 jpinft(1) = jpiglo - 1 406 ENDIF 407 IF( nbdysegs == -1 ) THEN 408 nbdysegs = 1 409 jpjsob(1) = 2 410 jpisdt(1) = 2 411 jpisft(1) = jpiglo - 1 412 ENDIF 413 414 ! ------------ T points ------------- 415 igrd = 1 416 icount = 0 417 DO ir = 1, nn_rimwidth(ib_bdy) 418 ! east 419 DO iseg = 1, nbdysege 420 DO ij = jpjedt(iseg), jpjeft(iseg) 421 icount = icount + 1 422 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 423 nbjdta(icount, igrd, ib_bdy) = ij 424 nbrdta(icount, igrd, ib_bdy) = ir 425 ENDDO 426 ENDDO 427 ! west 428 DO iseg = 1, nbdysegw 429 DO ij = jpjwdt(iseg), jpjwft(iseg) 430 icount = icount + 1 431 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 432 nbjdta(icount, igrd, ib_bdy) = ij 433 nbrdta(icount, igrd, ib_bdy) = ir 434 ENDDO 435 ENDDO 436 ! north 437 DO iseg = 1, nbdysegn 438 DO ii = jpindt(iseg), jpinft(iseg) 439 icount = icount + 1 440 nbidta(icount, igrd, ib_bdy) = ii 441 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 442 nbrdta(icount, igrd, ib_bdy) = ir 443 ENDDO 444 ENDDO 445 ! south 446 DO iseg = 1, nbdysegs 447 DO ii = jpisdt(iseg), jpisft(iseg) 448 icount = icount + 1 449 nbidta(icount, igrd, ib_bdy) = ii 450 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 451 nbrdta(icount, igrd, ib_bdy) = ir 452 ENDDO 453 ENDDO 454 ENDDO 455 456 ! ------------ U points ------------- 457 igrd = 2 458 icount = 0 459 DO ir = 1, nn_rimwidth(ib_bdy) 460 ! east 461 DO iseg = 1, nbdysege 462 DO ij = jpjedt(iseg), jpjeft(iseg) 463 icount = icount + 1 464 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 465 nbjdta(icount, igrd, ib_bdy) = ij 466 nbrdta(icount, igrd, ib_bdy) = ir 467 ENDDO 468 ENDDO 469 ! west 470 DO iseg = 1, nbdysegw 471 DO ij = jpjwdt(iseg), jpjwft(iseg) 472 icount = icount + 1 473 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 474 nbjdta(icount, igrd, ib_bdy) = ij 475 nbrdta(icount, igrd, ib_bdy) = ir 476 ENDDO 477 ENDDO 478 ! north 479 DO iseg = 1, nbdysegn 480 DO ii = jpindt(iseg), jpinft(iseg) - 1 481 icount = icount + 1 482 nbidta(icount, igrd, ib_bdy) = ii 483 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 484 nbrdta(icount, igrd, ib_bdy) = ir 485 ENDDO 486 ENDDO 487 ! south 488 DO iseg = 1, nbdysegs 489 DO ii = jpisdt(iseg), jpisft(iseg) - 1 490 icount = icount + 1 491 nbidta(icount, igrd, ib_bdy) = ii 492 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 493 nbrdta(icount, igrd, ib_bdy) = ir 494 ENDDO 495 ENDDO 496 ENDDO 497 498 ! ------------ V points ------------- 499 igrd = 3 500 icount = 0 501 DO ir = 1, nn_rimwidth(ib_bdy) 502 ! east 503 DO iseg = 1, nbdysege 504 DO ij = jpjedt(iseg), jpjeft(iseg) - 1 505 icount = icount + 1 506 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 507 nbjdta(icount, igrd, ib_bdy) = ij 508 nbrdta(icount, igrd, ib_bdy) = ir 509 ENDDO 510 ENDDO 511 ! west 512 DO iseg = 1, nbdysegw 513 DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 514 icount = icount + 1 515 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 516 nbjdta(icount, igrd, ib_bdy) = ij 517 nbrdta(icount, igrd, ib_bdy) = ir 518 ENDDO 519 ENDDO 520 ! north 521 DO iseg = 1, nbdysegn 522 DO ii = jpindt(iseg), jpinft(iseg) 523 icount = icount + 1 524 nbidta(icount, igrd, ib_bdy) = ii 525 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 526 nbrdta(icount, igrd, ib_bdy) = ir 527 ENDDO 528 ENDDO 529 ! south 530 DO iseg = 1, nbdysegs 531 DO ii = jpisdt(iseg), jpisft(iseg) 532 icount = icount + 1 533 nbidta(icount, igrd, ib_bdy) = ii 534 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 535 nbrdta(icount, igrd, ib_bdy) = ir 536 ENDDO 537 ENDDO 538 ENDDO 539 540 ELSE ! Read global index arrays from boundary coordinates file. 541 414 IF( ln_coords_file(ib_bdy) ) THEN 415 416 CALL iom_open( cn_coords_file(ib_bdy), inum ) 542 417 DO igrd = 1, jpbgrd 543 418 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) … … 560 435 IF (ibr_max < nn_rimwidth(ib_bdy)) & 561 436 CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 562 563 437 END DO 564 438 CALL iom_close( inum ) … … 566 440 ENDIF 567 441 568 ENDDO 442 ENDDO 443 444 ! 2. Now fill indices corresponding to straight open boundary arrays: 445 ! East 446 !----- 447 DO iseg = 1, nbdysege 448 ib_bdy = npckge(iseg) 449 ! 450 ! ------------ T points ------------- 451 igrd=1 452 icount=0 453 DO ir = 1, nn_rimwidth(ib_bdy) 454 DO ij = jpjedt(iseg), jpjeft(iseg) 455 icount = icount + 1 456 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 457 nbjdta(icount, igrd, ib_bdy) = ij 458 nbrdta(icount, igrd, ib_bdy) = ir 459 ENDDO 460 ENDDO 461 ! 462 ! ------------ U points ------------- 463 igrd=2 464 icount=0 465 DO ir = 1, nn_rimwidth(ib_bdy) 466 DO ij = jpjedt(iseg), jpjeft(iseg) 467 icount = icount + 1 468 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 469 nbjdta(icount, igrd, ib_bdy) = ij 470 nbrdta(icount, igrd, ib_bdy) = ir 471 ENDDO 472 ENDDO 473 ! 474 ! ------------ V points ------------- 475 igrd=3 476 icount=0 477 DO ir = 1, nn_rimwidth(ib_bdy) 478 ! DO ij = jpjedt(iseg), jpjeft(iseg) - 1 479 DO ij = jpjedt(iseg), jpjeft(iseg) 480 icount = icount + 1 481 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 482 nbjdta(icount, igrd, ib_bdy) = ij 483 nbrdta(icount, igrd, ib_bdy) = ir 484 ENDDO 485 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 486 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 487 ENDDO 488 ENDDO 489 ! 490 ! West 491 !----- 492 DO iseg = 1, nbdysegw 493 ib_bdy = npckgw(iseg) 494 ! 495 ! ------------ T points ------------- 496 igrd=1 497 icount=0 498 DO ir = 1, nn_rimwidth(ib_bdy) 499 DO ij = jpjwdt(iseg), jpjwft(iseg) 500 icount = icount + 1 501 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 502 nbjdta(icount, igrd, ib_bdy) = ij 503 nbrdta(icount, igrd, ib_bdy) = ir 504 ENDDO 505 ENDDO 506 ! 507 ! ------------ U points ------------- 508 igrd=2 509 icount=0 510 DO ir = 1, nn_rimwidth(ib_bdy) 511 DO ij = jpjwdt(iseg), jpjwft(iseg) 512 icount = icount + 1 513 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 514 nbjdta(icount, igrd, ib_bdy) = ij 515 nbrdta(icount, igrd, ib_bdy) = ir 516 ENDDO 517 ENDDO 518 ! 519 ! ------------ V points ------------- 520 igrd=3 521 icount=0 522 DO ir = 1, nn_rimwidth(ib_bdy) 523 ! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 524 DO ij = jpjwdt(iseg), jpjwft(iseg) 525 icount = icount + 1 526 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 527 nbjdta(icount, igrd, ib_bdy) = ij 528 nbrdta(icount, igrd, ib_bdy) = ir 529 ENDDO 530 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 531 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 532 ENDDO 533 ENDDO 534 ! 535 ! North 536 !----- 537 DO iseg = 1, nbdysegn 538 ib_bdy = npckgn(iseg) 539 ! 540 ! ------------ T points ------------- 541 igrd=1 542 icount=0 543 DO ir = 1, nn_rimwidth(ib_bdy) 544 DO ii = jpindt(iseg), jpinft(iseg) 545 icount = icount + 1 546 nbidta(icount, igrd, ib_bdy) = ii 547 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 548 nbrdta(icount, igrd, ib_bdy) = ir 549 ENDDO 550 ENDDO 551 ! 552 ! ------------ U points ------------- 553 igrd=2 554 icount=0 555 DO ir = 1, nn_rimwidth(ib_bdy) 556 ! DO ii = jpindt(iseg), jpinft(iseg) - 1 557 DO ii = jpindt(iseg), jpinft(iseg) 558 icount = icount + 1 559 nbidta(icount, igrd, ib_bdy) = ii 560 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 561 nbrdta(icount, igrd, ib_bdy) = ir 562 ENDDO 563 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 564 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 565 ENDDO 566 ! 567 ! ------------ V points ------------- 568 igrd=3 569 icount=0 570 DO ir = 1, nn_rimwidth(ib_bdy) 571 DO ii = jpindt(iseg), jpinft(iseg) 572 icount = icount + 1 573 nbidta(icount, igrd, ib_bdy) = ii 574 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 575 nbrdta(icount, igrd, ib_bdy) = ir 576 ENDDO 577 ENDDO 578 ENDDO 579 ! 580 ! South 581 !----- 582 DO iseg = 1, nbdysegs 583 ib_bdy = npckgs(iseg) 584 ! 585 ! ------------ T points ------------- 586 igrd=1 587 icount=0 588 DO ir = 1, nn_rimwidth(ib_bdy) 589 DO ii = jpisdt(iseg), jpisft(iseg) 590 icount = icount + 1 591 nbidta(icount, igrd, ib_bdy) = ii 592 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 593 nbrdta(icount, igrd, ib_bdy) = ir 594 ENDDO 595 ENDDO 596 ! 597 ! ------------ U points ------------- 598 igrd=2 599 icount=0 600 DO ir = 1, nn_rimwidth(ib_bdy) 601 ! DO ii = jpisdt(iseg), jpisft(iseg) - 1 602 DO ii = jpisdt(iseg), jpisft(iseg) 603 icount = icount + 1 604 nbidta(icount, igrd, ib_bdy) = ii 605 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 606 nbrdta(icount, igrd, ib_bdy) = ir 607 ENDDO 608 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 609 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 610 ENDDO 611 ! 612 ! ------------ V points ------------- 613 igrd=3 614 icount=0 615 DO ir = 1, nn_rimwidth(ib_bdy) 616 DO ii = jpisdt(iseg), jpisft(iseg) 617 icount = icount + 1 618 nbidta(icount, igrd, ib_bdy) = ii 619 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 620 nbrdta(icount, igrd, ib_bdy) = ir 621 ENDDO 622 ENDDO 623 ENDDO 624 625 ! Deal with duplicated points 626 !----------------------------- 627 ! We assign negative indices to duplicated points (to remove them from bdy points to be updated) 628 ! if their distance to the bdy is greater than the other 629 ! If their distance are the same, just keep only one to avoid updating a point twice 630 DO igrd = 1, jpbgrd 631 DO ib_bdy1 = 1, nb_bdy 632 DO ib_bdy2 = 1, nb_bdy 633 IF (ib_bdy1/=ib_bdy2) THEN 634 DO ib1 = 1, nblendta(igrd,ib_bdy1) 635 DO ib2 = 1, nblendta(igrd,ib_bdy2) 636 IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 637 & (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 638 ! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', & 639 ! & nbidta(ib1, igrd, ib_bdy1), & 640 ! & nbjdta(ib2, igrd, ib_bdy2) 641 ! keep only points with the lowest distance to boundary: 642 IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 643 nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 644 nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 645 ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 646 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 647 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 648 ! Arbitrary choice if distances are the same: 649 ELSE 650 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 651 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 652 ENDIF 653 END IF 654 END DO 655 END DO 656 ENDIF 657 END DO 658 END DO 659 END DO 569 660 570 661 ! Work out dimensions of boundary data on each processor 571 662 ! ------------------------------------------------------ 572 573 iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 574 ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 575 is = mjg(1) + 1 ! if monotasking and no zoom, is=2 576 in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 663 664 ! Rather assume that boundary data indices are given on global domain 665 ! TO BE DISCUSSED ? 666 ! iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 667 ! ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 668 ! is = mjg(1) + 1 ! if monotasking and no zoom, is=2 669 ! in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 670 iw = mig(1) - jpizoom + 2 ! if monotasking and no zoom, iw=2 671 ie = mig(1) + nlci - jpizoom - 1 ! if monotasking and no zoom, ie=jpim1 672 is = mjg(1) - jpjzoom + 2 ! if monotasking and no zoom, is=2 673 in = mjg(1) + nlcj - jpjzoom - 1 ! if monotasking and no zoom, in=jpjm1 577 674 578 675 DO ib_bdy = 1, nb_bdy … … 610 707 ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 611 708 ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 709 ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 612 710 ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 613 711 ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) … … 629 727 ! 630 728 icount = icount + 1 631 idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1 632 idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 729 730 ! Rather assume that boundary data indices are given on global domain 731 ! TO BE DISCUSSED ? 732 ! idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1 733 ! idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 734 idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom 735 idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom 633 736 idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy) 634 737 idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib … … 643 746 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 644 747 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 645 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation 646 idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic 647 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear 748 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation 749 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic 750 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)) ! linear 751 END DO 752 END DO 753 754 ! Compute damping coefficients 755 ! ---------------------------- 756 DO igrd = 1, jpbgrd 757 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 758 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 759 idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 760 & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic 648 761 END DO 649 762 END DO … … 660 773 ! = 0 elsewhere 661 774 662 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN ! EEL configuration at 5km resolution 663 zmask( : ,:) = 0.e0 664 zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0 665 ELSE IF( ln_mask_file ) THEN 775 IF( ln_mask_file ) THEN 666 776 CALL iom_open( cn_mask_file, inum ) 667 CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )777 CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 668 778 CALL iom_close( inum ) 669 ELSE 670 zmask(:,:) = 1.e0 671 ENDIF 672 673 DO ij = 1, nlcj ! Save mask over local domain 674 DO ii = 1, nlci 675 bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) ) 779 780 ! Derive mask on U and V grid from mask on T grid 781 bdyumask(:,:) = 0.e0 782 bdyvmask(:,:) = 0.e0 783 DO ij=1, jpjm1 784 DO ii=1, jpim1 785 bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 786 bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1) 787 END DO 676 788 END DO 677 END DO 678 679 ! Derive mask on U and V grid from mask on T grid 680 bdyumask(:,:) = 0.e0 681 bdyvmask(:,:) = 0.e0 682 DO ij=1, jpjm1 683 DO ii=1, jpim1 684 bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 685 bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1) 789 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond. 790 791 792 ! Mask corrections 793 ! ---------------- 794 DO ik = 1, jpkm1 795 DO ij = 1, jpj 796 DO ii = 1, jpi 797 tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 798 umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 799 vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 800 bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij) 801 END DO 802 END DO 686 803 END DO 687 END DO 688 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond. 689 690 691 ! Mask corrections 692 ! ---------------- 693 DO ik = 1, jpkm1 694 DO ij = 1, jpj 695 DO ii = 1, jpi 696 tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 697 umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 698 vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 699 bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij) 700 END DO 804 805 DO ik = 1, jpkm1 806 DO ij = 2, jpjm1 807 DO ii = 2, jpim1 808 fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) & 809 & * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 810 END DO 811 END DO 701 812 END DO 702 END DO 703 704 DO ik = 1, jpkm1 705 DO ij = 2, jpjm1 706 DO ii = 2, jpim1 707 fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) & 708 & * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 709 END DO 710 END DO 711 END DO 712 713 tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 813 814 tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 815 816 ENDIF ! ln_mask_file=.TRUE. 817 714 818 bdytmask(:,:) = tmask(:,:,1) 715 819 … … 800 904 ! Compute total lateral surface for volume correction: 801 905 ! ---------------------------------------------------- 906 ! JC: this must be done at each time step with key_vvl 802 907 bdysurftot = 0.e0 803 908 IF( ln_vol ) THEN … … 806 911 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 807 912 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 808 nbj => idx_bdy(ib_bdy)%nb i(ib,igrd)913 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 809 914 flagu => idx_bdy(ib_bdy)%flagu(ib) 810 915 bdysurftot = bdysurftot + hu (nbi , nbj) & … … 819 924 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 820 925 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 821 nbj => idx_bdy(ib_bdy)%nb i(ib,igrd)926 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 822 927 flagv => idx_bdy(ib_bdy)%flagv(ib) 823 928 bdysurftot = bdysurftot + hv (nbi, nbj ) & … … 833 938 ! Tidy up 834 939 !-------- 835 DEALLOCATE(nbidta, nbjdta, nbrdta) 940 IF (nb_bdy>0) THEN 941 DEALLOCATE(nbidta, nbjdta, nbrdta) 942 ENDIF 836 943 837 944 IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 838 945 839 946 END SUBROUTINE bdy_init 947 948 SUBROUTINE bdy_ctl_seg 949 !!---------------------------------------------------------------------- 950 !! *** ROUTINE bdy_ctl_seg *** 951 !! 952 !! ** Purpose : Check straight open boundary segments location 953 !! 954 !! ** Method : - Look for open boundary corners 955 !! - Check that segments start or end on land 956 !!---------------------------------------------------------------------- 957 INTEGER :: ib, ib1, ib2, ji ,jj, itest 958 INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 959 REAL(wp), DIMENSION(2) :: ztestmask 960 !!---------------------------------------------------------------------- 961 ! 962 IF (lwp) WRITE(numout,*) ' ' 963 IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments' 964 IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~' 965 ! 966 IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege 967 IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw 968 IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn 969 IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs 970 ! 1. Check bounds 971 !---------------- 972 DO ib = 1, nbdysegn 973 IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib) 974 IF ((jpjnob(ib).ge.jpjglo-1).or.& 975 &(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 976 IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 977 IF (jpindt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' ) 978 IF (jpinft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' ) 979 END DO 980 ! 981 DO ib = 1, nbdysegs 982 IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib) 983 IF ((jpjsob(ib).ge.jpjglo-1).or.& 984 &(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 985 IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 986 IF (jpisdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' ) 987 IF (jpisft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' ) 988 END DO 989 ! 990 DO ib = 1, nbdysege 991 IF (lwp) WRITE(numout,*) '**check east seg bounds pckg: ', npckge(ib) 992 IF ((jpieob(ib).ge.jpiglo-1).or.& 993 &(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 994 IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 995 IF (jpjedt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' ) 996 IF (jpjeft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' ) 997 END DO 998 ! 999 DO ib = 1, nbdysegw 1000 IF (lwp) WRITE(numout,*) '**check west seg bounds pckg: ', npckgw(ib) 1001 IF ((jpiwob(ib).ge.jpiglo-1).or.& 1002 &(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 1003 IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 1004 IF (jpjwdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' ) 1005 IF (jpjwft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' ) 1006 ENDDO 1007 ! 1008 ! 1009 ! 2. Look for segment crossings 1010 !------------------------------ 1011 IF (lwp) WRITE(numout,*) '**Look for segments corners :' 1012 ! 1013 itest = 0 ! corner number 1014 ! 1015 ! flag to detect if start or end of open boundary belongs to a corner 1016 ! if not (=0), it must be on land. 1017 ! if a corner is detected, save bdy package number for further tests 1018 icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0. 1019 ! South/West crossings 1020 IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN 1021 DO ib1 = 1, nbdysegw 1022 DO ib2 = 1, nbdysegs 1023 IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. & 1024 & ( jpisft(ib2)>=jpiwob(ib1)).AND. & 1025 & ( jpjwdt(ib1)<=jpjsob(ib2)).AND. & 1026 & ( jpjwft(ib1)>=jpjsob(ib2))) THEN 1027 IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 1028 ! We have a possible South-West corner 1029 ! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1) 1030 ! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2) 1031 icornw(ib1,1) = npckgs(ib2) 1032 icorns(ib2,1) = npckgw(ib1) 1033 ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 1034 IF(lwp) WRITE(numout,*) 1035 IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 1036 & jpisft(ib2), jpjwft(ib1) 1037 IF(lwp) WRITE(numout,*) ' ========== Not allowed yet' 1038 IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1039 & ' and South segment: ',npckgs(ib2) 1040 IF(lwp) WRITE(numout,*) 1041 nstop = nstop + 1 1042 ELSE 1043 IF(lwp) WRITE(numout,*) 1044 IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices' 1045 IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1) , & 1046 & ' and South segment: ',npckgs(ib2) 1047 IF(lwp) WRITE(numout,*) 1048 nstop = nstop+1 1049 END IF 1050 END IF 1051 END DO 1052 END DO 1053 END IF 1054 ! 1055 ! South/East crossings 1056 IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN 1057 DO ib1 = 1, nbdysege 1058 DO ib2 = 1, nbdysegs 1059 IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. & 1060 & ( jpisft(ib2)>=jpieob(ib1)+1).AND. & 1061 & ( jpjedt(ib1)<=jpjsob(ib2) ).AND. & 1062 & ( jpjeft(ib1)>=jpjsob(ib2) )) THEN 1063 IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN 1064 ! We have a possible South-East corner 1065 ! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1) 1066 ! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2) 1067 icorne(ib1,1) = npckgs(ib2) 1068 icorns(ib2,2) = npckge(ib1) 1069 ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 1070 IF(lwp) WRITE(numout,*) 1071 IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 1072 & jpisdt(ib2), jpjeft(ib1) 1073 IF(lwp) WRITE(numout,*) ' ========== Not allowed yet' 1074 IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), & 1075 & ' and South segment: ',npckgs(ib2) 1076 IF(lwp) WRITE(numout,*) 1077 nstop = nstop + 1 1078 ELSE 1079 IF(lwp) WRITE(numout,*) 1080 IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices' 1081 IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), & 1082 & ' and South segment: ',npckgs(ib2) 1083 IF(lwp) WRITE(numout,*) 1084 nstop = nstop + 1 1085 END IF 1086 END IF 1087 END DO 1088 END DO 1089 END IF 1090 ! 1091 ! North/West crossings 1092 IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN 1093 DO ib1 = 1, nbdysegw 1094 DO ib2 = 1, nbdysegn 1095 IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. & 1096 & ( jpinft(ib2)>=jpiwob(ib1) ).AND. & 1097 & ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. & 1098 & ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN 1099 IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN 1100 ! We have a possible North-West corner 1101 ! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1) 1102 ! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2) 1103 icornw(ib1,2) = npckgn(ib2) 1104 icornn(ib2,1) = npckgw(ib1) 1105 ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 1106 IF(lwp) WRITE(numout,*) 1107 IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 1108 & jpinft(ib2), jpjwdt(ib1) 1109 IF(lwp) WRITE(numout,*) ' ========== Not allowed yet' 1110 IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1111 & ' and North segment: ',npckgn(ib2) 1112 IF(lwp) WRITE(numout,*) 1113 nstop = nstop + 1 1114 ELSE 1115 IF(lwp) WRITE(numout,*) 1116 IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices' 1117 IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1), & 1118 & ' and North segment: ',npckgn(ib2) 1119 IF(lwp) WRITE(numout,*) 1120 nstop = nstop + 1 1121 END IF 1122 END IF 1123 END DO 1124 END DO 1125 END IF 1126 ! 1127 ! North/East crossings 1128 IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN 1129 DO ib1 = 1, nbdysege 1130 DO ib2 = 1, nbdysegn 1131 IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. & 1132 & ( jpinft(ib2)>=jpieob(ib1)+1).AND. & 1133 & ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. & 1134 & ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN 1135 IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN 1136 ! We have a possible North-East corner 1137 ! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1) 1138 ! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2) 1139 icorne(ib1,2) = npckgn(ib2) 1140 icornn(ib2,2) = npckge(ib1) 1141 ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 1142 IF(lwp) WRITE(numout,*) 1143 IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 1144 & jpindt(ib2), jpjedt(ib1) 1145 IF(lwp) WRITE(numout,*) ' ========== Not allowed yet' 1146 IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), & 1147 & ' and North segment: ',npckgn(ib2) 1148 IF(lwp) WRITE(numout,*) 1149 nstop = nstop + 1 1150 ELSE 1151 IF(lwp) WRITE(numout,*) 1152 IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices' 1153 IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), & 1154 & ' and North segment: ',npckgn(ib2) 1155 IF(lwp) WRITE(numout,*) 1156 nstop = nstop + 1 1157 END IF 1158 END IF 1159 END DO 1160 END DO 1161 END IF 1162 ! 1163 ! 3. Check if segment extremities are on land 1164 !-------------------------------------------- 1165 ! 1166 ! West segments 1167 DO ib = 1, nbdysegw 1168 ! get mask at boundary extremities: 1169 ztestmask(1:2)=0. 1170 DO ji = 1, jpi 1171 DO jj = 1, jpj 1172 IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 1173 & ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1) 1174 IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 1175 & ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 1176 END DO 1177 END DO 1178 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 1179 1180 IF (ztestmask(1)==1) THEN 1181 IF (icornw(ib,1)==0) THEN 1182 IF(lwp) WRITE(numout,*) 1183 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 1184 IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner' 1185 IF(lwp) WRITE(numout,*) 1186 nstop = nstop + 1 1187 ELSE 1188 ! This is a corner 1189 WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib) 1190 CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1)) 1191 itest=itest+1 1192 ENDIF 1193 ENDIF 1194 IF (ztestmask(2)==1) THEN 1195 IF (icornw(ib,2)==0) THEN 1196 IF(lwp) WRITE(numout,*) 1197 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 1198 IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner' 1199 IF(lwp) WRITE(numout,*) 1200 nstop = nstop + 1 1201 ELSE 1202 ! This is a corner 1203 WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib) 1204 CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2)) 1205 itest=itest+1 1206 ENDIF 1207 ENDIF 1208 END DO 1209 ! 1210 ! East segments 1211 DO ib = 1, nbdysege 1212 ! get mask at boundary extremities: 1213 ztestmask(1:2)=0. 1214 DO ji = 1, jpi 1215 DO jj = 1, jpj 1216 IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 1217 & ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1) 1218 IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 1219 & ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 1220 END DO 1221 END DO 1222 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 1223 1224 IF (ztestmask(1)==1) THEN 1225 IF (icorne(ib,1)==0) THEN 1226 IF(lwp) WRITE(numout,*) 1227 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 1228 IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner' 1229 IF(lwp) WRITE(numout,*) 1230 nstop = nstop + 1 1231 ELSE 1232 ! This is a corner 1233 WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib) 1234 CALL bdy_ctl_corn(npckge(ib), icorne(ib,1)) 1235 itest=itest+1 1236 ENDIF 1237 ENDIF 1238 IF (ztestmask(2)==1) THEN 1239 IF (icorne(ib,2)==0) THEN 1240 IF(lwp) WRITE(numout,*) 1241 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 1242 IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner' 1243 IF(lwp) WRITE(numout,*) 1244 nstop = nstop + 1 1245 ELSE 1246 ! This is a corner 1247 WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib) 1248 CALL bdy_ctl_corn(npckge(ib), icorne(ib,2)) 1249 itest=itest+1 1250 ENDIF 1251 ENDIF 1252 END DO 1253 ! 1254 ! South segments 1255 DO ib = 1, nbdysegs 1256 ! get mask at boundary extremities: 1257 ztestmask(1:2)=0. 1258 DO ji = 1, jpi 1259 DO jj = 1, jpj 1260 IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 1261 & ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1) 1262 IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 1263 & ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 1264 END DO 1265 END DO 1266 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 1267 1268 IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 1269 IF(lwp) WRITE(numout,*) 1270 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 1271 IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner' 1272 IF(lwp) WRITE(numout,*) 1273 nstop = nstop + 1 1274 ENDIF 1275 IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 1276 IF(lwp) WRITE(numout,*) 1277 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 1278 IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner' 1279 IF(lwp) WRITE(numout,*) 1280 nstop = nstop + 1 1281 ENDIF 1282 END DO 1283 ! 1284 ! North segments 1285 DO ib = 1, nbdysegn 1286 ! get mask at boundary extremities: 1287 ztestmask(1:2)=0. 1288 DO ji = 1, jpi 1289 DO jj = 1, jpj 1290 IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 1291 & ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1) 1292 IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 1293 & ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 1294 END DO 1295 END DO 1296 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 1297 1298 IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 1299 IF(lwp) WRITE(numout,*) 1300 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 1301 IF(lwp) WRITE(numout,*) ' ========== does not start on land' 1302 IF(lwp) WRITE(numout,*) 1303 nstop = nstop + 1 1304 ENDIF 1305 IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 1306 IF(lwp) WRITE(numout,*) 1307 IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 1308 IF(lwp) WRITE(numout,*) ' ========== does not end on land' 1309 IF(lwp) WRITE(numout,*) 1310 nstop = nstop + 1 1311 ENDIF 1312 END DO 1313 ! 1314 IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found' 1315 ! 1316 ! Other tests TBD: 1317 ! segments completly on land 1318 ! optimized open boundary array length according to landmask 1319 ! Nudging layers that overlap with interior domain 1320 ! 1321 END SUBROUTINE bdy_ctl_seg 1322 1323 SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 1324 !!---------------------------------------------------------------------- 1325 !! *** ROUTINE bdy_ctl_corn *** 1326 !! 1327 !! ** Purpose : Check numerical schemes consistency between 1328 !! segments having a common corner 1329 !! 1330 !! ** Method : 1331 !!---------------------------------------------------------------------- 1332 INTEGER, INTENT(in) :: ib1, ib2 1333 INTEGER :: itest 1334 !!---------------------------------------------------------------------- 1335 itest = 0 1336 1337 IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1 1338 IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1 1339 IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1 1340 ! 1341 IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1 1342 IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1 1343 IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1 1344 ! 1345 IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1 1346 ! 1347 IF ( itest>0 ) THEN 1348 IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2 1349 IF(lwp) WRITE(numout,*) ' ========== have different open bdy schemes' 1350 IF(lwp) WRITE(numout,*) 1351 nstop = nstop + 1 1352 ENDIF 1353 ! 1354 END SUBROUTINE bdy_ctl_corn 840 1355 841 1356 #else -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r3367 r3490 8 8 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 9 9 !! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes 10 !! 3.4 ! 2012-09 (G. Reffray and J. Chanut) New inputs + mods 10 11 !!---------------------------------------------------------------------- 11 12 #if defined key_bdy … … 27 28 USE bdy_oce ! ocean open boundary conditions 28 29 USE daymod ! calendar 30 USE wrk_nemo ! Memory allocation 29 31 USE tideini 30 USE tide_mod 32 ! USE tide_mod ! Useless ?? 31 33 USE fldread, ONLY: fld_map 32 34 … … 34 36 PRIVATE 35 37 36 PUBLIC bdytide_init ! routine called in nemo_init37 PUBLIC tide_update ! routine called in bdydyn38 PUBLIC bdytide_init ! routine called in bdy_init 39 PUBLIC bdytide_update ! routine called in bdy_dta 38 40 39 41 TYPE, PUBLIC :: TIDES_DATA !: Storage for external tidal harmonics data … … 47 49 48 50 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data 49 51 50 52 !!---------------------------------------------------------------------- 51 53 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 66 68 !!------------------- 67 69 CHARACTER(len=80) :: filtide !: Filename root for tidal input files 68 LOGICAL :: ln_conjug = .FALSE. !: F/T : tidal data in complex/complex conjugate 70 LOGICAL :: ln_bdytide_2ddta !: If true, read 2d harmonic data 71 LOGICAL :: ln_bdytide_conj !: If true, assume complex conjugate tidal data 69 72 !! 70 73 INTEGER :: ib_bdy, itide, ib !: dummy loop indices 74 INTEGER :: ii, ij !: dummy loop indices 71 75 INTEGER :: inum, igrd 72 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 76 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 77 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts 73 78 CHARACTER(len=80) :: clfile !: full file name for tidal input file 74 79 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data 80 REAL(wp), POINTER, DIMENSION(:,:) :: ztr, zti !: " " " " " " " " 75 81 !! 76 82 TYPE(TIDES_DATA), POINTER :: td !: local short cut 77 83 !! 78 NAMELIST/nambdy_tide/filtide, ln_ conjug84 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 79 85 !!---------------------------------------------------------------------- 80 86 81 87 IF( nn_timing == 1 ) CALL timing_start('bdytide_init') 82 88 83 IF(lwp) WRITE(numout,*) 84 IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 85 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 89 IF (nb_bdy>0) THEN 90 IF(lwp) WRITE(numout,*) 91 IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 92 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 93 ENDIF 94 95 ln_bdytide_2ddta = .FALSE. 96 ln_bdytide_conj = .FALSE. 86 97 87 98 REWIND(numnam) … … 90 101 91 102 td => tides(ib_bdy) 103 nblen => idx_bdy(ib_bdy)%nblen 104 nblenrim => idx_bdy(ib_bdy)%nblenrim 92 105 93 106 ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries … … 97 110 READ ( numnam, nambdy_tide ) 98 111 ! ! Parameter control and print 112 IF(lwp) WRITE(numout,*) ' ' 99 113 IF(lwp) WRITE(numout,*) ' Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 100 IF(lwp) WRITE(numout,*) ' tidal components specified ', nb_harmo 101 IF(lwp) WRITE(numout,*) ' ', Wave(ntide(1:nb_harmo))%cname_tide 102 IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : ' 103 IF(lwp) WRITE(numout,*) ' ', omega_tide(1:nb_harmo) 114 IF(lwp) WRITE(numout,*) ' read tidal data in 2d files: ', ln_bdytide_2ddta 115 IF(lwp) WRITE(numout,*) ' assume complex conjugate : ', ln_bdytide_conj 116 IF(lwp) WRITE(numout,*) ' Number of tidal components to read: ', nb_harmo 117 IF(lwp) THEN 118 WRITE(numout,*) ' Tidal cpt name - Phase speed (deg/hr)' 119 DO itide = 1, nb_harmo 120 WRITE(numout,*) ' ', Wave(ntide(itide))%cname_tide, omega_tide(itide) 121 END DO 122 ENDIF 123 IF(lwp) WRITE(numout,*) ' ' 104 124 105 125 ! Allocate space for tidal harmonics data - get size from OBC data arrays 106 126 ! ----------------------------------------------------------------------- 107 127 108 ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 128 ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 129 ! relaxation area 130 IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 131 ilen0(:)=nblen(:) 132 ELSE 133 ilen0(:)=nblenrim(:) 134 ENDIF 135 109 136 ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 110 137 ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 111 138 112 ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d )113 139 ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 114 140 ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 115 141 116 ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d )117 142 ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 118 143 ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 119 144 120 ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 121 122 ! Open files and read in tidal forcing data 123 ! ----------------------------------------- 124 125 DO itide = 1, nb_harmo 126 ! ! SSH fields 127 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 128 CALL iom_open( clfile, inum ) 129 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 130 td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 131 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 132 td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 145 td%ssh0(:,:,:) = 0.e0 146 td%ssh(:,:,:) = 0.e0 147 td%u0(:,:,:) = 0.e0 148 td%u(:,:,:) = 0.e0 149 td%v0(:,:,:) = 0.e0 150 td%v(:,:,:) = 0.e0 151 152 IF (ln_bdytide_2ddta) THEN 153 ! It is assumed that each data file contains all complex harmonic amplitudes 154 ! given on the data domain (ie global, jpidta x jpjdta) 155 ! 156 CALL wrk_alloc( jpi, jpj, zti, ztr ) 157 ! 158 ! SSH fields 159 clfile = TRIM(filtide)//'_grid_T.nc' 160 CALL iom_open (clfile , inum ) 161 igrd = 1 ! Everything is at T-points here 162 DO itide = 1, nb_harmo 163 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) ) 164 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 165 DO ib = 1, ilen0(igrd) 166 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 167 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 168 td%ssh0(ib,itide,1) = ztr(ii,ij) 169 td%ssh0(ib,itide,2) = zti(ii,ij) 170 END DO 171 END DO 133 172 CALL iom_close( inum ) 134 ! ! U fields 135 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 136 CALL iom_open( clfile, inum ) 137 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 138 td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 139 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 140 td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 173 ! 174 ! U fields 175 clfile = TRIM(filtide)//'_grid_U.nc' 176 CALL iom_open (clfile , inum ) 177 igrd = 2 ! Everything is at U-points here 178 DO itide = 1, nb_harmo 179 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) ) 180 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) ) 181 DO ib = 1, ilen0(igrd) 182 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 183 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 184 td%u0(ib,itide,1) = ztr(ii,ij) 185 td%u0(ib,itide,2) = zti(ii,ij) 186 END DO 187 END DO 141 188 CALL iom_close( inum ) 142 ! ! V fields 143 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 144 CALL iom_open( clfile, inum ) 145 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 146 td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 147 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 148 td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 189 ! 190 ! V fields 191 clfile = TRIM(filtide)//'_grid_V.nc' 192 CALL iom_open (clfile , inum ) 193 igrd = 3 ! Everything is at V-points here 194 DO itide = 1, nb_harmo 195 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) ) 196 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) ) 197 DO ib = 1, ilen0(igrd) 198 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 199 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 200 td%v0(ib,itide,1) = ztr(ii,ij) 201 td%v0(ib,itide,2) = zti(ii,ij) 202 END DO 203 END DO 149 204 CALL iom_close( inum ) 150 IF ( ln_conjug ) THEN 151 IF(lwp) WRITE(numout,*) ' The tidal input data are written in complex conjugate' 152 td%ssh0(:,itide,2) = - td%ssh0(:,itide,2) 153 td%u0 (:,itide,2) = - td%u0 (:,itide,2) 154 td%v0 (:,itide,2) = - td%v0 (:,itide,2) 155 ENDIF 156 ! 157 END DO ! end loop on tidal components 205 ! 206 CALL wrk_dealloc( jpi, jpj, ztr, zti ) 207 ! 208 ELSE 209 ! 210 ! Read tidal data only on bdy segments 211 ! 212 ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 213 214 ! Open files and read in tidal forcing data 215 ! ----------------------------------------- 216 217 DO itide = 1, nb_harmo 218 ! ! SSH fields 219 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 220 CALL iom_open( clfile, inum ) 221 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 222 td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 223 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 224 td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 225 CALL iom_close( inum ) 226 ! ! U fields 227 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 228 CALL iom_open( clfile, inum ) 229 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 230 td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 231 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 232 td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 233 CALL iom_close( inum ) 234 ! ! V fields 235 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 236 CALL iom_open( clfile, inum ) 237 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 238 td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 239 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 240 td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 241 CALL iom_close( inum ) 242 ! 243 END DO ! end loop on tidal components 244 ! 245 DEALLOCATE( dta_read ) 246 ENDIF ! ln_bdytide_2ddta=.true. 158 247 ! 159 DEALLOCATE( dta_read ) 248 IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files 249 td%ssh0(:,:,2) = - td%ssh0(:,:,2) 250 td%u0 (:,:,2) = - td%u0 (:,:,2) 251 td%v0 (:,:,2) = - td%v0 (:,:,2) 252 ENDIF 160 253 ! 161 254 ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 … … 167 260 END SUBROUTINE bdytide_init 168 261 169 SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset )170 !!---------------------------------------------------------------------- 171 !! *** SUBROUTINE tide_update ***262 SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset ) 263 !!---------------------------------------------------------------------- 264 !! *** SUBROUTINE bdytide_update *** 172 265 !! 173 266 !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. … … 186 279 ! etc. 187 280 !! 281 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 188 282 INTEGER :: itide, igrd, ib ! dummy loop indices 189 283 INTEGER :: time_add ! time offset in units of timesteps 190 REAL(wp) :: z_arg, z_sarg, zflag 284 REAL(wp) :: z_arg, z_sarg, zflag, zramp 191 285 REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost 192 286 !!---------------------------------------------------------------------- 193 287 194 IF( nn_timing == 1 ) CALL timing_start('tide_update') 288 IF( nn_timing == 1 ) CALL timing_start('bdytide_update') 289 290 ilen0(1) = SIZE(td%ssh(:,1,1)) 291 ilen0(2) = SIZE(td%u(:,1,1)) 292 ilen0(3) = SIZE(td%v(:,1,1)) 195 293 196 294 zflag=1 … … 205 303 IF(lwp) THEN 206 304 WRITE(numout,*) 207 WRITE(numout,*) 'bdy _tide : (re)Initialization of the tidal bdy forcing at kt=',kt208 WRITE(numout,*) '~~~~~~~ '305 WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 306 WRITE(numout,*) '~~~~~~~~~~~~~~ ' 209 307 ENDIF 210 308 ! … … 225 323 ENDIF 226 324 325 ! Linear ramp on tidal component at open boundaries 326 zramp = 1. 327 IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 328 227 329 DO itide = 1, nb_harmo 228 330 z_sarg = z_arg * omega_tide(itide) … … 233 335 DO itide = 1, nb_harmo 234 336 igrd=1 ! SSH on tracer grid 235 DO ib = 1, i dx%nblenrim(igrd)236 dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)337 DO ib = 1, ilen0(igrd) 338 dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) 237 339 END DO 238 340 igrd=2 ! U grid 239 DO ib =1, idx%nblenrim(igrd)240 dta%u2d(ib) = dta%u2d(ib) + td%u (ib,itide,1)*z_cost(itide) + td%u (ib,itide,2)*z_sist(itide)341 DO ib = 1, ilen0(igrd) 342 dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u (ib,itide,1)*z_cost(itide) + td%u (ib,itide,2)*z_sist(itide)) 241 343 END DO 242 344 igrd=3 ! V grid 243 DO ib =1, idx%nblenrim(igrd)244 dta%v2d(ib) = dta%v2d(ib) + td%v (ib,itide,1)*z_cost(itide) + td%v (ib,itide,2)*z_sist(itide)345 DO ib = 1, ilen0(igrd) 346 dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v (ib,itide,1)*z_cost(itide) + td%v (ib,itide,2)*z_sist(itide)) 245 347 END DO 246 348 END DO 247 349 ! 248 IF( nn_timing == 1 ) CALL timing_stop(' tide_update')350 IF( nn_timing == 1 ) CALL timing_stop('bdytide_update') 249 351 ! 250 END SUBROUTINE tide_update352 END SUBROUTINE bdytide_update 251 353 252 354 SUBROUTINE tide_init_elevation( idx, td ) … … 257 359 TYPE(TIDES_DATA),INTENT( inout ) :: td ! tidal harmonics data 258 360 !! * Local declarations 361 INTEGER, DIMENSION(1) :: ilen0 !: length of boundary data (from OBC arrays) 259 362 REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide 260 363 INTEGER :: itide, igrd, ib ! dummy loop indices 261 364 262 igrd=1 ! SSH on tracer grid. 263 264 ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd))) 365 igrd=1 366 ! SSH on tracer grid. 367 368 ilen0(1) = SIZE(td%ssh0(:,1,1)) 369 370 ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 265 371 266 372 DO itide = 1, nb_harmo 267 DO ib = 1, i dx%nblenrim(igrd)373 DO ib = 1, ilen0(igrd) 268 374 mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.) 269 375 phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 270 376 END DO 271 DO ib = 1 , idx%nblenrim(igrd)377 DO ib = 1 , ilen0(igrd) 272 378 mod_tide(ib)=mod_tide(ib)*ftide(itide) 273 379 phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 274 380 ENDDO 275 DO ib = 1 , idx%nblenrim(igrd)381 DO ib = 1 , ilen0(igrd) 276 382 td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 277 383 td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) … … 290 396 TYPE(TIDES_DATA),INTENT( inout ) :: td ! tidal harmonics data 291 397 !! * Local declarations 398 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 292 399 REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide 293 400 INTEGER :: itide, igrd, ib ! dummy loop indices 294 401 402 ilen0(2) = SIZE(td%u0(:,1,1)) 403 ilen0(3) = SIZE(td%v0(:,1,1)) 404 295 405 igrd=2 ! U grid. 296 406 297 ALLOCATE(mod_tide(i dx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))407 ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 298 408 299 409 DO itide = 1, nb_harmo 300 DO ib = 1, i dx%nblenrim(igrd)410 DO ib = 1, ilen0(igrd) 301 411 mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.) 302 412 phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 303 413 END DO 304 DO ib = 1, i dx%nblenrim(igrd)414 DO ib = 1, ilen0(igrd) 305 415 mod_tide(ib)=mod_tide(ib)*ftide(itide) 306 416 phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 307 417 ENDDO 308 DO ib = 1, i dx%nblenrim(igrd)418 DO ib = 1, ilen0(igrd) 309 419 td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 310 420 td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) … … 314 424 DEALLOCATE(mod_tide,phi_tide) 315 425 316 igrd=3 ! Ugrid.317 318 ALLOCATE(mod_tide(i dx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))426 igrd=3 ! V grid. 427 428 ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 319 429 320 430 DO itide = 1, nb_harmo 321 DO ib = 1, i dx%nblenrim(igrd)431 DO ib = 1, ilen0(igrd) 322 432 mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.) 323 433 phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 324 434 END DO 325 DO ib = 1, i dx%nblenrim(igrd)435 DO ib = 1, ilen0(igrd) 326 436 mod_tide(ib)=mod_tide(ib)*ftide(itide) 327 437 phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 328 438 ENDDO 329 DO ib = 1, i dx%nblenrim(igrd)439 DO ib = 1, ilen0(igrd) 330 440 td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 331 441 td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) … … 341 451 !!---------------------------------------------------------------------- 342 452 CONTAINS 343 SUBROUTINE bdytide_init ! Empty routine 453 SUBROUTINE bdytide_init ! Empty routine 454 WRITE(*,*) 'bdytide_init: You should not have seen this print! error?' 344 455 END SUBROUTINE bdytide_init 345 SUBROUTINE tide_data ! Empty routine 346 END SUBROUTINE tide_data 347 SUBROUTINE tide_update( kt, jit ) ! Empty routine 348 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit 349 END SUBROUTINE tide_update 456 SUBROUTINE bdytide_update( kt, jit ) ! Empty routine 457 WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 458 END SUBROUTINE bdytide_update 350 459 #endif 351 460 -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
r3367 r3490 22 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 23 USE in_out_manager ! I/O manager 24 USE phycst 24 25 25 26 26 IMPLICIT NONE … … 28 28 29 29 PUBLIC bdy_tra ! routine called in tranxt.F90 30 PUBLIC bdy_tra_dmp ! routine called in tranxt.F9030 PUBLIC bdy_tra_dmp ! routine called in step.F90 31 31 32 32 !!---------------------------------------------------------------------- … … 65 65 END SELECT 66 66 ENDDO 67 ! 68 ! Boundary points should be updated 69 IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) 70 IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 67 71 68 72 END SUBROUTINE bdy_tra … … 98 102 END DO 99 103 ! 100 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated101 !102 104 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 103 105 ! … … 133 135 END DO 134 136 END DO 135 !136 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated137 137 ! 138 138 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) … … 190 190 END DO 191 191 ! 192 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated193 !194 192 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 195 193 ! … … 229 227 END DO 230 228 END DO 231 !232 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated233 229 ! 234 230 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) … … 262 258 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 263 259 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 264 zwgt = idx_bdy(ib_bdy)%nb w(ib,igrd) / ( rn_time_dmp(ib_bdy) * rday)260 zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd) 265 261 DO ik = 1, jpkm1 266 tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik) 267 tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik) 262 zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik) 263 zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik) 264 tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta 265 tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa 268 266 END DO 269 267 END DO 270 268 ENDIF 271 269 ENDDO 272 !273 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated274 270 ! 275 271 IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp') -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
r3367 r3490 26 26 ftide 27 27 28 LOGICAL, PUBLIC :: ln_tide_pot = .false. 28 LOGICAL, PUBLIC :: ln_tide_pot = .false., ln_tide_ramp = .false. 29 REAL(wp), PUBLIC :: rdttideramp 29 30 INTEGER, PUBLIC :: nb_harmo 30 31 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide … … 46 47 CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 47 48 ! 48 NAMELIST/nam_tide/ln_tide_pot, clname49 NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname 49 50 !!---------------------------------------------------------------------- 50 51 … … 54 55 WRITE(numout,*) 55 56 WRITE(numout,*) 'tide_init : Initialization of the tidal components' 56 WRITE(numout,*) '~~~~~~~ '57 WRITE(numout,*) '~~~~~~~~~ ' 57 58 ENDIF 58 59 ! … … 77 78 WRITE(numout,*) ' Namelist nam_tide' 78 79 WRITE(numout,*) ' nb_harmo = ', nb_harmo 80 WRITE(numout,*) ' ln_tide_ramp = ', ln_tide_ramp 81 WRITE(numout,*) ' rdttideramp = ', rdttideramp 82 IF (ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp)) & 83 & CALL ctl_stop('rdttideramp must be lower than run duration') 84 IF (ln_tide_ramp.AND.(rdttideramp<0.)) & 85 & CALL ctl_stop('rdttideramp must be positive') 79 86 CALL flush(numout) 80 87 ENDIF -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90
r3294 r3490 13 13 USE sbctide 14 14 USE dynspg_oce 15 USE tideini, ONLY: ln_tide_ramp, rdttideramp 15 16 16 17 IMPLICIT NONE … … 33 34 INTEGER, INTENT( in ) :: kt,kit ! ocean time-step index 34 35 INTEGER :: ji,jj,jk 36 REAL (wp) :: zramp 35 37 REAL (wp), DIMENSION(nb_harmo) :: zwt 36 !...............................................................................37 ! Potentiel astronomique38 38 !............................................................................... 39 39 40 40 pot_astro(:,:)=0.e0 41 zramp = 1.e0 41 42 42 43 IF (lk_dynspg_ts) THEN 43 44 zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp))) 45 IF (ln_tide_ramp) THEN 46 zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.) 47 ENDIF 44 48 ELSE 45 49 zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt 50 IF (ln_tide_ramp) THEN 51 zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.) 52 ENDIF 46 53 ENDIF 47 54 … … 49 56 do ji=1,jpi 50 57 do jj=1,jpj 51 pot_astro(ji,jj)=pot_astro(ji,jj) + (amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))58 pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk))) 52 59 enddo 53 60 enddo -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r3367 r3490 76 76 USE mod_ioclient 77 77 #endif 78 USE sbctide, ONLY: lk_tide 79 78 80 79 81 IMPLICIT NONE … … 318 320 CALL istate_init ! ocean initial state (Dynamics and tracers) 319 321 320 322 IF( lk_tide ) CALL tide_init( nit000 ) ! Initialisation of the tidal harmonics 321 323 322 324 IF( lk_bdy ) CALL bdy_init ! Open boundaries initialisation -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step.F90
r3367 r3490 95 95 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 96 96 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 97 IF( kstp /= nit000) CALL tide_init ( kstp )97 IF( lk_tide.AND.(kstp /= nit000 )) CALL tide_init ( kstp ) 98 98 IF( lk_tide ) CALL sbc_tide( kstp ) 99 99 IF( lk_obc ) CALL obc_dta( kstp ) ! update dynamic and tracer data at open boundaries … … 188 188 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 189 189 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 190 IF( lk_bdy .AND. & 191 & ln_tra_dmp_tot) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 190 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 192 191 CALL tra_adv ( kstp ) ! horizontal & vertical advection 193 192 IF( lk_zdfkpp ) CALL tra_kpp ( kstp ) ! KPP non-local tracer fluxes … … 222 221 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 223 222 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 224 IF( lk_bdy .AND. & 225 & ln_dyn3d_dmp_tot) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 223 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 226 224 CALL dyn_adv( kstp ) ! advection (vector or flux form) 227 225 CALL dyn_vor( kstp ) ! vorticity term including Coriolis
Note: See TracChangeset
for help on using the changeset viewer.