Changeset 1125 for trunk/NEMO/OPA_SRC/BDY/bdyini.F90
- Timestamp:
- 2008-06-23T11:05:02+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/BDY/bdyini.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/BDY/bdyini.F90
r911 r1125 1 2 MODULE bdyini 3 !!================================================================================= 1 MODULE bdyini 2 !!====================================================================== 4 3 !! *** MODULE bdyini *** 5 !! Initialization of unstructured open boundaries 6 !!================================================================================= 7 #if defined key_bdy || defined key_bdy_tides 8 !!--------------------------------------------------------------------------------- 9 !! 'key_bdy' Unstructured Open Boundary Conditions 10 !!--------------------------------------------------------------------------------- 4 !! Unstructured open boundaries : initialisation 5 !!====================================================================== 6 !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code 7 !! - ! 2007-01 (D. Storkey) Update to use IOM module 8 !! - ! 2007-01 (D. Storkey) Tidal forcing 9 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 10 !!---------------------------------------------------------------------- 11 #if defined key_bdy 12 !!---------------------------------------------------------------------- 13 !! 'key_bdy' Unstructured Open Boundary Conditions 14 !!---------------------------------------------------------------------- 11 15 !! bdy_init : Initialization of unstructured open boundaries 12 !!--------------------------------------------------------------------------------- 13 !! * Modules used 16 !!---------------------------------------------------------------------- 14 17 USE oce ! ocean dynamics and tracers variables 15 18 USE dom_oce ! ocean space and time domain … … 24 27 PRIVATE 25 28 26 !! * Routine accessibility 27 PUBLIC bdy_init ! routine called by opa.F90 28 29 !! * Substitutions 30 31 !!--------------------------------------------------------------------------------- 32 !! OPA 9.0 , LODYC-IPSL (2003) 29 PUBLIC bdy_init ! routine called by opa.F90 30 31 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 33 !! $Id: $ 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 33 35 !!--------------------------------------------------------------------------------- 34 36 … … 47 49 !! ** Input : bdy_init.nc, input file for unstructured open boundaries 48 50 !! 49 !! History :50 !! OPA 9.0 ! 05-01 (J. Chanut, A. Sellar) Original code51 !! ! 07-01 (D. Storkey) Update to use IOM module.52 !! ! 07-01 (D. Storkey) Tidal forcing.53 51 !!---------------------------------------------------------------------- 54 !! * Local declarations 55 INTEGER :: ji, jj, jk, jgrd, & ! dummy loop indices 56 jb, jr, icount, & 57 icountr, nb_rim, nb_len, nbr_max 52 INTEGER :: ii, ij, ik, igrd, ib, ir ! dummy loop indices 53 INTEGER :: icount, icountr 54 INTEGER :: ib_len, ibr_max 58 55 INTEGER :: iw, ie, is, in 59 56 INTEGER :: inum ! temporary logical unit 60 INTEGER :: & ! temporary integers 61 dummy_id 62 INTEGER, DIMENSION (2) :: kdimsz 63 INTEGER, DIMENSION(jpbdta, jpbgrd) :: & !: Index arrays 64 nbidta, nbjdta, & !: i and j indices of bdy dta 65 nbrdta !: Discrete distance from rim points 66 REAL(wp) :: & 67 efl, wfl, nfl, sfl ! temporary scalars 68 REAL(wp) , DIMENSION(jpidta,jpjdta) :: & 69 tmpmsk ! global domain mask 70 REAL(wp) , DIMENSION(jpbdta,1) :: & 71 ndta ! temporary array 72 CHARACTER(LEN=80),DIMENSION(3) :: bdyfile 73 74 NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, & 75 ln_bdy_clim, ln_bdy_vol, ln_bdy_fla, ln_bdy_mask, & 76 nbdy_dta, nb_rimwidth, volbdy 77 57 INTEGER :: id_dummy ! temporary integers 58 INTEGER :: igrd_start, igrd_end ! start and end of loops on igrd 59 INTEGER, DIMENSION (2) :: kdimsz 60 INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 61 INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbrdta ! Discrete distance from rim points 62 REAL(wp) :: zefl, zwfl, znfl, zsfl ! temporary scalars 63 REAL(wp) , DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask 64 REAL(wp) , DIMENSION(jpbdta,1) :: zdta ! temporary array 65 CHARACTER(LEN=80),DIMENSION(3) :: clfile 66 !! 67 NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, & 68 & ln_bdy_tides, ln_bdy_clim, ln_bdy_vol, ln_bdy_mask, & 69 & ln_bdy_dyn_fla, ln_bdy_dyn_frs, ln_bdy_tra_frs, & 70 & nbdy_dta , nb_rimwidth , volbdy 78 71 !!---------------------------------------------------------------------- 79 72 … … 81 74 IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries' 82 75 IF(lwp) WRITE(numout,*) '~~~~~~~~' 83 84 IF( jperio /= 0 ) THEN 85 IF(lwp) WRITE(numout,*) 86 IF(lwp) WRITE(numout,*) ' E R R O R : Cyclic or symmetric,', & 87 ' and unstructured open boundary condition are not compatible' 88 IF(lwp) WRITE(numout,*) ' ========== ' 89 IF(lwp) WRITE(numout,*) 90 nstop = nstop + 1 91 END IF 76 ! 77 IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', & 78 ' and unstructured open boundary condition are not compatible' ) 92 79 93 80 #if defined key_obc 94 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) ' E R R O R : Straight open boundaries,', & 96 ' and unstructured open boundaries are not compatible' 97 IF(lwp) WRITE(numout,*) ' ========== ' 98 IF(lwp) WRITE(numout,*) 99 nstop = nstop + 1 81 CALL ctl_stop( 'Straight open boundaries,', & 82 ' and unstructured open boundaries are not compatible' ) 100 83 #endif 101 84 102 85 # if defined key_dynspg_rl 103 IF(lwp) WRITE(numout,*) 104 IF(lwp) WRITE(numout,*) ' E R R O R : Rigid lid,', & 105 ' and unstructured open boundaries are not compatible' 106 IF(lwp) WRITE(numout,*) ' ========== ' 107 IF(lwp) WRITE(numout,*) 108 nstop = nstop + 1 86 CALL ctl_stop( 'Rigid lid,', & 87 ' and unstructured open boundaries are not compatible' ) 109 88 #endif 110 89 111 ! 0.Read namelist parameters90 ! Read namelist parameters 112 91 ! --------------------------- 113 114 92 REWIND( numnam ) 115 93 READ ( numnam, nambdy ) … … 118 96 IF(lwp) WRITE(numout,*) ' nambdy' 119 97 120 IF ((nbdy_dta/=0).AND.(nbdy_dta/=1)) THEN ! Check nbdy_dta value 121 IF(lwp) WRITE(numout,*) 122 IF(lwp) WRITE(numout,*) ' E R R O R : nbdy_dta =',nbdy_dta, & 123 'but it should have been 0 or 1' 124 IF(lwp) WRITE(numout,*) ' ========== ' 125 IF(lwp) WRITE(numout,*) 126 nstop = nstop + 1 127 ELSE 128 IF(lwp) WRITE(numout,*) ' ' 129 IF(lwp) WRITE(numout,*) ' data in file (=1) or nbdy_dta = ', nbdy_dta 130 IF(lwp) WRITE(numout,*) ' initial state used (=0)' 131 END IF 98 ! Check nbdy_dta value 99 IF(lwp) WRITE(numout,*) 'nbdy_dta =', nbdy_dta 100 IF(lwp) WRITE(numout,*) ' ' 101 SELECT CASE( nbdy_dta ) 102 CASE( 0 ) 103 IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 104 CASE( 1 ) 105 IF(lwp) WRITE(numout,*) ' boundary data taken from file' 106 CASE DEFAULT 107 CALL ctl_stop( 'nbdy_dta must be 0 or 1' ) 108 END SELECT 132 109 133 110 IF(lwp) WRITE(numout,*) ' ' 134 111 IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth 135 112 113 IF(lwp) WRITE(numout,*) ' ' 114 IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy 115 136 116 IF (ln_bdy_vol) THEN 137 IF (volbdy==1) THEN ! Check volbdy value 138 IF(lwp) WRITE(numout,*) ' ' 139 IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy 117 SELECT CASE ( volbdy ) ! Check volbdy value 118 CASE( 1 ) 140 119 IF(lwp) WRITE(numout,*) ' The total volume will be constant' 141 142 ELSEIF (volbdy==0) THEN 143 IF(lwp) WRITE(numout,*) ' ' 144 IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy 145 IF(lwp) WRITE(numout,*) ' The total volume will vary according to & 146 &the surface E-P flux' 147 ELSE 148 IF(lwp) WRITE(numout,*) ' ' 149 IF(lwp) WRITE(numout,*) ' E R R O R : volbdy =',volbdy, & 150 'but it should have been 0 or 1' 151 IF(lwp) WRITE(numout,*) ' ========== ' 152 IF(lwp) WRITE(numout,*) 153 nstop = nstop + 1 154 END IF 120 CASE( 0 ) 121 IF(lwp) WRITE(numout,*) ' The total volume will vary according' 122 IF(lwp) WRITE(numout,*) ' to the surface E-P flux' 123 CASE DEFAULT 124 CALL ctl_stop( 'volbdy must be 0 or 1' ) 125 END SELECT 155 126 ELSE 156 IF(lwp) WRITE(numout,*) ' '157 127 IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries' 158 128 IF(lwp) WRITE(numout,*) ' ' 159 129 ENDIF 160 130 161 IF (ln_bdy_fla) THEN 162 IF(lwp) WRITE(numout,*) ' ' 163 IF(lwp) WRITE(numout,*) 'Flather bc with unstructured open boundaries' 164 IF(lwp) WRITE(numout,*) ' ' 165 ELSE 166 IF(lwp) WRITE(numout,*) ' ' 167 IF(lwp) WRITE(numout,*) 'NO Flather bc with unstructured open boundaries' 168 IF(lwp) WRITE(numout,*) ' ' 169 ENDIF 170 171 ! 0.5 Read tides namelist 131 IF (ln_bdy_tides) THEN 132 IF(lwp) WRITE(numout,*) ' ' 133 IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 134 IF(lwp) WRITE(numout,*) ' ' 135 ENDIF 136 137 IF (ln_bdy_dyn_fla) THEN 138 IF(lwp) WRITE(numout,*) ' ' 139 IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries' 140 IF(lwp) WRITE(numout,*) ' ' 141 ENDIF 142 143 IF (ln_bdy_dyn_frs) THEN 144 IF(lwp) WRITE(numout,*) ' ' 145 IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries' 146 IF(lwp) WRITE(numout,*) ' ' 147 ENDIF 148 149 IF (ln_bdy_tra_frs) THEN 150 IF(lwp) WRITE(numout,*) ' ' 151 IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries' 152 IF(lwp) WRITE(numout,*) ' ' 153 ENDIF 154 155 ! Read tides namelist 172 156 ! ------------------------ 173 174 IF ( lk_bdy_tides ) CALL tide_init 175 176 ! 1. Read arrays defining unstructured open boundaries 177 ! ---------------------------------------------------- 178 179 ! 1.1 Read global 2D mask at T-points: bdytmask 180 ! ********************************************* 181 ! bdytmask=1 on the computational domain AND on open boundaries 182 ! =0 elsewhere 157 IF( ln_bdy_tides ) CALL tide_init 158 159 ! Read arrays defining unstructured open boundaries 160 ! ------------------------------------------------- 161 162 ! Read global 2D mask at T-points: bdytmask 163 ! ***************************************** 164 ! bdytmask = 1 on the computational domain AND on open boundaries 165 ! = 0 elsewhere 183 166 184 167 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 185 tmpmsk(: , :) = 0.e0186 tmpmsk(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0168 zmask( : ,:) = 0.e0 169 zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0 187 170 ELSE IF ( ln_bdy_mask ) THEN 188 CALL iom_open( filbdy_mask, inum )189 CALL iom_get ( inum, jpdom_data, 'bdy_msk', tmpmsk(:,:) )190 CALL iom_close( inum )171 CALL iom_open( filbdy_mask, inum ) 172 CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 173 CALL iom_close( inum ) 191 174 ELSE 192 tmpmsk(:,:) = 1.0175 zmask(:,:) = 1.e0 193 176 ENDIF 194 177 195 178 ! Save mask over local domain 196 DO jj = 1, nlcj197 DO ji = 1, nlci198 bdytmask(ji,jj) = tmpmsk( mig(ji), mjg(jj))199 END DO179 DO ij = 1, nlcj 180 DO ii = 1, nlci 181 bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) ) 182 END DO 200 183 END DO 201 184 202 185 ! Derive mask on U and V grid from mask on T grid 203 bdyumask(:,:) =0.e0204 bdyvmask(:,:) =0.e0205 DO jj=1, jpjm1206 DO ji=1, jpim1207 bdyumask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji+1, jj )208 bdyvmask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji ,jj+1)209 END DO186 bdyumask(:,:) = 0.e0 187 bdyvmask(:,:) = 0.e0 188 DO ij=1, jpjm1 189 DO ii=1, jpim1 190 bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 191 bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1) 192 END DO 210 193 END DO 211 194 … … 214 197 CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 215 198 216 ! 1.2Read discrete distance and mapping indices217 ! ****************************************** ****218 nbidta(:,:)=0.219 nbjdta(:,:)=0.220 nbrdta(:,:)=0.199 ! Read discrete distance and mapping indices 200 ! ****************************************** 201 nbidta(:,:) = 0.e0 202 nbjdta(:,:) = 0.e0 203 nbrdta(:,:) = 0.e0 221 204 222 205 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 223 224 icount = 0 225 ! Define west boundary (from ji=2 to ji=1+nb_rimwidth): 226 DO jr=1,nb_rimwidth 227 DO jj=3,jpjglo-2 228 icount=icount+1 229 nbidta(icount,:) = jr + 1 + (jpizoom-1) 230 nbjdta(icount,:) = jj + (jpjzoom-1) 231 nbrdta(icount,:) = jr 232 END DO 233 END DO 234 235 ! Define east boundary (from ji=jpiglo-1 to ji=jpiglo-nb_rimwidth): 236 DO jr=1,nb_rimwidth 237 DO jj=3,jpjglo-2 238 icount=icount+1 239 nbidta(icount,:) = jpiglo-jr + (jpizoom-1) 240 nbidta(icount,2) = jpiglo-jr-1 + (jpizoom-1) ! special case for u points 241 nbjdta(icount,:) = jj + (jpjzoom-1) 242 nbrdta(icount,:) = jr 243 END DO 244 END DO 206 icount = 0 207 ! Define west boundary (from ii=2 to ii=1+nb_rimwidth): 208 DO ir = 1, nb_rimwidth 209 DO ij = 3, jpjglo-2 210 icount=icount+1 211 nbidta(icount,:) = ir + 1 + (jpizoom-1) 212 nbjdta(icount,:) = ij + (jpjzoom-1) 213 nbrdta(icount,:) = ir 214 END DO 215 END DO 216 217 ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nb_rimwidth): 218 DO ir=1,nb_rimwidth 219 DO ij=3,jpjglo-2 220 icount=icount+1 221 nbidta(icount,:) = jpiglo-ir + (jpizoom-1) 222 nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points 223 nbjdta(icount,:) = ij + (jpjzoom-1) 224 nbrdta(icount,:) = ir 225 END DO 226 END DO 245 227 246 ELSE 247 248 ! Read indices and distances in unstructured boundary data files 249 250 IF ( lk_bdy ) THEN 251 bdyfile(1) = filbdy_data_T 252 bdyfile(2) = filbdy_data_U 253 bdyfile(3) = filbdy_data_V 254 ELSE 255 ! In this case we have tides only at the boundaries 256 ! so read index arrays from tides files for first tidal component 257 bdyfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' 258 bdyfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' 259 bdyfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' 260 ENDIF 261 262 DO jgrd = 1,3 263 CALL iom_open( bdyfile(jgrd), inum ) 264 dummy_id = iom_varid( inum, 'nbidta', kdimsz=kdimsz ) 265 WRITE(numout,*) 'kdimsz : ',kdimsz 266 nb_len = kdimsz(1) 267 IF (nb_len > jpbdta) THEN 268 IF(lwp) WRITE(numout,*) 269 IF(lwp) WRITE(numout,*) ' E R R O R : jpbdta is too small:' 270 IF(lwp) WRITE(numout,*) ' ========== Boundary array length in file is ', nb_len 271 IF(lwp) WRITE(numout,*) ' But jpbdta is ', jpbdta 272 IF(lwp) WRITE(numout,*) ' File : ', bdyfile(jgrd) 273 IF(lwp) WRITE(numout,*) 274 nstop = nstop + 1 275 ENDIF 276 CALL iom_get ( inum, jpdom_unknown, 'nbidta', ndta(1:nb_len,:) ) 277 DO ji=1,nb_len 278 nbidta(ji,jgrd) = INT( ndta(ji,1) ) 279 ENDDO 280 CALL iom_get ( inum, jpdom_unknown, 'nbjdta', ndta(1:nb_len,:) ) 281 DO ji=1,nb_len 282 nbjdta(ji,jgrd) = INT( ndta(ji,1) ) 283 ENDDO 284 CALL iom_get ( inum, jpdom_unknown, 'nbrdta', ndta(1:nb_len,:) ) 285 DO ji=1,nb_len 286 nbrdta(ji,jgrd) = INT( ndta(ji,1) ) 287 ENDDO 288 CALL iom_close( inum ) 289 290 ! Check that rimwidth in file is big enough: 291 nbr_max = MAXVAL(nbrdta(:,jgrd)) 292 IF (nbr_max < nb_rimwidth) THEN 293 IF(lwp) WRITE(numout,*) 294 IF(lwp) WRITE(numout,*) ' E R R O R : Maximum rimwidth in file is ', nbr_max 295 IF(lwp) WRITE(numout,*) ' ========== but nb_rimwidth is ', nb_rimwidth 296 IF(lwp) WRITE(numout,*) 297 nstop = nstop + 1 298 ELSE 299 IF(lwp) WRITE(numout,*) 300 IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', nbr_max 301 IF(lwp) WRITE(numout,*) 302 END IF 303 304 ENDDO 305 306 END IF 307 308 ! 1.3 Dispatch mapping indices and discrete distances on each processor 309 ! ********************************************************************* 228 ELSE ! Read indices and distances in unstructured boundary data files 229 230 IF( ln_bdy_tides ) THEN 231 ! Read tides input files for preference in case there are 232 ! no bdydata files. 233 clfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' 234 clfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' 235 clfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' 236 ELSE 237 clfile(1) = filbdy_data_T 238 clfile(2) = filbdy_data_U 239 clfile(3) = filbdy_data_V 240 ENDIF 241 242 ! how many files are we to read in? 243 igrd_start = 1 244 igrd_end = 3 245 IF(.NOT. ln_bdy_tides ) THEN 246 IF(.NOT. (ln_bdy_dyn_fla) .AND..NOT. (ln_bdy_tra_frs)) THEN 247 ! No T-grid file. 248 igrd_start = 2 249 ELSEIF ( .NOT. ln_bdy_dyn_frs .AND..NOT. ln_bdy_dyn_fla ) THEN 250 ! No U-grid or V-grid file. 251 igrd_end = 1 252 ENDIF 253 ENDIF 254 255 DO igrd = igrd_start, igrd_end 256 CALL iom_open( clfile(igrd), inum ) 257 id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz ) 258 WRITE(numout,*) 'kdimsz : ',kdimsz 259 ib_len = kdimsz(1) 260 IF( ib_len > jpbdta) CALL ctl_stop( & 261 'Boundary data array in file too long.', & 262 'File :', TRIM(clfile(igrd)), & 263 'increase parameter jpbdta.' ) 264 265 CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) ) 266 DO ii = 1,ib_len 267 nbidta(ii,igrd) = INT( zdta(ii,1) ) 268 END DO 269 CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) ) 270 DO ii = 1,ib_len 271 nbjdta(ii,igrd) = INT( zdta(ii,1) ) 272 END DO 273 CALL iom_get ( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) ) 274 DO ii = 1,ib_len 275 nbrdta(ii,igrd) = INT( zdta(ii,1) ) 276 END DO 277 CALL iom_close( inum ) 278 279 ! Check that rimwidth in file is big enough: 280 ibr_max = MAXVAL( nbrdta(:,igrd) ) 281 IF(lwp) WRITE(numout,*) 282 IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 283 IF(lwp) WRITE(numout,*) ' nb_rimwidth from namelist is ', nb_rimwidth 284 IF (ibr_max < nb_rimwidth) CALL ctl_stop( & 285 'nb_rimwidth is larger than maximum rimwidth in file' ) 286 ! 287 END DO 288 ! 289 ENDIF 290 291 ! Dispatch mapping indices and discrete distances on each processor 292 ! ***************************************************************** 310 293 311 iw = mig(1) +1 ! if monotasking and no zoom, iw=2312 ie = mig(1) + nlci-1 -1 ! if monotasking and no zoom, ie=jpim1313 is = mjg(1) +1 ! if monotasking and no zoom, is=2314 in = mjg(1) + nlcj-1 -1 ! if monotasking and no zoom, in=jpjm1315 316 DO jgrd = 1, jpbgrd294 iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 295 ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 296 is = mjg(1) + 1 ! if monotasking and no zoom, is=2 297 in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 298 299 DO igrd = igrd_start, igrd_end 317 300 icount = 0 318 301 icountr = 0 319 DO nb_rim=1, nb_rimwidth320 DO jb = 1, jpbdta321 ! check if point is in local domain and equals nb_rim322 IF ( (nbidta(jb,jgrd) >= iw ).AND. &323 (nbidta(jb,jgrd) <= ie ).AND. &324 (nbjdta(jb,jgrd) >= is ).AND. &325 (nbjdta(jb,jgrd) <= in ).AND.&326 (nbrdta(jb,jgrd) == nb_rim ) ) THEN327 328 icount = icount + 1329 330 IF (nb_rim==1) icountr = icountr+1331 332 IF (icount > jpbdim) THEN333 IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'334 nstop = nstop + 1335 ELSE336 nbi(icount, jgrd) = nbidta(jb,jgrd)- mig(1)+1337 nbj(icount, jgrd) = nbjdta(jb,jgrd)- mjg(1)+1338 nbr(icount, jgrd) = nbrdta(jb,jgrd)339 nbmap(icount,jgrd) = jb340 ENDIF341 ENDIF342 END DO343 END DO344 nblenrim(jgrd) = icountr !: length of rim boundary data on each proc345 nblen (jgrd) = icount !: length of boundary data on each proc302 nblen(igrd) = 0 303 nblenrim(igrd) = 0 304 nblendta(igrd) = 0 305 DO ir=1, nb_rimwidth 306 DO ib = 1, jpbdta 307 ! check if point is in local domain and equals ir 308 IF( nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND. & 309 & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND. & 310 & nbrdta(ib,igrd) == ir ) THEN 311 ! 312 icount = icount + 1 313 ! 314 IF( ir == 1 ) icountr = icountr+1 315 IF (icount > jpbdim) THEN 316 IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small' 317 nstop = nstop + 1 318 ELSE 319 nbi(icount, igrd) = nbidta(ib,igrd)- mig(1)+1 320 nbj(icount, igrd) = nbjdta(ib,igrd)- mjg(1)+1 321 nbr(icount, igrd) = nbrdta(ib,igrd) 322 nbmap(icount,igrd) = ib 323 ENDIF 324 ENDIF 325 END DO 326 END DO 327 nblenrim(igrd) = icountr !: length of rim boundary data on each proc 328 nblen (igrd) = icount !: length of boundary data on each proc 346 329 END DO 347 330 348 ! 2. Compute rim weights 349 ! ---------------------- 350 351 DO jgrd = 1, jpbgrd 352 DO jb = 1, nblen(jgrd) 353 ! tanh formulation 354 nbw(jb,jgrd) = 1.-TANH((FLOAT(nbr(jb,jgrd)-1))/2.) 355 ! quadratic 356 ! nbw(jb,jgrd) = (FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth))**2 357 ! linear 358 ! nbw(jb,jgrd) = FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth) 359 END DO 331 ! Compute rim weights 332 ! ------------------- 333 DO igrd = igrd_start, igrd_end 334 DO ib = 1, nblen(igrd) 335 ! tanh formulation 336 nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 ) 337 ! quadratic 338 ! nbw(ib,igrd) = (FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth))**2 339 ! linear 340 ! nbw(ib,igrd) = FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth) 341 END DO 360 342 END DO 361 343 362 ! 3. Mask corrections 363 ! ------------------- 364 365 DO jk=1, jpkm1 366 DO jj=1, jpj 367 DO ji=1, jpi 368 tmask(ji,jj,jk)=tmask(ji,jj,jk)*bdytmask(ji,jj) 369 umask(ji,jj,jk)=umask(ji,jj,jk)*bdyumask(ji,jj) 370 vmask(ji,jj,jk)=vmask(ji,jj,jk)*bdyvmask(ji,jj) 371 bmask(ji,jj)=bmask(ji,jj)*bdytmask(ji,jj) 372 END DO 373 END DO 374 END DO 375 376 ! I am not sure that it is useful: 377 DO jk=1, jpkm1 378 DO jj=2, jpjm1 379 DO ji=2, jpim1 380 fmask(ji,jj,jk) = fmask(ji,jj,jk) & 381 & * bdytmask(ji, jj ) * bdytmask(ji+1, jj ) & 382 & * bdytmask(ji,jj+1) * bdytmask(ji+1,jj+1) 383 END DO 384 END DO 385 END DO 386 387 tmask_i(:,:) = tmask(:,:,1)*tmask_i(:,:) 388 389 bdytmask(:,:)=tmask(:,:,1) 344 ! Mask corrections 345 ! ---------------- 346 DO ik = 1, jpkm1 347 DO ij = 1, jpj 348 DO ii = 1, jpi 349 tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 350 umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 351 vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 352 bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij) 353 END DO 354 END DO 355 END DO 356 357 DO ik = 1, jpkm1 358 DO ij = 2, jpjm1 359 DO ii = 2, jpim1 360 fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) & 361 & * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 362 END DO 363 END DO 364 END DO 365 366 tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 367 bdytmask(:,:) = tmask(:,:,1) 390 368 391 369 ! bdy masks and bmask are now set to zero on boundary points: 392 393 jgrd=1 ! In the free surface case, bmask is at T-points 394 DO jb=1, nblenrim(jgrd) 395 bmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 396 END DO 397 398 jgrd=1 399 DO jb=1, nblenrim(jgrd) 400 bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 401 END DO 402 403 jgrd=2 404 DO jb=1, nblenrim(jgrd) 405 bdyumask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 406 END DO 407 408 jgrd=3 409 DO jb=1, nblenrim(jgrd) 410 bdyvmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 411 END DO 412 413 ! Lateral boundary conditions 414 415 CALL lbc_lnk( fmask, 'F', 1. ) 370 igrd = 1 ! In the free surface case, bmask is at T-points 371 DO ib = 1, nblenrim(igrd) 372 bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 373 END DO 374 ! 375 igrd = 1 376 DO ib = 1, nblenrim(igrd) 377 bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 378 END DO 379 ! 380 igrd = 2 381 DO ib = 1, nblenrim(igrd) 382 bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 383 END DO 384 ! 385 igrd = 3 386 DO ib = 1, nblenrim(igrd) 387 bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 388 END DO 389 390 ! Lateral boundary conditions 391 CALL lbc_lnk( fmask , 'F', 1. ) 416 392 CALL lbc_lnk( bdytmask(:,:), 'T', 1. ) 417 393 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) 418 394 CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 419 395 420 IF ((ln_bdy_vol).OR.(ln_bdy_fla)) THEN 421 422 ! 4 Indices and directions of rim velocity components 423 ! --------------------------------------------------- 424 !flagu = -1 : u component is normal to the dynamical boundary 425 ! but its direction is outward 426 ! 427 !flagu = 0 : u is tangential 428 ! 429 !flagu = 1 : u is normal to the boundary 430 ! and is direction is inward 396 IF( ln_bdy_vol .OR. ln_bdy_dyn_fla ) THEN ! Indices and directions of rim velocity components 397 ! 398 !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 399 !flagu = 0 : u is tangential 400 !flagu = 1 : u is normal to the boundary and is direction is inward 401 icount = 0 402 flagu(:) = 0.e0 431 403 432 icount = 0 433 434 flagu(:)=0.e0 404 igrd = 2 ! u-component 405 DO ib = 1, nblenrim(igrd) 406 zefl=bdytmask(nbi(ib,igrd) , nbj(ib,igrd)) 407 zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd)) 408 IF( zefl + zwfl ==2 ) THEN 409 icount = icount +1 410 ELSE 411 flagu(ib)=-zefl+zwfl 412 ENDIF 413 END DO 414 415 !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward 416 !flagv = 0 : u is tangential 417 !flagv = 1 : u is normal to the boundary and is direction is inward 418 flagv(:) = 0.e0 419 420 igrd = 3 ! v-component 421 DO ib = 1, nblenrim(igrd) 422 znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)) 423 zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1) 424 IF( znfl + zsfl ==2 ) THEN 425 icount = icount + 1 426 ELSE 427 flagv(ib) = -znfl + zsfl 428 END IF 429 END DO 435 430 436 jgrd=2 ! u-component 437 DO jb=1, nblenrim(jgrd) 438 efl=bdytmask(nbi(jb,jgrd) , nbj(jb,jgrd)) 439 wfl=bdytmask(nbi(jb,jgrd)+1, nbj(jb,jgrd)) 440 IF ((efl+wfl)==2) THEN 441 icount = icount +1 442 ELSE 443 flagu(jb)=-efl+wfl 444 END IF 445 END DO 446 447 !flagv = -1 : u component is normal to the dynamical boundary 448 ! but its direction is outward 449 ! 450 !flagv = 0 : u is tangential 451 ! 452 !flagv = 1 : u is normal to the boundary 453 ! and is direction is inward 454 455 flagv(:)=0.e0 456 457 jgrd=3 ! v-component 458 DO jb=1, nblenrim(jgrd) 459 nfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) 460 sfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)+1) 461 IF ((nfl+sfl)==2) THEN 462 icount = icount +1 463 ELSE 464 flagv(jb)=-nfl+sfl 465 END IF 466 END DO 467 468 IF( icount /= 0 ) THEN 469 IF(lwp) WRITE(numout,*) 470 IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & 471 ' are not boundary points. Check nbi, nbj, indices.' 472 IF(lwp) WRITE(numout,*) ' ========== ' 473 IF(lwp) WRITE(numout,*) 474 nstop = nstop + 1 475 END IF 431 IF( icount /= 0 ) THEN 432 IF(lwp) WRITE(numout,*) 433 IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & 434 ' are not boundary points. Check nbi, nbj, indices.' 435 IF(lwp) WRITE(numout,*) ' ========== ' 436 IF(lwp) WRITE(numout,*) 437 nstop = nstop + 1 438 ENDIF 476 439 477 END 478 479 ! 5Compute total lateral surface for volume correction:480 ! ---------------------------------------------------- --440 ENDIF 441 442 ! Compute total lateral surface for volume correction: 443 ! ---------------------------------------------------- 481 444 482 445 bdysurftot = 0.e0 483 484 IF (ln_bdy_vol) THEN 485 jgrd=2 ! Lateral surface at U-points 486 DO jb=1, nblenrim(jgrd) 487 bdysurftot = bdysurftot + & 488 hu(nbi(jb,jgrd), nbj(jb,jgrd)) * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & 489 * ABS(flagu(jb))*tmask_i(nbi(jb,jgrd) , nbj(jb,jgrd)) & 490 *tmask_i(nbi(jb,jgrd)+1, nbj(jb,jgrd)) 491 END DO 492 493 jgrd=3 ! Add lateral surface at V-points 494 DO jb=1, nblenrim(jgrd) 495 bdysurftot = bdysurftot + & 496 hv(nbi(jb,jgrd), nbj(jb,jgrd)) * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 497 * ABS(flagv(jb))*tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)) & 498 *tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)+1) 499 END DO 500 501 IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain 446 IF( ln_bdy_vol ) THEN 447 igrd = 2 ! Lateral surface at U-points 448 DO ib = 1, nblenrim(igrd) 449 bdysurftot = bdysurftot + hu (nbi(ib,igrd) ,nbj(ib,igrd)) & 450 & * e2u (nbi(ib,igrd) ,nbj(ib,igrd)) * ABS( flagu(ib) ) & 451 & * tmask_i(nbi(ib,igrd) ,nbj(ib,igrd)) & 452 & * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd)) 453 END DO 454 455 igrd=3 ! Add lateral surface at V-points 456 DO ib = 1, nblenrim(igrd) 457 bdysurftot = bdysurftot + hv (nbi(ib,igrd),nbj(ib,igrd) ) & 458 & * e1v (nbi(ib,igrd),nbj(ib,igrd) ) * ABS( flagv(ib) ) & 459 & * tmask_i(nbi(ib,igrd),nbj(ib,igrd) ) & 460 & * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1) 461 END DO 462 463 IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain 502 464 END IF 503 465 504 505 ! 6. Initialise bdy data arrays 506 ! ----------------------------- 507 466 ! Initialise bdy data arrays 467 ! -------------------------- 508 468 tbdy(:,:) = 0.e0 509 469 sbdy(:,:) = 0.e0 … … 514 474 vbtbdy(:) = 0.e0 515 475 516 ! 7. Read in tidal constituents and adjust for model start time 517 ! ------------------------------------------------------------- 518 519 IF ( lk_bdy_tides ) CALL tide_data 520 476 ! Read in tidal constituents and adjust for model start time 477 ! ---------------------------------------------------------- 478 IF( ln_bdy_tides ) CALL tide_data 479 ! 521 480 END SUBROUTINE bdy_init 522 481
Note: See TracChangeset
for help on using the changeset viewer.
