New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
fldread.F90 in NEMO/trunk/src/OCE/SBC – NEMO

source: NEMO/trunk/src/OCE/SBC/fldread.F90 @ 10522

Last change on this file since 10522 was 10425, checked in by smasson, 5 years ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 108.6 KB
RevLine 
[888]1MODULE fldread
2   !!======================================================================
3   !!                       ***  MODULE  fldread  ***
4   !! Ocean forcing:  read input field for surface boundary condition
5   !!=====================================================================
[7646]6   !! History :  2.0  !  2006-06  (S. Masson, G. Madec)  Original code
7   !!            3.0  !  2008-05  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid
8   !!            3.4  !  2013-10  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation
9   !!                 !  12-2015  (J. Harle) Adding BDY on-the-fly interpolation
[888]10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
[7646]13   !!   fld_read      : read input fields used for the computation of the surface boundary condition
14   !!   fld_init      : initialization of field read
15   !!   fld_rec       : determined the record(s) to be read
16   !!   fld_get       : read the data
17   !!   fld_map       : read global data from file and map onto local data using a general mapping (use for open boundaries)
18   !!   fld_rot       : rotate the vector fields onto the local grid direction
19   !!   fld_clopn     : update the data file name and close/open the files
20   !!   fld_fill      : fill the data structure with the associated information read in namelist
21   !!   wgt_list      : manage the weights used for interpolation
22   !!   wgt_print     : print the list of known weights
23   !!   fld_weight    : create a WGT structure and fill in data from file, restructuring as required
24   !!   apply_seaoverland : fill land with ocean values
25   !!   seaoverland   : create shifted matrices for seaoverland application
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
[888]28   !!----------------------------------------------------------------------
[6140]29   USE oce            ! ocean dynamics and tracers
30   USE dom_oce        ! ocean space and time domain
31   USE phycst         ! physical constant
32   USE sbc_oce        ! surface boundary conditions : fields
33   USE geo2ocean      ! for vector rotation on to model grid
34   !
35   USE in_out_manager ! I/O manager
36   USE iom            ! I/O manager library
37   USE ioipsl  , ONLY : ymds2ju, ju2ymds   ! for calendar
38   USE lib_mpp        ! MPP library
39   USE lbclnk         ! ocean lateral boundary conditions (C1D case)
[4230]40   
[888]41   IMPLICIT NONE
42   PRIVATE   
[3294]43 
44   PUBLIC   fld_map    ! routine called by tides_init
[3851]45   PUBLIC   fld_read, fld_fill   ! called by sbc... modules
[5768]46   PUBLIC   fld_clopn
[888]47
48   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
[1730]49      CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file
[4245]50      REAL(wp)             ::   nfreqh      ! frequency of each flux file
[1730]51      CHARACTER(len = 34)  ::   clvar       ! generic name of the variable in the NetCDF flux file
52      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F)
53      LOGICAL              ::   ln_clim     ! climatology or not (T/F)
[2528]54      CHARACTER(len = 8)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly'
[4663]55      CHARACTER(len = 256) ::   wname       ! generic name of a NetCDF weights file to be used, blank if not
[1730]56      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation
[2715]57      !                                     ! a string starting with "U" or "V" for each component   
58      !                                     ! chars 2 onwards identify which components go together 
[4230]59      CHARACTER(len = 34)  ::   lname       ! generic name of a NetCDF land/sea mask file to be used, blank if not
60      !                                     ! 0=sea 1=land
[888]61   END TYPE FLD_N
62
63   TYPE, PUBLIC ::   FLD        !: Input field related variables
64      CHARACTER(len = 256)            ::   clrootname   ! generic name of the NetCDF file
65      CHARACTER(len = 256)            ::   clname       ! current name of the NetCDF file
[4245]66      REAL(wp)                        ::   nfreqh       ! frequency of each flux file
[888]67      CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file
68      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F)
[1132]69      LOGICAL                         ::   ln_clim      ! climatology or not (T/F)
[2528]70      CHARACTER(len = 8)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly'
[1132]71      INTEGER                         ::   num          ! iom id of the jpfld files to be read
[1730]72      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
73      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year)
[6140]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
[1275]76      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
[2715]77      !                                                 ! into the WGTLIST structure
[1275]78      CHARACTER(len = 34)             ::   vcomp        ! symbolic name for a vector component that needs rotation
[3851]79      LOGICAL, DIMENSION(2)           ::   rotn         ! flag to indicate whether before/after field has been rotated
80      INTEGER                         ::   nreclast     ! last record to be read in the current file
[4230]81      CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key
[7646]82      INTEGER                         ::   igrd         ! grid type for bdy data
83      INTEGER                         ::   ibdy         ! bdy set id number
[888]84   END TYPE FLD
85
[5132]86   TYPE, PUBLIC ::   MAP_POINTER      !: Map from input data file to local domain
87      INTEGER, POINTER, DIMENSION(:)  ::  ptr           ! Array of integer pointers to 1D arrays
88      LOGICAL                         ::  ll_unstruc    ! Unstructured (T) or structured (F) boundary data file
[3294]89   END TYPE MAP_POINTER
90
[1275]91!$AGRIF_DO_NOT_TREAT
92
93   !! keep list of all weights variables so they're only read in once
94   !! need to add AGRIF directives not to process this structure
95   !! also need to force wgtname to include AGRIF nest number
96   TYPE         ::   WGT        !: Input weights related variables
97      CHARACTER(len = 256)                    ::   wgtname      ! current name of the NetCDF weight file
98      INTEGER , DIMENSION(2)                  ::   ddims        ! shape of input grid
99      INTEGER , DIMENSION(2)                  ::   botleft      ! top left corner of box in input grid containing
[2715]100      !                                                         ! current processor grid
[1275]101      INTEGER , DIMENSION(2)                  ::   topright     ! top right corner of box
102      INTEGER                                 ::   jpiwgt       ! width of box on input grid
103      INTEGER                                 ::   jpjwgt       ! height of box on input grid
104      INTEGER                                 ::   numwgt       ! number of weights (4=bilinear, 16=bicubic)
105      INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in
[2528]106      INTEGER                                 ::   overlap      ! =0 when cyclic grid has no overlapping EW columns
[2715]107      !                                                         ! =>1 when they have one or more overlapping columns     
108      !                                                         ! =-1 not cyclic
[1275]109      LOGICAL                                 ::   cyclic       ! east-west cyclic or not
[2528]110      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpi     ! array of source integers
111      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpj     ! array of source integers
[1275]112      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid
[2528]113      REAL(wp), DIMENSION(:,:,:), POINTER     ::   fly_dta      ! array of values on input grid
114      REAL(wp), DIMENSION(:,:,:), POINTER     ::   col          ! temporary array for reading in columns
[1275]115   END TYPE WGT
116
[9019]117   INTEGER,     PARAMETER             ::   tot_wgts = 20
[1275]118   TYPE( WGT ), DIMENSION(tot_wgts)   ::   ref_wgts     ! array of wgts
119   INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array
[4230]120   REAL(wp), PARAMETER                ::   undeff_lsm = -999.00_wp
[1275]121
122!$AGRIF_END_DO_NOT_TREAT
123
[888]124   !!----------------------------------------------------------------------
[9598]125   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1156]126   !! $Id$
[10068]127   !! Software governed by the CeCILL license (see ./LICENSE)
[888]128   !!----------------------------------------------------------------------
129CONTAINS
130
[7646]131   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl )
[888]132      !!---------------------------------------------------------------------
133      !!                    ***  ROUTINE fld_read  ***
134      !!                   
135      !! ** Purpose :   provide at each time step the surface ocean fluxes
136      !!                (momentum, heat, freshwater and runoff)
137      !!
138      !! ** Method  :   READ each input fields in NetCDF files using IOM
139      !!      and intepolate it to the model time-step.
140      !!         Several assumptions are made on the input file:
141      !!      blahblahblah....
142      !!----------------------------------------------------------------------
143      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
[1132]144      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)
[888]145      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
[3851]146      TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping indices
147      INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option
148      INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now"
[6140]149      !                                                     !   kt_offset = -1 => fields at "before" time level
150      !                                                     !   kt_offset = +1 => fields at "after"  time level
151      !                                                     !   etc.
[7646]152      INTEGER  , INTENT(in   ), OPTIONAL     ::   jpk_bdy   ! number of vertical levels in the BDY data
153      LOGICAL  , INTENT(in   ), OPTIONAL     ::   fvl   ! number of vertical levels in the BDY data
154      !!
[6140]155      INTEGER  ::   itmp         ! local variable
156      INTEGER  ::   imf          ! size of the structure sd
157      INTEGER  ::   jf           ! dummy indices
158      INTEGER  ::   isecend      ! number of second since Jan. 1st 00h of nit000 year at nitend
159      INTEGER  ::   isecsbc      ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
160      INTEGER  ::   it_offset    ! local time offset variable
161      LOGICAL  ::   llnxtyr      ! open next year  file?
162      LOGICAL  ::   llnxtmth     ! open next month file?
163      LOGICAL  ::   llstop       ! stop is the file does not exist
[3294]164      LOGICAL  ::   ll_firstcall ! true if this is the first call to fld_read for this set of fields
[6140]165      REAL(wp) ::   ztinta       ! ratio applied to after  records when doing time interpolation
166      REAL(wp) ::   ztintb       ! ratio applied to before records when doing time interpolation
167      CHARACTER(LEN=1000) ::   clfmt  ! write format
168      TYPE(MAP_POINTER)   ::   imap   ! global-to-local mapping indices
[888]169      !!---------------------------------------------------------------------
[3851]170      ll_firstcall = kt == nit000
171      IF( PRESENT(kit) )   ll_firstcall = ll_firstcall .and. kit == 1
[3294]172
[5407]173      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
174      ELSE                                      ;   it_offset = 0
175      ENDIF
[3851]176      IF( PRESENT(kt_offset) )   it_offset = kt_offset
177
178      imap%ptr => NULL()
179
[2528]180      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
[3851]181      IF( present(kit) ) THEN   ! ignore kn_fsbc in this case
182         isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) )
183      ELSE                      ! middle of sbc time step
[6140]184         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdt) + it_offset * NINT(rdt)
[3294]185      ENDIF
[1275]186      imf = SIZE( sd )
[2323]187      !
[3294]188      IF( ll_firstcall ) THEN                      ! initialization
[3851]189         DO jf = 1, imf 
190            IF( PRESENT(map) ) imap = map(jf)
[7646]191               IF( PRESENT(jpk_bdy) ) THEN
192                  CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl )  ! read each before field (put them in after as they will be swapped)
193               ELSE
194                  CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped)
195               ENDIF
[3851]196         END DO
[2528]197         IF( lwp ) CALL wgt_print()                ! control print
198      ENDIF
199      !                                            ! ====================================== !
200      IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! update field at each kn_fsbc time-step !
201         !                                         ! ====================================== !
[888]202         !
[2528]203         DO jf = 1, imf                            ! ---   loop over field   --- !
204           
[3851]205            IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN    ! read/update the after data?
[888]206
[3851]207               IF( PRESENT(map) )   imap = map(jf)   ! temporary definition of map
[1132]208
[3851]209               sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)                                  ! swap before record informations
210               sd(jf)%rotn(1) = sd(jf)%rotn(2)                                      ! swap before rotate informations
211               IF( sd(jf)%ln_tint )   sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! swap before record field
212
213               CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit )    ! update after record informations
214
[6140]215               ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd),
[3851]216               ! it is possible that the before value is no more the good one... we have to re-read it
217               ! if before is not the last record of the file currently opened and after is the first record to be read
218               ! in a new file which means after = 1 (the file to be opened corresponds to the current time)
219               ! or after = nreclast + 1 (the file to be opened corresponds to a future time step)
220               IF( .NOT. ll_firstcall .AND. sd(jf)%ln_tint .AND. sd(jf)%nrec_b(1) /= sd(jf)%nreclast &
221                  &                   .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) == 1 ) THEN
222                  itmp = sd(jf)%nrec_a(1)                       ! temporary storage
223                  sd(jf)%nrec_a(1) = sd(jf)%nreclast            ! read the last record of the file currently opened
224                  CALL fld_get( sd(jf), imap )                  ! read after data
225                  sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
226                  sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
[4245]227                  sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case
[3851]228                  sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations
229                  sd(jf)%nrec_a(1) = itmp                       ! move back to after record
[3294]230               ENDIF
[1132]231
[3851]232               CALL fld_clopn( sd(jf) )   ! Do we need to open a new year/month/week/day file?
233               
[2528]234               IF( sd(jf)%ln_tint ) THEN
[3851]235                 
[6140]236                  ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd),
[3851]237                  ! it is possible that the before value is no more the good one... we have to re-read it
238                  ! if before record is not just just before the after record...
239                  IF( .NOT. ll_firstcall .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) /= 1 &
240                     &                   .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN   
241                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1       ! move back to before record
242                     CALL fld_get( sd(jf), imap )                  ! read after data
243                     sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
244                     sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
[4245]245                     sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case
[3851]246                     sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations
247                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) + 1       ! move back to after record
248                  ENDIF
[9807]249               ENDIF ! temporal interpolation?
[3851]250
[9807]251               ! do we have to change the year/month/week/day of the forcing field??
252               ! if we do time interpolation we will need to open next year/month/week/day file before the end of the current
253               ! 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)
254               ! will be larger than the record number that should be read for current year/month/week/day
255               ! do we need next file data?
256               ! This applies to both cases with or without time interpolation
257               IF( sd(jf)%nrec_a(1) > sd(jf)%nreclast ) THEN
258                 
259                  sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - sd(jf)%nreclast   !
260                 
261                  IF( .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) ) THEN   ! close/open the current/new file
[3851]262                     
[9807]263                     llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth)      ! open next month file?
264                     llnxtyr  = sd(jf)%cltype == 'yearly'  .OR. (nmonth == 12 .AND. llnxtmth)   ! open next year  file?
[1132]265
[9807]266                     ! if the run finishes at the end of the current year/month/week/day, we will allow next
267                     ! year/month/week/day file to be not present. If the run continue further than the current
268                     ! year/month/week/day, next year/month/week/day file must exist
269                     isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt)   ! second at the end of the run
270                     llstop = isecend > sd(jf)%nrec_a(2)                                   ! read more than 1 record of next year
271                     ! we suppose that the date of next file is next day (should be ok even for weekly files...)
272                     CALL fld_clopn( sd(jf), nyear  + COUNT((/llnxtyr /))                                           ,         &
273                        &                    nmonth + COUNT((/llnxtmth/)) - 12                 * COUNT((/llnxtyr /)),         &
274                        &                    nday   + 1                   - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
[2528]275
[9807]276                     IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN    ! next year file does not exist
277                        CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)//     &
278                           &     ' not present -> back to current year/month/day')
279                        CALL fld_clopn( sd(jf) )       ! back to the current year/month/day
280                        sd(jf)%nrec_a(1) = sd(jf)%nreclast     ! force to read the last record in the current year file
[1132]281                     ENDIF
[9807]282                     
283                  ENDIF
284               ENDIF   ! open need next file?
[3851]285                 
[2528]286               ! read after data
[7646]287               IF( PRESENT(jpk_bdy) ) THEN
288                  CALL fld_get( sd(jf), imap, jpk_bdy, fvl)
289               ELSE
290                  CALL fld_get( sd(jf), imap )
291               ENDIF
[3851]292            ENDIF   ! read new data?
[2528]293         END DO                                    ! --- end loop over field --- !
[1132]294
[3851]295         CALL fld_rot( kt, sd )                    ! rotate vector before/now/after fields if needed
[888]296
[2528]297         DO jf = 1, imf                            ! ---   loop over field   --- !
[888]298            !
[2528]299            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
[1191]300               IF(lwp .AND. kt - nit000 <= 100 ) THEN
[7646]301                  clfmt = "('   fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
[4245]302                     &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
[3294]303                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &           
[1730]304                     & 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
[7646]305                  WRITE(numout, *) '      it_offset is : ',it_offset
[1191]306               ENDIF
[2528]307               ! temporal interpolation weights
[2323]308               ztinta =  REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
[1132]309               ztintb =  1. - ztinta
[2528]310               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
311            ELSE   ! nothing to do...
[1191]312               IF(lwp .AND. kt - nit000 <= 100 ) THEN
[7646]313                  clfmt = "('   fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
[4245]314                     &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
[2528]315                  WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
316                     &                 sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
[1191]317               ENDIF
[888]318            ENDIF
319            !
[2528]320            IF( kt == nitend - kn_fsbc + 1 )   CALL iom_close( sd(jf)%num )   ! Close the input files
[1132]321
[2528]322         END DO                                    ! --- end loop over field --- !
323         !
[6140]324      ENDIF
[2528]325      !
[888]326   END SUBROUTINE fld_read
327
328
[7646]329   SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl)
[888]330      !!---------------------------------------------------------------------
[1132]331      !!                    ***  ROUTINE fld_init  ***
332      !!
[4371]333      !! ** Purpose :  - first call to fld_rec to define before values
334      !!               - if time interpolation, read before data
[1132]335      !!----------------------------------------------------------------------
[7646]336      INTEGER  , INTENT(in   ) ::   kn_fsbc      ! sbc computation period (in time step)
337      TYPE(FLD), INTENT(inout) ::   sdjf         ! input field related variables
338      TYPE(MAP_POINTER),INTENT(in) ::   map      ! global-to-local mapping indices
339      INTEGER  , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data
340      LOGICAL  , INTENT(in), OPTIONAL :: fvl     ! number of vertical levels in the BDY data
[1132]341      !!
[2528]342      LOGICAL :: llprevyr              ! are we reading previous year  file?
343      LOGICAL :: llprevmth             ! are we reading previous month file?
344      LOGICAL :: llprevweek            ! are we reading previous week  file?
345      LOGICAL :: llprevday             ! are we reading previous day   file?
346      LOGICAL :: llprev                ! llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
347      INTEGER :: idvar                 ! variable id
348      INTEGER :: inrec                 ! number of record existing for this variable
349      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
350      INTEGER :: isec_week             ! number of seconds since start of the weekly file
[1191]351      CHARACTER(LEN=1000) ::   clfmt   ! write format
[1132]352      !!---------------------------------------------------------------------
[2528]353      llprevyr   = .FALSE.
354      llprevmth  = .FALSE.
355      llprevweek = .FALSE.
356      llprevday  = .FALSE.
357      isec_week  = 0
[6140]358      !
[1132]359      ! define record informations
[2528]360      CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. )  ! return before values in sdjf%nrec_a (as we will swap it later)
[6140]361      !
[2528]362      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
[6140]363      !
[1132]364      IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
[6140]365         !
[2528]366         IF( sdjf%nrec_a(1) == 0  ) THEN   ! we redefine record sdjf%nrec_a(1) with the last record of previous year file
367            IF    ( sdjf%nfreqh == -12 ) THEN   ! yearly mean
368               IF( sdjf%cltype == 'yearly' ) THEN             ! yearly file
369                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
370                  llprevyr  = .NOT. sdjf%ln_clim                                           ! use previous year  file?
371               ELSE
[3851]372                  CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clrootname) )
[2528]373               ENDIF
374            ELSEIF( sdjf%nfreqh ==  -1 ) THEN   ! monthly mean
375               IF( sdjf%cltype == 'monthly' ) THEN            ! monthly file
376                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
377                  llprevmth = .TRUE.                                                       ! use previous month file?
[1628]378                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
[2528]379               ELSE                                           ! yearly file
380                  sdjf%nrec_a(1) = 12                                                      ! force to read december mean
[1628]381                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
382               ENDIF
[2528]383            ELSE                                ! higher frequency mean (in hours)
384               IF    ( sdjf%cltype      == 'monthly' ) THEN   ! monthly file
[4245]385                  sdjf%nrec_a(1) = NINT( 24 * nmonth_len(nmonth-1) / sdjf%nfreqh )         ! last record of previous month
[2528]386                  llprevmth = .TRUE.                                                       ! use previous month file?
[1628]387                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
[2528]388               ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ! weekly file
389                  llprevweek = .TRUE.                                                      ! use previous week  file?
[4245]390                  sdjf%nrec_a(1) = NINT( 24 * 7 / sdjf%nfreqh )                            ! last record of previous week
[2528]391                  isec_week = NINT(rday) * 7                                               ! add a shift toward previous week
392               ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ! daily file
[4245]393                  sdjf%nrec_a(1) = NINT( 24 / sdjf%nfreqh )                                ! last record of previous day
[2528]394                  llprevday = .TRUE.                                                       ! use previous day   file?
[1628]395                  llprevmth = llprevday .AND. nday   == 1                                  ! use previous month file?
396                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
[2528]397               ELSE                                           ! yearly file
[4245]398                  sdjf%nrec_a(1) = NINT( 24 * nyear_len(0) / sdjf%nfreqh )                 ! last record of previous year
[1628]399                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
[1132]400               ENDIF
401            ENDIF
402         ENDIF
[3851]403         !
[2528]404         IF ( sdjf%cltype(1:4) == 'week' ) THEN
405            isec_week = isec_week + ksec_week( sdjf%cltype(6:8) )   ! second since the beginning of the week
406            llprevmth = isec_week > nsec_month                      ! longer time since the beginning of the week than the month
407            llprevyr  = llprevmth .AND. nmonth == 1
408         ENDIF
409         llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
410         !
411         iyear  = nyear  - COUNT((/llprevyr /))
412         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
413         iday   = nday   - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
414         !
415         CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev )
[6140]416         !
[1628]417         ! if previous year/month/day file does not exist, we switch to the current year/month/day
[1818]418         IF( llprev .AND. sdjf%num <= 0 ) THEN
[3851]419            CALL ctl_warn( 'previous year/month/week/day file: '//TRIM(sdjf%clrootname)//   &
[2528]420               &           ' not present -> back to current year/month/week/day' )
[1628]421            ! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
[2528]422            llprev = .FALSE.
423            sdjf%nrec_a(1) = 1
[3851]424            CALL fld_clopn( sdjf )
[1132]425         ENDIF
[6140]426         !
[3851]427         IF( llprev ) THEN   ! check if the record sdjf%nrec_a(1) exists in the file
[1132]428            idvar = iom_varid( sdjf%num, sdjf%clvar )                                        ! id of the variable sdjf%clvar
429            IF( idvar <= 0 )   RETURN
430            inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar )   ! size of the last dim of idvar
[2528]431            sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec )   ! make sure we select an existing record
[1132]432         ENDIF
[6140]433         !
[7646]434         ! read before data in after arrays(as we will swap it later)
435         IF( PRESENT(jpk_bdy) ) THEN
436            CALL fld_get( sdjf, map, jpk_bdy, fvl )
437         ELSE
438            CALL fld_get( sdjf, map )
439         ENDIF
[6140]440         !
[7646]441         clfmt = "('   fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
[2528]442         IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
[6140]443         !
[1132]444      ENDIF
[2715]445      !
[1132]446   END SUBROUTINE fld_init
447
448
[3851]449   SUBROUTINE fld_rec( kn_fsbc, sdjf, ldbefore, kit, kt_offset )
[1132]450      !!---------------------------------------------------------------------
[888]451      !!                    ***  ROUTINE fld_rec  ***
452      !!
[2528]453      !! ** Purpose : Compute
454      !!              if sdjf%ln_tint = .TRUE.
455      !!                  nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping)
456      !!              if sdjf%ln_tint = .FALSE.
457      !!                  nrec_a(1): record number
458      !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record (for print only)
[888]459      !!----------------------------------------------------------------------
[2528]460      INTEGER  , INTENT(in   )           ::   kn_fsbc   ! sbc computation period (in time step)
461      TYPE(FLD), INTENT(inout)           ::   sdjf      ! input field related variables
462      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldbefore  ! sent back before record values (default = .FALSE.)
[3851]463      INTEGER  , INTENT(in   ), OPTIONAL ::   kit       ! index of barotropic subcycle
[6140]464      !                                                 ! used only if sdjf%ln_tint = .TRUE.
[3851]465      INTEGER  , INTENT(in   ), OPTIONAL ::   kt_offset ! Offset of required time level compared to "now"
[6140]466      !                                                 !   time level in units of time steps.
467      !
[2528]468      LOGICAL  ::   llbefore    ! local definition of ldbefore
469      INTEGER  ::   iendrec     ! end of this record (in seconds)
470      INTEGER  ::   imth        ! month number
471      INTEGER  ::   ifreq_sec   ! frequency mean (in seconds)
472      INTEGER  ::   isec_week   ! number of seconds since the start of the weekly file
[3851]473      INTEGER  ::   it_offset   ! local time offset variable
[1132]474      REAL(wp) ::   ztmp        ! temporary variable
[888]475      !!----------------------------------------------------------------------
476      !
[2528]477      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
[2323]478      !
[2528]479      IF( PRESENT(ldbefore) ) THEN   ;   llbefore = ldbefore .AND. sdjf%ln_tint   ! needed only if sdjf%ln_tint = .TRUE.
480      ELSE                           ;   llbefore = .FALSE.
481      ENDIF
482      !
[5407]483      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
484      ELSE                                      ;   it_offset = 0
485      ENDIF
[3851]486      IF( PRESENT(kt_offset) )   it_offset = kt_offset
487      IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) )
[6140]488      ELSE                      ;   it_offset =         it_offset   * NINT(       rdt            )
[3851]489      ENDIF
[3294]490      !
[2528]491      !                                      ! =========== !
492      IF    ( sdjf%nfreqh == -12 ) THEN      ! yearly mean
493         !                                   ! =========== !
[888]494         !
[1132]495         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
496            !
497            !                  INT( ztmp )
498            !                     /|\
499            !                    1 |    *----
500            !                    0 |----(             
501            !                      |----+----|--> time
[2528]502            !                      0   /|\   1   (nday/nyear_len(1))
503            !                           |   
504            !                           |   
505            !       forcing record :    1
506            !                           
[4784]507            ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &
508           &       + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )
[2528]509            sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
510            ! swap at the middle of the year
[4784]511            IF( llbefore ) THEN   ;   sdjf%nrec_a(2) = nsec1jan000 - (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(0) + &
512                                    & INT(ztmp) * NINT( 0.5 * rday) * nyear_len(1) 
513            ELSE                  ;   sdjf%nrec_a(2) = nsec1jan000 + (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(1) + &
514                                    & INT(ztmp) * INT(rday) * nyear_len(1) + INT(ztmp) * NINT( 0.5 * rday) * nyear_len(2) 
[2528]515            ENDIF
516         ELSE                                    ! no time interpolation
517            sdjf%nrec_a(1) = 1
518            sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000   ! swap at the end    of the year
519            sdjf%nrec_b(2) = nsec1jan000                               ! beginning of the year (only for print)
520         ENDIF
521         !
522         !                                   ! ============ !
523      ELSEIF( sdjf%nfreqh ==  -1 ) THEN      ! monthly mean !
524         !                                   ! ============ !
525         !
526         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
527            !
528            !                  INT( ztmp )
529            !                     /|\
530            !                    1 |    *----
531            !                    0 |----(             
532            !                      |----+----|--> time
[1132]533            !                      0   /|\   1   (nday/nmonth_len(nmonth))
534            !                           |   
535            !                           |   
536            !       forcing record :  nmonth
537            !                           
[4784]538            ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &
539           &       + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )
[2528]540            imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
541            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
542            ELSE                                  ;   sdjf%nrec_a(1) = imth
543            ENDIF
544            sdjf%nrec_a(2) = nmonth_half(   imth ) + nsec1jan000   ! swap at the middle of the month
545         ELSE                                    ! no time interpolation
546            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1
547            ELSE                                  ;   sdjf%nrec_a(1) = nmonth
548            ENDIF
549            sdjf%nrec_a(2) =  nmonth_end(nmonth  ) + nsec1jan000   ! swap at the end    of the month
550            sdjf%nrec_b(2) =  nmonth_end(nmonth-1) + nsec1jan000   ! beginning of the month (only for print)
[888]551         ENDIF
552         !
[2528]553         !                                   ! ================================ !
554      ELSE                                   ! higher frequency mean (in hours)
555         !                                   ! ================================ !
[888]556         !
[4245]557         ifreq_sec = NINT( sdjf%nfreqh * 3600 )                                         ! frequency mean (in seconds)
[2528]558         IF( sdjf%cltype(1:4) == 'week' )   isec_week = ksec_week( sdjf%cltype(6:8) )   ! since the first day of the current week
[1132]559         ! number of second since the beginning of the file
[2528]560         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month,wp)  ! since the first day of the current month
561         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ztmp = REAL(isec_week ,wp)  ! since the first day of the current week
562         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day  ,wp)  ! since 00h of the current day
563         ELSE                                           ;   ztmp = REAL(nsec_year ,wp)  ! since 00h on Jan 1 of the current year
[1132]564         ENDIF
[6140]565         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt + REAL( it_offset, wp )        ! centrered in the middle of sbc time step
566         ztmp = ztmp + 0.01 * rdt                                                       ! avoid truncation error
[1132]567         IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record
568            !
[3851]569            !          INT( ztmp/ifreq_sec + 0.5 )
[1132]570            !                     /|\
571            !                    2 |        *-----(
572            !                    1 |  *-----(
573            !                    0 |--(             
574            !                      |--+--|--+--|--+--|--> time
[3851]575            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
[1132]576            !                         |     |     |
577            !                         |     |     |
578            !       forcing record :  1     2     3
579            !                   
[2528]580            ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
581         ELSE                                   ! no time interpolation
[1132]582            !
[3851]583            !           INT( ztmp/ifreq_sec )
[1132]584            !                     /|\
585            !                    2 |           *-----(
586            !                    1 |     *-----(
587            !                    0 |-----(             
588            !                      |--+--|--+--|--+--|--> time
[3851]589            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
[1132]590            !                         |     |     |
591            !                         |     |     |
592            !       forcing record :  1     2     3
593            !                           
[2528]594            ztmp= ztmp / REAL(ifreq_sec, wp)
[1132]595         ENDIF
[3851]596         sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))   ! record number to be read
[1132]597
[2528]598         iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000       ! end of this record (in second)
599         ! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
600         IF( sdjf%cltype      == 'monthly' )   iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
601         IF( sdjf%cltype(1:4) == 'week'    )   iendrec = iendrec + ( nsec_year - isec_week )
602         IF( sdjf%cltype      == 'daily'   )   iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
603         IF( sdjf%ln_tint ) THEN
604             sdjf%nrec_a(2) = iendrec - ifreq_sec / 2        ! swap at the middle of the record
605         ELSE
606             sdjf%nrec_a(2) = iendrec                        ! swap at the end    of the record
607             sdjf%nrec_b(2) = iendrec - ifreq_sec            ! beginning of the record (only for print)
608         ENDIF
[888]609         !
610      ENDIF
611      !
[1132]612   END SUBROUTINE fld_rec
613
614
[7646]615   SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl )
[2528]616      !!---------------------------------------------------------------------
[3294]617      !!                    ***  ROUTINE fld_get  ***
[2528]618      !!
619      !! ** Purpose :   read the data
620      !!----------------------------------------------------------------------
[6140]621      TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables
622      TYPE(MAP_POINTER), INTENT(in   ) ::   map    ! global-to-local mapping indices
[7646]623      INTEGER  , INTENT(in), OPTIONAL  ::   jpk_bdy ! number of vertical levels in the bdy data
624      LOGICAL  , INTENT(in), OPTIONAL  ::   fvl     ! number of vertical levels in the bdy data
[6140]625      !
626      INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
627      INTEGER ::   iw       ! index into wgts array
628      INTEGER ::   ipdom    ! index of the domain
629      INTEGER ::   idvar    ! variable ID
630      INTEGER ::   idmspc   ! number of spatial dimensions
631      LOGICAL ::   lmoor    ! C1D case: point data
[2528]632      !!---------------------------------------------------------------------
[3851]633      !
[2528]634      ipk = SIZE( sdjf%fnow, 3 )
[3680]635      !
[3851]636      IF( ASSOCIATED(map%ptr) ) THEN
[7646]637         IF( PRESENT(jpk_bdy) ) THEN
638            IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2),                &
639                                                        sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )
640            ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ),                &
641                                                        sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )
642            ENDIF
643         ELSE
644            IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map )
645            ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map )
646            ENDIF
647         ENDIF       
[3294]648      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
[2528]649         CALL wgt_list( sdjf, iw )
[6140]650         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2),          & 
651            &                                                                          sdjf%nrec_a(1), sdjf%lsmname )
652         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,:  ),          &
653            &                                                                          sdjf%nrec_a(1), sdjf%lsmname )
[2528]654         ENDIF
655      ELSE
[6140]656         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN   ;   ipdom = jpdom_data
657         ELSE                                   ;   ipdom = jpdom_unknown
[3680]658         ENDIF
[4245]659         ! C1D case: If product of spatial dimensions == ipk, then x,y are of
660         ! size 1 (point/mooring data): this must be read onto the central grid point
661         idvar  = iom_varid( sdjf%num, sdjf%clvar )
[6140]662         idmspc = iom_file ( sdjf%num )%ndims( idvar )
[4245]663         IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1
[6140]664         lmoor  = (  idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk  )
[4245]665         !
[2528]666         SELECT CASE( ipk )
[3680]667         CASE(1)
[4245]668            IF( lk_c1d .AND. lmoor ) THEN
669               IF( sdjf%ln_tint ) THEN
670                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) )
[10425]671                  CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,1,2),'Z',1. )
[4245]672               ELSE
673                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1  ), sdjf%nrec_a(1) )
[10425]674                  CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,1  ),'Z',1. )
[4245]675               ENDIF
676            ELSE
677               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
678               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1  ), sdjf%nrec_a(1) )
679               ENDIF
[2528]680            ENDIF
681         CASE DEFAULT
[4245]682            IF (lk_c1d .AND. lmoor ) THEN
683               IF( sdjf%ln_tint ) THEN
684                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
[10425]685                  CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,:,2),'Z',1. )
[4245]686               ELSE
687                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,:  ), sdjf%nrec_a(1) )
[10425]688                  CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,:  ),'Z',1. )
[4245]689               ENDIF
690            ELSE
691               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
692               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) )
693               ENDIF
[2528]694            ENDIF
695         END SELECT
696      ENDIF
697      !
[3851]698      sdjf%rotn(2) = .false.   ! vector not yet rotated
[6140]699      !
[2528]700   END SUBROUTINE fld_get
701
[7646]702   SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl )
[3294]703      !!---------------------------------------------------------------------
[3851]704      !!                    ***  ROUTINE fld_map  ***
[3294]705      !!
706      !! ** Purpose :   read global data from file and map onto local data
707      !!                using a general mapping (for open boundaries)
708      !!----------------------------------------------------------------------
[7646]709
710      USE bdy_oce, ONLY: ln_bdy, idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz                 ! workspace to read in global data arrays
711
[3294]712      INTEGER                   , INTENT(in ) ::   num     ! stream number
713      CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name
[7646]714      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta     ! output field on model grid (2 dimensional)
[3294]715      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
[5132]716      TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices
[7646]717      INTEGER  , INTENT(in), OPTIONAL         ::   igrd, ibdy, jpk_bdy  ! grid type, set number and number of vertical levels in the bdy data
718      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl     ! grid type, set number and number of vertical levels in the bdy data
719      INTEGER                                 ::   jpkm1_bdy! number of vertical levels in the bdy data minus 1
[3294]720      !!
721      INTEGER                                 ::   ipi      ! length of boundary data on local process
722      INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 )
723      INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
724      INTEGER                                 ::   ilendta  ! length of data in file
725      INTEGER                                 ::   idvar    ! variable ID
[3651]726      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
[3294]727      INTEGER                                 ::   ierr
[7646]728      REAL(wp)                                ::   fv          ! fillvalue
729      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read    ! work space for global data
730      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_z  ! work space for global data
731      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_dz ! work space for global data
[3294]732      !!---------------------------------------------------------------------
[6140]733      !
[3294]734      ipi = SIZE( dta, 1 )
735      ipj = 1
736      ipk = SIZE( dta, 3 )
[6140]737      !
[3294]738      idvar   = iom_varid( num, clvar )
739      ilendta = iom_file(num)%dimsz(1,idvar)
[3651]740
[7646]741      IF ( ln_bdy ) THEN
742         ipj = iom_file(num)%dimsz(2,idvar)
743         IF( map%ll_unstruc) THEN   ! unstructured open boundary data file
744            dta_read => dta_global
745            IF( PRESENT(jpk_bdy) ) THEN
746               IF( jpk_bdy>0 ) THEN
747                  dta_read_z => dta_global_z
748                  dta_read_dz => dta_global_dz
749                  jpkm1_bdy = jpk_bdy-1
750               ENDIF
751            ENDIF
752         ELSE                       ! structured open boundary file
753            dta_read => dta_global2
754            IF( PRESENT(jpk_bdy) ) THEN
755               IF( jpk_bdy>0 ) THEN
756                  dta_read_z => dta_global2_z
757                  dta_read_dz => dta_global2_dz
758                  jpkm1_bdy = jpk_bdy-1
759               ENDIF
760            ENDIF
761         ENDIF
[3651]762      ENDIF
763
[6140]764      IF(lwp) WRITE(numout,*) 'Dim size for ',        TRIM(clvar),' is ', ilendta
[3294]765      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
[6140]766      !
[3294]767      SELECT CASE( ipk )
[7646]768      CASE(1)        ;   
769      CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
770         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
771            DO ib = 1, ipi
772              DO ik = 1, ipk
773                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
774              END DO
775            END DO
776         ELSE ! we assume that this is a structured open boundary file
777            DO ib = 1, ipi
778               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
779               ji=map%ptr(ib)-(jj-1)*ilendta
780               DO ik = 1, ipk
781                  dta(ib,1,ik) =  dta_read(ji,jj,ik)
782               END DO
783            END DO
784         ENDIF
785
786      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
787      ! Do we include something here to adjust barotropic velocities !
788      ! in case of a depth difference between bdy files and          !
789      ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0?  !
790      ! [as the enveloping and parital cells could change H]         !
791      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
792
793      CASE DEFAULT   ;
794
795      IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation
796         CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar )
797         dta_read(:,:,:) = -ABS(fv)
798         dta_read_z(:,:,:) = 0._wp
799         dta_read_dz(:,:,:) = 0._wp
800         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec )
801         SELECT CASE( igrd )                 
802            CASE(1)
803               CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
804               CALL iom_get ( num, jpdom_unknown, 'e3t',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
805            CASE(2) 
806               CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
807               CALL iom_get ( num, jpdom_unknown, 'e3u',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
808            CASE(3)
809               CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
810               CALL iom_get ( num, jpdom_unknown, 'e3v',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
811         END SELECT
812
813      IF ( ln_bdy ) & 
814         CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)
815
816      ELSE ! boundary data assumed to be on model grid
817         
818         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )                   
819         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
820            DO ib = 1, ipi
821              DO ik = 1, ipk
822                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
823              END DO
824            END DO
825         ELSE ! we assume that this is a structured open boundary file
826            DO ib = 1, ipi
827               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
828               ji=map%ptr(ib)-(jj-1)*ilendta
829               DO ik = 1, ipk
830                  dta(ib,1,ik) =  dta_read(ji,jj,ik)
831               END DO
832            END DO
833         ENDIF
834      ENDIF ! PRESENT(jpk_bdy)
[3294]835      END SELECT
[7646]836
837   END SUBROUTINE fld_map
838   
839   SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)
840
841      !!---------------------------------------------------------------------
842      !!                    ***  ROUTINE fld_bdy_interp  ***
843      !!
844      !! ** Purpose :   on the fly vertical interpolation to allow the use of
845      !!                boundary data from non-native vertical grid
846      !!----------------------------------------------------------------------
847      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation
848
849      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read      ! work space for global data
850      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z    ! work space for global data
851      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz   ! work space for global data
852      REAL(wp) , INTENT(in)                                ::   fv            ! fillvalue and alternative -ABS(fv)
853      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta                        ! output field on model grid (2 dimensional)
854      TYPE(MAP_POINTER)         , INTENT(in ) ::   map                        ! global-to-local mapping indices
855      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl                        ! grid type, set number and number of vertical levels in the bdy data
856      INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy        ! number of levels in bdy data
857      INTEGER  , INTENT(in)                   ::   ilendta                    ! length of data in file
858      !!
859      INTEGER                                 ::   ipi                        ! length of boundary data on local process
860      INTEGER                                 ::   ipj                        ! length of dummy dimension ( = 1 )
861      INTEGER                                 ::   ipk                        ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
862      INTEGER                                 ::   jpkm1_bdy                  ! number of levels in bdy data minus 1
863      INTEGER                                 ::   ib, ik, ikk                ! loop counters
864      INTEGER                                 ::   ji, jj, zij, zjj           ! temporary indices
865      REAL(wp)                                ::   zl, zi, zh                 ! tmp variable for current depth and interpolation factor
866      REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv)
867      CHARACTER (LEN=10)                      ::   ibstr
868      !!---------------------------------------------------------------------
869     
870
871      ipi       = SIZE( dta, 1 )
872      ipj       = SIZE( dta_read, 2 )
873      ipk       = SIZE( dta, 3 )
874      jpkm1_bdy = jpk_bdy-1
875     
876      fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
877      DO ib = 1, ipi
878            zij = idx_bdy(ibdy)%nbi(ib,igrd)
879            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
880         IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj
881      ENDDO
[3294]882      !
[7646]883      IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file
884
[3651]885         DO ib = 1, ipi
[7646]886            DO ik = 1, jpk_bdy
887               IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN
888                  dta_read_z(map%ptr(ib),1,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data
889                  dta_read_dz(map%ptr(ib),1,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct
890               ENDIF
891            ENDDO
892         ENDDO 
893
894         DO ib = 1, ipi
895            zij = idx_bdy(ibdy)%nbi(ib,igrd)
896            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
897            zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
898            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
899            SELECT CASE( igrd )                         
900               CASE(1)
901                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN
902                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
903                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
904                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
905                  ENDIF
906               CASE(2)
907                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN
908                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
909                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
910                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), &
911                       &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , &
912                        &                dta_read(map%ptr(ib),1,:)
913                  ENDIF
914               CASE(3)
915                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN
916                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
917                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
918                  ENDIF
919            END SELECT
920            DO ik = 1, ipk                     
921               SELECT CASE( igrd )                       
922                  CASE(1)
923                     zl =  gdept_n(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_n?
924                  CASE(2)
925                     IF(ln_sco) THEN
926                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
927                     ELSE
928                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 
929                     ENDIF
930                  CASE(3)
931                     IF(ln_sco) THEN
932                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
933                     ELSE
934                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
935                     ENDIF
936               END SELECT
937               IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN                                         ! above the first level of external data
938                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,1)
939               ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN                           ! below the last level of external data
940                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1))
941               ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1
942                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
943                     IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) &
944                    &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN
945                        zi           = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / &
946                       &               ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) )
947                        dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + &
948                       &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi
949                     ENDIF
950                  END DO
951               ENDIF
[3651]952            END DO
[3294]953         END DO
[7646]954
955         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
956            DO ib = 1, ipi
957              zij = idx_bdy(ibdy)%nbi(ib,igrd)
958              zjj = idx_bdy(ibdy)%nbj(ib,igrd)
959              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
960              ztrans = 0._wp
961              ztrans_new = 0._wp
962              DO ik = 1, jpk_bdy                            ! calculate transport on input grid
963                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
964              ENDDO
965              DO ik = 1, ipk                                ! calculate transport on model grid
966                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
967              ENDDO
968              DO ik = 1, ipk                                ! make transport correction
969                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
970                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
971                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
972                    IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) &
973                   &   CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at')
974                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik)
975                 ENDIF
976              ENDDO
977            ENDDO
978         ENDIF
979
980         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
981            DO ib = 1, ipi
982              zij = idx_bdy(ibdy)%nbi(ib,igrd)
983              zjj = idx_bdy(ibdy)%nbj(ib,igrd)
984              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
985              ztrans = 0._wp
986              ztrans_new = 0._wp
987              DO ik = 1, jpk_bdy                            ! calculate transport on input grid
988                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
989              ENDDO
990              DO ik = 1, ipk                                ! calculate transport on model grid
991                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
992              ENDDO
993              DO ik = 1, ipk                                ! make transport correction
994                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
995                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
996                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
997                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik)
998                 ENDIF
999              ENDDO
1000            ENDDO
1001         ENDIF
1002 
1003      ELSE ! structured open boundary file
1004
[3651]1005         DO ib = 1, ipi
[5132]1006            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
1007            ji=map%ptr(ib)-(jj-1)*ilendta
[7646]1008            DO ik = 1, jpk_bdy                     
1009               IF( ( dta_read(ji,jj,ik) == fv ) ) THEN
1010                  dta_read_z(ji,jj,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data
1011                  dta_read_dz(ji,jj,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct
1012               ENDIF
1013            ENDDO
1014         ENDDO 
1015       
1016
1017         DO ib = 1, ipi
1018            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
1019            ji=map%ptr(ib)-(jj-1)*ilendta
1020            zij = idx_bdy(ibdy)%nbi(ib,igrd)
1021            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
1022            zh  = SUM(dta_read_dz(ji,jj,:) )
1023            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
1024            SELECT CASE( igrd )                         
1025               CASE(1)
1026                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN
1027                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
1028                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
1029                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
1030                  ENDIF
1031               CASE(2)
1032                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN
1033                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
1034                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
1035                  ENDIF
1036               CASE(3)
1037                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN
1038                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
1039                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
1040                  ENDIF
1041            END SELECT
1042            DO ik = 1, ipk                     
1043               SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min
1044                  CASE(1)
1045                     zl =  gdept_n(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_n?
1046                  CASE(2)
1047                     IF(ln_sco) THEN
1048                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
1049                     ELSE
1050                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) )
1051                     ENDIF
1052                  CASE(3)
1053                     IF(ln_sco) THEN
1054                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
1055                     ELSE
1056                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
1057                     ENDIF
1058               END SELECT
1059               IF( zl < dta_read_z(ji,jj,1) ) THEN                                      ! above the first level of external data
1060                  dta(ib,1,ik) =  dta_read(ji,jj,1)
1061               ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                        ! below the last level of external data
1062                  dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1))
1063               ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1
1064                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
1065                     IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) &
1066                    &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN
1067                        zi           = ( zl - dta_read_z(ji,jj,ikk) ) / &
1068                       &               ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) )
1069                        dta(ib,1,ik) = dta_read(ji,jj,ikk) + &
1070                       &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi
1071                     ENDIF
1072                  END DO
1073               ENDIF
[3651]1074            END DO
1075         END DO
[3294]1076
[7646]1077         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
1078            DO ib = 1, ipi
1079               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
1080               ji=map%ptr(ib)-(jj-1)*ilendta
1081               zij = idx_bdy(ibdy)%nbi(ib,igrd)
1082               zjj = idx_bdy(ibdy)%nbj(ib,igrd)
1083               zh = SUM(dta_read_dz(ji,jj,:) )
1084               ztrans = 0._wp
1085               ztrans_new = 0._wp
1086               DO ik = 1, jpk_bdy                            ! calculate transport on input grid
1087                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
1088               ENDDO
1089               DO ik = 1, ipk                                ! calculate transport on model grid
1090                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
1091               ENDDO
1092               DO ik = 1, ipk                                ! make transport correction
1093                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
1094                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
1095                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
1096                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
1097                  ENDIF
1098               ENDDO
1099            ENDDO
1100         ENDIF
[3294]1101
[7646]1102         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
1103            DO ib = 1, ipi
1104               jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
1105               ji  = map%ptr(ib)-(jj-1)*ilendta
1106               zij = idx_bdy(ibdy)%nbi(ib,igrd)
1107               zjj = idx_bdy(ibdy)%nbj(ib,igrd)
1108               zh  = SUM(dta_read_dz(ji,jj,:) )
1109               ztrans = 0._wp
1110               ztrans_new = 0._wp
1111               DO ik = 1, jpk_bdy                            ! calculate transport on input grid
1112                  ztrans     = ztrans     + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
1113               ENDDO
1114               DO ik = 1, ipk                                ! calculate transport on model grid
1115                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
1116               ENDDO
1117               DO ik = 1, ipk                                ! make transport correction
1118                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
1119                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
1120                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
1121                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
1122                  ENDIF
1123               ENDDO
1124            ENDDO
1125         ENDIF
1126
1127      ENDIF ! endif unstructured or structured
1128
1129   END SUBROUTINE fld_bdy_interp
1130
1131
[2528]1132   SUBROUTINE fld_rot( kt, sd )
1133      !!---------------------------------------------------------------------
[3294]1134      !!                    ***  ROUTINE fld_rot  ***
[2528]1135      !!
1136      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
[2715]1137      !!----------------------------------------------------------------------
[6140]1138      INTEGER                , INTENT(in   ) ::   kt   ! ocean time step
1139      TYPE(FLD), DIMENSION(:), INTENT(inout) ::   sd   ! input field related variables
1140      !
1141      INTEGER ::   ju, jv, jk, jn  ! loop indices
1142      INTEGER ::   imf             ! size of the structure sd
1143      INTEGER ::   ill             ! character length
1144      INTEGER ::   iv              ! indice of V component
[9125]1145      CHARACTER (LEN=100)          ::   clcomp       ! dummy weight name
1146      REAL(wp), DIMENSION(jpi,jpj) ::   utmp, vtmp   ! temporary arrays for vector rotation
[2528]1147      !!---------------------------------------------------------------------
[6140]1148      !
[2528]1149      !! (sga: following code should be modified so that pairs arent searched for each time
1150      !
1151      imf = SIZE( sd )
1152      DO ju = 1, imf
1153         ill = LEN_TRIM( sd(ju)%vcomp )
[3851]1154         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
1155            IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN   ! find vector rotations required             
1156               IF( sd(ju)%vcomp(1:1) == 'U' ) THEN      ! east-west component has symbolic name starting with 'U'
1157                  ! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
1158                  clcomp = 'V' // sd(ju)%vcomp(2:ill)   ! works even if ill == 1
1159                  iv = -1
1160                  DO jv = 1, imf
1161                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
1162                  END DO
1163                  IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together
1164                     DO jk = 1, SIZE( sd(ju)%fnow, 3 )
1165                        IF( sd(ju)%ln_tint )THEN
1166                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
1167                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
1168                           sd(ju)%fdta(:,:,jk,jn) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
1169                        ELSE
1170                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->i', utmp(:,:) )
1171                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->j', vtmp(:,:) )
1172                           sd(ju)%fnow(:,:,jk   ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk   ) = vtmp(:,:)
1173                        ENDIF
1174                     END DO
1175                     sd(ju)%rotn(jn) = .TRUE.               ! vector was rotated
1176                     IF( lwp .AND. kt == nit000 )   WRITE(numout,*)   &
1177                        &   'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
1178                  ENDIF
1179               ENDIF
1180            ENDIF
1181         END DO
[2528]1182       END DO
[2715]1183      !
[2528]1184   END SUBROUTINE fld_rot
1185
1186
[1628]1187   SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
[1132]1188      !!---------------------------------------------------------------------
1189      !!                    ***  ROUTINE fld_clopn  ***
1190      !!
[7646]1191      !! ** Purpose :   update the file name and close/open the files
[1132]1192      !!----------------------------------------------------------------------
[2715]1193      TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
[3851]1194      INTEGER, OPTIONAL, INTENT(in   ) ::   kyear    ! year value
1195      INTEGER, OPTIONAL, INTENT(in   ) ::   kmonth   ! month value
1196      INTEGER, OPTIONAL, INTENT(in   ) ::   kday     ! day value
[2715]1197      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
[6140]1198      !
[3851]1199      LOGICAL :: llprevyr              ! are we reading previous year  file?
1200      LOGICAL :: llprevmth             ! are we reading previous month file?
1201      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
1202      INTEGER :: isec_week             ! number of seconds since start of the weekly file
1203      INTEGER :: indexyr               ! year undex (O/1/2: previous/current/next)
1204      INTEGER :: iyear_len, imonth_len ! length (days) of iyear and imonth             !
1205      CHARACTER(len = 256)::   clname  ! temporary file name
[2715]1206      !!----------------------------------------------------------------------
[3851]1207      IF( PRESENT(kyear) ) THEN                             ! use given values
1208         iyear = kyear
1209         imonth = kmonth
1210         iday = kday
[5751]1211         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
[5749]1212            isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 
[5751]1213            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
[5749]1214            llprevyr   = llprevmth .AND. nmonth == 1
1215            iyear  = nyear  - COUNT((/llprevyr /))
1216            imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
[5751]1217            iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
[5749]1218         ENDIF
[3851]1219      ELSE                                                  ! use current day values
1220         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
1221            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
1222            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
1223            llprevyr   = llprevmth .AND. nmonth == 1
1224         ELSE
1225            isec_week  = 0
1226            llprevmth  = .FALSE.
1227            llprevyr   = .FALSE.
1228         ENDIF
1229         iyear  = nyear  - COUNT((/llprevyr /))
1230         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
1231         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
1232      ENDIF
[1132]1233
1234      ! build the new filename if not climatological data
[3851]1235      clname=TRIM(sdjf%clrootname)
[2528]1236      !
[3851]1237      ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
[2528]1238      IF( .NOT. sdjf%ln_clim ) THEN   
[3851]1239                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear    ! add year
1240         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname          ), imonth   ! add month
[2528]1241      ELSE
1242         ! build the new filename if climatological data
[3851]1243         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth   ! add month
[888]1244      ENDIF
[2528]1245      IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
[3851]1246            &                            WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), iday     ! add day
[2528]1247      !
[3851]1248      IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open
[6140]1249         !
[3851]1250         sdjf%clname = TRIM(clname)
1251         IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
1252         CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
[6140]1253         !
[3851]1254         ! find the last record to be read -> update sdjf%nreclast
1255         indexyr = iyear - nyear + 1
1256         iyear_len = nyear_len( indexyr )
1257         SELECT CASE ( indexyr )
1258         CASE ( 0 )   ;   imonth_len = 31   ! previous year -> imonth = 12
1259         CASE ( 1 )   ;   imonth_len = nmonth_len(imonth) 
1260         CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1
1261         END SELECT
[6140]1262         !
[3851]1263         ! last record to be read in the current file
1264         IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean
1265         ELSEIF( sdjf%nfreqh ==  -1 ) THEN                                          ! monthly mean
1266            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = 1
1267            ELSE                                           ;   sdjf%nreclast = 12
1268            ENDIF
1269         ELSE                                                                       ! higher frequency mean (in hours)
[4245]1270            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
1271            ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24 * 7          / sdjf%nfreqh )
1272            ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24              / sdjf%nfreqh )
1273            ELSE                                           ;   sdjf%nreclast = NINT( 24 * iyear_len  / sdjf%nfreqh )
[3851]1274            ENDIF
1275         ENDIF
[6140]1276         !
[3851]1277      ENDIF
1278      !
[1132]1279   END SUBROUTINE fld_clopn
1280
1281
[7646]1282   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint )
[1132]1283      !!---------------------------------------------------------------------
1284      !!                    ***  ROUTINE fld_fill  ***
1285      !!
[7646]1286      !! ** Purpose :   fill the data structure (sdf) with the associated information
1287      !!              read in namelist (sdf_n) and control print
[1132]1288      !!----------------------------------------------------------------------
[7646]1289      TYPE(FLD)  , DIMENSION(:)          , INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
1290      TYPE(FLD_N), DIMENSION(:)          , INTENT(in   ) ::   sdf_n      ! array of namelist information structures
1291      CHARACTER(len=*)                   , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
1292      CHARACTER(len=*)                   , INTENT(in   ) ::   cdcaller   ! name of the calling routine
1293      CHARACTER(len=*)                   , INTENT(in   ) ::   cdtitle    ! description of the calling routine
1294      CHARACTER(len=*)                   , INTENT(in   ) ::   cdnam      ! name of the namelist from which sdf_n comes
1295      INTEGER                  , OPTIONAL, INTENT(in   ) ::   knoprint   ! no calling routine information printed
[888]1296      !
[7646]1297      INTEGER  ::   jf   ! dummy indices
[1132]1298      !!---------------------------------------------------------------------
[6140]1299      !
[1132]1300      DO jf = 1, SIZE(sdf)
1301         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
[3851]1302         sdf(jf)%clname     = "not yet defined"
[1730]1303         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
[1132]1304         sdf(jf)%clvar      = sdf_n(jf)%clvar
1305         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
1306         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
[2528]1307         sdf(jf)%cltype     = sdf_n(jf)%cltype
[3851]1308         sdf(jf)%num        = -1
1309         sdf(jf)%wgtname    = " "
[1730]1310         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
[4230]1311         sdf(jf)%lsmname = " "
1312         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
[3851]1313         sdf(jf)%vcomp      = sdf_n(jf)%vcomp
1314         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
1315         IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   &
1316            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
1317         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
1318            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
[4371]1319         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
[1132]1320      END DO
[6140]1321      !
[1132]1322      IF(lwp) THEN      ! control print
1323         WRITE(numout,*)
[7646]1324         IF( .NOT.PRESENT( knoprint) ) THEN
1325            WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
1326            WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
1327         ENDIF
1328         WRITE(numout,*) '   fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
1329         WRITE(numout,*) '   ~~~~~~~~'
1330         WRITE(numout,*) '      list of files and frequency (>0: in hours ; <0 in months)'
[1132]1331         DO jf = 1, SIZE(sdf)
[7646]1332            WRITE(numout,*) '      root filename: '  , TRIM( sdf(jf)%clrootname ), '   variable name: ', TRIM( sdf(jf)%clvar )
1333            WRITE(numout,*) '         frequency: '      ,       sdf(jf)%nfreqh      ,   &
1334               &                  '   time interp: '    ,       sdf(jf)%ln_tint     ,   &
1335               &                  '   climatology: '    ,       sdf(jf)%ln_clim
1336            WRITE(numout,*) '         weights: '        , TRIM( sdf(jf)%wgtname    ),   &
1337               &                  '   pairing: '        , TRIM( sdf(jf)%vcomp      ),   &
1338               &                  '   data type: '      ,       sdf(jf)%cltype      ,   &
1339               &                  '   land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
[2528]1340            call flush(numout)
[1132]1341         END DO
1342      ENDIF
[6140]1343      !
[1132]1344   END SUBROUTINE fld_fill
1345
1346
[1275]1347   SUBROUTINE wgt_list( sd, kwgt )
1348      !!---------------------------------------------------------------------
1349      !!                    ***  ROUTINE wgt_list  ***
1350      !!
[7646]1351      !! ** Purpose :   search array of WGTs and find a weights file entry,
1352      !!                or return a new one adding it to the end if new entry.
1353      !!                the weights data is read in and restructured (fld_weight)
[1275]1354      !!----------------------------------------------------------------------
[2715]1355      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
1356      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
[6140]1357      !
[2715]1358      INTEGER ::   kw, nestid   ! local integer
1359      LOGICAL ::   found        ! local logical
[1275]1360      !!----------------------------------------------------------------------
1361      !
1362      !! search down linked list
1363      !! weights filename is either present or we hit the end of the list
1364      found = .FALSE.
[6140]1365      !
[1275]1366      !! because agrif nest part of filenames are now added in iom_open
1367      !! to distinguish between weights files on the different grids, need to track
1368      !! nest number explicitly
1369      nestid = 0
1370#if defined key_agrif
1371      nestid = Agrif_Fixed()
1372#endif
1373      DO kw = 1, nxt_wgt-1
1374         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
1375             ref_wgts(kw)%nestid == nestid) THEN
1376            kwgt = kw
1377            found = .TRUE.
1378            EXIT
1379         ENDIF
1380      END DO
1381      IF( .NOT.found ) THEN
1382         kwgt = nxt_wgt
1383         CALL fld_weight( sd )
1384      ENDIF
[2715]1385      !
[1275]1386   END SUBROUTINE wgt_list
1387
[2715]1388
[1275]1389   SUBROUTINE wgt_print( )
1390      !!---------------------------------------------------------------------
1391      !!                    ***  ROUTINE wgt_print  ***
1392      !!
1393      !! ** Purpose :   print the list of known weights
1394      !!----------------------------------------------------------------------
[2715]1395      INTEGER ::   kw   !
[1275]1396      !!----------------------------------------------------------------------
1397      !
1398      DO kw = 1, nxt_wgt-1
1399         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
1400         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1401         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
1402         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
1403         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
1404         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
1405         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
1406         IF( ref_wgts(kw)%cyclic ) THEN
1407            WRITE(numout,*) '       cyclical'
[2528]1408            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap
[1275]1409         ELSE
1410            WRITE(numout,*) '       not cyclical'
1411         ENDIF
1412         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
1413      END DO
[2715]1414      !
[1275]1415   END SUBROUTINE wgt_print
1416
[2715]1417
[1275]1418   SUBROUTINE fld_weight( sd )
1419      !!---------------------------------------------------------------------
1420      !!                    ***  ROUTINE fld_weight  ***
1421      !!
[7646]1422      !! ** Purpose :   create a new WGT structure and fill in data from file,
1423      !!              restructuring as required
[1275]1424      !!----------------------------------------------------------------------
[2715]1425      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
1426      !!
[6140]1427      INTEGER ::   jn         ! dummy loop indices
1428      INTEGER ::   inum       ! local logical unit
1429      INTEGER ::   id         ! local variable id
1430      INTEGER ::   ipk        ! local vertical dimension
1431      INTEGER ::   zwrap      ! local integer
1432      LOGICAL ::   cyclical   !
1433      CHARACTER (len=5) ::   aname   !
[5399]1434      INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims
[9125]1435      INTEGER,  DIMENSION(jpi,jpj) ::   data_src
1436      REAL(wp), DIMENSION(jpi,jpj) ::   data_tmp
[1275]1437      !!----------------------------------------------------------------------
1438      !
1439      IF( nxt_wgt > tot_wgts ) THEN
[2777]1440        CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
[1275]1441      ENDIF
1442      !
1443      !! new weights file entry, add in extra information
1444      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
1445      !! input data file is representative of all other files to be opened and processed with the
1446      !! current weights file
1447
1448      !! open input data file (non-model grid)
[1319]1449      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
[1275]1450
1451      !! get dimensions
[5399]1452      IF ( SIZE(sd%fnow, 3) > 1 ) THEN
1453         ALLOCATE( ddims(4) )
1454      ELSE
1455         ALLOCATE( ddims(3) )
1456      ENDIF
[1275]1457      id = iom_varid( inum, sd%clvar, ddims )
1458
[2528]1459      !! close it
1460      CALL iom_close( inum )
[1275]1461
[2528]1462      !! now open the weights file
[1275]1463
[2528]1464      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
1465      IF ( inum > 0 ) THEN
[1275]1466
[2528]1467         !! determine whether we have an east-west cyclic grid
1468         !! from global attribute called "ew_wrap" in the weights file
1469         !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1470         !! since this is the most common forcing configuration
[1275]1471
[2528]1472         CALL iom_getatt(inum, 'ew_wrap', zwrap)
1473         IF( zwrap >= 0 ) THEN
[1275]1474            cyclical = .TRUE.
[2528]1475         ELSE IF( zwrap == -999 ) THEN
[1275]1476            cyclical = .TRUE.
[2528]1477            zwrap = 0
1478         ELSE
1479            cyclical = .FALSE.
[1275]1480         ENDIF
1481
1482         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1483         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1484         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
[2528]1485         ref_wgts(nxt_wgt)%overlap = zwrap
1486         ref_wgts(nxt_wgt)%cyclic = cyclical
[1275]1487         ref_wgts(nxt_wgt)%nestid = 0
1488#if defined key_agrif
1489         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1490#endif
1491         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1492         !! for each weight wgtNN there is an integer array srcNN which gives the point in
1493         !! the input data grid which is to be multiplied by the weight
1494         !! they are both arrays on the model grid so the result of the multiplication is
1495         !! added into an output array on the model grid as a running sum
1496
1497         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1498         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1499         IF( id <= 0) THEN
1500            ref_wgts(nxt_wgt)%numwgt = 4
1501         ELSE
1502            ref_wgts(nxt_wgt)%numwgt = 16
1503         ENDIF
1504
1505         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1506         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1507         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1508
1509         DO jn = 1,4
1510            aname = ' '
1511            WRITE(aname,'(a3,i2.2)') 'src',jn
1512            data_tmp(:,:) = 0
[1955]1513            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
[1275]1514            data_src(:,:) = INT(data_tmp(:,:))
1515            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1516            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1517         END DO
1518
1519         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1520            aname = ' '
1521            WRITE(aname,'(a3,i2.2)') 'wgt',jn
[1955]1522            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1523            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
[1275]1524         END DO
1525         CALL iom_close (inum)
1526 
1527         ! find min and max indices in grid
[1955]1528         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1529         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1530         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1531         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
[1275]1532
1533         ! and therefore dimensions of the input box
1534         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1535         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1536
1537         ! shift indexing of source grid
1538         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1539         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1540
1541         ! create input grid, give it a halo to allow gradient calculations
[1702]1542         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1543         ! a more robust solution will be given in next release
[2528]1544         ipk =  SIZE(sd%fnow, 3)
1545         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1546         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
[6140]1547         !
[1275]1548         nxt_wgt = nxt_wgt + 1
[6140]1549         !
[1275]1550      ELSE
1551         CALL ctl_stop( '    fld_weight : unable to read the file ' )
1552      ENDIF
1553
[5399]1554      DEALLOCATE (ddims )
[2715]1555      !
[1275]1556   END SUBROUTINE fld_weight
1557
[2715]1558
[6140]1559   SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm,   &
[7646]1560      &                          jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
[1275]1561      !!---------------------------------------------------------------------
[4230]1562      !!                    ***  ROUTINE apply_seaoverland  ***
1563      !!
1564      !! ** Purpose :   avoid spurious fluxes in coastal or near-coastal areas
1565      !!                due to the wrong usage of "land" values from the coarse
1566      !!                atmospheric model when spatial interpolation is required
1567      !!      D. Delrosso INGV         
1568      !!----------------------------------------------------------------------
[6140]1569      INTEGER,                   INTENT(in   ) :: itmpi,itmpj,itmpz                    ! lengths
1570      INTEGER,                   INTENT(in   ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
1571      INTEGER, DIMENSION(3),     INTENT(in   ) :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
1572      REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo                              ! input/output array for seaoverland application
1573      CHARACTER (len=100),       INTENT(in   ) :: clmaskfile                           ! land/sea mask file name
1574      !
1575      INTEGER :: inum,jni,jnj,jnz,jc   ! local indices
1576      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1             ! local array for land point detection
1577      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfieldn   ! array of forcing field with undeff for land points
1578      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfield    ! array of forcing field
[4230]1579      !!---------------------------------------------------------------------
[6140]1580      !
1581      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) )
1582      !
[4230]1583      ! Retrieve the land sea mask data
1584      CALL iom_open( clmaskfile, inum )
1585      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1586      CASE(1)
[6140]1587         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
[4230]1588      CASE DEFAULT
[6140]1589         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
[4230]1590      END SELECT
1591      CALL iom_close( inum )
[6140]1592      !
1593      DO jnz=1,rec1_lsm(3)             !! Loop over k dimension
1594         !
1595         DO jni = 1, itmpi                               !! copy the original field into a tmp array
1596            DO jnj = 1, itmpj                            !! substituting undeff over land points
1597               zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1598               IF( zslmec1(jni,jnj,jnz) == 1. )   zfieldn(jni,jnj) = undeff_lsm
[4230]1599            END DO
1600         END DO
[6140]1601         !
1602         CALL seaoverland( zfieldn, itmpi, itmpj, zfield )
1603         DO jc = 1, nn_lsm
1604            CALL seaoverland( zfield, itmpi, itmpj, zfield )
1605         END DO
1606         !
1607         !   Check for Undeff and substitute original values
1608         IF( ANY(zfield==undeff_lsm) ) THEN
1609            DO jni = 1, itmpi
1610               DO jnj = 1, itmpj
1611                  IF( zfield(jni,jnj)==undeff_lsm )   zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1612               END DO
1613            END DO
1614         ENDIF
1615         !
1616         zfieldo(:,:,jnz) = zfield(:,:)
1617         !
1618      END DO                           !! End Loop over k dimension
1619      !
1620      DEALLOCATE ( zslmec1, zfieldn, zfield )
1621      !
[4230]1622   END SUBROUTINE apply_seaoverland 
1623
1624
[6140]1625   SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield )
[4230]1626      !!---------------------------------------------------------------------
1627      !!                    ***  ROUTINE seaoverland  ***
1628      !!
1629      !! ** Purpose :   create shifted matrices for seaoverland application 
1630      !!      D. Delrosso INGV
1631      !!----------------------------------------------------------------------
[6140]1632      INTEGER                      , INTENT(in   ) :: ileni,ilenj   ! lengths
1633      REAL, DIMENSION (ileni,ilenj), INTENT(in   ) :: zfieldn       ! array of forcing field with undeff for land points
1634      REAL, DIMENSION (ileni,ilenj), INTENT(  out) :: zfield        ! array of forcing field
1635      !
1636      REAL   , DIMENSION (ileni,ilenj)   :: zmat1, zmat2, zmat3, zmat4  ! local arrays
1637      REAL   , DIMENSION (ileni,ilenj)   :: zmat5, zmat6, zmat7, zmat8  !   -     -
1638      REAL   , DIMENSION (ileni,ilenj)   :: zlsm2d                      !   -     -
1639      REAL   , DIMENSION (ileni,ilenj,8) :: zlsm3d                      !   -     -
1640      LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d                 ! logical mask for undeff detection
1641      LOGICAL, DIMENSION (ileni,ilenj)   :: ll_msknan2d                 ! logical mask for undeff detection
[4230]1642      !!----------------------------------------------------------------------
[6140]1643      zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/)     , DIM=2 )
1644      zmat1 = eoshift( zmat8   , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/)       , DIM=1 )
1645      zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/)     , DIM=1 )
1646      zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 )
1647      zmat3 = eoshift( zmat4   , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/)       , DIM=1 )
1648      zmat5 = eoshift( zmat4   , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/)   , DIM=1 )
1649      zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 )
1650      zmat7 = eoshift( zmat8   , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/)   , DIM=1 )
1651      !
[4230]1652      zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
[6140]1653      ll_msknan3d = .NOT.( zlsm3d  == undeff_lsm )
1654      ll_msknan2d = .NOT.( zfieldn == undeff_lsm )  ! FALSE where is Undeff (land)
1655      zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) )
1656      WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp )   zlsm2d = undeff_lsm
1657      zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d )
1658      !
[4230]1659   END SUBROUTINE seaoverland
1660
1661
1662   SUBROUTINE fld_interp( num, clvar, kw, kk, dta,  &
1663                          &         nrec, lsmfile)     
1664      !!---------------------------------------------------------------------
[1275]1665      !!                    ***  ROUTINE fld_interp  ***
1666      !!
1667      !! ** Purpose :   apply weights to input gridded data to create data
1668      !!                on model grid
1669      !!----------------------------------------------------------------------
[2715]1670      INTEGER                   , INTENT(in   ) ::   num     ! stream number
1671      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
1672      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
1673      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
1674      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
1675      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
[4230]1676      CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
[6140]1677      !
1678      INTEGER, DIMENSION(3) ::   rec1, recn           ! temporary arrays for start and length
1679      INTEGER, DIMENSION(3) ::   rec1_lsm, recn_lsm   ! temporary arrays for start and length in case of seaoverland
1680      INTEGER ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2    ! temporary indices
1681      INTEGER ::   jk, jn, jm, jir, jjr               ! loop counters
1682      INTEGER ::   ni, nj                             ! lengths
1683      INTEGER ::   jpimin,jpiwid                      ! temporary indices
1684      INTEGER ::   jpimin_lsm,jpiwid_lsm              ! temporary indices
1685      INTEGER ::   jpjmin,jpjwid                      ! temporary indices
1686      INTEGER ::   jpjmin_lsm,jpjwid_lsm              ! temporary indices
1687      INTEGER ::   jpi1,jpi2,jpj1,jpj2                ! temporary indices
1688      INTEGER ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
1689      INTEGER ::   itmpi,itmpj,itmpz                     ! lengths
1690      REAL(wp),DIMENSION(:,:,:), ALLOCATABLE ::   ztmp_fly_dta                 ! local array of values on input grid     
[1275]1691      !!----------------------------------------------------------------------
1692      !
1693      !! for weighted interpolation we have weights at four corners of a box surrounding
1694      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
1695      !! or by a grid value and gradients at the corner point (bicubic case)
1696      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1697
[2528]1698      !! sub grid from non-model input grid which encloses all grid points in this nemo process
[1275]1699      jpimin = ref_wgts(kw)%botleft(1)
1700      jpjmin = ref_wgts(kw)%botleft(2)
1701      jpiwid = ref_wgts(kw)%jpiwgt
1702      jpjwid = ref_wgts(kw)%jpjwgt
1703
[2528]1704      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
[1275]1705      rec1(1) = MAX( jpimin-1, 1 )
1706      rec1(2) = MAX( jpjmin-1, 1 )
[2528]1707      rec1(3) = 1
[1275]1708      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1709      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
[2528]1710      recn(3) = kk
[1275]1711
[2528]1712      !! where we need to put it in the non-nemo grid fly_dta
1713      !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1714      !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
[1275]1715      jpi1 = 2 + rec1(1) - jpimin
1716      jpj1 = 2 + rec1(2) - jpjmin
1717      jpi2 = jpi1 + recn(1) - 1
1718      jpj2 = jpj1 + recn(2) - 1
1719
1720
[4230]1721      IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1722      !! indeces for ztmp_fly_dta
1723      ! --------------------------
1724         rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1)  ! starting index for enlarged external data, x direction
1725         rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1)  ! starting index for enlarged external data, y direction
1726         rec1_lsm(3) = 1                    ! vertical dimension
1727         recn_lsm(1)=MIN(rec1(1)-rec1_lsm(1)+recn(1)+nn_lsm,ref_wgts(kw)%ddims(1)-rec1_lsm(1)) ! n points in x direction
1728         recn_lsm(2)=MIN(rec1(2)-rec1_lsm(2)+recn(2)+nn_lsm,ref_wgts(kw)%ddims(2)-rec1_lsm(2)) ! n points in y direction
1729         recn_lsm(3) = kk                   ! number of vertical levels in the input file
1730
1731      !  Avoid out of bound
1732         jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1733         jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1734         jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1735         jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1736
1737         jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1738         jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1739         jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1740         jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1741
1742
[6140]1743         itmpi=jpi2_lsm-jpi1_lsm+1
1744         itmpj=jpj2_lsm-jpj1_lsm+1
[4230]1745         itmpz=kk
1746         ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1747         ztmp_fly_dta(:,:,:) = 0.0
1748         SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1749         CASE(1)
1750              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   &
1751                 &                                                                nrec, rec1_lsm, recn_lsm)
1752         CASE DEFAULT
1753              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   &
1754                 &                                                                nrec, rec1_lsm, recn_lsm)
1755         END SELECT
1756         CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),                  &
1757                 &                                      jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm,                  &
1758                 &                                      itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1759
1760
1761         ! Relative indeces for remapping
1762         ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1763         ii_lsm2 = (ii_lsm1+recn(1))-1
1764         ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1765         ij_lsm2 = (ij_lsm1+recn(2))-1
1766
1767         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1768         ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1769         DEALLOCATE(ztmp_fly_dta)
1770
1771      ELSE
1772         
1773         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1774         SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1775         CASE(1)
1776              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1777         CASE DEFAULT
1778              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1779         END SELECT
1780      ENDIF
1781     
1782
[1275]1783      !! first four weights common to both bilinear and bicubic
[2528]1784      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
[1275]1785      !! note that we have to offset by 1 into fly_dta array because of halo
[2528]1786      dta(:,:,:) = 0.0
[1275]1787      DO jk = 1,4
1788        DO jn = 1, jpj
1789          DO jm = 1,jpi
1790            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1791            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]1792            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
[1275]1793          END DO
1794        END DO
1795      END DO
1796
1797      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1798
1799        !! fix up halo points that we couldnt read from file
1800        IF( jpi1 == 2 ) THEN
[2528]1801           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
[1275]1802        ENDIF
1803        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]1804           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
[1275]1805        ENDIF
1806        IF( jpj1 == 2 ) THEN
[2528]1807           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
[1275]1808        ENDIF
1809        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
[2528]1810           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
[1275]1811        ENDIF
1812
1813        !! if data grid is cyclic we can do better on east-west edges
1814        !! but have to allow for whether first and last columns are coincident
1815        IF( ref_wgts(kw)%cyclic ) THEN
1816           rec1(2) = MAX( jpjmin-1, 1 )
[2528]1817           recn(1) = 1
[1275]1818           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1819           jpj1 = 2 + rec1(2) - jpjmin
1820           jpj2 = jpj1 + recn(2) - 1
1821           IF( jpi1 == 2 ) THEN
[2528]1822              rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1823              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1824              CASE(1)
1825                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1826              CASE DEFAULT
1827                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1828              END SELECT     
1829              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]1830           ENDIF
1831           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]1832              rec1(1) = 1 + ref_wgts(kw)%overlap
1833              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1834              CASE(1)
1835                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1836              CASE DEFAULT
1837                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1838              END SELECT
1839              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]1840           ENDIF
1841        ENDIF
1842
1843        ! gradient in the i direction
1844        DO jk = 1,4
1845          DO jn = 1, jpj
1846            DO jm = 1,jpi
1847              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1848              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]1849              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1850                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
[1275]1851            END DO
1852          END DO
1853        END DO
1854
1855        ! gradient in the j direction
1856        DO jk = 1,4
1857          DO jn = 1, jpj
1858            DO jm = 1,jpi
1859              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1860              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]1861              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         &
1862                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
[1275]1863            END DO
1864          END DO
1865        END DO
1866
[2715]1867         ! gradient in the ij direction
1868         DO jk = 1,4
1869            DO jn = 1, jpj
1870               DO jm = 1,jpi
1871                  ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1872                  nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1873                  dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
[2528]1874                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1875                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
[2715]1876               END DO
[1275]1877            END DO
[2715]1878         END DO
1879         !
[1275]1880      END IF
[2715]1881      !
[1275]1882   END SUBROUTINE fld_interp
[2528]1883
1884
1885   FUNCTION ksec_week( cdday )
1886      !!---------------------------------------------------------------------
1887      !!                    ***  FUNCTION kshift_week ***
1888      !!
[7646]1889      !! ** Purpose :   return the first 3 letters of the first day of the weekly file
[2528]1890      !!---------------------------------------------------------------------
[7646]1891      CHARACTER(len=*), INTENT(in)   ::   cdday   ! first 3 letters of the first day of the weekly file
[2528]1892      !!
[7646]1893      INTEGER                        ::   ksec_week      ! output variable
1894      INTEGER                        ::   ijul, ishift   ! local integer
[2528]1895      CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
1896      !!----------------------------------------------------------------------
1897      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1898      DO ijul = 1, 7
1899         IF( cl_week(ijul) == TRIM(cdday) ) EXIT
[2715]1900      END DO
[2528]1901      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1902      !
1903      ishift = ijul * NINT(rday)
1904      !
1905      ksec_week = nsec_week + ishift
1906      ksec_week = MOD( ksec_week, 7*NINT(rday) )
1907      !
1908   END FUNCTION ksec_week
1909
[2715]1910   !!======================================================================
[888]1911END MODULE fldread
Note: See TracBrowser for help on using the repository browser.