Changeset 12377 for NEMO/trunk/src/OCE/SBC/fldread.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/SBC/fldread.F90
r12367 r12377 13 13 !! fld_read : read input fields used for the computation of the surface boundary condition 14 14 !! fld_init : initialization of field read 15 !! fld_ rec : determined the record(s) to be read15 !! fld_def : define the record(s) of the file and its name 16 16 !! fld_get : read the data 17 17 !! fld_map : read global data from file and map onto local data using a general mapping (use for open boundaries) 18 18 !! fld_rot : rotate the vector fields onto the local grid direction 19 !! fld_clopn : update the data file name andclose/open the files19 !! fld_clopn : close/open the files 20 20 !! fld_fill : fill the data structure with the associated information read in namelist 21 21 !! wgt_list : manage the weights used for interpolation … … 25 25 !! seaoverland : create shifted matrices for seaoverland application 26 26 !! fld_interp : apply weights to input gridded data to create data on model grid 27 !! ksec_week : function returning the first 3 letters of the first day of the weekly file 27 !! fld_filename : define the filename according to a given date 28 !! ksec_week : function returning seconds between 00h of the beginning of the week and half of the current time step 28 29 !!---------------------------------------------------------------------- 29 30 USE oce ! ocean dynamics and tracers … … 44 45 PUBLIC fld_map ! routine called by tides_init 45 46 PUBLIC fld_read, fld_fill ! called by sbc... modules 46 PUBLIC fld_ clopn47 PUBLIC fld_def 47 48 48 49 TYPE, PUBLIC :: FLD_N !: Namelist field informations … … 72 73 INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 73 74 INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year) 74 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,: ) :: fnow ! input fields interpolated to now time step 75 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta ! 2 consecutive record of input fields 75 INTEGER , ALLOCATABLE, DIMENSION(: ) :: nrecsec ! 76 REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) :: fnow ! input fields interpolated to now time step 77 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta ! 2 consecutive record of input fields 76 78 CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key 77 79 ! ! into the WGTLIST structure … … 118 120 TYPE( WGT ), DIMENSION(tot_wgts) :: ref_wgts ! array of wgts 119 121 INTEGER :: nxt_wgt = 1 ! point to next available space in ref_wgts array 122 INTEGER :: nflag = 0 120 123 REAL(wp), PARAMETER :: undeff_lsm = -999.00_wp 121 124 122 125 !$AGRIF_END_DO_NOT_TREAT 123 126 127 !! * Substitutions 128 # include "do_loop_substitute.h90" 124 129 !!---------------------------------------------------------------------- 125 130 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 129 134 CONTAINS 130 135 131 SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, kt_offset)136 SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, pt_offset, Kmm ) 132 137 !!--------------------------------------------------------------------- 133 138 !! *** ROUTINE fld_read *** … … 145 150 TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables 146 151 INTEGER , INTENT(in ), OPTIONAL :: kit ! subcycle timestep for timesplitting option 147 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! provide fields at time other than "now" 148 ! ! kt_offset = -1 => fields at "before" time level 149 ! ! kt_offset = +1 => fields at "after" time level 150 ! ! etc. 151 !! 152 INTEGER :: itmp ! local variable 152 REAL(wp) , INTENT(in ), OPTIONAL :: pt_offset ! provide fields at time other than "now" 153 INTEGER , INTENT(in ), OPTIONAL :: Kmm ! ocean time level index 154 !! 153 155 INTEGER :: imf ! size of the structure sd 154 156 INTEGER :: jf ! dummy indices 155 INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend156 157 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 157 INTEGER :: it_offset ! local time offset variable158 LOGICAL :: llnxtyr ! open next year file?159 LOGICAL :: llnxtmth ! open next month file?160 LOGICAL :: llstop ! stop is the file does not exist161 158 LOGICAL :: ll_firstcall ! true if this is the first call to fld_read for this set of fields 159 REAL(wp) :: zt_offset ! local time offset variable 162 160 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 163 161 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation … … 167 165 IF( PRESENT(kit) ) ll_firstcall = ll_firstcall .and. kit == 1 168 166 169 IF ( nn_components == jp_iam_sas ) THEN ; it_offset = nn_fsbc170 ELSE ; it_offset = 0171 ENDIF 172 IF( PRESENT( kt_offset) ) it_offset = kt_offset173 174 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar175 IF( present(kit) ) THEN ! ignore kn_fsbc in this case176 isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) )167 IF( nn_components == jp_iam_sas ) THEN ; zt_offset = REAL( nn_fsbc, wp ) 168 ELSE ; zt_offset = 0. 169 ENDIF 170 IF( PRESENT(pt_offset) ) zt_offset = pt_offset 171 172 ! Note that all varibles starting by nsec_* are shifted time by +1/2 time step to be centrered 173 IF( PRESENT(kit) ) THEN ! ignore kn_fsbc in this case 174 isecsbc = nsec_year + nsec1jan000 + NINT( ( REAL( kit,wp) + zt_offset ) * rdt / REAL(nn_baro,wp) ) 177 175 ELSE ! middle of sbc time step 178 isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdt) + it_offset * NINT(rdt) 176 ! note: we use kn_fsbc-1 because nsec_year is defined at the middle of the current time step 177 isecsbc = nsec_year + nsec1jan000 + NINT( ( 0.5*REAL(kn_fsbc-1,wp) + zt_offset ) * rdt ) 179 178 ENDIF 180 179 imf = SIZE( sd ) … … 183 182 DO jf = 1, imf 184 183 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 185 CALL fld_init( kn_fsbc, sd(jf) ) ! read each before field (put them in after as they will be swapped)184 CALL fld_init( isecsbc, sd(jf) ) ! read each before field (put them in after as they will be swapped) 186 185 END DO 187 186 IF( lwp ) CALL wgt_print() ! control print … … 192 191 ! 193 192 DO jf = 1, imf ! --- loop over field --- ! 194 193 ! 195 194 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 196 197 IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN ! read/update the after data? 198 199 sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:) ! swap before record informations 200 sd(jf)%rotn(1) = sd(jf)%rotn(2) ! swap before rotate informations 201 IF( sd(jf)%ln_tint ) sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! swap before record field 202 203 CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit ) ! update after record informations 204 205 ! if kn_fsbc*rdt is larger than freqh (which is kind of odd), 206 ! it is possible that the before value is no more the good one... we have to re-read it 207 ! if before is not the last record of the file currently opened and after is the first record to be read 208 ! in a new file which means after = 1 (the file to be opened corresponds to the current time) 209 ! or after = nreclast + 1 (the file to be opened corresponds to a future time step) 210 IF( .NOT. ll_firstcall .AND. sd(jf)%ln_tint .AND. sd(jf)%nrec_b(1) /= sd(jf)%nreclast & 211 & .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) == 1 ) THEN 212 itmp = sd(jf)%nrec_a(1) ! temporary storage 213 sd(jf)%nrec_a(1) = sd(jf)%nreclast ! read the last record of the file currently opened 214 CALL fld_get( sd(jf) ) ! read after data 215 sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field 216 sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations 217 sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%freqh * 3600. ) ! assume freq to be in hours in this case 218 sd(jf)%rotn(1) = sd(jf)%rotn(2) ! update before rotate informations 219 sd(jf)%nrec_a(1) = itmp ! move back to after record 220 ENDIF 221 222 CALL fld_clopn( sd(jf) ) ! Do we need to open a new year/month/week/day file? 223 224 IF( sd(jf)%ln_tint ) THEN 225 226 ! if kn_fsbc*rdt is larger than freqh (which is kind of odd), 227 ! it is possible that the before value is no more the good one... we have to re-read it 228 ! if before record is not just just before the after record... 229 IF( .NOT. ll_firstcall .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) /= 1 & 230 & .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN 231 sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1 ! move back to before record 232 CALL fld_get( sd(jf) ) ! read after data 233 sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field 234 sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations 235 sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%freqh * 3600. ) ! assume freq to be in hours in this case 236 sd(jf)%rotn(1) = sd(jf)%rotn(2) ! update before rotate informations 237 sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) + 1 ! move back to after record 238 ENDIF 239 ENDIF ! temporal interpolation? 240 241 ! do we have to change the year/month/week/day of the forcing field?? 242 ! if we do time interpolation we will need to open next year/month/week/day file before the end of the current 243 ! one. If so, we are still before the end of the year/month/week/day when calling fld_rec so sd(jf)%nrec_a(1) 244 ! will be larger than the record number that should be read for current year/month/week/day 245 ! do we need next file data? 246 ! This applies to both cases with or without time interpolation 247 IF( sd(jf)%nrec_a(1) > sd(jf)%nreclast ) THEN 248 249 sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - sd(jf)%nreclast ! 250 251 IF( .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) ) THEN ! close/open the current/new file 252 253 llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth) ! open next month file? 254 llnxtyr = sd(jf)%cltype == 'yearly' .OR. (nmonth == 12 .AND. llnxtmth) ! open next year file? 255 256 ! if the run finishes at the end of the current year/month/week/day, we will allow next 257 ! year/month/week/day file to be not present. If the run continue further than the current 258 ! year/month/week/day, next year/month/week/day file must exist 259 isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt) ! second at the end of the run 260 llstop = isecend > sd(jf)%nrec_a(2) ! read more than 1 record of next year 261 ! we suppose that the date of next file is next day (should be ok even for weekly files...) 262 CALL fld_clopn( sd(jf), nyear + COUNT((/llnxtyr /)) , & 263 & nmonth + COUNT((/llnxtmth/)) - 12 * COUNT((/llnxtyr /)), & 264 & nday + 1 - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop ) 265 266 IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN ! next year file does not exist 267 CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)// & 268 & ' not present -> back to current year/month/day') 269 CALL fld_clopn( sd(jf) ) ! back to the current year/month/day 270 sd(jf)%nrec_a(1) = sd(jf)%nreclast ! force to read the last record in the current year file 271 ENDIF 272 273 ENDIF 274 ENDIF ! open need next file? 275 276 ! read after data 277 CALL fld_get( sd(jf) ) 278 279 ENDIF ! read new data? 195 CALL fld_update( isecsbc, sd(jf), Kmm ) 196 ! 280 197 END DO ! --- end loop over field --- ! 281 198 … … 292 209 WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, & 293 210 & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday 294 WRITE(numout, *) ' it_offset is : ',it_offset211 WRITE(numout, *) ' zt_offset is : ',zt_offset 295 212 ENDIF 296 213 ! temporal interpolation weights … … 316 233 317 234 318 SUBROUTINE fld_init( k n_fsbc, sdjf )235 SUBROUTINE fld_init( ksecsbc, sdjf ) 319 236 !!--------------------------------------------------------------------- 320 237 !! *** ROUTINE fld_init *** 321 238 !! 322 !! ** Purpose : - first call to fld_recto define before values323 !! - if time interpolation, read before data324 !!---------------------------------------------------------------------- 325 INTEGER , INTENT(in ) :: k n_fsbc ! sbc computation period (in time step)239 !! ** Purpose : - first call(s) to fld_def to define before values 240 !! - open file 241 !!---------------------------------------------------------------------- 242 INTEGER , INTENT(in ) :: ksecsbc ! 326 243 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 327 !! 328 LOGICAL :: llprevyr ! are we reading previous year file? 329 LOGICAL :: llprevmth ! are we reading previous month file? 330 LOGICAL :: llprevweek ! are we reading previous week file? 331 LOGICAL :: llprevday ! are we reading previous day file? 332 LOGICAL :: llprev ! llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday 333 INTEGER :: idvar ! variable id 334 INTEGER :: inrec ! number of record existing for this variable 335 INTEGER :: iyear, imonth, iday ! first day of the current file in yyyy mm dd 336 INTEGER :: isec_week ! number of seconds since start of the weekly file 337 CHARACTER(LEN=1000) :: clfmt ! write format 338 !!--------------------------------------------------------------------- 339 ! 340 llprevyr = .FALSE. 341 llprevmth = .FALSE. 342 llprevweek = .FALSE. 343 llprevday = .FALSE. 344 isec_week = 0 345 ! 346 ! define record informations 347 CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. ) ! return before values in sdjf%nrec_a (as we will swap it later) 348 ! 349 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 350 ! 351 IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure 352 ! 353 IF( sdjf%nrec_a(1) == 0 ) THEN ! we redefine record sdjf%nrec_a(1) with the last record of previous year file 354 IF ( NINT(sdjf%freqh) == -12 ) THEN ! yearly mean 355 IF( sdjf%cltype == 'yearly' ) THEN ! yearly file 356 sdjf%nrec_a(1) = 1 ! force to read the unique record 357 llprevyr = .NOT. sdjf%ln_clim ! use previous year file? 358 ELSE 359 CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clrootname) ) 360 ENDIF 361 ELSEIF( NINT(sdjf%freqh) == -1 ) THEN ! monthly mean 362 IF( sdjf%cltype == 'monthly' ) THEN ! monthly file 363 sdjf%nrec_a(1) = 1 ! force to read the unique record 364 llprevmth = .TRUE. ! use previous month file? 365 llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file? 366 ELSE ! yearly file 367 sdjf%nrec_a(1) = 12 ! force to read december mean 368 llprevyr = .NOT. sdjf%ln_clim ! use previous year file? 369 ENDIF 370 ELSE ! higher frequency mean (in hours) 371 IF ( sdjf%cltype == 'monthly' ) THEN ! monthly file 372 sdjf%nrec_a(1) = NINT( 24. * REAL(nmonth_len(nmonth-1),wp) / sdjf%freqh )! last record of previous month 373 llprevmth = .TRUE. ! use previous month file? 374 llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file? 375 ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ! weekly file 376 llprevweek = .TRUE. ! use previous week file? 377 sdjf%nrec_a(1) = NINT( 24. * 7. / sdjf%freqh ) ! last record of previous week 378 isec_week = NINT(rday) * 7 ! add a shift toward previous week 379 ELSEIF( sdjf%cltype == 'daily' ) THEN ! daily file 380 sdjf%nrec_a(1) = NINT( 24. / sdjf%freqh ) ! last record of previous day 381 llprevday = .TRUE. ! use previous day file? 382 llprevmth = llprevday .AND. nday == 1 ! use previous month file? 383 llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file? 384 ELSE ! yearly file 385 sdjf%nrec_a(1) = NINT( 24. * REAL(nyear_len(0),wp) / sdjf%freqh ) ! last record of previous year 386 llprevyr = .NOT. sdjf%ln_clim ! use previous year file? 387 ENDIF 388 ENDIF 389 ENDIF 390 ! 391 IF ( sdjf%cltype(1:4) == 'week' ) THEN 392 isec_week = isec_week + ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week 393 llprevmth = isec_week > nsec_month ! longer time since the beginning of the week than the month 394 llprevyr = llprevmth .AND. nmonth == 1 395 ENDIF 396 llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday 397 ! 398 iyear = nyear - COUNT((/llprevyr /)) 399 imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)) 400 iday = nday - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday) 401 ! 402 CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev ) 403 ! 404 ! if previous year/month/day file does not exist, we switch to the current year/month/day 405 IF( llprev .AND. sdjf%num <= 0 ) THEN 406 CALL ctl_warn( 'previous year/month/week/day file: '//TRIM(sdjf%clrootname)// & 407 & ' not present -> back to current year/month/week/day' ) 408 ! we force to read the first record of the current year/month/day instead of last record of previous year/month/day 409 llprev = .FALSE. 410 sdjf%nrec_a(1) = 1 411 CALL fld_clopn( sdjf ) 412 ENDIF 413 ! 414 IF( llprev ) THEN ! check if the record sdjf%nrec_a(1) exists in the file 415 idvar = iom_varid( sdjf%num, sdjf%clvar ) ! id of the variable sdjf%clvar 416 IF( idvar <= 0 ) RETURN 417 inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar ) ! size of the last dim of idvar 418 sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec ) ! make sure we select an existing record 419 ENDIF 420 ! 421 ! read before data in after arrays(as we will swap it later) 422 CALL fld_get( sdjf ) 423 ! 424 clfmt = "(' fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 425 IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday 426 ! 427 ENDIF 244 !!--------------------------------------------------------------------- 245 ! 246 IF( nflag == 0 ) nflag = -( HUGE(0) - 10 ) 247 ! 248 CALL fld_def( sdjf ) 249 IF( sdjf%ln_tint .AND. ksecsbc < sdjf%nrecsec(1) ) CALL fld_def( sdjf, ldprev = .TRUE. ) 250 ! 251 CALL fld_clopn( sdjf ) 252 sdjf%nrec_a(:) = (/ 1, nflag /) ! default definition to force flp_update to read the file. 428 253 ! 429 254 END SUBROUTINE fld_init 430 255 431 256 432 SUBROUTINE fld_ rec( kn_fsbc, sdjf, ldbefore, kit, kt_offset)433 !!--------------------------------------------------------------------- 434 !! *** ROUTINE fld_ rec***257 SUBROUTINE fld_update( ksecsbc, sdjf, Kmm ) 258 !!--------------------------------------------------------------------- 259 !! *** ROUTINE fld_update *** 435 260 !! 436 261 !! ** Purpose : Compute … … 441 266 !! nrec_b(2) and nrec_a(2): time of the beginning and end of the record 442 267 !!---------------------------------------------------------------------- 443 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 444 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 445 LOGICAL , INTENT(in ), OPTIONAL :: ldbefore ! sent back before record values (default = .FALSE.) 446 INTEGER , INTENT(in ), OPTIONAL :: kit ! index of barotropic subcycle 447 ! ! used only if sdjf%ln_tint = .TRUE. 448 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! Offset of required time level compared to "now" 449 ! ! time level in units of time steps. 450 ! 451 LOGICAL :: llbefore ! local definition of ldbefore 452 INTEGER :: iendrec ! end of this record (in seconds) 453 INTEGER :: imth ! month number 454 INTEGER :: ifreq_sec ! frequency mean (in seconds) 455 INTEGER :: isec_week ! number of seconds since the start of the weekly file 456 INTEGER :: it_offset ! local time offset variable 457 REAL(wp) :: ztmp ! temporary variable 458 !!---------------------------------------------------------------------- 459 ! 460 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 461 ! 462 IF( PRESENT(ldbefore) ) THEN ; llbefore = ldbefore .AND. sdjf%ln_tint ! needed only if sdjf%ln_tint = .TRUE. 463 ELSE ; llbefore = .FALSE. 464 ENDIF 465 ! 466 IF ( nn_components == jp_iam_sas ) THEN ; it_offset = nn_fsbc 467 ELSE ; it_offset = 0 468 ENDIF 469 IF( PRESENT(kt_offset) ) it_offset = kt_offset 470 IF( PRESENT(kit) ) THEN ; it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) ) 471 ELSE ; it_offset = it_offset * NINT( rdt ) 472 ENDIF 473 ! 474 ! ! =========== ! 475 IF ( NINT(sdjf%freqh) == -12 ) THEN ! yearly mean 476 ! ! =========== ! 477 ! 478 IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record 479 ! 480 ! INT( ztmp ) 481 ! /|\ 482 ! 1 | *---- 483 ! 0 |----( 484 ! |----+----|--> time 485 ! 0 /|\ 1 (nday/nyear_len(1)) 486 ! | 487 ! | 488 ! forcing record : 1 489 ! 490 ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 & 491 & + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday ) 492 sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) 493 ! swap at the middle of the year 494 IF( llbefore ) THEN ; sdjf%nrec_a(2) = nsec1jan000 - (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(0) + & 495 & INT(ztmp) * NINT( 0.5 * rday) * nyear_len(1) 496 ELSE ; sdjf%nrec_a(2) = nsec1jan000 + (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(1) + & 497 & INT(ztmp) * INT(rday) * nyear_len(1) + INT(ztmp) * NINT( 0.5 * rday) * nyear_len(2) 268 INTEGER , INTENT(in ) :: ksecsbc ! 269 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 270 INTEGER , OPTIONAL, INTENT(in ) :: Kmm ! ocean time level index 271 ! 272 INTEGER :: ja ! end of this record (in seconds) 273 !!---------------------------------------------------------------------- 274 ! 275 IF( ksecsbc > sdjf%nrec_a(2) ) THEN ! --> we need to update after data 276 277 ! find where is the new after record... (it is not necessary sdjf%nrec_a(1)+1 ) 278 ja = sdjf%nrec_a(1) 279 DO WHILE ( ksecsbc >= sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast ) ! Warning: make sure ja <= sdjf%nreclast in this test 280 ja = ja + 1 281 END DO 282 IF( ksecsbc > sdjf%nrecsec(ja) ) ja = ja + 1 ! in case ksecsbc > sdjf%nrecsec(sdjf%nreclast) 283 284 ! if ln_tint and if the new after is not ja+1, we need also to update after data before the swap 285 ! so, after the swap, sdjf%nrec_b(2) will still be the closest value located just before ksecsbc 286 IF( sdjf%ln_tint .AND. ( ja > sdjf%nrec_a(1) + 1 .OR. sdjf%nrec_a(2) == nflag ) ) THEN 287 sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec_a with before information 288 CALL fld_get( sdjf, Kmm ) ! read after data that will be used as before data 289 ENDIF 290 291 ! if after is in the next file... 292 IF( ja > sdjf%nreclast ) THEN 293 294 CALL fld_def( sdjf ) 295 IF( ksecsbc > sdjf%nrecsec(sdjf%nreclast) ) CALL fld_def( sdjf, ldnext = .TRUE. ) 296 CALL fld_clopn( sdjf ) ! open next file 297 298 ! find where is after in this new file 299 ja = 1 300 DO WHILE ( ksecsbc > sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast ) 301 ja = ja + 1 302 END DO 303 IF( ksecsbc > sdjf%nrecsec(ja) ) ja = ja + 1 ! in case ksecsbc > sdjf%nrecsec(sdjf%nreclast) 304 305 IF( ja > sdjf%nreclast ) THEN 306 CALL ctl_stop( "STOP", "fld_def: need next-next file? we should not be there... file: "//TRIM(sdjf%clrootname) ) 498 307 ENDIF 499 ELSE ! no time interpolation 500 sdjf%nrec_a(1) = 1 501 sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000 ! swap at the end of the year 502 sdjf%nrec_b(2) = nsec1jan000 ! beginning of the year (only for print) 503 ENDIF 504 ! 505 ! ! ============ ! 506 ELSEIF( NINT(sdjf%freqh) == -1 ) THEN ! monthly mean ! 507 ! ! ============ ! 508 ! 509 IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record 510 ! 511 ! INT( ztmp ) 512 ! /|\ 513 ! 1 | *---- 514 ! 0 |----( 515 ! |----+----|--> time 516 ! 0 /|\ 1 (nday/nmonth_len(nmonth)) 517 ! | 518 ! | 519 ! forcing record : nmonth 520 ! 521 ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 & 522 & + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) 523 imth = nmonth + INT( ztmp ) - COUNT((/llbefore/)) 524 IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) 525 ELSE ; sdjf%nrec_a(1) = imth 308 309 ! if ln_tint and if after is not the first record, we must (potentially again) update after data before the swap 310 IF( sdjf%ln_tint .AND. ja > 1 ) THEN 311 IF( sdjf%nrecsec(0) /= nflag ) THEN ! no trick used: after file is not the current file 312 sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec_a with before information 313 CALL fld_get( sdjf, Kmm ) ! read after data that will be used as before data 314 ENDIF 526 315 ENDIF 527 sdjf%nrec_a(2) = nmonth_half( imth ) + nsec1jan000 ! swap at the middle of the month 528 ELSE ! no time interpolation 529 IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1 530 ELSE ; sdjf%nrec_a(1) = nmonth 531 ENDIF 532 sdjf%nrec_a(2) = nmonth_end(nmonth ) + nsec1jan000 ! swap at the end of the month 533 sdjf%nrec_b(2) = nmonth_end(nmonth-1) + nsec1jan000 ! beginning of the month (only for print) 534 ENDIF 535 ! 536 ! ! ================================ ! 537 ELSE ! higher frequency mean (in hours) 538 ! ! ================================ ! 539 ! 540 ifreq_sec = NINT( sdjf%freqh * 3600. ) ! frequency mean (in seconds) 541 IF( sdjf%cltype(1:4) == 'week' ) isec_week = ksec_week( sdjf%cltype(6:8) ) ! since the first day of the current week 542 ! number of second since the beginning of the file 543 IF( sdjf%cltype == 'monthly' ) THEN ; ztmp = REAL(nsec_month,wp) ! since the first day of the current month 544 ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; ztmp = REAL(isec_week ,wp) ! since the first day of the current week 545 ELSEIF( sdjf%cltype == 'daily' ) THEN ; ztmp = REAL(nsec_day ,wp) ! since 00h of the current day 546 ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year 547 ENDIF 548 ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt + REAL( it_offset, wp ) ! centrered in the middle of sbc time step 549 ztmp = ztmp + 0.01 * rdt ! avoid truncation error 550 IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record 551 ! 552 ! INT( ztmp/ifreq_sec + 0.5 ) 553 ! /|\ 554 ! 2 | *-----( 555 ! 1 | *-----( 556 ! 0 |--( 557 ! |--+--|--+--|--+--|--> time 558 ! 0 /|\ 1 /|\ 2 /|\ 3 (ztmp/ifreq_sec) 559 ! | | | 560 ! | | | 561 ! forcing record : 1 2 3 562 ! 563 ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5 564 ELSE ! no time interpolation 565 ! 566 ! INT( ztmp/ifreq_sec ) 567 ! /|\ 568 ! 2 | *-----( 569 ! 1 | *-----( 570 ! 0 |-----( 571 ! |--+--|--+--|--+--|--> time 572 ! 0 /|\ 1 /|\ 2 /|\ 3 (ztmp/ifreq_sec) 573 ! | | | 574 ! | | | 575 ! forcing record : 1 2 3 576 ! 577 ztmp= ztmp / REAL(ifreq_sec, wp) 578 ENDIF 579 sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) ! record number to be read 580 581 iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000 ! end of this record (in second) 582 ! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1) 583 IF( sdjf%cltype == 'monthly' ) iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1)) 584 IF( sdjf%cltype(1:4) == 'week' ) iendrec = iendrec + ( nsec_year - isec_week ) 585 IF( sdjf%cltype == 'daily' ) iendrec = iendrec + NINT(rday) * ( nday_year - 1 ) 586 IF( sdjf%ln_tint ) THEN 587 sdjf%nrec_a(2) = iendrec - ifreq_sec / 2 ! swap at the middle of the record 316 317 ENDIF 318 319 IF( sdjf%ln_tint ) THEN 320 ! Swap data 321 sdjf%nrec_b(:) = sdjf%nrec_a(:) ! swap before record informations 322 sdjf%rotn(1) = sdjf%rotn(2) ! swap before rotate informations 323 sdjf%fdta(:,:,:,1) = sdjf%fdta(:,:,:,2) ! swap before record field 588 324 ELSE 589 sdjf%nrec_a(2) = iendrec ! swap at the end of the record 590 sdjf%nrec_b(2) = iendrec - ifreq_sec ! beginning of the record (only for print) 591 ENDIF 592 ! 593 ENDIF 594 ! 595 IF( .NOT. sdjf%ln_tint ) sdjf%nrec_a(2) = sdjf%nrec_a(2) - 1 ! last second belongs to bext record : *----( 596 ! 597 END SUBROUTINE fld_rec 598 599 600 SUBROUTINE fld_get( sdjf ) 325 sdjf%nrec_b(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! only for print 326 ENDIF 327 328 ! read new after data 329 sdjf%nrec_a(:) = (/ ja, sdjf%nrecsec(ja) /) ! update nrec_a as it is used by fld_get 330 CALL fld_get( sdjf, Kmm ) ! read after data (with nrec_a informations) 331 332 ENDIF 333 ! 334 END SUBROUTINE fld_update 335 336 337 SUBROUTINE fld_get( sdjf, Kmm ) 601 338 !!--------------------------------------------------------------------- 602 339 !! *** ROUTINE fld_get *** … … 604 341 !! ** Purpose : read the data 605 342 !!---------------------------------------------------------------------- 606 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 343 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 344 INTEGER , OPTIONAL, INTENT(in ) :: Kmm ! ocean time level index 607 345 ! 608 346 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 618 356 IF( ASSOCIATED(sdjf%imap) ) THEN 619 357 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), & 620 & sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint )358 & sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm ) 621 359 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), & 622 & sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint )360 & sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm ) 623 361 ENDIF 624 362 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN … … 656 394 ENDIF 657 395 CASE DEFAULT 658 IF 396 IF(lk_c1d .AND. lmoor ) THEN 659 397 IF( sdjf%ln_tint ) THEN 660 398 CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) ) … … 677 415 678 416 679 SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel, ldzint )417 SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel, ldzint, Kmm ) 680 418 !!--------------------------------------------------------------------- 681 419 !! *** ROUTINE fld_map *** … … 694 432 LOGICAL, OPTIONAL , INTENT(in ) :: ldtotvel ! true if total ( = barotrop + barocline) velocity 695 433 LOGICAL, OPTIONAL , INTENT(in ) :: ldzint ! true if 3D variable requires a vertical interpolation 434 INTEGER, OPTIONAL , INTENT(in ) :: Kmm ! ocean time level index 696 435 !! 697 436 INTEGER :: ipi ! length of boundary data on local process … … 758 497 759 498 CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar ) 760 CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel )499 CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel, Kmm) 761 500 DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz ) 762 501 … … 822 561 END SUBROUTINE fld_map_core 823 562 824 825 SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel) 563 SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel, Kmm ) 826 564 !!--------------------------------------------------------------------- 827 565 !! *** ROUTINE fld_bdy_interp *** … … 840 578 INTEGER , INTENT(in ) :: kgrd ! grid type (t, u, v) 841 579 INTEGER , INTENT(in ) :: kbdy ! bdy number 580 INTEGER, OPTIONAL , INTENT(in ) :: Kmm ! ocean time level index 842 581 !! 843 582 INTEGER :: ipi ! length of boundary data on local process … … 866 605 SELECT CASE( kgrd ) 867 606 CASE(1) ! depth of T points: 868 zdepth(:) = gdept _n(ji,jj,:)607 zdepth(:) = gdept(ji,jj,:,Kmm) 869 608 CASE(2) ! depth of U points: we must not use gdept_n as we don't want to do a communication 870 609 ! --> copy what is done for gdept_n in domvvl... 871 610 zdhalf(1) = 0.0_wp 872 zdepth(1) = 0.5_wp * e3uw _n(ji,jj,1)611 zdepth(1) = 0.5_wp * e3uw(ji,jj,1,Kmm) 873 612 DO jk = 2, jpk ! vertical sum 874 613 ! zcoef = umask - wumask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt … … 877 616 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 878 617 zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 879 zdhalf(jk) = zdhalf(jk-1) + e3u _n(ji,jj,jk-1)880 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3uw_n(ji,jj,jk)) &881 & + (1. -zcoef) * ( zdepth(jk-1) + e3uw_n(ji,jj,jk))618 zdhalf(jk) = zdhalf(jk-1) + e3u(ji,jj,jk-1,Kmm) 619 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5_wp * e3uw(ji,jj,jk,Kmm)) & 620 & + (1._wp-zcoef) * ( zdepth(jk-1) + e3uw(ji,jj,jk,Kmm)) 882 621 END DO 883 622 CASE(3) ! depth of V points: we must not use gdept_n as we don't want to do a communication 884 623 ! --> copy what is done for gdept_n in domvvl... 885 624 zdhalf(1) = 0.0_wp 886 zdepth(1) = 0.5_wp * e3vw _n(ji,jj,1)625 zdepth(1) = 0.5_wp * e3vw(ji,jj,1,Kmm) 887 626 DO jk = 2, jpk ! vertical sum 888 627 ! zcoef = vmask - wvmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt … … 891 630 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 892 631 zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 893 zdhalf(jk) = zdhalf(jk-1) + e3v _n(ji,jj,jk-1)894 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3vw_n(ji,jj,jk)) &895 & + (1. -zcoef) * ( zdepth(jk-1) + e3vw_n(ji,jj,jk))632 zdhalf(jk) = zdhalf(jk-1) + e3v(ji,jj,jk-1,Kmm) 633 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5_wp * e3vw(ji,jj,jk,Kmm)) & 634 & + (1._wp-zcoef) * ( zdepth(jk-1) + e3vw(ji,jj,jk,Kmm)) 896 635 END DO 897 636 END SELECT … … 911 650 END DO 912 651 ENDIF 913 END DO 652 END DO ! jpk 914 653 ! 915 654 END DO ! ipi … … 937 676 ztrans_new = 0._wp 938 677 DO jk = 1, jpk ! calculate transport on model grid 939 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)678 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u(ji,jj,jk,Kmm ) * umask(ji,jj,jk) 940 679 ENDDO 941 680 DO jk = 1, jpk ! make transport correction 942 681 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 943 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu _n(ji,jj) ) * umask(ji,jj,jk)682 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu(ji,jj,Kmm) ) * umask(ji,jj,jk) 944 683 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 945 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu_n(ji,jj) )* umask(ji,jj,jk)684 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu(ji,jj,Kmm) * umask(ji,jj,jk) 946 685 ENDIF 947 686 ENDDO … … 958 697 ztrans_new = 0._wp 959 698 DO jk = 1, jpk ! calculate transport on model grid 960 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)699 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v(ji,jj,jk,Kmm ) * vmask(ji,jj,jk) 961 700 ENDDO 962 701 DO jk = 1, jpk ! make transport correction 963 702 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 964 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv _n(ji,jj) ) * vmask(ji,jj,jk)703 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv(ji,jj,Kmm) ) * vmask(ji,jj,jk) 965 704 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 966 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv_n(ji,jj) )* vmask(ji,jj,jk)705 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv(ji,jj,Kmm) * vmask(ji,jj,jk) 967 706 ENDIF 968 707 ENDDO 969 708 ENDDO 970 709 END SELECT 971 710 972 711 END SUBROUTINE fld_bdy_interp 973 712 974 713 975 714 SUBROUTINE fld_rot( kt, sd ) 976 715 !!--------------------------------------------------------------------- … … 1013 752 sd(ju)%fdta(:,:,jk,jn) = utmp(:,:) ; sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:) 1014 753 ELSE 1015 CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk), 'T', 'en->i', utmp(:,:) )1016 CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk), 'T', 'en->j', vtmp(:,:) )754 CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->i', utmp(:,:) ) 755 CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->j', vtmp(:,:) ) 1017 756 sd(ju)%fnow(:,:,jk ) = utmp(:,:) ; sd(iv)%fnow(:,:,jk ) = vtmp(:,:) 1018 757 ENDIF … … 1030 769 1031 770 1032 SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop ) 771 SUBROUTINE fld_def( sdjf, ldprev, ldnext ) 772 !!--------------------------------------------------------------------- 773 !! *** ROUTINE fld_def *** 774 !! 775 !! ** Purpose : define the record(s) of the file and its name 776 !!---------------------------------------------------------------------- 777 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 778 LOGICAL, OPTIONAL, INTENT(in ) :: ldprev ! 779 LOGICAL, OPTIONAL, INTENT(in ) :: ldnext ! 780 ! 781 INTEGER :: jt 782 INTEGER :: idaysec ! number of seconds in 1 day = NINT(rday) 783 INTEGER :: iyr, imt, idy, isecwk 784 INTEGER :: indexyr, indexmt 785 INTEGER :: ireclast 786 INTEGER :: ishift, istart 787 INTEGER, DIMENSION(2) :: isave 788 REAL(wp) :: zfreqs 789 LOGICAL :: llprev, llnext, llstop 790 LOGICAL :: llprevmt, llprevyr 791 LOGICAL :: llnextmt, llnextyr 792 !!---------------------------------------------------------------------- 793 idaysec = NINT(rday) 794 ! 795 IF( PRESENT(ldprev) ) THEN ; llprev = ldprev 796 ELSE ; llprev = .FALSE. 797 ENDIF 798 IF( PRESENT(ldnext) ) THEN ; llnext = ldnext 799 ELSE ; llnext = .FALSE. 800 ENDIF 801 802 ! current file parameters 803 IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the current week 804 isecwk = ksec_week( sdjf%cltype(6:8) ) ! seconds between the beginning of the week and half of current time step 805 llprevmt = isecwk > nsec_month ! longer time since beginning of the current week than the current month 806 llprevyr = llprevmt .AND. nmonth == 1 807 iyr = nyear - COUNT((/llprevyr/)) 808 imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/)) 809 idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec 810 isecwk = nsec_year - isecwk ! seconds between 00h jan 1st of current year and current week beginning 811 ELSE 812 iyr = nyear 813 imt = nmonth 814 idy = nday 815 isecwk = 0 816 ENDIF 817 818 ! previous file parameters 819 IF( llprev ) THEN 820 IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of previous week 821 isecwk = isecwk + 7 * idaysec ! seconds between the beginning of previous week and half of the time step 822 llprevmt = isecwk > nsec_month ! longer time since beginning of the previous week than the current month 823 llprevyr = llprevmt .AND. nmonth == 1 824 iyr = nyear - COUNT((/llprevyr/)) 825 imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/)) 826 idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec 827 isecwk = nsec_year - isecwk ! seconds between 00h jan 1st of current year and previous week beginning 828 ELSE 829 idy = nday - COUNT((/ sdjf%cltype == 'daily' /)) 830 imt = nmonth - COUNT((/ sdjf%cltype == 'monthly' .OR. idy == 0 /)) 831 iyr = nyear - COUNT((/ sdjf%cltype == 'yearly' .OR. imt == 0 /)) 832 IF( idy == 0 ) idy = nmonth_len(imt) 833 IF( imt == 0 ) imt = 12 834 isecwk = 0 835 ENDIF 836 ENDIF 837 838 ! next file parameters 839 IF( llnext ) THEN 840 IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of next week 841 isecwk = 7 * idaysec - isecwk ! seconds between half of the time step and the beginning of next week 842 llnextmt = isecwk > ( nmonth_len(nmonth)*idaysec - nsec_month ) ! larger than the seconds to the end of the month 843 llnextyr = llnextmt .AND. nmonth == 12 844 iyr = nyear + COUNT((/llnextyr/)) 845 imt = nmonth + COUNT((/llnextmt/)) - 12 * COUNT((/llnextyr/)) 846 idy = nday - nmonth_len(nmonth) * COUNT((/llnextmt/)) + isecwk / idaysec + 1 847 isecwk = nsec_year + isecwk ! seconds between 00h jan 1st of current year and next week beginning 848 ELSE 849 idy = nday + COUNT((/ sdjf%cltype == 'daily' /)) 850 imt = nmonth + COUNT((/ sdjf%cltype == 'monthly' .OR. idy > nmonth_len(nmonth) /)) 851 iyr = nyear + COUNT((/ sdjf%cltype == 'yearly' .OR. imt == 13 /)) 852 IF( idy > nmonth_len(nmonth) ) idy = 1 853 IF( imt == 13 ) imt = 1 854 isecwk = 0 855 ENDIF 856 ENDIF 857 ! 858 ! find the last record to be read -> update sdjf%nreclast 859 indexyr = iyr - nyear + 1 ! which year are we looking for? previous(0), current(1) or next(2)? 860 indexmt = imt + 12 * ( indexyr - 1 ) ! which month are we looking for (relatively to current year)? 861 ! 862 ! Last record to be read in the current file 863 ! Predefine the number of record in the file according of its type. 864 ! We could compare this number with the number of records in the file and make a stop if the 2 numbers do not match... 865 ! However this would be much less fexible (e.g. for tests) and will force to rewite input files according to nleapy... 866 IF ( NINT(sdjf%freqh) == -12 ) THEN ; ireclast = 1 ! yearly mean: consider only 1 record 867 ELSEIF( NINT(sdjf%freqh) == -1 ) THEN ! monthly mean: 868 IF( sdjf%cltype == 'monthly' ) THEN ; ireclast = 1 ! consider that the file has 1 record 869 ELSE ; ireclast = 12 ! consider that the file has 12 record 870 ENDIF 871 ELSE ! higher frequency mean (in hours) 872 IF( sdjf%cltype == 'monthly' ) THEN ; ireclast = NINT( 24. * REAL(nmonth_len(indexmt), wp) / sdjf%freqh ) 873 ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; ireclast = NINT( 24. * 7. / sdjf%freqh ) 874 ELSEIF( sdjf%cltype == 'daily' ) THEN ; ireclast = NINT( 24. / sdjf%freqh ) 875 ELSE ; ireclast = NINT( 24. * REAL( nyear_len(indexyr), wp) / sdjf%freqh ) 876 ENDIF 877 ENDIF 878 879 sdjf%nreclast = ireclast 880 ! Allocate arrays for beginning/middle/end of each record (seconds since Jan. 1st 00h of nit000 year) 881 IF( ALLOCATED(sdjf%nrecsec) ) DEALLOCATE( sdjf%nrecsec ) 882 ALLOCATE( sdjf%nrecsec( 0:ireclast ) ) 883 ! 884 IF ( NINT(sdjf%freqh) == -12 ) THEN ! yearly mean and yearly file 885 SELECT CASE( indexyr ) 886 CASE(0) ; sdjf%nrecsec(0) = nsec1jan000 - nyear_len( 0 ) * idaysec 887 CASE(1) ; sdjf%nrecsec(0) = nsec1jan000 888 CASE(2) ; sdjf%nrecsec(0) = nsec1jan000 + nyear_len( 1 ) * idaysec 889 ENDSELECT 890 sdjf%nrecsec(1) = sdjf%nrecsec(0) + nyear_len( indexyr ) * idaysec 891 ELSEIF( NINT(sdjf%freqh) == -1 ) THEN ! monthly mean: 892 IF( sdjf%cltype == 'monthly' ) THEN ! monthly file 893 sdjf%nrecsec(0 ) = nsec1jan000 + nmonth_beg(indexmt ) 894 sdjf%nrecsec(1 ) = nsec1jan000 + nmonth_beg(indexmt+1) 895 ELSE ! yearly file 896 ishift = 12 * ( indexyr - 1 ) 897 sdjf%nrecsec(0:12) = nsec1jan000 + nmonth_beg(1+ishift:13+ishift) 898 ENDIF 899 ELSE ! higher frequency mean (in hours) 900 IF( sdjf%cltype == 'monthly' ) THEN ; istart = nsec1jan000 + nmonth_beg(indexmt) 901 ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; istart = nsec1jan000 + isecwk 902 ELSEIF( sdjf%cltype == 'daily' ) THEN ; istart = nsec1jan000 + nmonth_beg(indexmt) + ( idy - 1 ) * idaysec 903 ELSEIF( indexyr == 0 ) THEN ; istart = nsec1jan000 - nyear_len( 0 ) * idaysec 904 ELSEIF( indexyr == 2 ) THEN ; istart = nsec1jan000 + nyear_len( 1 ) * idaysec 905 ELSE ; istart = nsec1jan000 906 ENDIF 907 zfreqs = sdjf%freqh * rhhmm * rmmss 908 DO jt = 0, sdjf%nreclast 909 sdjf%nrecsec(jt) = istart + NINT( zfreqs * REAL(jt,wp) ) 910 END DO 911 ENDIF 912 ! 913 IF( sdjf%ln_tint ) THEN ! record time defined in the middle of the record, computed using an implementation 914 ! of the rounded average that is valid over the full integer range 915 sdjf%nrecsec(1:sdjf%nreclast) = sdjf%nrecsec(0:sdjf%nreclast-1) / 2 + sdjf%nrecsec(1:sdjf%nreclast) / 2 + & 916 & MAX( MOD( sdjf%nrecsec(0:sdjf%nreclast-1), 2 ), MOD( sdjf%nrecsec(1:sdjf%nreclast), 2 ) ) 917 END IF 918 ! 919 sdjf%clname = fld_filename( sdjf, idy, imt, iyr ) 920 ! 921 END SUBROUTINE fld_def 922 923 924 SUBROUTINE fld_clopn( sdjf ) 1033 925 !!--------------------------------------------------------------------- 1034 926 !! *** ROUTINE fld_clopn *** 1035 927 !! 1036 !! ** Purpose : update the file name and close/open the files 1037 !!---------------------------------------------------------------------- 1038 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 1039 INTEGER, OPTIONAL, INTENT(in ) :: kyear ! year value 1040 INTEGER, OPTIONAL, INTENT(in ) :: kmonth ! month value 1041 INTEGER, OPTIONAL, INTENT(in ) :: kday ! day value 1042 LOGICAL, OPTIONAL, INTENT(in ) :: ldstop ! stop if open to read a non-existing file (default = .TRUE.) 1043 ! 1044 LOGICAL :: llprevyr ! are we reading previous year file? 1045 LOGICAL :: llprevmth ! are we reading previous month file? 1046 INTEGER :: iyear, imonth, iday ! first day of the current file in yyyy mm dd 1047 INTEGER :: isec_week ! number of seconds since start of the weekly file 1048 INTEGER :: indexyr ! year undex (O/1/2: previous/current/next) 1049 REAL(wp) :: zyear_len, zmonth_len ! length (days) of iyear and imonth ! 1050 CHARACTER(len = 256) :: clname ! temporary file name 1051 !!---------------------------------------------------------------------- 1052 IF( PRESENT(kyear) ) THEN ! use given values 1053 iyear = kyear 1054 imonth = kmonth 1055 iday = kday 1056 IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week 1057 isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 1058 llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month 1059 llprevyr = llprevmth .AND. nmonth == 1 1060 iyear = nyear - COUNT((/llprevyr /)) 1061 imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)) 1062 iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday) 1063 ENDIF 1064 ELSE ! use current day values 1065 IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week 1066 isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week 1067 llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month 1068 llprevyr = llprevmth .AND. nmonth == 1 1069 ELSE 1070 isec_week = 0 1071 llprevmth = .FALSE. 1072 llprevyr = .FALSE. 1073 ENDIF 1074 iyear = nyear - COUNT((/llprevyr /)) 1075 imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)) 1076 iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday) 1077 ENDIF 1078 1079 ! build the new filename if not climatological data 1080 clname=TRIM(sdjf%clrootname) 1081 ! 1082 ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name 1083 IF( .NOT. sdjf%ln_clim ) THEN 1084 WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear ! add year 1085 IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), imonth ! add month 1086 ELSE 1087 ! build the new filename if climatological data 1088 IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth ! add month 1089 ENDIF 1090 IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) & 1091 & WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), iday ! add day 1092 ! 1093 IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN ! new file to be open 1094 ! 1095 sdjf%clname = TRIM(clname) 1096 IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open 1097 CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 ) 1098 ! 1099 ! find the last record to be read -> update sdjf%nreclast 1100 indexyr = iyear - nyear + 1 1101 zyear_len = REAL(nyear_len( indexyr ), wp) 1102 SELECT CASE ( indexyr ) 1103 CASE ( 0 ) ; zmonth_len = 31. ! previous year -> imonth = 12 1104 CASE ( 1 ) ; zmonth_len = REAL(nmonth_len(imonth), wp) 1105 CASE ( 2 ) ; zmonth_len = 31. ! next year -> imonth = 1 1106 END SELECT 1107 ! 1108 ! last record to be read in the current file 1109 IF ( sdjf%freqh == -12. ) THEN ; sdjf%nreclast = 1 ! yearly mean 1110 ELSEIF( sdjf%freqh == -1. ) THEN ! monthly mean 1111 IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = 1 1112 ELSE ; sdjf%nreclast = 12 1113 ENDIF 1114 ELSE ! higher frequency mean (in hours) 1115 IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = NINT( 24. * zmonth_len / sdjf%freqh ) 1116 ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; sdjf%nreclast = NINT( 24. * 7. / sdjf%freqh ) 1117 ELSEIF( sdjf%cltype == 'daily' ) THEN ; sdjf%nreclast = NINT( 24. / sdjf%freqh ) 1118 ELSE ; sdjf%nreclast = NINT( 24. * zyear_len / sdjf%freqh ) 1119 ENDIF 1120 ENDIF 928 !! ** Purpose : close/open the files 929 !!---------------------------------------------------------------------- 930 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 931 ! 932 INTEGER, DIMENSION(2) :: isave 933 LOGICAL :: llprev, llnext, llstop 934 !!---------------------------------------------------------------------- 935 ! 936 llprev = sdjf%nrecsec(sdjf%nreclast) < nsec000_1jan000 ! file ends before the beginning of the job -> file may not exist 937 llnext = sdjf%nrecsec( 0 ) > nsecend_1jan000 ! file begins after the end of the job -> file may not exist 938 939 llstop = sdjf%ln_clim .OR. .NOT. ( llprev .OR. llnext ) 940 941 IF( sdjf%num <= 0 .OR. .NOT. sdjf%ln_clim ) THEN 942 IF( sdjf%num > 0 ) CALL iom_close( sdjf%num ) ! close file if already open 943 CALL iom_open( sdjf%clname, sdjf%num, ldstop = llstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 ) 944 ENDIF 945 ! 946 IF( sdjf%num <= 0 .AND. .NOT. llstop ) THEN ! file not found but we do accept this... 947 ! 948 IF( llprev ) THEN ! previous file does not exist : go back to current and accept to read only the first record 949 CALL ctl_warn('previous file: '//TRIM(sdjf%clname)//' not found -> go back to current year/month/week/day file') 950 isave(1:2) = sdjf%nrecsec(sdjf%nreclast-1:sdjf%nreclast) ! save previous file info 951 CALL fld_def( sdjf ) ! go back to current file 952 sdjf%nreclast = 1 ! force to use only the first record (do as if other were not existing...) 953 sdjf%nrecsec(0:1) = isave(1:2) 954 ENDIF 955 ! 956 IF( llnext ) THEN ! next file does not exist : go back to current and accept to read only the last record 957 CALL ctl_warn('next file: '//TRIM(sdjf%clname)//' not found -> go back to current year/month/week/day file') 958 isave(1:2) = sdjf%nrecsec(0:1) ! save next file info 959 CALL fld_def( sdjf ) ! go back to current file 960 ! -> read last record but keep record info from the first record of next file 961 sdjf%nrecsec(sdjf%nreclast-1:sdjf%nreclast) = isave(1:2) 962 sdjf%nrecsec(0:sdjf%nreclast-2) = nflag 963 ENDIF 964 ! 965 CALL iom_open( sdjf%clname, sdjf%num, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 ) 1121 966 ! 1122 967 ENDIF … … 1300 1145 CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 ) 1301 1146 1302 !! get dimensions 1303 IF ( SIZE(sd%fnow, 3) > 1) THEN1147 !! get dimensions: we consider 2D data as 3D data with vertical dim size = 1 1148 IF( SIZE(sd%fnow, 3) > 0 ) THEN 1304 1149 ALLOCATE( ddims(4) ) 1305 1150 ELSE … … 1314 1159 1315 1160 CALL iom_open ( sd%wgtname, inum ) ! interpolation weights 1316 IF 1161 IF( inum > 0 ) THEN 1317 1162 1318 1163 !! determine whether we have an east-west cyclic grid … … 1623 1468 1624 1469 ref_wgts(kw)%fly_dta(:,:,:) = 0.0 1625 SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) ) 1626 CASE(1) 1627 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn) 1628 CASE DEFAULT 1629 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 1630 END SELECT 1470 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 1631 1471 ENDIF 1632 1472 … … 1646 1486 END DO 1647 1487 1648 IF 1488 IF(ref_wgts(kw)%numwgt .EQ. 16) THEN 1649 1489 1650 1490 !! fix up halo points that we couldnt read from file … … 1672 1512 IF( jpi1 == 2 ) THEN 1673 1513 rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap 1674 SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1675 CASE(1) 1676 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1677 CASE DEFAULT 1678 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1679 END SELECT 1514 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1680 1515 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1681 1516 ENDIF 1682 1517 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 1683 1518 rec1(1) = 1 + ref_wgts(kw)%overlap 1684 SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1685 CASE(1) 1686 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1687 CASE DEFAULT 1688 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1689 END SELECT 1519 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1690 1520 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1691 1521 ENDIF … … 1729 1559 END DO 1730 1560 ! 1731 END 1561 ENDIF 1732 1562 ! 1733 1563 END SUBROUTINE fld_interp 1734 1564 1735 1565 1566 FUNCTION fld_filename( sdjf, kday, kmonth, kyear ) 1567 !!--------------------------------------------------------------------- 1568 !! *** FUNCTION fld_filename *** 1569 !! 1570 !! ** Purpose : define the filename according to a given date 1571 !!--------------------------------------------------------------------- 1572 TYPE(FLD), INTENT(in) :: sdjf ! input field related variables 1573 INTEGER , INTENT(in) :: kday, kmonth, kyear 1574 ! 1575 CHARACTER(len = 256) :: clname, fld_filename 1576 !!--------------------------------------------------------------------- 1577 1578 1579 ! build the new filename if not climatological data 1580 clname=TRIM(sdjf%clrootname) 1581 ! 1582 ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name 1583 IF( .NOT. sdjf%ln_clim ) THEN 1584 WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year 1585 IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a, "m",i2.2)' ) TRIM( clname ), kmonth ! add month 1586 ELSE 1587 ! build the new filename if climatological data 1588 IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth ! add month 1589 ENDIF 1590 IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) & 1591 & WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), kday ! add day 1592 1593 fld_filename = clname 1594 1595 END FUNCTION fld_filename 1596 1597 1736 1598 FUNCTION ksec_week( cdday ) 1737 1599 !!--------------------------------------------------------------------- 1738 !! *** FUNCTION ks hift_week ***1739 !! 1740 !! ** Purpose : return the first 3 letters of the first day of the weekly file1600 !! *** FUNCTION ksec_week *** 1601 !! 1602 !! ** Purpose : seconds between 00h of the beginning of the week and half of the current time step 1741 1603 !!--------------------------------------------------------------------- 1742 1604 CHARACTER(len=*), INTENT(in) :: cdday ! first 3 letters of the first day of the weekly file … … 1754 1616 ishift = ijul * NINT(rday) 1755 1617 ! 1756 ksec_week = nsec_ week+ ishift1618 ksec_week = nsec_monday + ishift 1757 1619 ksec_week = MOD( ksec_week, 7*NINT(rday) ) 1758 1620 !
Note: See TracChangeset
for help on using the changeset viewer.