Changeset 1125 for trunk/NEMO/OPA_SRC/BDY/bdydta.F90
- Timestamp:
- 2008-06-23T11:05:02+02:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/BDY/bdydta.F90
r1058 r1125 1 1 MODULE bdydta 2 !!====================================================================== ========3 !! 2 !!====================================================================== 3 !! *** MODULE bdydta *** 4 4 !! Open boundary data : read the data for the unstructured open boundaries. 5 !!============================================================================== 6 #if defined key_bdy || defined key_bdy_tides 7 !!------------------------------------------------------------------------------ 8 !! 'key_bdy' : Unstructured Open Boundary Conditions 9 !!------------------------------------------------------------------------------ 10 !! bdy_dta : read u, v, t, s data along each open boundary 11 !!------------------------------------------------------------------------------ 12 !! * Modules used 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-07 (D. Storkey) add bdy_dta_bt 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 !!---------------------------------------------------------------------- 15 !! bdy_dta : read u, v, t, s data along open boundaries 16 !! bdy_dta_bt : read depth-mean velocities and elevation along open 17 !! boundaries 18 !!---------------------------------------------------------------------- 13 19 USE oce ! ocean dynamics and tracers 14 20 USE dom_oce ! ocean space and time domain … … 24 30 PRIVATE 25 31 26 !! * Accessibility 27 PUBLIC bdy_dta ! routines called by step.F90 28 PUBLIC bdy_dta_bt 32 PUBLIC bdy_dta ! routines called by step.F90 33 PUBLIC bdy_dta_bt 34 35 INTEGER :: numbdyt, numbdyu, numbdyv !: logical units for T-, U-, & V-points data file, resp. 36 INTEGER :: ntimes_bdy !: exact number of time dumps in data files 37 INTEGER :: nbdy_b, nbdy_a !: record of bdy data file for before and after model time step 38 INTEGER :: numbdyt_bt, numbdyu_bt, numbdyv_bt !: logical unit for T-, U- & V-points data file, resp. 39 INTEGER :: ntimes_bdy_bt !: exact number of time dumps in data files 40 INTEGER :: nbdy_b_bt, nbdy_a_bt !: record of bdy data file for before and after model time step 41 42 INTEGER, DIMENSION (jpbtime) :: istep, istep_bt !: time array in seconds in each data file 43 44 REAL(wp) :: zoffset !: time offset between time origin in file & start time of model run 45 46 REAL(wp), DIMENSION(jpbdim,jpk,2) :: tbdydta, sbdydta !: time interpolated values of T and S bdy data 47 REAL(wp), DIMENSION(jpbdim,jpk,2) :: ubdydta, vbdydta !: time interpolated values of U and V bdy data 48 REAL(wp), DIMENSION(jpbdim,2) :: ubtbdydta, vbtbdydta !: Arrays used for time interpolation of bdy data 49 REAL(wp), DIMENSION(jpbdim,2) :: sshbdydta !: bdy data of ssh 50 51 !!---------------------------------------------------------------------- 52 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 53 !! $Id: $ 54 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 55 !!---------------------------------------------------------------------- 56 57 CONTAINS 58 59 SUBROUTINE bdy_dta( kt ) 60 !!---------------------------------------------------------------------- 61 !! *** SUBROUTINE bdy_dta *** 62 !! 63 !! ** Purpose : Read unstructured boundary data for FRS condition. 64 !! 65 !! ** Method : At the first timestep, read in boundary data for two 66 !! times from the file and time-interpolate. At other 67 !! timesteps, check to see if we need another time from 68 !! the file. If so read it in. Time interpolate. 69 !!---------------------------------------------------------------------- 70 INTEGER, INTENT( in ) :: kt ! ocean time-step index (for timesplitting option, otherwise zero) 71 !! 72 CHARACTER(LEN=80), DIMENSION(3) :: clfile ! names of input files 73 CHARACTER(LEN=70 ) :: clunits ! units attribute of time coordinate 74 LOGICAL :: lect ! flag for reading 75 INTEGER :: it, ib, ik, igrd ! dummy loop indices 76 INTEGER :: igrd_start, igrd_end ! start and end of loops on igrd 77 INTEGER :: idvar ! netcdf var ID 78 INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing 79 INTEGER :: ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 80 INTEGER :: itimer, totime 81 INTEGER :: ii, ij ! array addresses 82 INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read) 83 INTEGER :: iyear0, imonth0, iday0 84 INTEGER :: ihours0, iminutes0, isec0 85 INTEGER :: iyear, imonth, iday, isecs 86 INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files 87 REAL(wp) :: dayfrac, zxy, zoffsett 88 REAL(wp) :: zoffsetu, zoffsetv 89 REAL(wp) :: dayjul0, zdayjulini 90 REAL(wp), DIMENSION(jpbtime) :: zstepr ! REAL time array from data files 91 REAL(wp), DIMENSION(jpbdta,1,jpk) :: zdta ! temporary array for data fields 92 !!--------------------------------------------------------------------------- 93 94 IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_frs ) THEN ! If these are both false then this routine 95 ! does nothing. 96 97 ! -------------------- ! 98 ! Initialization ! 99 ! -------------------- ! 100 101 lect = .false. ! If true, read a time record 102 103 ! Some time variables for monthly climatological forcing: 104 ! ******************************************************* 105 !!gm here use directely daymod variables 106 107 iman = INT( raamo ) ! Number of months in a year 108 109 i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 110 ! i15=0 if the current day is in the first half of the month, else i15=1 111 112 imois = nmonth + i15 - 1 ! imois is the first month record 113 IF( imois == 0 ) imois = iman 114 115 ! Time variable for non-climatological forcing: 116 ! ********************************************* 117 itimer = (kt-nit000+1)*rdt ! current time in seconds for interpolation 118 119 120 ! !-------------------! 121 IF( kt == nit000 ) THEN ! First call only ! 122 ! !-------------------! 123 istep(:) = 0 124 nbdy_b = 0 125 nbdy_a = 0 126 127 ! Get time information from bdy data file 128 ! *************************************** 129 130 IF(lwp) WRITE(numout,*) 131 IF(lwp) WRITE(numout,*) 'bdy_dta : Initialize unstructured boundary data' 132 IF(lwp) WRITE(numout,*) '~~~~~~~' 133 134 IF ( nbdy_dta == 0 ) THEN 135 ! 136 IF(lwp) WRITE(numout,*) ' Bdy data are taken from initial conditions' 137 ! 138 ELSEIF (nbdy_dta == 1) THEN 139 ! 140 IF(lwp) WRITE(numout,*) ' Bdy data are read in netcdf files' 141 ! 142 dayfrac = adatrj - FLOAT( itimer ) / 86400. ! day fraction at time step kt-1 143 dayfrac = dayfrac - INT ( dayfrac ) ! 144 totime = ( nitend - nit000 + 1 ) * rdt ! Total time of the run to verify that all the 145 ! ! necessary time dumps in file are included 146 ! 147 clfile(1) = filbdy_data_T 148 clfile(2) = filbdy_data_U 149 clfile(3) = filbdy_data_V 150 ! 151 ! how many files are we to read in? 152 igrd_start = 1 153 igrd_end = 3 154 IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN 155 ! No T-grid file. 156 igrd_start = 2 157 ELSEIF ( .NOT. ln_bdy_dyn_frs ) THEN 158 ! No U-grid or V-grid file. 159 igrd_end = 1 160 ENDIF 161 162 DO igrd = igrd_start, igrd_end ! loop over T, U & V grid ! 163 ! !---------------------------! 164 CALL iom_open( clfile(igrd), inum ) 165 CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits ) 166 167 SELECT CASE( igrd ) 168 CASE (1) 169 numbdyt = inum 170 CASE (2) 171 numbdyu = inum 172 CASE (3) 173 numbdyv = inum 174 END SELECT 175 176 ! Calculate time offset 177 READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 178 ! Convert time origin in file to julian days 179 isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 180 CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 181 ! Compute model initialization time 182 iyear = ndastp / 10000 183 imonth = ( ndastp - iyear * 10000 ) / 100 184 iday = ndastp - iyear * 10000 - imonth * 100 185 isecs = dayfrac * 86400 186 CALL ymds2ju(iyear, imonth, iday, real(isecs) , zdayjulini) 187 ! offset from initialization date: 188 zoffset = (dayjul0-zdayjulini)*86400 189 ! 190 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 191 192 !! TO BE DONE... Check consistency between calendar from file 193 !! (available optionally from iom_gettime) and calendar in model 194 !! when calendar in model available outside of IOIPSL. 195 196 IF(lwp) WRITE(numout,*) 'number of times: ',ntimes_bdy 197 IF(lwp) WRITE(numout,*) 'offset: ',zoffset 198 IF(lwp) WRITE(numout,*) 'totime: ',totime 199 IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr 200 201 ! Check that there are not too many times in the file. 202 IF( ntimes_bdy > jpbtime ) THEN 203 WRITE(ctmp1,*) 'Check file: ', clfile(igrd), 'jpbtime= ', jpbtime, ' ntimes_bdy= ', ntimes_bdy 204 CALL ctl_stop( 'Number of time dumps in files exceed jpbtime parameter', ctmp1 ) 205 ENDIF 206 207 ! Check that time array increases: 208 209 it = 1 210 DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 ) 211 it = it + 1 212 END DO 213 214 IF( it.NE.ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN 215 WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 216 CALL ctl_stop( 'Time array in unstructured boundary data files', & 217 & 'does not continuously increase.' , ctmp1 ) 218 ENDIF 219 ! 220 ! Check that times in file span model run time: 221 IF( zstepr(1) + zoffset > 0 ) THEN 222 WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 223 CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 ) 224 END IF 225 IF( zstepr(ntimes_bdy) + zoffset < totime ) THEN 226 WRITE(ctmp1,*) 'Check file: ', clfile(igrd) 227 CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 ) 228 END IF 229 ! 230 IF ( igrd == 1 ) THEN 231 ntimes_bdyt = ntimes_bdy 232 zoffsett = zoffset 233 istept(:) = INT( zstepr(:) + zoffset ) 234 ELSEIF(igrd == 2 ) THEN 235 ntimes_bdyu = ntimes_bdy 236 zoffsetu = zoffset 237 istepu(:) = INT( zstepr(:) + zoffset ) 238 ELSEIF(igrd == 3 ) THEN 239 ntimes_bdyv = ntimes_bdy 240 zoffsetv = zoffset 241 istepv(:) = INT( zstepr(:) + zoffset ) 242 ENDIF 243 ! 244 END DO ! end loop over T, U & V grid 245 246 IF (igrd_start == 1 .and. igrd_end == 3) THEN 247 ! Only test differences if we are reading in 3 files 248 ! Verify time consistency between files 249 IF( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 250 CALL ctl_stop( 'Bdy data files must have the same number of time dumps', & 251 & 'Multiple time frequencies not implemented yet' ) 252 ENDIF 253 ntimes_bdy = ntimes_bdyt 254 ! 255 IF( zoffsetu /= zoffsett .OR. zoffsetv /= zoffsett ) THEN 256 CALL ctl_stop( 'Bdy data files must have the same time origin', & 257 & 'Multiple time frequencies not implemented yet' ) 258 ENDIF 259 zoffset = zoffsett 260 ENDIF 261 262 IF( igrd_start == 1 ) THEN 263 istep(:) = istept(:) 264 ELSE 265 istep(:) = istepu(:) 266 ENDIF 267 268 ! Check number of time dumps: 269 IF( ntimes_bdy == 1 .AND. .NOT. ln_bdy_clim ) THEN 270 CALL ctl_stop( 'There is only one time dump in data files', & 271 & 'Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' ) 272 ENDIF 273 274 IF( ln_bdy_clim ) THEN 275 IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN 276 CALL ctl_stop( 'For climatological boundary forcing (ln_bdy_clim=.true.),', & 277 & 'bdy data files must contain 1 or 12 time dumps.' ) 278 ELSEIF( ntimes_bdy == 1 ) THEN 279 IF(lwp) WRITE(numout,*) 280 IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 281 ELSEIF( ntimes_bdy == 12 ) THEN 282 IF(lwp) WRITE(numout,*) 283 IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 284 ENDIF 285 ENDIF 286 287 ! Find index of first record to read (before first model time). 288 it = 1 289 DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 ) 290 it = it + 1 291 END DO 292 nbdy_b = it 293 ! 294 WRITE(numout,*) 'Time offset is ',zoffset 295 WRITE(numout,*) 'First record to read is ',nbdy_b 296 297 ENDIF ! endif (nbdy_dta == 1) 298 299 300 ! 1.2 Read first record in file if necessary (ie if nbdy_dta == 1) 301 ! ***************************************************************** 302 303 IF( nbdy_dta == 0) THEN ! boundary data arrays are filled with initial conditions 304 ! 305 IF (ln_bdy_tra_frs) THEN 306 igrd = 1 ! T-points data 307 DO ib = 1, nblen(igrd) 308 ii = nbi(ib,igrd) 309 ij = nbj(ib,igrd) 310 DO ik = 1, jpkm1 311 tbdy(ib,ik) = tn(ii, ij, ik) 312 sbdy(ib,ik) = sn(ii, ij, ik) 313 ENDDO 314 END DO 315 ENDIF 316 317 IF(ln_bdy_dyn_frs) THEN 318 igrd = 2 ! U-points data 319 DO ib = 1, nblen(igrd) 320 ii = nbi(ib,igrd) 321 ij = nbj(ib,igrd) 322 DO ik = 1, jpkm1 323 ubdy(ib,ik) = un(ii, ij, ik) 324 ENDDO 325 END DO 326 327 igrd = 3 ! V-points data 328 DO ib = 1, nblen(igrd) 329 ii = nbi(ib,igrd) 330 ij = nbj(ib,igrd) 331 DO ik = 1, jpkm1 332 vbdy(ib,ik) = vn(ii, ij, ik) 333 ENDDO 334 END DO 335 ENDIF 336 ! 337 ELSEIF( nbdy_dta == 1 ) THEN ! Set first record in the climatological case: 338 ! 339 IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN 340 nbdy_a = 1 341 ELSEIF( ln_bdy_clim .AND. ntimes_bdy == iman ) THEN 342 nbdy_b = 0 343 nbdy_a = imois 344 ELSE 345 nbdy_a = nbdy_b 346 ENDIF 29 347 30 !! * local modules variables 31 INTEGER :: & 32 bdynumt, & ! logical unit for T-points data file 33 bdynumu, & ! logical unit for U-points data file 34 bdynumv, & ! logical unit for V-points data file 35 kbdy, & ! Exact number of time dumps in data files 36 nbdy1, & ! Time record in bdy data file BEFORE the model current time 37 nbdy2 ! Time record in bdy data file AFTER the model current time 38 39 !! * local modules variables for barotropic variables 40 INTEGER :: & 41 bdynumt_bt, & ! logical unit for T-points data file 42 bdynumu_bt, & ! logical unit for U-points data file 43 bdynumv_bt, & ! logical unit for V-points data file 44 kbdy_bt, & ! Exact number of time dumps in data files 45 nbdy1_bt, & ! Time record in bdy data file BEFORE the model current time 46 nbdy2_bt ! Time record in bdy data file AFTER the model current time 47 48 INTEGER, DIMENSION (jpbtime) :: & 49 istep, istep_bt ! time array in seconds in each data file 50 51 REAL(wp) :: & 52 offset ! time offset between time origin in file and start time of model run 53 54 REAL(wp), DIMENSION(jpbdim,jpk,2) :: & 55 tbdydta, sbdydta, & !: Arrays used for time interpolation of bdy data 56 ubdydta, vbdydta 57 58 REAL(wp), DIMENSION(jpbdim,2) :: & 59 ubtbdydta, vbtbdydta, & !: Arrays used for time interpolation of bdy data 60 sshbdydta !: bdy data of ssh 61 62 !!--------------------------------------------------------------------------------- 63 !! OPA 9.0 , LODYC-IPSL (2003) 64 !!--------------------------------------------------------------------------------- 65 66 CONTAINS 67 68 SUBROUTINE bdy_dta( kt ) 348 ! Read first record: 349 ipj = 1 350 ipk = jpk 351 igrd = 1 352 ipi = nblendta(igrd) 353 354 IF(ln_bdy_tra_frs) THEN 355 igrd = 1 ! Temperature 356 IF( nblendta(igrd) <= 0 ) THEN 357 idvar = iom_varid( numbdyt, 'votemper' ) 358 nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 359 ENDIF 360 WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd) 361 ipi = nblendta(igrd) 362 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 363 364 DO ib = 1, nblen(igrd) 365 DO ik = 1, jpkm1 366 tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 367 END DO 368 END DO 369 ! 370 igrd = 1 ! salinity 371 IF( nblendta(igrd) .le. 0 ) THEN 372 idvar = iom_varid( numbdyt, 'vosaline' ) 373 nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar) 374 ENDIF 375 WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd) 376 ipi = nblendta(igrd) 377 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 378 379 DO ib = 1, nblen(igrd) 380 DO ik = 1, jpkm1 381 sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 382 END DO 383 END DO 384 ENDIF ! ln_bdy_tra_frs 385 386 IF(ln_bdy_dyn_frs) THEN 387 388 igrd = 2 ! u-velocity 389 IF ( nblendta(igrd) .le. 0 ) THEN 390 idvar = iom_varid( numbdyu,'vozocrtx' ) 391 nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar) 392 ENDIF 393 WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd) 394 ipi = nblendta(igrd) 395 CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 396 DO ib = 1, nblen(igrd) 397 DO ik = 1, jpkm1 398 ubdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 399 END DO 400 END DO 401 ! 402 igrd = 3 ! v-velocity 403 IF ( nblendta(igrd) .le. 0 ) THEN 404 idvar = iom_varid( numbdyv,'vomecrty' ) 405 nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar) 406 ENDIF 407 WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd) 408 ipi = nblendta(igrd) 409 CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 410 DO ib = 1, nblen(igrd) 411 DO ik = 1, jpkm1 412 vbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 413 END DO 414 END DO 415 ENDIF ! ln_bdy_dyn_frs 416 417 418 IF ((.NOT.ln_bdy_clim) .AND. (istep(1) > 0)) THEN 419 ! First data time is after start of run 420 ! Put first value in both time levels 421 nbdy_b = nbdy_a 422 IF(ln_bdy_tra_frs) THEN 423 tbdydta(:,:,1) = tbdydta(:,:,2) 424 sbdydta(:,:,1) = sbdydta(:,:,2) 425 ENDIF 426 IF(ln_bdy_dyn_frs) THEN 427 ubdydta(:,:,1) = ubdydta(:,:,2) 428 vbdydta(:,:,1) = vbdydta(:,:,2) 429 ENDIF 430 END IF 431 432 END IF ! nbdy_dta == 0/1 433 434 ! In the case of constant boundary forcing fill bdy arrays once for all 435 IF ((ln_bdy_clim).AND.(ntimes_bdy==1)) THEN 436 IF(ln_bdy_tra_frs) THEN 437 tbdy (:,:) = tbdydta (:,:,2) 438 sbdy (:,:) = sbdydta (:,:,2) 439 ENDIF 440 IF(ln_bdy_dyn_frs) THEN 441 ubdy (:,:) = ubdydta (:,:,2) 442 vbdy (:,:) = vbdydta (:,:,2) 443 ENDIF 444 445 IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt ) 446 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu ) 447 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv ) 448 END IF 449 450 ENDIF ! End if nit000 451 452 453 ! !---------------------! 454 ! ! at each time step ! 455 ! !---------------------! 456 457 IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN 458 ! 459 ! Read one more record if necessary 460 !********************************** 461 462 IF( ln_bdy_clim .AND. imois /= nbdy_b ) THEN ! remember that nbdy_b=0 for kt=nit000 463 nbdy_b = imois 464 nbdy_a = imois + 1 465 nbdy_b = MOD( nbdy_b, iman ) ; IF( nbdy_b == 0 ) nbdy_b = iman 466 nbdy_a = MOD( nbdy_a, iman ) ; IF( nbdy_a == 0 ) nbdy_a = iman 467 lect=.true. 468 ELSEIF( .NOT.ln_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN 469 470 IF ( nbdy_a < ntimes_bdy ) THEN 471 nbdy_b = nbdy_a 472 nbdy_a = nbdy_a + 1 473 lect =.true. 474 ELSE 475 ! We have reached the end of the file 476 ! put the last data time into both time levels 477 nbdy_b = nbdy_a 478 IF(ln_bdy_tra_frs) THEN 479 tbdydta(:,:,1) = tbdydta(:,:,2) 480 sbdydta(:,:,1) = sbdydta(:,:,2) 481 ENDIF 482 IF(ln_bdy_dyn_frs) THEN 483 ubdydta(:,:,1) = ubdydta(:,:,2) 484 vbdydta(:,:,1) = vbdydta(:,:,2) 485 ENDIF 486 END IF ! nbdy_a < ntimes_bdy 487 488 END IF 489 490 IF( lect ) THEN 491 ! Swap arrays 492 IF(ln_bdy_tra_frs) THEN 493 tbdydta(:,:,1) = tbdydta(:,:,2) 494 sbdydta(:,:,1) = sbdydta(:,:,2) 495 ENDIF 496 IF(ln_bdy_dyn_frs) THEN 497 ubdydta(:,:,1) = ubdydta(:,:,2) 498 vbdydta(:,:,1) = vbdydta(:,:,2) 499 ENDIF 500 501 ! read another set 502 ipj = 1 503 ipk = jpk 504 505 IF(ln_bdy_tra_frs) THEN 506 ! 507 igrd = 1 ! temperature 508 ipi = nblendta(igrd) 509 CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 510 DO ib = 1, nblen(igrd) 511 DO ik = 1, jpkm1 512 tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 513 END DO 514 END DO 515 ! 516 igrd = 1 ! salinity 517 ipi = nblendta(igrd) 518 CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a ) 519 DO ib = 1, nblen(igrd) 520 DO ik = 1, jpkm1 521 sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 522 END DO 523 END DO 524 ENDIF ! ln_bdy_tra_frs 525 526 IF(ln_bdy_dyn_frs) THEN 527 ! 528 igrd = 2 ! u-velocity 529 ipi = nblendta(igrd) 530 CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 531 DO ib = 1, nblen(igrd) 532 DO ik = 1, jpkm1 533 ubdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 534 END DO 535 END DO 536 ! 537 igrd = 3 ! v-velocity 538 ipi = nblendta(igrd) 539 CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a ) 540 DO ib = 1, nblen(igrd) 541 DO ik = 1, jpkm1 542 vbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik) 543 END DO 544 END DO 545 ENDIF ! ln_bdy_dyn_frs 546 547 ! 548 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b ',nbdy_b 549 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy_a ',nbdy_a 550 IF (.NOT.ln_bdy_clim) THEN 551 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep(nbdy_b) 552 IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer 553 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy_a) 554 ENDIF 555 ! 556 ENDIF ! end lect=.true. 557 558 559 ! Interpolate linearly 560 ! ******************** 561 ! 562 IF( ln_bdy_clim ) THEN ; zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b) ) + 0.5 - i15 563 ELSE ; zxy = FLOAT( istep(nbdy_b) - itimer ) / FLOAT( istep(nbdy_b) - istep(nbdy_a) ) 564 END IF 565 566 IF(ln_bdy_tra_frs) THEN 567 igrd = 1 ! temperature & salinity 568 DO ib = 1, nblen(igrd) 569 DO ik = 1, jpkm1 570 tbdy(ib,ik) = zxy * tbdydta(ib,ik,2) + (1.-zxy) * tbdydta(ib,ik,1) 571 sbdy(ib,ik) = zxy * sbdydta(ib,ik,2) + (1.-zxy) * sbdydta(ib,ik,1) 572 END DO 573 END DO 574 ENDIF 575 576 IF(ln_bdy_dyn_frs) THEN 577 igrd = 2 ! u-velocity 578 DO ib = 1, nblen(igrd) 579 DO ik = 1, jpkm1 580 ubdy(ib,ik) = zxy * ubdydta(ib,ik,2) + (1.-zxy) * ubdydta(ib,ik,1) 581 END DO 582 END DO 583 ! 584 igrd = 3 ! v-velocity 585 DO ib = 1, nblen(igrd) 586 DO ik = 1, jpkm1 587 vbdy(ib,ik) = zxy * vbdydta(ib,ik,2) + (1.-zxy) * vbdydta(ib,ik,1) 588 END DO 589 END DO 590 ENDIF 591 592 END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy>1)) 593 594 595 ! !---------------------! 596 ! ! last call ! 597 ! !---------------------! 598 IF( kt == nitend ) THEN 599 IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt ) ! Closing of the 3 files 600 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu ) 601 IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv ) 602 ENDIF 603 ! 604 ENDIF ! ln_bdy_dyn_frs .OR. ln_bdy_tra_frs 605 606 END SUBROUTINE bdy_dta 607 608 609 SUBROUTINE bdy_dta_bt( kt, jit ) 69 610 !!--------------------------------------------------------------------------- 70 !! *** SUBROUTINE bdy_dta ***611 !! *** SUBROUTINE bdy_dta_bt *** 71 612 !! 72 !! ** Purpose : Read unstructured boundary data 613 !! ** Purpose : Read unstructured boundary data for Flather condition 73 614 !! 74 !! ** Method : 615 !! ** Method : At the first timestep, read in boundary data for two 616 !! times from the file and time-interpolate. At other 617 !! timesteps, check to see if we need another time from 618 !! the file. If so read it in. Time interpolate. 619 !!--------------------------------------------------------------------------- 620 !!gm DOCTOR names : argument integer : start with "k" 621 INTEGER, INTENT( in ) :: kt ! ocean time-step index 622 INTEGER, INTENT( in ) :: jit ! barotropic time step index 623 ! ! (for timesplitting option, otherwise zero) 75 624 !! 76 !! History : OPA 9.0 ! 05-01 (J. Chanut, A. Sellar) Original77 !! " ! 07-01 (D. Storkey) Update to use IOM module.78 !!---------------------------------------------------------------------------79 !! * Arguments80 INTEGER, INTENT( in ) :: kt ! ocean time-step index81 ! (for timesplitting option, otherwise zero)82 83 !! * Local declarations84 625 LOGICAL :: lect ! flag for reading 85 INTEGER :: jt, jb, jk, jgrd! dummy loop indices626 INTEGER :: it, ib, igrd ! dummy loop indices 86 627 INTEGER :: idvar ! netcdf var ID 87 628 INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing 88 INTEGER :: kbdyt, kbdyu, kbdyv629 INTEGER :: ntimes_bdyt, ntimes_bdyu, ntimes_bdyv 89 630 INTEGER :: itimer, totime 90 631 INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read) 91 632 INTEGER :: iyear0, imonth0, iday0 92 633 INTEGER :: ihours0, iminutes0, isec0 93 INTEGER :: kyear, kmonth, kday, ksecs634 INTEGER :: iyear, imonth, iday, isecs 94 635 INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files 95 REAL(wp) :: dayfrac, zxy, offsett 96 REAL(wp) :: offsetu, offsetv 97 REAL(wp) :: dayjul0, zdayjulini 98 REAL(wp), DIMENSION(jpbtime) :: istepr ! REAL time array from data files 99 REAL(wp), DIMENSION(jpbdta,1,jpk) :: pdta3 ! temporary array for data fields 100 CHARACTER(LEN=80), DIMENSION(3) :: bdyfile 636 REAL(wp) :: dayfrac, zxy, zoffsett 637 REAL(wp) :: zoffsetu, zoffsetv 638 REAL(wp) :: dayjul0, zdayjulini, zntotime 639 REAL(wp) :: zinterval_s, zinterval_e ! First and last interval in time axis 640 REAL(wp), DIMENSION(jpbtime) :: zstepr ! REAL time array from data files 641 REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array for data fields 642 CHARACTER(LEN=80), DIMENSION(3) :: clfile 101 643 CHARACTER(LEN=70 ) :: clunits ! units attribute of time coordinate 102 103 644 !!--------------------------------------------------------------------------- 645 646 !!gm add here the same style as in bdy_dta 647 !!gm clearly bdy_dta_bt and bdy_dta can be combined... 648 !!gm too many things duplicated in the read of data... simplification can be done 104 649 105 650 ! -------------------- ! … … 111 656 ! Some time variables for monthly climatological forcing: 112 657 ! ******************************************************* 658 !!gm here use directely daymod variables 113 659 114 660 iman = INT( raamo ) ! Number of months in a year … … 123 669 ! ********************************************* 124 670 125 itimer = (kt-nit000+1)*rdt ! current time in seconds for interpolation 126 127 ! -------------------- ! 128 ! 1. First call only ! 129 ! -------------------- ! 130 131 IF ( kt == nit000 ) THEN 132 133 istep(:) = 0 134 135 nbdy1=0 136 nbdy2=0 137 138 ! 1.1 Get time information from bdy data file 139 ! ******************************************** 671 itimer = ((kt-1)-nit000+1)*rdt ! current time in seconds for interpolation 672 itimer = itimer + jit*rdtbt ! in non-climatological case 673 674 IF ( ln_bdy_tides ) THEN 675 676 ! -------------------------------------! 677 ! Update BDY fields with tidal forcing ! 678 ! -------------------------------------! 679 680 CALL tide_update( kt, jit ) 681 682 ENDIF 683 684 IF ( ln_bdy_dyn_fla ) THEN 685 686 ! -------------------------------------! 687 ! Update BDY fields with model data ! 688 ! -------------------------------------! 689 690 ! !-------------------! 691 IF( kt == nit000 ) THEN ! First call only ! 692 ! !-------------------! 693 istep_bt(:) = 0 694 nbdy_b_bt = 0 695 nbdy_a_bt = 0 696 697 ! Get time information from bdy data file 698 ! *************************************** 140 699 141 700 IF(lwp) WRITE(numout,*) 142 IF(lwp) WRITE(numout,*) 'bdy_dta :Initialize unstructured boundary data'701 IF(lwp) WRITE(numout,*) 'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.' 143 702 IF(lwp) WRITE(numout,*) '~~~~~~~' 144 703 145 IF (nbdy_dta == 0) THEN704 IF( nbdy_dta == 0 ) THEN 146 705 IF(lwp) WRITE(numout,*) 'Bdy data are taken from initial conditions' 147 706 … … 154 713 ! necessary time dumps in file are included 155 714 156 bdyfile(1) = filbdy_data_T 157 bdyfile(2) = filbdy_data_U 158 bdyfile(3) = filbdy_data_V 159 160 DO jgrd = 1,3 161 162 CALL iom_open( bdyfile(jgrd), inum ) 163 CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits ) 164 CALL iom_close( inum ) 715 clfile(1) = filbdy_data_bt_T 716 clfile(2) = filbdy_data_bt_U 717 clfile(3) = filbdy_data_bt_V 718 719 DO igrd = 1,3 720 721 CALL iom_open( clfile(igrd), inum ) 722 CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits ) 723 724 SELECT CASE( igrd ) 725 CASE (1) 726 numbdyt = inum 727 CASE (2) 728 numbdyu = inum 729 CASE (3) 730 numbdyv = inum 731 END SELECT 165 732 166 733 ! Calculate time offset … … 170 737 CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 171 738 ! Compute model initialization time 172 kyear = ndastp / 10000173 kmonth = ( ndastp - kyear * 10000 ) / 100174 kday = ndastp - kyear * 10000 - kmonth * 100175 ksecs = dayfrac * 86400176 CALL ymds2ju( kyear, kmonth, kday, real(ksecs) , zdayjulini)177 ! offset from initialization date:178 offset = (dayjul0-zdayjulini)*86400739 iyear = ndastp / 10000 740 imonth = ( ndastp - iyear * 10000 ) / 100 741 iday = ndastp - iyear * 10000 - imonth * 100 742 isecs = dayfrac * 86400 743 CALL ymds2ju(iyear, imonth, iday, real(isecs) , zdayjulini) 744 ! zoffset from initialization date: 745 zoffset = (dayjul0-zdayjulini)*86400 179 746 ! 180 747 181 748 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 182 183 749 184 750 !! TO BE DONE... Check consistency between calendar from file … … 186 752 !! when calendar in model available outside of IOIPSL. 187 753 188 IF(lwp) WRITE(numout,*) 'kdby (number of times): ',kbdy189 IF(lwp) WRITE(numout,*) 'offset: ',offset190 IF(lwp) WRITE(numout,*) 'totime: ',totime191 IF(lwp) WRITE(numout,*) 'istepr: ',istepr192 193 754 ! Check that there are not too many times in the file. 194 IF (kbdy > jpbtime) THEN 195 IF(lwp) WRITE(numout,*) 196 IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter' 197 IF(lwp) WRITE(numout,*) ' ========== jpbtime= ',jpbtime,' kbdy= ', kbdy 198 IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd) 199 nstop = nstop + 1 755 IF (ntimes_bdy_bt > jpbtime) CALL ctl_stop( & 756 'Number of time dumps in bdy file exceed jpbtime parameter', & 757 'Check file:' // TRIM(clfile(igrd)) ) 758 759 ! Check that time array increases (or interp will fail): 760 DO it = 2, ntimes_bdy 761 IF ( zstepr(it-1) >= zstepr(it) ) THEN 762 CALL ctl_stop('Time array in unstructured boundary data file', & 763 'does not continuously increase.', & 764 'Check file:' // TRIM(clfile(igrd)) ) 765 EXIT 766 END IF 767 END DO 768 769 IF ( .NOT. ln_bdy_clim ) THEN 770 ! Check that times in file span model run time: 771 772 ! Note: the fields may be time means, so we allow nit000 to be before 773 ! first time in the file, provided that it falls inside the meaning 774 ! period of the first field. Until we can get the meaning period 775 ! from the file, use the interval between fields as a proxy. 776 ! If nit000 is before the first time, use the value at first time 777 ! instead of extrapolating. This is done by putting time 1 into 778 ! both time levels. 779 ! The same applies to the last time level: see setting of lect below. 780 781 IF ( ntimes_bdy == 1 ) CALL ctl_stop( & 782 'There is only one time dump in data files', & 783 'Set ln_bdy_clim=.true. in namelist for constant bdy forcing.' ) 784 785 zinterval_s = zstepr(2) - zstepr(1) 786 zinterval_e = zstepr(ntimes_bdy) - zstepr(ntimes_bdy-1) 787 788 IF ( zstepr(1) - zinterval_s / 2.0 > 0 ) THEN 789 IF(lwp) WRITE(numout,*) 'First bdy time relative to nit000:', zstepr(1) 790 IF(lwp) WRITE(numout,*) 'Interval between first two times: ', zinterval_s 791 CALL ctl_stop( 'First data time is after start of run', & 792 'by more than half a meaning period', & 793 'Check file: ' // TRIM(clfile(igrd)) ) 794 END IF 795 796 IF ( zstepr(ntimes_bdy) + zinterval_e / 2.0 < zntotime ) THEN 797 IF(lwp) WRITE(numout,*) 'Last bdy time relative to nit000:', zstepr(ntimes_bdy) 798 IF(lwp) WRITE(numout,*) 'Interval between last two times: ', zinterval_e 799 CALL ctl_stop( 'Last data time is before end of run', & 800 'by more than half a meaning period', & 801 'Check file: ' // TRIM(clfile(igrd)) ) 802 END IF 803 804 END IF ! .NOT. ln_bdy_clim 805 806 IF ( igrd .EQ. 1) THEN 807 ntimes_bdyt = ntimes_bdy_bt 808 zoffsett = zoffset 809 istept(:) = INT( zstepr(:) + zoffset ) 810 ELSE IF (igrd .EQ. 2) THEN 811 ntimes_bdyu = ntimes_bdy_bt 812 zoffsetu = zoffset 813 istepu(:) = INT( zstepr(:) + zoffset ) 814 ELSE IF (igrd .EQ. 3) THEN 815 ntimes_bdyv = ntimes_bdy_bt 816 zoffsetv = zoffset 817 istepv(:) = INT( zstepr(:) + zoffset ) 200 818 ENDIF 201 819 202 ! Check that time array increases:203 204 jt=1205 DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy-1)))206 jt=jt+1207 END DO208 209 IF ((jt.NE.kbdy).AND.(kbdy.GT.1)) THEN210 IF(lwp) WRITE(numout,*)211 IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files'212 IF(lwp) WRITE(numout,*) ' ========== does not continuously increase. Check file: '213 IF(lwp) WRITE(numout,*) bdyfile(jgrd)214 IF(lwp) WRITE(numout,*)215 nstop = nstop + 1216 ENDIF217 218 ! Check that times in file span model run time:219 220 IF (istepr(1)+offset > 0) THEN221 IF(lwp) WRITE(numout,*)222 IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time'223 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd)224 IF(lwp) WRITE(numout,*)225 nstop = nstop + 1226 END IF227 228 IF (istepr(kbdy)+offset < totime) THEN229 IF(lwp) WRITE(numout,*)230 IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time'231 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd)232 IF(lwp) WRITE(numout,*)233 nstop = nstop + 1234 END IF235 236 IF ( jgrd .EQ. 1) THEN237 kbdyt = kbdy238 offsett = offset239 istept(:) = INT( istepr(:) + offset )240 ELSE IF (jgrd .EQ. 2) THEN241 kbdyu = kbdy242 offsetu = offset243 istepu(:) = INT( istepr(:) + offset )244 ELSE IF (jgrd .EQ. 3) THEN245 kbdyv = kbdy246 offsetv = offset247 istepv(:) = INT( istepr(:) + offset )248 ENDIF249 250 820 ENDDO 251 821 252 822 ! Verify time consistency between files 253 823 254 IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN 255 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps' 256 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 257 nstop = nstop + 1 824 IF ( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN 825 CALL ctl_stop( & 826 'Time axis lengths differ between bdy data files', & 827 'Multiple time frequencies not implemented yet' ) 828 ELSE 829 ntimes_bdy_bt = ntimes_bdyt 258 830 ENDIF 259 kbdy = kbdyt 260 261 IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN 262 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin' 263 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 264 nstop = nstop + 1 831 832 IF (zoffsetu.NE.zoffsett .OR. zoffsetv.NE.zoffsett) THEN 833 CALL ctl_stop( & 834 'Bdy data files must have the same time origin', & 835 'Multiple time frequencies not implemented yet' ) 265 836 ENDIF 266 offset =offsett837 zoffset = zoffsett 267 838 268 839 !! Check that times are the same in the three files... HERE. 269 istep (:) = istept(:)840 istep_bt(:) = istept(:) 270 841 271 842 ! Check number of time dumps: 272 IF ((kbdy.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN273 IF(lwp) WRITE(numout,*)274 IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files'275 IF(lwp) WRITE(numout,*) ' ========== Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.'276 IF(lwp) WRITE(numout,*)277 nstop = nstop + 1278 ENDIF279 280 843 IF (ln_bdy_clim) THEN 281 IF ((kbdy.NE.1).AND.(kbdy.NE.12)) THEN 282 IF(lwp) WRITE(numout,*) 283 IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),' 284 IF(lwp) WRITE(numout,*) ' ========== bdy data files must contain 1 or 12 time dumps.' 285 IF(lwp) WRITE(numout,*) 286 nstop = nstop + 1 287 ELSEIF (kbdy.EQ.1 ) THEN 844 SELECT CASE ( ntimes_bdy_bt ) 845 CASE( 1 ) 288 846 IF(lwp) WRITE(numout,*) 289 847 IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 290 848 IF(lwp) WRITE(numout,*) 291 ELSEIF (kbdy.EQ.12) THEN849 CASE( 12 ) 292 850 IF(lwp) WRITE(numout,*) 293 851 IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 294 852 IF(lwp) WRITE(numout,*) 295 ENDIF 853 CASE DEFAULT 854 CALL ctl_stop( & 855 'For climatological boundary forcing (ln_bdy_clim=.true.),',& 856 'bdy data files must contain 1 or 12 time dumps.' ) 857 END SELECT 296 858 ENDIF 297 859 298 860 ! Find index of first record to read (before first model time). 299 861 300 jt=1301 DO WHILE ( ((istep (jt+1)) <= 0 ).AND.(jt.LE.(kbdy-1)))302 jt=jt+1303 END DO 304 nbdy 1 = jt305 306 WRITE(numout,*) 'Time offset is ', offset307 WRITE(numout,*) 'First record to read is ',nbdy 1862 it=1 863 DO WHILE ( ((istep_bt(it+1)) <= 0 ).AND.(it.LE.(ntimes_bdy_bt-1))) 864 it=it+1 865 END DO 866 nbdy_b_bt = it 867 868 WRITE(numout,*) 'Time offset is ',zoffset 869 WRITE(numout,*) 'First record to read is ',nbdy_b_bt 308 870 309 871 ENDIF ! endif (nbdy_dta == 1) … … 314 876 IF ( nbdy_dta == 0) THEN 315 877 ! boundary data arrays are filled with initial conditions 316 DO jk = 1, jpkm1 317 318 #if ! defined key_barotropic 319 jgrd = 1 ! T-points data 320 DO jb = 1, nblen(jgrd) 321 tbdy(jb,jk) = tn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 322 sbdy(jb,jk) = sn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 323 END DO 324 #endif 325 326 jgrd = 2 ! U-points data 327 DO jb = 1, nblen(jgrd) 328 ubdy(jb,jk) = un(nbi(jb,jgrd), nbj(jb,jgrd), jk) 329 END DO 330 331 jgrd = 3 ! V-points data 332 DO jb = 1, nblen(jgrd) 333 vbdy(jb,jk) = vn(nbi(jb,jgrd), nbj(jb,jgrd), jk) 334 END DO 335 END DO 878 igrd = 2 ! U-points data 879 DO ib = 1, nblen(igrd) 880 ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1) 881 END DO 882 883 igrd = 3 ! V-points data 884 DO ib = 1, nblen(igrd) 885 vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1) 886 END DO 887 888 igrd = 1 ! T-points data 889 DO ib = 1, nblen(igrd) 890 sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd)) 891 END DO 336 892 337 893 ELSEIF (nbdy_dta == 1) THEN 338 894 339 895 ! Set first record in the climatological case: 340 IF ((ln_bdy_clim).AND.( kbdy==1)) THEN341 nbdy 2= 1342 ELSEIF ((ln_bdy_clim).AND.( kbdy==iman)) THEN343 nbdy 1= 0344 nbdy 2= imois896 IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN 897 nbdy_a_bt = 1 898 ELSEIF ((ln_bdy_clim).AND.(ntimes_bdy_bt==iman)) THEN 899 nbdy_b_bt = 0 900 nbdy_a_bt = imois 345 901 ELSE 346 nbdy 2 = nbdy1902 nbdy_a_bt = nbdy_b_bt 347 903 END IF 348 904 349 905 ! Open Netcdf files: 350 906 351 CALL iom_open ( filbdy_data_ T, bdynumt )352 CALL iom_open ( filbdy_data_ U, bdynumu)353 CALL iom_open ( filbdy_data_ V, bdynumv)907 CALL iom_open ( filbdy_data_bt_T, numbdyt_bt ) 908 CALL iom_open ( filbdy_data_bt_U, numbdyu_bt ) 909 CALL iom_open ( filbdy_data_bt_V, numbdyv_bt ) 354 910 355 911 ! Read first record: 356 912 ipj=1 357 ipk=jpk 358 jgrd=1 359 ipi=nblendta(jgrd) 360 361 #if ! defined key_barotropic 362 !temperature 363 364 jgrd=1 365 IF ( nblendta(jgrd) .le. 0 ) THEN 366 idvar = iom_varid( bdynumt,'votemper' ) 367 nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar) 913 igrd=1 914 ipi=nblendta(igrd) 915 916 ! ssh 917 igrd=1 918 IF ( nblendta(igrd) .le. 0 ) THEN 919 idvar = iom_varid( numbdyt_bt,'sossheig' ) 920 nblendta(igrd) = iom_file(numbdyt_bt)%dimsz(1,idvar) 368 921 ENDIF 369 WRITE(numout,*) 'Dim size for votemper is ',nblendta(jgrd) 370 ipi=nblendta(jgrd) 371 372 CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 373 374 DO jb=1, nblen(jgrd) 375 DO jk=1,jpkm1 376 tbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 377 END DO 378 END DO 379 380 ! salinity 381 382 jgrd=1 383 IF ( nblendta(jgrd) .le. 0 ) THEN 384 idvar = iom_varid( bdynumt,'vosaline' ) 385 nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar) 922 WRITE(numout,*) 'Dim size for sossheig is ',nblendta(igrd) 923 ipi=nblendta(igrd) 924 925 CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 926 927 DO ib=1, nblen(igrd) 928 sshbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 929 END DO 930 931 ! u-velocity 932 igrd=2 933 IF ( nblendta(igrd) .le. 0 ) THEN 934 idvar = iom_varid( numbdyu_bt,'vobtcrtx' ) 935 nblendta(igrd) = iom_file(numbdyu_bt)%dimsz(1,idvar) 386 936 ENDIF 387 WRITE(numout,*) 'Dim size for vosaline is ',nblendta(jgrd) 388 ipi=nblendta(jgrd) 389 390 CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 391 392 DO jb=1, nblen(jgrd) 393 DO jk=1,jpkm1 394 sbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 395 END DO 396 END DO 397 #endif 398 399 ! u-velocity 400 401 jgrd=2 402 IF ( nblendta(jgrd) .le. 0 ) THEN 403 idvar = iom_varid( bdynumu,'vozocrtx' ) 404 nblendta(jgrd) = iom_file(bdynumu)%dimsz(1,idvar) 937 WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(igrd) 938 ipi=nblendta(igrd) 939 940 CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 941 942 DO ib=1, nblen(igrd) 943 ubtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 944 END DO 945 946 ! v-velocity 947 igrd=3 948 IF ( nblendta(igrd) .le. 0 ) THEN 949 idvar = iom_varid( numbdyv_bt,'vobtcrty' ) 950 nblendta(igrd) = iom_file(numbdyv_bt)%dimsz(1,idvar) 405 951 ENDIF 406 WRITE(numout,*) 'Dim size for vozocrtx is ',nblendta(jgrd) 407 ipi=nblendta(jgrd) 408 409 CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 410 411 DO jb=1, nblen(jgrd) 412 DO jk=1,jpkm1 413 ubdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 414 END DO 415 END DO 416 417 ! v-velocity 418 jgrd=3 419 IF ( nblendta(jgrd) .le. 0 ) THEN 420 idvar = iom_varid( bdynumv,'vomecrty' ) 421 nblendta(jgrd) = iom_file(bdynumv)%dimsz(1,idvar) 422 ENDIF 423 WRITE(numout,*) 'Dim size for vomecrty is ',nblendta(jgrd) 424 ipi=nblendta(jgrd) 425 426 CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 427 428 DO jb=1, nblen(jgrd) 429 DO jk=1,jpkm1 430 vbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 431 END DO 952 WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(igrd) 953 ipi=nblendta(igrd) 954 955 CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 956 957 DO ib=1, nblen(igrd) 958 vbtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 432 959 END DO 433 960 … … 435 962 436 963 ! In the case of constant boundary forcing fill bdy arrays once for all 437 IF ((ln_bdy_clim).AND.(kbdy==1)) THEN 438 #if ! defined key_barotropic 439 tbdy (:,:) = tbdydta (:,:,2) 440 sbdy (:,:) = sbdydta (:,:,2) 441 #endif 442 ubdy (:,:) = ubdydta (:,:,2) 443 vbdy (:,:) = vbdydta (:,:,2) 444 445 CALL iom_close( bdynumt ) 446 CALL iom_close( bdynumu ) 447 CALL iom_close( bdynumv ) 964 IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN 965 966 ubtbdy (:) = ubtbdydta (:,2) 967 vbtbdy (:) = vbtbdydta (:,2) 968 sshbdy (:) = sshbdydta (:,2) 969 970 CALL iom_close( numbdyt_bt ) 971 CALL iom_close( numbdyu_bt ) 972 CALL iom_close( numbdyv_bt ) 973 448 974 END IF 449 975 … … 454 980 ! -------------------- ! 455 981 456 IF ((nbdy_dta==1).AND.( kbdy>1)) THEN982 IF ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) THEN 457 983 458 984 ! 2.1 Read one more record if necessary 459 985 !************************************** 460 986 461 IF ( (ln_bdy_clim).AND.(imois/=nbdy 1) ) THEN ! remember that nbdy1=0 for kt=nit000462 nbdy 1= imois463 nbdy 2= imois+1464 nbdy 1 = MOD( nbdy1, iman )465 IF( nbdy 1 == 0 ) nbdy1= iman466 nbdy 2 = MOD( nbdy2, iman )467 IF( nbdy 2 == 0 ) nbdy2= iman987 IF ( (ln_bdy_clim).AND.(imois/=nbdy_b_bt) ) THEN ! remember that nbdy_b_bt=0 for kt=nit000 988 nbdy_b_bt = imois 989 nbdy_a_bt = imois+1 990 nbdy_b_bt = MOD( nbdy_b_bt, iman ) 991 IF( nbdy_b_bt == 0 ) nbdy_b_bt = iman 992 nbdy_a_bt = MOD( nbdy_a_bt, iman ) 993 IF( nbdy_a_bt == 0 ) nbdy_a_bt = iman 468 994 lect=.true. 469 995 470 ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep(nbdy2))) THEN 471 nbdy1=nbdy2 472 nbdy2=nbdy2+1 473 lect=.true. 474 END IF 475 476 IF (lect) THEN 477 478 ! Swap arrays 479 #if ! defined key_barotropic 480 tbdydta(:,:,1) = tbdydta(:,:,2) 481 sbdydta(:,:,1) = sbdydta(:,:,2) 482 #endif 483 ubdydta(:,:,1) = ubdydta(:,:,2) 484 vbdydta(:,:,1) = vbdydta(:,:,2) 485 486 ! read another set 487 488 ipj=1 489 ipk=jpk 490 jgrd=1 491 ipi=nblendta(jgrd) 492 493 #if ! defined key_barotropic 494 !temperature 495 496 jgrd=1 497 ipi=nblendta(jgrd) 498 499 CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 500 501 DO jb=1, nblen(jgrd) 502 DO jk=1,jpkm1 503 tbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 504 END DO 505 END DO 506 507 ! salinity 508 509 jgrd=1 510 ipi=nblendta(jgrd) 511 512 CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 513 514 DO jb=1, nblen(jgrd) 515 DO jk=1,jpkm1 516 sbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 517 END DO 518 END DO 519 #endif 520 521 ! u-velocity 522 523 jgrd=2 524 ipi=nblendta(jgrd) 525 526 CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 527 528 DO jb=1, nblen(jgrd) 529 DO jk=1,jpkm1 530 ubdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 531 END DO 532 END DO 533 534 ! v-velocity 535 jgrd=3 536 ipi=nblendta(jgrd) 537 538 CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 ) 539 540 DO jb=1, nblen(jgrd) 541 DO jk=1,jpkm1 542 vbdydta(jb,jk,2) = pdta3(nbmap(jb,jgrd),1,jk) 543 END DO 544 END DO 545 546 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy1 ',nbdy1 547 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy2 ',nbdy2 548 IF (.NOT.ln_bdy_clim) THEN 549 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep(nbdy1) 550 IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer 551 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy2) 552 ENDIF 553 END IF ! end lect=.true. 554 555 556 ! 2.2 Interpolate linearly: 557 ! *************************** 558 559 IF (ln_bdy_clim) THEN 560 zxy = FLOAT( nday ) / FLOAT( nobis(nbdy1) ) + 0.5 - i15 561 ELSE 562 zxy = FLOAT(istep(nbdy1)-itimer) / FLOAT(istep(nbdy1)-istep(nbdy2)) 563 END IF 564 565 #if ! defined key_barotropic 566 jgrd=1 567 DO jb=1, nblen(jgrd) 568 DO jk=1, jpkm1 569 tbdy(jb,jk) = zxy * tbdydta(jb,jk,2) + & 570 (1.-zxy) * tbdydta(jb,jk,1) 571 sbdy(jb,jk) = zxy * sbdydta(jb,jk,2) + & 572 (1.-zxy) * sbdydta(jb,jk,1) 573 END DO 574 END DO 575 #endif 576 577 jgrd=2 578 DO jb=1, nblen(jgrd) 579 DO jk=1, jpkm1 580 ubdy(jb,jk) = zxy * ubdydta(jb,jk,2) + & 581 (1.-zxy) * ubdydta(jb,jk,1) 582 END DO 583 END DO 584 585 jgrd=3 586 DO jb=1, nblen(jgrd) 587 DO jk=1, jpkm1 588 vbdy(jb,jk) = zxy * vbdydta(jb,jk,2) + & 589 (1.-zxy) * vbdydta(jb,jk,1) 590 END DO 591 END DO 592 593 END IF !end if ((nbdy_dta==1).AND.(kbdy>1)) 594 595 ! ------------------- ! 596 ! Last call kt=nitend ! 597 ! ------------------- ! 598 599 ! Closing of the 3 files 600 IF( kt == nitend ) THEN 601 CALL iom_close( bdynumt ) 602 CALL iom_close( bdynumu ) 603 CALL iom_close( bdynumv ) 604 ENDIF 605 606 END SUBROUTINE bdy_dta 607 608 SUBROUTINE bdy_dta_bt( kt, jit ) 609 !!--------------------------------------------------------------------------- 610 !! *** SUBROUTINE bdy_dta_bt *** 611 !! 612 !! ** Purpose : Read unstructured boundary data for barotropic variables 613 !! 614 !! ** Method : 615 !! 616 !! History : OPA 9.0 ! 07-07 (D. Storkey) Original 617 !!--------------------------------------------------------------------------- 618 !! * Arguments 619 INTEGER, INTENT( in ) :: kt ! ocean time-step index 620 INTEGER, INTENT( in ) :: jit ! barotropic time step index 621 ! (for timesplitting option, otherwise zero) 622 623 !! * Local declarations 624 LOGICAL :: lect ! flag for reading 625 INTEGER :: jt, jb, jgrd ! dummy loop indices 626 INTEGER :: idvar ! netcdf var ID 627 INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing 628 INTEGER :: kbdyt, kbdyu, kbdyv 629 INTEGER :: itimer, totime 630 INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read) 631 INTEGER :: iyear0, imonth0, iday0 632 INTEGER :: ihours0, iminutes0, isec0 633 INTEGER :: kyear, kmonth, kday, ksecs 634 INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files 635 REAL(wp) :: dayfrac, zxy, offsett 636 REAL(wp) :: offsetu, offsetv 637 REAL(wp) :: dayjul0, zdayjulini 638 REAL(wp), DIMENSION(jpbtime) :: istepr ! REAL time array from data files 639 REAL(wp), DIMENSION(jpbdta,1) :: pdta2 ! temporary array for data fields 640 REAL(wp), DIMENSION(jpbdta,1,jpk) :: pdta3 ! temporary array for data fields 641 CHARACTER(LEN=80), DIMENSION(3) :: bdyfile 642 CHARACTER(LEN=70 ) :: clunits ! units attribute of time coordinate 643 644 !!--------------------------------------------------------------------------- 645 646 ! -------------------- ! 647 ! Initialization ! 648 ! -------------------- ! 649 650 lect = .false. ! If true, read a time record 651 652 ! Some time variables for monthly climatological forcing: 653 ! ******************************************************* 654 655 iman = INT( raamo ) ! Number of months in a year 656 657 i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 658 ! i15=0 if the current day is in the first half of the month, else i15=1 659 660 imois = nmonth + i15 - 1 ! imois is the first month record 661 IF( imois == 0 ) imois = iman 662 663 ! Time variable for non-climatological forcing: 664 ! ********************************************* 665 666 itimer = ((kt-1)-nit000+1)*rdt ! current time in seconds for interpolation 667 itimer = itimer + jit*rdtbt ! in non-climatological case 668 669 ! -------------------------------------! 670 ! Update BDY fields with tidal forcing ! 671 ! -------------------------------------! 672 673 IF ( lk_bdy_tides ) CALL tide_update( kt, jit ) 674 675 IF ( lk_bdy ) THEN 676 677 ! -------------------- ! 678 ! 1. First call only ! 679 ! -------------------- ! 680 681 IF ( kt == nit000 ) THEN 682 683 istep_bt(:) = 0 684 685 nbdy1_bt=0 686 nbdy2_bt=0 687 688 ! 1.1 Get time information from bdy data file 689 ! ******************************************** 690 691 IF(lwp) WRITE(numout,*) 692 IF(lwp) WRITE(numout,*) 'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.' 693 IF(lwp) WRITE(numout,*) '~~~~~~~' 694 695 IF (nbdy_dta == 0) THEN 696 IF(lwp) WRITE(numout,*) 'Bdy data are taken from initial conditions' 697 698 ELSEIF (nbdy_dta == 1) THEN 699 IF(lwp) WRITE(numout,*) 'Bdy data are read in netcdf files' 700 701 dayfrac = adatrj-float(itimer)/86400. ! day fraction at time step kt-1 702 dayfrac = dayfrac - INT(dayfrac) ! 703 totime = (nitend-nit000+1)*rdt ! Total time of the run to verify that all the 704 ! necessary time dumps in file are included 705 706 bdyfile(1) = filbdy_data_bt_T 707 bdyfile(2) = filbdy_data_bt_U 708 bdyfile(3) = filbdy_data_bt_V 709 710 DO jgrd = 1,3 711 712 CALL iom_open( bdyfile(jgrd), inum ) 713 CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits ) 714 CALL iom_close( inum ) 715 716 ! Calculate time offset 717 READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0 718 ! Convert time origin in file to julian days 719 isec0 = isec0 + ihours0*60.*60. + iminutes0*60. 720 CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0) 721 ! Compute model initialization time 722 kyear = ndastp / 10000 723 kmonth = ( ndastp - kyear * 10000 ) / 100 724 kday = ndastp - kyear * 10000 - kmonth * 100 725 ksecs = dayfrac * 86400 726 CALL ymds2ju(kyear, kmonth, kday, real(ksecs) , zdayjulini) 727 ! offset from initialization date: 728 offset = (dayjul0-zdayjulini)*86400 729 ! 730 731 7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2) 732 733 !! TO BE DONE... Check consistency between calendar from file 734 !! (available optionally from iom_gettime) and calendar in model 735 !! when calendar in model available outside of IOIPSL. 736 737 IF(lwp) WRITE(numout,*) 'kdby (number of times): ',kbdy_bt 738 IF(lwp) WRITE(numout,*) 'offset: ',offset 739 IF(lwp) WRITE(numout,*) 'totime: ',totime 740 IF(lwp) WRITE(numout,*) 'istepr: ',istepr 741 742 ! Check that there are not too many times in the file. 743 IF (kbdy_bt > jpbtime) THEN 744 IF(lwp) WRITE(numout,*) 745 IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter' 746 IF(lwp) WRITE(numout,*) ' ========== jpbtime= ',jpbtime,' kbdy_bt= ', kbdy_bt 747 IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd) 748 nstop = nstop + 1 749 ENDIF 750 751 ! Check that time array increases: 752 753 jt=1 754 DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy_bt-1))) 755 jt=jt+1 756 END DO 757 758 IF ((jt.NE.kbdy_bt).AND.(kbdy_bt.GT.1)) THEN 759 IF(lwp) WRITE(numout,*) 760 IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files' 761 IF(lwp) WRITE(numout,*) ' ========== does not continuously increase. Check file: ' 762 IF(lwp) WRITE(numout,*) bdyfile(jgrd) 763 IF(lwp) WRITE(numout,*) 764 nstop = nstop + 1 765 ENDIF 766 767 ! Check that times in file span model run time: 768 769 IF (istepr(1)+offset > 0) THEN 770 IF(lwp) WRITE(numout,*) 771 IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time' 772 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd) 773 IF(lwp) WRITE(numout,*) 774 nstop = nstop + 1 775 END IF 776 777 IF (istepr(kbdy_bt)+offset < totime) THEN 778 IF(lwp) WRITE(numout,*) 779 IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time' 780 IF(lwp) WRITE(numout,*) ' ========== Check file: ',bdyfile(jgrd) 781 IF(lwp) WRITE(numout,*) 782 nstop = nstop + 1 783 END IF 784 785 IF ( jgrd .EQ. 1) THEN 786 kbdyt = kbdy_bt 787 offsett = offset 788 istept(:) = INT( istepr(:) + offset ) 789 ELSE IF (jgrd .EQ. 2) THEN 790 kbdyu = kbdy_bt 791 offsetu = offset 792 istepu(:) = INT( istepr(:) + offset ) 793 ELSE IF (jgrd .EQ. 3) THEN 794 kbdyv = kbdy_bt 795 offsetv = offset 796 istepv(:) = INT( istepr(:) + offset ) 797 ENDIF 798 799 ENDDO 800 801 ! Verify time consistency between files 802 803 IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN 804 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps' 805 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 806 nstop = nstop + 1 807 ENDIF 808 kbdy_bt = kbdyt 809 810 IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN 811 IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin' 812 IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet' 813 nstop = nstop + 1 814 ENDIF 815 offset = offsett 816 817 !! Check that times are the same in the three files... HERE. 818 istep_bt(:) = istept(:) 819 820 ! Check number of time dumps: 821 IF ((kbdy_bt.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN 822 IF(lwp) WRITE(numout,*) 823 IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files' 824 IF(lwp) WRITE(numout,*) ' ========== Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' 825 IF(lwp) WRITE(numout,*) 826 nstop = nstop + 1 827 ENDIF 828 829 IF (ln_bdy_clim) THEN 830 IF ((kbdy_bt.NE.1).AND.(kbdy_bt.NE.12)) THEN 831 IF(lwp) WRITE(numout,*) 832 IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),' 833 IF(lwp) WRITE(numout,*) ' ========== bdy data files must contain 1 or 12 time dumps.' 834 IF(lwp) WRITE(numout,*) 835 nstop = nstop + 1 836 ELSEIF (kbdy_bt.EQ.1 ) THEN 837 IF(lwp) WRITE(numout,*) 838 IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files' 839 IF(lwp) WRITE(numout,*) 840 ELSEIF (kbdy_bt.EQ.12) THEN 841 IF(lwp) WRITE(numout,*) 842 IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files' 843 IF(lwp) WRITE(numout,*) 844 ENDIF 845 ENDIF 846 847 ! Find index of first record to read (before first model time). 848 849 jt=1 850 DO WHILE ( ((istep_bt(jt+1)) <= 0 ).AND.(jt.LE.(kbdy_bt-1))) 851 jt=jt+1 852 END DO 853 nbdy1_bt = jt 854 855 WRITE(numout,*) 'Time offset is ',offset 856 WRITE(numout,*) 'First record to read is ',nbdy1_bt 857 858 ENDIF ! endif (nbdy_dta == 1) 859 860 ! 1.2 Read first record in file if necessary (ie if nbdy_dta == 1) 861 ! ***************************************************************** 862 863 IF ( nbdy_dta == 0) THEN 864 ! boundary data arrays are filled with initial conditions 865 jgrd = 2 ! U-points data 866 DO jb = 1, nblen(jgrd) 867 ubtbdy(jb) = un(nbi(jb,jgrd), nbj(jb,jgrd), 1) 868 END DO 869 870 jgrd = 3 ! V-points data 871 DO jb = 1, nblen(jgrd) 872 vbtbdy(jb) = vn(nbi(jb,jgrd), nbj(jb,jgrd), 1) 873 END DO 874 875 jgrd = 1 ! T-points data 876 DO jb = 1, nblen(jgrd) 877 sshbdy(jb) = sshn(nbi(jb,jgrd), nbj(jb,jgrd)) 878 END DO 879 880 ELSEIF (nbdy_dta == 1) THEN 881 882 ! Set first record in the climatological case: 883 IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN 884 nbdy2_bt = 1 885 ELSEIF ((ln_bdy_clim).AND.(kbdy_bt==iman)) THEN 886 nbdy1_bt = 0 887 nbdy2_bt = imois 888 ELSE 889 nbdy2_bt = nbdy1_bt 890 END IF 891 892 ! Open Netcdf files: 893 894 CALL iom_open ( filbdy_data_bt_T, bdynumt_bt ) 895 CALL iom_open ( filbdy_data_bt_U, bdynumu_bt ) 896 CALL iom_open ( filbdy_data_bt_V, bdynumv_bt ) 897 898 ! Read first record: 899 ipj=1 900 jgrd=1 901 ipi=nblendta(jgrd) 902 903 ! ssh 904 jgrd=1 905 IF ( nblendta(jgrd) .le. 0 ) THEN 906 idvar = iom_varid( bdynumt_bt,'sossheig' ) 907 nblendta(jgrd) = iom_file(bdynumt_bt)%dimsz(1,idvar) 908 ENDIF 909 WRITE(numout,*) 'Dim size for sossheig is ',nblendta(jgrd) 910 ipi=nblendta(jgrd) 911 912 CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt ) 913 914 DO jb=1, nblen(jgrd) 915 sshbdydta(jb,2) = pdta2(nbmap(jb,jgrd),1) 916 END DO 917 918 ! u-velocity 919 jgrd=2 920 IF ( nblendta(jgrd) .le. 0 ) THEN 921 idvar = iom_varid( bdynumu_bt,'vobtcrtx' ) 922 nblendta(jgrd) = iom_file(bdynumu_bt)%dimsz(1,idvar) 923 ENDIF 924 WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(jgrd) 925 ipi=nblendta(jgrd) 926 927 CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt ) 928 929 DO jb=1, nblen(jgrd) 930 ubtbdydta(jb,2) = pdta2(nbmap(jb,jgrd),1) 931 END DO 932 933 ! v-velocity 934 jgrd=3 935 IF ( nblendta(jgrd) .le. 0 ) THEN 936 idvar = iom_varid( bdynumv_bt,'vobtcrty' ) 937 nblendta(jgrd) = iom_file(bdynumv_bt)%dimsz(1,idvar) 938 ENDIF 939 WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(jgrd) 940 ipi=nblendta(jgrd) 941 942 CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt ) 943 944 DO jb=1, nblen(jgrd) 945 vbtbdydta(jb,2) = pdta2(nbmap(jb,jgrd),1) 946 END DO 947 948 END IF 949 950 ! In the case of constant boundary forcing fill bdy arrays once for all 951 IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN 952 953 ubtbdy (:) = ubtbdydta (:,2) 954 vbtbdy (:) = vbtbdydta (:,2) 955 sshbdy (:) = sshbdydta (:,2) 956 957 CALL iom_close( bdynumt_bt ) 958 CALL iom_close( bdynumu_bt ) 959 CALL iom_close( bdynumv_bt ) 960 961 END IF 962 963 ENDIF ! End if nit000 964 965 ! -------------------- ! 966 ! 2. At each time step ! 967 ! -------------------- ! 968 969 IF ((nbdy_dta==1).AND.(kbdy_bt>1)) THEN 970 971 ! 2.1 Read one more record if necessary 972 !************************************** 973 974 IF ( (ln_bdy_clim).AND.(imois/=nbdy1_bt) ) THEN ! remember that nbdy1_bt=0 for kt=nit000 975 nbdy1_bt = imois 976 nbdy2_bt = imois+1 977 nbdy1_bt = MOD( nbdy1_bt, iman ) 978 IF( nbdy1_bt == 0 ) nbdy1_bt = iman 979 nbdy2_bt = MOD( nbdy2_bt, iman ) 980 IF( nbdy2_bt == 0 ) nbdy2_bt = iman 981 lect=.true. 982 983 ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy2_bt))) THEN 984 nbdy1_bt=nbdy2_bt 985 nbdy2_bt=nbdy2_bt+1 996 ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy_a_bt))) THEN 997 nbdy_b_bt=nbdy_a_bt 998 nbdy_a_bt=nbdy_a_bt+1 986 999 lect=.true. 987 1000 END IF … … 998 1011 ipj=1 999 1012 ipk=jpk 1000 jgrd=11001 ipi=nblendta( jgrd)1013 igrd=1 1014 ipi=nblendta(igrd) 1002 1015 1003 1016 1004 1017 ! ssh 1005 jgrd=11006 ipi=nblendta( jgrd)1007 1008 CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt )1009 1010 DO jb=1, nblen(jgrd)1011 sshbdydta( jb,2) = pdta2(nbmap(jb,jgrd),1)1018 igrd=1 1019 ipi=nblendta(igrd) 1020 1021 CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt ) 1022 1023 DO ib=1, nblen(igrd) 1024 sshbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 1012 1025 END DO 1013 1026 1014 1027 ! u-velocity 1015 jgrd=21016 ipi=nblendta( jgrd)1017 1018 CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt )1019 1020 DO jb=1, nblen(jgrd)1021 ubtbdydta( jb,2) = pdta3(nbmap(jb,jgrd),1)1028 igrd=2 1029 ipi=nblendta(igrd) 1030 1031 CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt ) 1032 1033 DO ib=1, nblen(igrd) 1034 ubtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 1022 1035 END DO 1023 1036 1024 1037 ! v-velocity 1025 jgrd=31026 ipi=nblendta( jgrd)1027 1028 CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt )1029 1030 DO jb=1, nblen(jgrd)1031 vbtbdydta( jb,2) = pdta3(nbmap(jb,jgrd),1)1032 END DO 1033 1034 1035 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy 1_bt ',nbdy1_bt1036 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy 2_bt ',nbdy2_bt1038 igrd=3 1039 ipi=nblendta(igrd) 1040 1041 CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt ) 1042 1043 DO ib=1, nblen(igrd) 1044 vbtbdydta(ib,2) = zdta(nbmap(ib,igrd),1) 1045 END DO 1046 1047 1048 IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b_bt ',nbdy_b_bt 1049 IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy_a_bt ',nbdy_a_bt 1037 1050 IF (.NOT.ln_bdy_clim) THEN 1038 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep_bt(nbdy 1_bt)1051 IF(lwp) WRITE(numout,*) 'first record time (s): ', istep_bt(nbdy_b_bt) 1039 1052 IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer 1040 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy 2_bt)1053 IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt) 1041 1054 ENDIF 1042 1055 END IF ! end lect=.true. … … 1047 1060 1048 1061 IF (ln_bdy_clim) THEN 1049 zxy = FLOAT( nday ) / FLOAT( nobis(nbdy 1_bt) ) + 0.5 - i151062 zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b_bt) ) + 0.5 - i15 1050 1063 ELSE 1051 zxy = FLOAT(istep_bt(nbdy 1_bt)-itimer) / FLOAT(istep_bt(nbdy1_bt)-istep_bt(nbdy2_bt))1064 zxy = FLOAT(istep_bt(nbdy_b_bt)-itimer) / FLOAT(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt)) 1052 1065 END IF 1053 1066 1054 jgrd=11055 DO jb=1, nblen(jgrd)1056 sshbdy( jb) = zxy * sshbdydta(jb,2) + &1057 (1.-zxy) * sshbdydta( jb,1)1058 END DO 1059 1060 jgrd=21061 DO jb=1, nblen(jgrd)1062 ubtbdy( jb) = zxy * ubtbdydta(jb,2) + &1063 (1.-zxy) * ubtbdydta( jb,1)1064 END DO 1065 1066 jgrd=31067 DO jb=1, nblen(jgrd)1068 vbtbdy( jb) = zxy * vbtbdydta(jb,2) + &1069 (1.-zxy) * vbtbdydta( jb,1)1070 END DO 1071 1072 1073 END IF !end if ((nbdy_dta==1).AND.( kbdy_bt>1))1067 igrd=1 1068 DO ib=1, nblen(igrd) 1069 sshbdy(ib) = zxy * sshbdydta(ib,2) + & 1070 (1.-zxy) * sshbdydta(ib,1) 1071 END DO 1072 1073 igrd=2 1074 DO ib=1, nblen(igrd) 1075 ubtbdy(ib) = zxy * ubtbdydta(ib,2) + & 1076 (1.-zxy) * ubtbdydta(ib,1) 1077 END DO 1078 1079 igrd=3 1080 DO ib=1, nblen(igrd) 1081 vbtbdy(ib) = zxy * vbtbdydta(ib,2) + & 1082 (1.-zxy) * vbtbdydta(ib,1) 1083 END DO 1084 1085 1086 END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) 1074 1087 1075 1088 ! ------------------- ! … … 1079 1092 ! Closing of the 3 files 1080 1093 IF( kt == nitend ) THEN 1081 CALL iom_close( bdynumt_bt )1082 CALL iom_close( bdynumu_bt )1083 CALL iom_close( bdynumv_bt )1094 CALL iom_close( numbdyt_bt ) 1095 CALL iom_close( numbdyu_bt ) 1096 CALL iom_close( numbdyv_bt ) 1084 1097 ENDIF 1085 1098 1086 ENDIF ! l k_bdy1099 ENDIF ! ln_bdy_dyn_frs 1087 1100 1088 1101 END SUBROUTINE bdy_dta_bt … … 1090 1103 1091 1104 #else 1092 !!---------------------------------------------------------------------- --------1093 !! default option: Dummy moduleNO Unstruct Open Boundary Conditions1094 !!---------------------------------------------------------------------- --------1105 !!---------------------------------------------------------------------- 1106 !! Dummy module NO Unstruct Open Boundary Conditions 1107 !!---------------------------------------------------------------------- 1095 1108 CONTAINS 1096 SUBROUTINE bdy_dta( kt ) ! Dummy routine 1097 INTEGER, INTENT (in) :: kt 1109 SUBROUTINE bdy_dta( kt ) ! Empty routine 1098 1110 WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt 1099 1111 END SUBROUTINE bdy_dta 1100 1101 SUBROUTINE bdy_dta_bt( kt, jit ) ! Dummy routine 1102 INTEGER, INTENT (in) :: kt, jit 1103 WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, jit 1112 SUBROUTINE bdy_dta_bt( kt, kit ) ! Empty routine 1113 WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, kit 1104 1114 END SUBROUTINE bdy_dta_bt 1105 1115 #endif
Note: See TracChangeset
for help on using the changeset viewer.