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

Last change on this file since 12489 was 12489, checked in by davestorkey, 9 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

  • Property svn:keywords set to Id
File size: 88.1 KB
Line 
1MODULE fldread
2   !!======================================================================
3   !!                       ***  MODULE  fldread  ***
4   !! Ocean forcing:  read input field for surface boundary condition
5   !!=====================================================================
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
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   fld_read      : read input fields used for the computation of the surface boundary condition
14   !!   fld_init      : initialization of field read
15   !!   fld_def       : define the record(s) of the file and its name
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     : 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   !!   fld_filename  : define the filename according to a given date
28   !!   ksec_week     : function returning seconds between 00h of the beginning of the week and half of the current time step
29   !!----------------------------------------------------------------------
30   USE oce            ! ocean dynamics and tracers
31   USE dom_oce        ! ocean space and time domain
32   USE phycst         ! physical constant
33   USE sbc_oce        ! surface boundary conditions : fields
34   USE geo2ocean      ! for vector rotation on to model grid
35   !
36   USE in_out_manager ! I/O manager
37   USE iom            ! I/O manager library
38   USE ioipsl  , ONLY : ymds2ju, ju2ymds   ! for calendar
39   USE lib_mpp        ! MPP library
40   USE lbclnk         ! ocean lateral boundary conditions (C1D case)
41   
42   IMPLICIT NONE
43   PRIVATE   
44 
45   PUBLIC   fld_map    ! routine called by tides_init
46   PUBLIC   fld_read, fld_fill   ! called by sbc... modules
47   PUBLIC   fld_def
48
49   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
50      CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file
51      REAL(wp)             ::   freqh       ! frequency of each flux file
52      CHARACTER(len = 34)  ::   clvar       ! generic name of the variable in the NetCDF flux file
53      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F)
54      LOGICAL              ::   ln_clim     ! climatology or not (T/F)
55      CHARACTER(len = 8)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly'
56      CHARACTER(len = 256) ::   wname       ! generic name of a NetCDF weights file to be used, blank if not
57      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation
58      !                                     ! a string starting with "U" or "V" for each component   
59      !                                     ! chars 2 onwards identify which components go together 
60      CHARACTER(len = 34)  ::   lname       ! generic name of a NetCDF land/sea mask file to be used, blank if not
61      !                                     ! 0=sea 1=land
62   END TYPE FLD_N
63
64   TYPE, PUBLIC ::   FLD        !: Input field related variables
65      CHARACTER(len = 256)            ::   clrootname   ! generic name of the NetCDF file
66      CHARACTER(len = 256)            ::   clname       ! current name of the NetCDF file
67      REAL(wp)                        ::   freqh        ! frequency of each flux file
68      CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file
69      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F)
70      LOGICAL                         ::   ln_clim      ! climatology or not (T/F)
71      CHARACTER(len = 8)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly'
72      INTEGER                         ::   num          ! iom id of the jpfld files to be read
73      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
74      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year)
75      INTEGER , ALLOCATABLE, DIMENSION(:      ) ::   nrecsec   !
76      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow   ! input fields interpolated to now time step
77      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta   ! 2 consecutive record of input fields
78      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
79      !                                                 ! into the WGTLIST structure
80      CHARACTER(len = 34)             ::   vcomp        ! symbolic name for a vector component that needs rotation
81      LOGICAL, DIMENSION(2)           ::   rotn         ! flag to indicate whether before/after field has been rotated
82      INTEGER                         ::   nreclast     ! last record to be read in the current file
83      CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key
84      !                                                 !
85      !                                                 ! Variables related to BDY
86      INTEGER                         ::   igrd         !   grid type for bdy data
87      INTEGER                         ::   ibdy         !   bdy set id number
88      INTEGER, POINTER, DIMENSION(:)  ::   imap         !   Array of integer pointers to 1D arrays
89      LOGICAL                         ::   ltotvel      !   total velocity or not (T/F)
90      LOGICAL                         ::   lzint        !   T if it requires a vertical interpolation
91   END TYPE FLD
92
93!$AGRIF_DO_NOT_TREAT
94
95   !! keep list of all weights variables so they're only read in once
96   !! need to add AGRIF directives not to process this structure
97   !! also need to force wgtname to include AGRIF nest number
98   TYPE         ::   WGT        !: Input weights related variables
99      CHARACTER(len = 256)                    ::   wgtname      ! current name of the NetCDF weight file
100      INTEGER , DIMENSION(2)                  ::   ddims        ! shape of input grid
101      INTEGER , DIMENSION(2)                  ::   botleft      ! top left corner of box in input grid containing
102      !                                                         ! current processor grid
103      INTEGER , DIMENSION(2)                  ::   topright     ! top right corner of box
104      INTEGER                                 ::   jpiwgt       ! width of box on input grid
105      INTEGER                                 ::   jpjwgt       ! height of box on input grid
106      INTEGER                                 ::   numwgt       ! number of weights (4=bilinear, 16=bicubic)
107      INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in
108      INTEGER                                 ::   overlap      ! =0 when cyclic grid has no overlapping EW columns
109      !                                                         ! =>1 when they have one or more overlapping columns     
110      !                                                         ! =-1 not cyclic
111      LOGICAL                                 ::   cyclic       ! east-west cyclic or not
112      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpi     ! array of source integers
113      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpj     ! array of source integers
114      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid
115      REAL(wp), DIMENSION(:,:,:), POINTER     ::   fly_dta      ! array of values on input grid
116      REAL(wp), DIMENSION(:,:,:), POINTER     ::   col          ! temporary array for reading in columns
117   END TYPE WGT
118
119   INTEGER,     PARAMETER             ::   tot_wgts = 20
120   TYPE( WGT ), DIMENSION(tot_wgts)   ::   ref_wgts     ! array of wgts
121   INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array
122   INTEGER                            ::   nflag = 0
123   REAL(wp), PARAMETER                ::   undeff_lsm = -999.00_wp
124
125!$AGRIF_END_DO_NOT_TREAT
126
127   !! * Substitutions
128#  include "do_loop_substitute.h90"
129   !!----------------------------------------------------------------------
130   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
131   !! $Id$
132   !! Software governed by the CeCILL license (see ./LICENSE)
133   !!----------------------------------------------------------------------
134CONTAINS
135
136   SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, pt_offset, Kmm )
137      !!---------------------------------------------------------------------
138      !!                    ***  ROUTINE fld_read  ***
139      !!                   
140      !! ** Purpose :   provide at each time step the surface ocean fluxes
141      !!                (momentum, heat, freshwater and runoff)
142      !!
143      !! ** Method  :   READ each input fields in NetCDF files using IOM
144      !!      and intepolate it to the model time-step.
145      !!         Several assumptions are made on the input file:
146      !!      blahblahblah....
147      !!----------------------------------------------------------------------
148      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
149      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)
150      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
151      INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option
152      REAL(wp) , INTENT(in   ), OPTIONAL     ::   pt_offset ! provide fields at time other than "now"
153      INTEGER  , INTENT(in   ), OPTIONAL     ::   Kmm       ! ocean time level index
154      !!
155      INTEGER  ::   imf          ! size of the structure sd
156      INTEGER  ::   jf           ! dummy indices
157      INTEGER  ::   isecsbc      ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
158      LOGICAL  ::   ll_firstcall ! true if this is the first call to fld_read for this set of fields
159      REAL(wp) ::   zt_offset    ! local time offset variable
160      REAL(wp) ::   ztinta       ! ratio applied to after  records when doing time interpolation
161      REAL(wp) ::   ztintb       ! ratio applied to before records when doing time interpolation
162      CHARACTER(LEN=1000) ::   clfmt  ! write format
163      !!---------------------------------------------------------------------
164      ll_firstcall = kt == nit000
165      IF( PRESENT(kit) )   ll_firstcall = ll_firstcall .and. kit == 1
166
167      IF( nn_components == jp_iam_sas ) THEN   ;   zt_offset = REAL( nn_fsbc, wp )
168      ELSE                                     ;   zt_offset = 0.
169      ENDIF
170      IF( PRESENT(pt_offset) )   zt_offset = pt_offset
171
172      ! Note that all varibles starting by nsec_* are shifted time by +1/2 time step to be centrered
173      IF( PRESENT(kit) ) THEN   ! ignore kn_fsbc in this case
174         isecsbc = nsec_year + nsec1jan000 + NINT( (     REAL(      kit,wp) + zt_offset ) * rn_Dt / REAL(nn_e,wp) )
175      ELSE                      ! middle of sbc time step
176         ! note: we use kn_fsbc-1 because nsec_year is defined at the middle of the current time step
177         isecsbc = nsec_year + nsec1jan000 + NINT( ( 0.5*REAL(kn_fsbc-1,wp) + zt_offset ) * rn_Dt )
178      ENDIF
179      imf = SIZE( sd )
180      !
181      IF( ll_firstcall ) THEN                      ! initialization
182         DO jf = 1, imf 
183            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE
184            CALL fld_init( isecsbc, sd(jf) )       ! read each before field (put them in after as they will be swapped)
185         END DO
186         IF( lwp ) CALL wgt_print()                ! control print
187      ENDIF
188      !                                            ! ====================================== !
189      IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! update field at each kn_fsbc time-step !
190         !                                         ! ====================================== !
191         !
192         DO jf = 1, imf                            ! ---   loop over field   --- !
193            !
194            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE
195            CALL fld_update( isecsbc, sd(jf), Kmm )
196            !
197         END DO                                    ! --- end loop over field --- !
198
199         CALL fld_rot( kt, sd )                    ! rotate vector before/now/after fields if needed
200
201         DO jf = 1, imf                            ! ---   loop over field   --- !
202            !
203            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE
204            !
205            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
206               IF(lwp .AND. kt - nit000 <= 100 ) THEN
207                  clfmt = "('   fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
208                     &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
209                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &           
210                     & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
211                  WRITE(numout, *) '      zt_offset is : ',zt_offset
212               ENDIF
213               ! temporal interpolation weights
214               ztinta =  REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
215               ztintb =  1. - ztinta
216               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
217            ELSE   ! nothing to do...
218               IF(lwp .AND. kt - nit000 <= 100 ) THEN
219                  clfmt = "('   fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
220                     &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
221                  WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
222                     &                 sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
223               ENDIF
224            ENDIF
225            !
226            IF( kt == nitend - kn_fsbc + 1 )   CALL iom_close( sd(jf)%num )   ! Close the input files
227
228         END DO                                    ! --- end loop over field --- !
229         !
230      ENDIF
231      !
232   END SUBROUTINE fld_read
233
234
235   SUBROUTINE fld_init( ksecsbc, sdjf )
236      !!---------------------------------------------------------------------
237      !!                    ***  ROUTINE fld_init  ***
238      !!
239      !! ** Purpose :  - first call(s) to fld_def to define before values
240      !!               - open file
241      !!----------------------------------------------------------------------
242      INTEGER  , INTENT(in   ) ::   ksecsbc   !
243      TYPE(FLD), INTENT(inout) ::   sdjf         ! input field related variables
244      !!---------------------------------------------------------------------
245      !
246      IF( nflag == 0 )   nflag = -( HUGE(0) - 10 )
247      !
248      CALL fld_def( sdjf )
249      IF( sdjf%ln_tint .AND. ksecsbc < sdjf%nrecsec(1) )   CALL fld_def( sdjf, ldprev = .TRUE. )
250      !
251      CALL fld_clopn( sdjf )
252      sdjf%nrec_a(:) = (/ 1, nflag /)  ! default definition to force flp_update to read the file.
253      !
254   END SUBROUTINE fld_init
255
256
257   SUBROUTINE fld_update( ksecsbc, sdjf, Kmm )
258      !!---------------------------------------------------------------------
259      !!                    ***  ROUTINE fld_update  ***
260      !!
261      !! ** Purpose : Compute
262      !!              if sdjf%ln_tint = .TRUE.
263      !!                  nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping)
264      !!              if sdjf%ln_tint = .FALSE.
265      !!                  nrec_a(1): record number
266      !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record
267      !!----------------------------------------------------------------------
268      INTEGER  ,           INTENT(in   ) ::   ksecsbc   !
269      TYPE(FLD),           INTENT(inout) ::   sdjf      ! input field related variables
270      INTEGER  , OPTIONAL, INTENT(in   ) ::   Kmm    ! ocean time level index
271      !
272      INTEGER  ::   ja     ! end of this record (in seconds)
273      !!----------------------------------------------------------------------
274      !
275      IF( ksecsbc > sdjf%nrec_a(2) ) THEN     ! --> we need to update after data
276       
277         ! find where is the new after record... (it is not necessary sdjf%nrec_a(1)+1 )
278         ja = sdjf%nrec_a(1)
279         DO WHILE ( ksecsbc >= sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast )   ! Warning: make sure ja <= sdjf%nreclast in this test
280            ja = ja + 1
281         END DO
282         IF( ksecsbc > sdjf%nrecsec(ja) )   ja = ja + 1   ! in case ksecsbc > sdjf%nrecsec(sdjf%nreclast)
283
284         ! if ln_tint and if the new after is not ja+1, we need also to update after data before the swap
285         ! so, after the swap, sdjf%nrec_b(2) will still be the closest value located just before ksecsbc
286         IF( sdjf%ln_tint .AND. ( ja > sdjf%nrec_a(1) + 1 .OR. sdjf%nrec_a(2) == nflag ) ) THEN
287            sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /)   ! update nrec_a with before information
288            CALL fld_get( sdjf, Kmm )                         ! read after data that will be used as before data
289         ENDIF
290           
291         ! if after is in the next file...
292         IF( ja > sdjf%nreclast ) THEN
293           
294            CALL fld_def( sdjf )
295            IF( ksecsbc > sdjf%nrecsec(sdjf%nreclast) )   CALL fld_def( sdjf, ldnext = .TRUE. )
296            CALL fld_clopn( sdjf )           ! open next file
297           
298            ! find where is after in this new file
299            ja = 1
300            DO WHILE ( ksecsbc > sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast )
301               ja = ja + 1
302            END DO
303            IF( ksecsbc > sdjf%nrecsec(ja) )   ja = ja + 1   ! in case ksecsbc > sdjf%nrecsec(sdjf%nreclast)
304           
305            IF( ja > sdjf%nreclast ) THEN
306               CALL ctl_stop( "STOP", "fld_def: need next-next file? we should not be there... file: "//TRIM(sdjf%clrootname) )
307            ENDIF
308           
309            ! if ln_tint and if after is not the first record, we must (potentially again) update after data before the swap
310            IF( sdjf%ln_tint .AND. ja > 1 ) THEN
311               IF( sdjf%nrecsec(0) /= nflag ) THEN                  ! no trick used: after file is not the current file
312                  sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /)   ! update nrec_a with before information
313                  CALL fld_get( sdjf, Kmm )                         ! read after data that will be used as before data
314               ENDIF
315            ENDIF
316           
317         ENDIF
318
319         IF( sdjf%ln_tint ) THEN 
320            ! Swap data
321            sdjf%nrec_b(:)     = sdjf%nrec_a(:)                     ! swap before record informations
322            sdjf%rotn(1)       = sdjf%rotn(2)                       ! swap before rotate informations
323            sdjf%fdta(:,:,:,1) = sdjf%fdta(:,:,:,2)                 ! swap before record field
324         ELSE
325            sdjf%nrec_b(:) = (/ ja-1, sdjf%nrecsec(ja-1) /)         ! only for print
326         ENDIF
327           
328         ! read new after data
329         sdjf%nrec_a(:) = (/ ja, sdjf%nrecsec(ja) /)                ! update nrec_a as it is used by fld_get
330         CALL fld_get( sdjf, Kmm )                                  ! read after data (with nrec_a informations)
331       
332      ENDIF
333      !
334   END SUBROUTINE fld_update
335
336
337   SUBROUTINE fld_get( sdjf, Kmm )
338      !!---------------------------------------------------------------------
339      !!                    ***  ROUTINE fld_get  ***
340      !!
341      !! ** Purpose :   read the data
342      !!----------------------------------------------------------------------
343      TYPE(FLD),           INTENT(inout) ::   sdjf   ! input field related variables
344      INTEGER  , OPTIONAL, INTENT(in   ) ::   Kmm    ! ocean time level index
345      !
346      INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
347      INTEGER ::   iw       ! index into wgts array
348      INTEGER ::   ipdom    ! index of the domain
349      INTEGER ::   idvar    ! variable ID
350      INTEGER ::   idmspc   ! number of spatial dimensions
351      LOGICAL ::   lmoor    ! C1D case: point data
352      !!---------------------------------------------------------------------
353      !
354      ipk = SIZE( sdjf%fnow, 3 )
355      !
356      IF( ASSOCIATED(sdjf%imap) ) THEN
357         IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1),   &
358            &                                        sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm )
359         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1),   &
360            &                                        sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm )
361         ENDIF
362      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
363         CALL wgt_list( sdjf, iw )
364         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2),          & 
365            &                                                                          sdjf%nrec_a(1), sdjf%lsmname )
366         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,:  ),          &
367            &                                                                          sdjf%nrec_a(1), sdjf%lsmname )
368         ENDIF
369      ELSE
370         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN   ;   ipdom = jpdom_data
371         ELSE                                   ;   ipdom = jpdom_unknown
372         ENDIF
373         ! C1D case: If product of spatial dimensions == ipk, then x,y are of
374         ! size 1 (point/mooring data): this must be read onto the central grid point
375         idvar  = iom_varid( sdjf%num, sdjf%clvar )
376         idmspc = iom_file ( sdjf%num )%ndims( idvar )
377         IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1
378         lmoor  = (  idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk  )
379         !
380         SELECT CASE( ipk )
381         CASE(1)
382            IF( lk_c1d .AND. lmoor ) THEN
383               IF( sdjf%ln_tint ) THEN
384                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) )
385                  CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,1,2),'Z',1. )
386               ELSE
387                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1  ), sdjf%nrec_a(1) )
388                  CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,1  ),'Z',1. )
389               ENDIF
390            ELSE
391               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
392               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1  ), sdjf%nrec_a(1) )
393               ENDIF
394            ENDIF
395         CASE DEFAULT
396            IF(lk_c1d .AND. lmoor ) THEN
397               IF( sdjf%ln_tint ) THEN
398                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
399                  CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,:,2),'Z',1. )
400               ELSE
401                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,:  ), sdjf%nrec_a(1) )
402                  CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,:  ),'Z',1. )
403               ENDIF
404            ELSE
405               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
406               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) )
407               ENDIF
408            ENDIF
409         END SELECT
410      ENDIF
411      !
412      sdjf%rotn(2) = .false.   ! vector not yet rotated
413      !
414   END SUBROUTINE fld_get
415
416   
417   SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel, ldzint, Kmm )
418      !!---------------------------------------------------------------------
419      !!                    ***  ROUTINE fld_map  ***
420      !!
421      !! ** Purpose :   read global data from file and map onto local data
422      !!                using a general mapping (for open boundaries)
423      !!----------------------------------------------------------------------
424      INTEGER                   , INTENT(in   ) ::   knum         ! stream number
425      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdvar        ! variable name
426      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta         ! bdy output field on model grid
427      INTEGER                   , INTENT(in   ) ::   krec         ! record number to read (ie time slice)
428      INTEGER , DIMENSION(:)    , INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices
429      ! optional variables used for vertical interpolation:
430      INTEGER, OPTIONAL         , INTENT(in   ) ::   kgrd         ! grid type (t, u, v)
431      INTEGER, OPTIONAL         , INTENT(in   ) ::   kbdy         ! bdy number
432      LOGICAL, OPTIONAL         , INTENT(in   ) ::   ldtotvel     ! true if total ( = barotrop + barocline) velocity
433      LOGICAL, OPTIONAL         , INTENT(in   ) ::   ldzint       ! true if 3D variable requires a vertical interpolation
434      INTEGER, OPTIONAL         , INTENT(in   ) ::   Kmm          ! ocean time level index
435      !!
436      INTEGER                                   ::   ipi          ! length of boundary data on local process
437      INTEGER                                   ::   ipj          ! length of dummy dimension ( = 1 )
438      INTEGER                                   ::   ipk          ! number of vertical levels of pdta ( 2D: ipk=1 ; 3D: ipk=jpk )
439      INTEGER                                   ::   ipkb         ! number of vertical levels in boundary data file
440      INTEGER                                   ::   idvar        ! variable ID
441      INTEGER                                   ::   indims       ! number of dimensions of the variable
442      INTEGER, DIMENSION(4)                     ::   idimsz       ! size of variable dimensions
443      REAL(wp)                                  ::   zfv          ! fillvalue
444      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zz_read      ! work space for global boundary data
445      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read    ! work space local data requiring vertical interpolation
446      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_z  ! work space local data requiring vertical interpolation
447      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_dz ! work space local data requiring vertical interpolation
448      CHARACTER(LEN=1),DIMENSION(3)             ::   clgrid
449      LOGICAL                                   ::   lluld        ! is the variable using the unlimited dimension
450      LOGICAL                                   ::   llzint       ! local value of ldzint
451      !!---------------------------------------------------------------------
452      !
453      clgrid = (/'t','u','v'/)
454      !
455      ipi = SIZE( pdta, 1 )
456      ipj = SIZE( pdta, 2 )   ! must be equal to 1
457      ipk = SIZE( pdta, 3 )
458      !
459      llzint = .FALSE.
460      IF( PRESENT(ldzint) )   llzint = ldzint
461      !
462      idvar = iom_varid( knum, cdvar, kndims = indims, kdimsz = idimsz, lduld = lluld  )
463      IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipkb = idimsz(3)   ! xy(zl)t or xy(zl)
464      ELSE                                                            ;   ipkb = 1           ! xy or xyt
465      ENDIF
466      !
467      ALLOCATE( zz_read( idimsz(1), idimsz(2), ipkb ) )  ! ++++++++ !!! this can be very big...         
468      !
469      IF( ipk == 1 ) THEN
470
471         IF( ipkb /= 1 ) CALL ctl_stop( 'fld_map : we must have ipkb = 1 to read surface data' )
472         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,1), krec )   ! call iom_get with a 2D file
473         CALL fld_map_core( zz_read, kmap, pdta )
474
475      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
476      ! Do we include something here to adjust barotropic velocities !
477      ! in case of a depth difference between bdy files and          !
478      ! bathymetry in the case ln_totvel = .false. and ipkb>0?       !
479      ! [as the enveloping and parital cells could change H]         !
480      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481
482      ELSE
483         !
484         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,:), krec )   ! call iom_get with a 3D file
485         !
486         IF( ipkb /= ipk .OR. llzint ) THEN   ! boundary data not on model vertical grid : vertical interpolation
487            !
488            IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN
489               
490               ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) )
491               
492               CALL fld_map_core( zz_read, kmap, zdta_read )
493               CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution?
494               CALL fld_map_core( zz_read, kmap, zdta_read_z )
495               CALL iom_get ( knum, jpdom_unknown,   'e3'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution?
496               CALL fld_map_core( zz_read, kmap, zdta_read_dz )
497               
498               CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar )
499               CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel, Kmm)
500               DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz )
501               
502            ELSE
503               IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' )
504               WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires ' 
505               IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' )
506               IF( iom_varid(knum,   'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//  'e3'//clgrid(kgrd)//' variable' )
507
508            ENDIF
509            !
510         ELSE                            ! bdy data assumed to be the same levels as bdy variables
511            !
512            CALL fld_map_core( zz_read, kmap, pdta )
513            !
514         ENDIF   ! ipkb /= ipk
515      ENDIF   ! ipk == 1
516     
517      DEALLOCATE( zz_read )
518
519   END SUBROUTINE fld_map
520
521     
522   SUBROUTINE fld_map_core( pdta_read, kmap, pdta_bdy )
523      !!---------------------------------------------------------------------
524      !!                    ***  ROUTINE fld_map_core  ***
525      !!
526      !! ** Purpose :  inner core of fld_map
527      !!----------------------------------------------------------------------
528      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read    ! global boundary data
529      INTEGER,  DIMENSION(:    ), INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices
530      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta_bdy     ! bdy output field on model grid
531      !!
532      INTEGER,  DIMENSION(3) ::   idim_read,  idim_bdy            ! arrays dimensions
533      INTEGER                ::   ji, jj, jk, jb                  ! loop counters
534      INTEGER                ::   im1
535      !!---------------------------------------------------------------------
536      !
537      idim_read = SHAPE( pdta_read )
538      idim_bdy  = SHAPE( pdta_bdy  )
539      !
540      ! in all cases: idim_bdy(2) == 1 .AND. idim_read(1) * idim_read(2) == idim_bdy(1)
541      ! structured BDY with rimwidth > 1                     : idim_read(2) == rimwidth /= 1
542      ! structured BDY with rimwidth == 1 or unstructured BDY: idim_read(2) == 1
543      !
544      IF( idim_read(2) > 1 ) THEN    ! structured BDY with rimwidth > 1 
545         DO jk = 1, idim_bdy(3)
546            DO jb = 1, idim_bdy(1)
547               im1 = kmap(jb) - 1
548               jj = im1 / idim_read(1) + 1
549               ji = MOD( im1, idim_read(1) ) + 1
550               pdta_bdy(jb,1,jk) =  pdta_read(ji,jj,jk)
551            END DO
552         END DO
553      ELSE
554         DO jk = 1, idim_bdy(3)
555            DO jb = 1, idim_bdy(1)   ! horizontal remap of bdy data on the local bdy
556               pdta_bdy(jb,1,jk) = pdta_read(kmap(jb),1,jk)
557            END DO
558         END DO
559      ENDIF
560     
561   END SUBROUTINE fld_map_core
562   
563   SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel, Kmm )
564      !!---------------------------------------------------------------------
565      !!                    ***  ROUTINE fld_bdy_interp  ***
566      !!
567      !! ** Purpose :   on the fly vertical interpolation to allow the use of
568      !!                boundary data from non-native vertical grid
569      !!----------------------------------------------------------------------
570      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation
571
572      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file
573      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_z     ! depth of the data read in bdy file
574      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_dz    ! thickness of the levels in bdy file
575      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional)
576      REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file
577      LOGICAL                   , INTENT(in   ) ::   ldtotvel        ! true if toal ( = barotrop + barocline) velocity
578      INTEGER                   , INTENT(in   ) ::   kgrd            ! grid type (t, u, v)
579      INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number
580      INTEGER, OPTIONAL         , INTENT(in   ) ::   Kmm             ! ocean time level index
581      !!
582      INTEGER                  ::   ipi                 ! length of boundary data on local process
583      INTEGER                  ::   ipkb                ! number of vertical levels in boundary data file
584      INTEGER                  ::   ipkmax              ! number of vertical levels in boundary data file where no mask
585      INTEGER                  ::   jb, ji, jj, jk, jkb ! loop counters
586      REAL(wp)                 ::   zcoef, zi           !
587      REAL(wp)                 ::   ztrans, ztrans_new  ! transports
588      REAL(wp), DIMENSION(jpk) ::   zdepth, zdhalf      ! level and half-level depth
589      !!---------------------------------------------------------------------
590     
591      ipi  = SIZE( pdta, 1 )
592      ipkb = SIZE( pdta_read, 3 )
593     
594      DO jb = 1, ipi
595         ji = idx_bdy(kbdy)%nbi(jb,kgrd)
596         jj = idx_bdy(kbdy)%nbj(jb,kgrd)
597         !
598         ! --- max jk where input data /= FillValue --- !
599         ipkmax = 1
600         DO jkb = 2, ipkb
601            IF( pdta_read(jb,1,jkb) /= pfv )   ipkmax = MAX( ipkmax, jkb )
602         END DO
603         !
604         ! --- calculate depth at t,u,v points --- !
605         SELECT CASE( kgrd )                         
606         CASE(1)            ! depth of T points:
607            zdepth(:) = gdept(ji,jj,:,Kmm)
608         CASE(2)            ! depth of U points: we must not use gdept_n as we don't want to do a communication
609            !                 --> copy what is done for gdept_n in domvvl...
610            zdhalf(1) = 0.0_wp
611            zdepth(1) = 0.5_wp * e3uw(ji,jj,1,Kmm)
612            DO jk = 2, jpk                               ! vertical sum
613               !    zcoef = umask - wumask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
614               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
615               !                              ! 0.5 where jk = mikt     
616               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
617               zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
618               zdhalf(jk) = zdhalf(jk-1) + e3u(ji,jj,jk-1,Kmm)
619               zdepth(jk) =          zcoef  * ( zdhalf(jk  ) + 0.5_wp * e3uw(ji,jj,jk,Kmm))  &
620                  &         + (1._wp-zcoef) * ( zdepth(jk-1) +          e3uw(ji,jj,jk,Kmm))
621            END DO
622         CASE(3)            ! depth of V points: we must not use gdept_n as we don't want to do a communication
623            !                 --> copy what is done for gdept_n in domvvl...
624            zdhalf(1) = 0.0_wp
625            zdepth(1) = 0.5_wp * e3vw(ji,jj,1,Kmm)
626            DO jk = 2, jpk                               ! vertical sum
627               !    zcoef = vmask - wvmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
628               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
629               !                              ! 0.5 where jk = mikt     
630               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
631               zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
632               zdhalf(jk) = zdhalf(jk-1) + e3v(ji,jj,jk-1,Kmm)
633               zdepth(jk) =          zcoef  * ( zdhalf(jk  ) + 0.5_wp * e3vw(ji,jj,jk,Kmm))  &
634                  &         + (1._wp-zcoef) * ( zdepth(jk-1) +          e3vw(ji,jj,jk,Kmm))
635            END DO
636         END SELECT
637         !         
638         ! --- interpolate bdy data on the model grid --- !
639         DO jk = 1, jpk
640            IF(     zdepth(jk) <= pdta_read_z(jb,1,1)      ) THEN   ! above the first level of external data
641               pdta(jb,1,jk) = pdta_read(jb,1,1)
642            ELSEIF( zdepth(jk) >  pdta_read_z(jb,1,ipkmax) ) THEN   ! below the last level of external data /= FillValue
643               pdta(jb,1,jk) = pdta_read(jb,1,ipkmax)
644            ELSE                                                    ! inbetween: vertical interpolation between jkb & jkb+1
645               DO jkb = 1, ipkmax-1
646                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) ) <= 0._wp ) THEN ! linear interpolation between 2 levels
647                     zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
648                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) )
649                  ENDIF
650               END DO
651            ENDIF
652         END DO   ! jpk
653         !
654      END DO   ! ipi
655
656      ! --- mask data and adjust transport --- !
657      SELECT CASE( kgrd )                         
658
659      CASE(1)                                 ! mask data (probably unecessary)
660         DO jb = 1, ipi
661            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
662            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
663            DO jk = 1, jpk                     
664               pdta(jb,1,jk) = pdta(jb,1,jk) * tmask(ji,jj,jk)
665            END DO
666         END DO
667         
668      CASE(2)                                 ! adjust the U-transport term
669         DO jb = 1, ipi
670            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
671            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
672            ztrans = 0._wp
673            DO jkb = 1, ipkb                              ! calculate transport on input grid
674               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
675            ENDDO
676            ztrans_new = 0._wp
677            DO jk = 1, jpk                                ! calculate transport on model grid
678               ztrans_new = ztrans_new +      pdta(jb,1,jk ) * e3u(ji,jj,jk,Kmm ) * umask(ji,jj,jk)
679            ENDDO
680            DO jk = 1, jpk                                ! make transport correction
681               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
682                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu(ji,jj,Kmm) ) * umask(ji,jj,jk)
683               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
684                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu(ji,jj,Kmm)   * umask(ji,jj,jk)
685               ENDIF
686            ENDDO
687         ENDDO
688
689      CASE(3)                                 ! adjust the V-transport term
690         DO jb = 1, ipi
691            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
692            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
693            ztrans = 0._wp
694            DO jkb = 1, ipkb                              ! calculate transport on input grid
695               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
696            ENDDO
697            ztrans_new = 0._wp
698            DO jk = 1, jpk                                ! calculate transport on model grid
699               ztrans_new = ztrans_new +      pdta(jb,1,jk ) * e3v(ji,jj,jk,Kmm ) * vmask(ji,jj,jk)
700            ENDDO
701            DO jk = 1, jpk                                ! make transport correction
702               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
703                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv(ji,jj,Kmm) ) * vmask(ji,jj,jk)
704               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
705                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv(ji,jj,Kmm)   * vmask(ji,jj,jk)
706               ENDIF
707            ENDDO
708         ENDDO
709      END SELECT
710     
711   END SUBROUTINE fld_bdy_interp
712
713
714   SUBROUTINE fld_rot( kt, sd )
715      !!---------------------------------------------------------------------
716      !!                    ***  ROUTINE fld_rot  ***
717      !!
718      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
719      !!----------------------------------------------------------------------
720      INTEGER                , INTENT(in   ) ::   kt   ! ocean time step
721      TYPE(FLD), DIMENSION(:), INTENT(inout) ::   sd   ! input field related variables
722      !
723      INTEGER ::   ju, jv, jk, jn  ! loop indices
724      INTEGER ::   imf             ! size of the structure sd
725      INTEGER ::   ill             ! character length
726      INTEGER ::   iv              ! indice of V component
727      CHARACTER (LEN=100)          ::   clcomp       ! dummy weight name
728      REAL(wp), DIMENSION(jpi,jpj) ::   utmp, vtmp   ! temporary arrays for vector rotation
729      !!---------------------------------------------------------------------
730      !
731      !! (sga: following code should be modified so that pairs arent searched for each time
732      !
733      imf = SIZE( sd )
734      DO ju = 1, imf
735         IF( TRIM(sd(ju)%clrootname) == 'NOT USED' )   CYCLE
736         ill = LEN_TRIM( sd(ju)%vcomp )
737         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
738            IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN   ! find vector rotations required             
739               IF( sd(ju)%vcomp(1:1) == 'U' ) THEN      ! east-west component has symbolic name starting with 'U'
740                  ! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
741                  clcomp = 'V' // sd(ju)%vcomp(2:ill)   ! works even if ill == 1
742                  iv = -1
743                  DO jv = 1, imf
744                     IF( TRIM(sd(jv)%clrootname) == 'NOT USED' )   CYCLE
745                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
746                  END DO
747                  IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together
748                     DO jk = 1, SIZE( sd(ju)%fnow, 3 )
749                        IF( sd(ju)%ln_tint )THEN
750                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
751                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
752                           sd(ju)%fdta(:,:,jk,jn) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
753                        ELSE
754                           CALL rot_rep( sd(ju)%fnow(:,:,jk   ), sd(iv)%fnow(:,:,jk   ), 'T', 'en->i', utmp(:,:) )
755                           CALL rot_rep( sd(ju)%fnow(:,:,jk   ), sd(iv)%fnow(:,:,jk   ), 'T', 'en->j', vtmp(:,:) )
756                           sd(ju)%fnow(:,:,jk   ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk   ) = vtmp(:,:)
757                        ENDIF
758                     END DO
759                     sd(ju)%rotn(jn) = .TRUE.               ! vector was rotated
760                     IF( lwp .AND. kt == nit000 )   WRITE(numout,*)   &
761                        &   'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
762                  ENDIF
763               ENDIF
764            ENDIF
765         END DO
766       END DO
767      !
768   END SUBROUTINE fld_rot
769
770
771   SUBROUTINE fld_def( sdjf, ldprev, ldnext )
772      !!---------------------------------------------------------------------
773      !!                    ***  ROUTINE fld_def  ***
774      !!
775      !! ** Purpose :   define the record(s) of the file and its name
776      !!----------------------------------------------------------------------
777      TYPE(FLD)        , INTENT(inout) ::   sdjf       ! input field related variables
778      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldprev     !
779      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldnext     !
780      !
781      INTEGER  :: jt
782      INTEGER  :: idaysec               ! number of seconds in 1 day = NINT(rday)
783      INTEGER  :: iyr, imt, idy, isecwk
784      INTEGER  :: indexyr, indexmt
785      INTEGER  :: ireclast
786      INTEGER  :: ishift, istart
787      INTEGER, DIMENSION(2)  :: isave
788      REAL(wp) :: zfreqs
789      LOGICAL  :: llprev, llnext, llstop
790      LOGICAL  :: llprevmt, llprevyr
791      LOGICAL  :: llnextmt, llnextyr
792      !!----------------------------------------------------------------------
793      idaysec = NINT(rday)
794      !
795      IF( PRESENT(ldprev) ) THEN   ;   llprev = ldprev
796      ELSE                         ;   llprev = .FALSE.
797      ENDIF
798      IF( PRESENT(ldnext) ) THEN   ;   llnext = ldnext
799      ELSE                         ;   llnext = .FALSE.
800      ENDIF
801
802      ! current file parameters
803      IF( sdjf%cltype(1:4) == 'week' ) THEN          ! find the day of the beginning of the current week
804         isecwk = ksec_week( sdjf%cltype(6:8) )     ! seconds between the beginning of the week and half of current time step
805         llprevmt = isecwk > nsec_month               ! longer time since beginning of the current week than the current month
806         llprevyr = llprevmt .AND. nmonth == 1
807         iyr = nyear  - COUNT((/llprevyr/))
808         imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/))
809         idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec
810         isecwk = nsec_year - isecwk              ! seconds between 00h jan 1st of current year and current week beginning
811      ELSE
812         iyr = nyear
813         imt = nmonth
814         idy = nday
815         isecwk  = 0
816      ENDIF
817
818      ! previous file parameters
819      IF( llprev ) THEN
820         IF( sdjf%cltype(1:4) == 'week'    ) THEN     ! find the day of the beginning of previous week
821            isecwk = isecwk + 7 * idaysec         ! seconds between the beginning of previous week and half of the time step
822            llprevmt = isecwk > nsec_month            ! longer time since beginning of the previous week than the current month
823            llprevyr = llprevmt .AND. nmonth == 1
824            iyr = nyear  - COUNT((/llprevyr/))
825            imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/))
826            idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec
827            isecwk = nsec_year - isecwk           ! seconds between 00h jan 1st of current year and previous week beginning
828         ELSE
829            idy = nday   - COUNT((/ sdjf%cltype == 'daily'                 /))
830            imt = nmonth - COUNT((/ sdjf%cltype == 'monthly' .OR. idy == 0 /))
831            iyr = nyear  - COUNT((/ sdjf%cltype == 'yearly'  .OR. imt == 0 /))
832            IF( idy == 0 ) idy = nmonth_len(imt)
833            IF( imt == 0 ) imt = 12
834            isecwk = 0
835         ENDIF
836      ENDIF
837
838      ! next file parameters
839      IF( llnext ) THEN
840         IF( sdjf%cltype(1:4) == 'week'    ) THEN     ! find the day of the beginning of next week
841            isecwk = 7 * idaysec - isecwk         ! seconds between half of the time step and the beginning of next week
842            llnextmt = isecwk > ( nmonth_len(nmonth)*idaysec - nsec_month )   ! larger than the seconds to the end of the month
843            llnextyr = llnextmt .AND. nmonth == 12
844            iyr = nyear  + COUNT((/llnextyr/))
845            imt = nmonth + COUNT((/llnextmt/)) - 12 * COUNT((/llnextyr/))
846            idy = nday - nmonth_len(nmonth) * COUNT((/llnextmt/)) + isecwk / idaysec + 1
847            isecwk = nsec_year + isecwk           ! seconds between 00h jan 1st of current year and next week beginning
848         ELSE
849            idy = nday   + COUNT((/ sdjf%cltype == 'daily'                                 /))
850            imt = nmonth + COUNT((/ sdjf%cltype == 'monthly' .OR. idy > nmonth_len(nmonth) /))
851            iyr = nyear  + COUNT((/ sdjf%cltype == 'yearly'  .OR. imt == 13                /))
852            IF( idy > nmonth_len(nmonth) )   idy = 1
853            IF( imt == 13                )   imt = 1
854            isecwk = 0
855         ENDIF
856      ENDIF
857      !
858      ! find the last record to be read -> update sdjf%nreclast
859      indexyr = iyr - nyear + 1                 ! which  year are we looking for? previous(0), current(1) or next(2)?
860      indexmt = imt + 12 * ( indexyr - 1 )      ! which month are we looking for (relatively to current year)?
861      !
862      ! Last record to be read in the current file
863      ! Predefine the number of record in the file according of its type.
864      ! We could compare this number with the number of records in the file and make a stop if the 2 numbers do not match...
865      ! However this would be much less fexible (e.g. for tests) and will force to rewite input files according to nleapy...
866      IF    ( NINT(sdjf%freqh) == -12 ) THEN            ;   ireclast = 1    ! yearly mean: consider only 1 record
867      ELSEIF( NINT(sdjf%freqh) ==  -1 ) THEN                                ! monthly mean:
868         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ireclast = 1    !  consider that the file has  1 record
869         ELSE                                           ;   ireclast = 12   !  consider that the file has 12 record
870         ENDIF
871      ELSE                                                                  ! higher frequency mean (in hours)
872         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ireclast = NINT( 24. * REAL(nmonth_len(indexmt), wp) / sdjf%freqh )
873         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ireclast = NINT( 24. * 7.                            / sdjf%freqh )
874         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ireclast = NINT( 24.                                 / sdjf%freqh )
875         ELSE                                           ;   ireclast = NINT( 24. * REAL( nyear_len(indexyr), wp) / sdjf%freqh )
876         ENDIF
877      ENDIF
878
879      sdjf%nreclast = ireclast
880      ! Allocate arrays for beginning/middle/end of each record (seconds since Jan. 1st 00h of nit000 year)
881      IF( ALLOCATED(sdjf%nrecsec) )   DEALLOCATE( sdjf%nrecsec )
882      ALLOCATE( sdjf%nrecsec( 0:ireclast ) )
883      !
884      IF    ( NINT(sdjf%freqh) == -12 ) THEN                                     ! yearly mean and yearly file
885         SELECT CASE( indexyr )
886         CASE(0)   ;   sdjf%nrecsec(0) = nsec1jan000 - nyear_len( 0 ) * idaysec
887         CASE(1)   ;   sdjf%nrecsec(0) = nsec1jan000
888         CASE(2)   ;   sdjf%nrecsec(0) = nsec1jan000 + nyear_len( 1 ) * idaysec
889         ENDSELECT
890         sdjf%nrecsec(1) = sdjf%nrecsec(0) + nyear_len( indexyr ) * idaysec
891      ELSEIF( NINT(sdjf%freqh) ==  -1 ) THEN                                     ! monthly mean:
892         IF(     sdjf%cltype      == 'monthly' ) THEN                            !    monthly file
893            sdjf%nrecsec(0   ) = nsec1jan000 + nmonth_beg(indexmt  )
894            sdjf%nrecsec(1   ) = nsec1jan000 + nmonth_beg(indexmt+1)
895         ELSE                                                                    !    yearly  file
896            ishift = 12 * ( indexyr - 1 )
897            sdjf%nrecsec(0:12) = nsec1jan000 + nmonth_beg(1+ishift:13+ishift)
898         ENDIF
899      ELSE                                                                       ! higher frequency mean (in hours)
900         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   istart = nsec1jan000 + nmonth_beg(indexmt)
901         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   istart = nsec1jan000 + isecwk
902         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   istart = nsec1jan000 + nmonth_beg(indexmt) + ( idy - 1 ) * idaysec
903         ELSEIF( indexyr          == 0         ) THEN   ;   istart = nsec1jan000 - nyear_len( 0 ) * idaysec
904         ELSEIF( indexyr          == 2         ) THEN   ;   istart = nsec1jan000 + nyear_len( 1 ) * idaysec
905         ELSE                                           ;   istart = nsec1jan000
906         ENDIF
907         zfreqs = sdjf%freqh * rhhmm * rmmss
908         DO jt = 0, sdjf%nreclast
909            sdjf%nrecsec(jt) = istart + NINT( zfreqs * REAL(jt,wp) )
910         END DO
911      ENDIF
912      !
913      IF( sdjf%ln_tint ) THEN   ! record time defined in the middle of the record, computed using an implementation
914                                ! of the rounded average that is valid over the full integer range
915         sdjf%nrecsec(1:sdjf%nreclast) = sdjf%nrecsec(0:sdjf%nreclast-1) / 2 + sdjf%nrecsec(1:sdjf%nreclast) / 2 + &
916            & MAX( MOD( sdjf%nrecsec(0:sdjf%nreclast-1), 2 ), MOD( sdjf%nrecsec(1:sdjf%nreclast), 2 ) )
917      END IF
918      !
919      sdjf%clname = fld_filename( sdjf, idy, imt, iyr )
920      !
921   END SUBROUTINE fld_def
922
923   
924   SUBROUTINE fld_clopn( sdjf )
925      !!---------------------------------------------------------------------
926      !!                    ***  ROUTINE fld_clopn  ***
927      !!
928      !! ** Purpose :   close/open the files
929      !!----------------------------------------------------------------------
930      TYPE(FLD)        , INTENT(inout) ::   sdjf       ! input field related variables
931      !
932      INTEGER, DIMENSION(2)  :: isave
933      LOGICAL  :: llprev, llnext, llstop
934      !!----------------------------------------------------------------------
935      !
936      llprev = sdjf%nrecsec(sdjf%nreclast) < nsec000_1jan000   ! file ends before the beginning of the job -> file may not exist
937      llnext = sdjf%nrecsec(       0     ) > nsecend_1jan000   ! file begins after the end of the job -> file may not exist
938
939      llstop = sdjf%ln_clim .OR. .NOT. ( llprev .OR. llnext )
940
941      IF( sdjf%num <= 0 .OR. .NOT. sdjf%ln_clim  ) THEN
942         IF( sdjf%num > 0 )   CALL iom_close( sdjf%num )   ! close file if already open
943         CALL iom_open( sdjf%clname, sdjf%num, ldstop = llstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
944      ENDIF
945      !
946      IF( sdjf%num <= 0 .AND. .NOT. llstop ) THEN   ! file not found but we do accept this...
947         !
948         IF( llprev ) THEN   ! previous file does not exist : go back to current and accept to read only the first record
949            CALL ctl_warn('previous file: '//TRIM(sdjf%clname)//' not found -> go back to current year/month/week/day file')
950            isave(1:2) = sdjf%nrecsec(sdjf%nreclast-1:sdjf%nreclast)   ! save previous file info
951            CALL fld_def( sdjf )   ! go back to current file
952            sdjf%nreclast = 1   ! force to use only the first record (do as if other were not existing...)
953            sdjf%nrecsec(0:1) = isave(1:2)
954         ENDIF
955         !
956         IF( llnext ) THEN   ! next     file does not exist : go back to current and accept to read only the last  record
957            CALL ctl_warn('next file: '//TRIM(sdjf%clname)//' not found -> go back to current year/month/week/day file')
958            isave(1:2) = sdjf%nrecsec(0:1)    ! save next file info
959            CALL fld_def( sdjf )   ! go back to current file
960            ! -> read last record but keep record info from the first record of next file
961            sdjf%nrecsec(sdjf%nreclast-1:sdjf%nreclast) = isave(1:2)
962            sdjf%nrecsec(0:sdjf%nreclast-2) = nflag
963         ENDIF
964         !
965         CALL iom_open( sdjf%clname, sdjf%num, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )   
966         !
967      ENDIF
968      !
969   END SUBROUTINE fld_clopn
970
971
972   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint )
973      !!---------------------------------------------------------------------
974      !!                    ***  ROUTINE fld_fill  ***
975      !!
976      !! ** Purpose :   fill the data structure (sdf) with the associated information
977      !!              read in namelist (sdf_n) and control print
978      !!----------------------------------------------------------------------
979      TYPE(FLD)  , DIMENSION(:)          , INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
980      TYPE(FLD_N), DIMENSION(:)          , INTENT(in   ) ::   sdf_n      ! array of namelist information structures
981      CHARACTER(len=*)                   , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
982      CHARACTER(len=*)                   , INTENT(in   ) ::   cdcaller   ! name of the calling routine
983      CHARACTER(len=*)                   , INTENT(in   ) ::   cdtitle    ! description of the calling routine
984      CHARACTER(len=*)                   , INTENT(in   ) ::   cdnam      ! name of the namelist from which sdf_n comes
985      INTEGER                  , OPTIONAL, INTENT(in   ) ::   knoprint   ! no calling routine information printed
986      !
987      INTEGER  ::   jf   ! dummy indices
988      !!---------------------------------------------------------------------
989      !
990      DO jf = 1, SIZE(sdf)
991         sdf(jf)%clrootname = sdf_n(jf)%clname
992         IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' )   sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname
993         sdf(jf)%clname     = "not yet defined"
994         sdf(jf)%freqh      = sdf_n(jf)%freqh
995         sdf(jf)%clvar      = sdf_n(jf)%clvar
996         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
997         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
998         sdf(jf)%cltype     = sdf_n(jf)%cltype
999         sdf(jf)%num        = -1
1000         sdf(jf)%wgtname    = " "
1001         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname
1002         sdf(jf)%lsmname = " "
1003         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname
1004         sdf(jf)%vcomp      = sdf_n(jf)%vcomp
1005         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
1006         IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   &
1007            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
1008         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
1009            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
1010         sdf(jf)%nreclast   = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
1011         sdf(jf)%igrd       = 0
1012         sdf(jf)%ibdy       = 0
1013         sdf(jf)%imap       => NULL()
1014         sdf(jf)%ltotvel    = .FALSE.
1015         sdf(jf)%lzint      = .FALSE.
1016      END DO
1017      !
1018      IF(lwp) THEN      ! control print
1019         WRITE(numout,*)
1020         IF( .NOT.PRESENT( knoprint) ) THEN
1021            WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
1022            WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
1023         ENDIF
1024         WRITE(numout,*) '   fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
1025         WRITE(numout,*) '   ~~~~~~~~'
1026         WRITE(numout,*) '      list of files and frequency (>0: in hours ; <0 in months)'
1027         DO jf = 1, SIZE(sdf)
1028            WRITE(numout,*) '      root filename: '  , TRIM( sdf(jf)%clrootname ), '   variable name: ', TRIM( sdf(jf)%clvar )
1029            WRITE(numout,*) '         frequency: '      ,       sdf(jf)%freqh       ,   &
1030               &                  '   time interp: '    ,       sdf(jf)%ln_tint     ,   &
1031               &                  '   climatology: '    ,       sdf(jf)%ln_clim
1032            WRITE(numout,*) '         weights: '        , TRIM( sdf(jf)%wgtname    ),   &
1033               &                  '   pairing: '        , TRIM( sdf(jf)%vcomp      ),   &
1034               &                  '   data type: '      ,       sdf(jf)%cltype      ,   &
1035               &                  '   land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
1036            call flush(numout)
1037         END DO
1038      ENDIF
1039      !
1040   END SUBROUTINE fld_fill
1041
1042
1043   SUBROUTINE wgt_list( sd, kwgt )
1044      !!---------------------------------------------------------------------
1045      !!                    ***  ROUTINE wgt_list  ***
1046      !!
1047      !! ** Purpose :   search array of WGTs and find a weights file entry,
1048      !!                or return a new one adding it to the end if new entry.
1049      !!                the weights data is read in and restructured (fld_weight)
1050      !!----------------------------------------------------------------------
1051      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
1052      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
1053      !
1054      INTEGER ::   kw, nestid   ! local integer
1055      LOGICAL ::   found        ! local logical
1056      !!----------------------------------------------------------------------
1057      !
1058      !! search down linked list
1059      !! weights filename is either present or we hit the end of the list
1060      found = .FALSE.
1061      !
1062      !! because agrif nest part of filenames are now added in iom_open
1063      !! to distinguish between weights files on the different grids, need to track
1064      !! nest number explicitly
1065      nestid = 0
1066#if defined key_agrif
1067      nestid = Agrif_Fixed()
1068#endif
1069      DO kw = 1, nxt_wgt-1
1070         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
1071             ref_wgts(kw)%nestid == nestid) THEN
1072            kwgt = kw
1073            found = .TRUE.
1074            EXIT
1075         ENDIF
1076      END DO
1077      IF( .NOT.found ) THEN
1078         kwgt = nxt_wgt
1079         CALL fld_weight( sd )
1080      ENDIF
1081      !
1082   END SUBROUTINE wgt_list
1083
1084
1085   SUBROUTINE wgt_print( )
1086      !!---------------------------------------------------------------------
1087      !!                    ***  ROUTINE wgt_print  ***
1088      !!
1089      !! ** Purpose :   print the list of known weights
1090      !!----------------------------------------------------------------------
1091      INTEGER ::   kw   !
1092      !!----------------------------------------------------------------------
1093      !
1094      DO kw = 1, nxt_wgt-1
1095         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
1096         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1097         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
1098         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
1099         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
1100         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
1101         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
1102         IF( ref_wgts(kw)%cyclic ) THEN
1103            WRITE(numout,*) '       cyclical'
1104            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap
1105         ELSE
1106            WRITE(numout,*) '       not cyclical'
1107         ENDIF
1108         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
1109      END DO
1110      !
1111   END SUBROUTINE wgt_print
1112
1113
1114   SUBROUTINE fld_weight( sd )
1115      !!---------------------------------------------------------------------
1116      !!                    ***  ROUTINE fld_weight  ***
1117      !!
1118      !! ** Purpose :   create a new WGT structure and fill in data from file,
1119      !!              restructuring as required
1120      !!----------------------------------------------------------------------
1121      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
1122      !!
1123      INTEGER ::   jn         ! dummy loop indices
1124      INTEGER ::   inum       ! local logical unit
1125      INTEGER ::   id         ! local variable id
1126      INTEGER ::   ipk        ! local vertical dimension
1127      INTEGER ::   zwrap      ! local integer
1128      LOGICAL ::   cyclical   !
1129      CHARACTER (len=5) ::   aname   !
1130      INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims
1131      INTEGER,  DIMENSION(jpi,jpj) ::   data_src
1132      REAL(wp), DIMENSION(jpi,jpj) ::   data_tmp
1133      !!----------------------------------------------------------------------
1134      !
1135      IF( nxt_wgt > tot_wgts ) THEN
1136        CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1137      ENDIF
1138      !
1139      !! new weights file entry, add in extra information
1140      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
1141      !! input data file is representative of all other files to be opened and processed with the
1142      !! current weights file
1143
1144      !! open input data file (non-model grid)
1145      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
1146
1147      !! get dimensions: we consider 2D data as 3D data with vertical dim size = 1
1148      IF( SIZE(sd%fnow, 3) > 0 ) THEN
1149         ALLOCATE( ddims(4) )
1150      ELSE
1151         ALLOCATE( ddims(3) )
1152      ENDIF
1153      id = iom_varid( inum, sd%clvar, ddims )
1154
1155      !! close it
1156      CALL iom_close( inum )
1157
1158      !! now open the weights file
1159
1160      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
1161      IF( inum > 0 ) THEN
1162
1163         !! determine whether we have an east-west cyclic grid
1164         !! from global attribute called "ew_wrap" in the weights file
1165         !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1166         !! since this is the most common forcing configuration
1167
1168         CALL iom_getatt(inum, 'ew_wrap', zwrap)
1169         IF( zwrap >= 0 ) THEN
1170            cyclical = .TRUE.
1171         ELSE IF( zwrap == -999 ) THEN
1172            cyclical = .TRUE.
1173            zwrap = 0
1174         ELSE
1175            cyclical = .FALSE.
1176         ENDIF
1177
1178         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1179         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1180         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1181         ref_wgts(nxt_wgt)%overlap = zwrap
1182         ref_wgts(nxt_wgt)%cyclic = cyclical
1183         ref_wgts(nxt_wgt)%nestid = 0
1184#if defined key_agrif
1185         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1186#endif
1187         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1188         !! for each weight wgtNN there is an integer array srcNN which gives the point in
1189         !! the input data grid which is to be multiplied by the weight
1190         !! they are both arrays on the model grid so the result of the multiplication is
1191         !! added into an output array on the model grid as a running sum
1192
1193         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1194         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1195         IF( id <= 0) THEN
1196            ref_wgts(nxt_wgt)%numwgt = 4
1197         ELSE
1198            ref_wgts(nxt_wgt)%numwgt = 16
1199         ENDIF
1200
1201         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1202         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1203         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1204
1205         DO jn = 1,4
1206            aname = ' '
1207            WRITE(aname,'(a3,i2.2)') 'src',jn
1208            data_tmp(:,:) = 0
1209            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1210            data_src(:,:) = INT(data_tmp(:,:))
1211            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1212            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1213         END DO
1214
1215         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1216            aname = ' '
1217            WRITE(aname,'(a3,i2.2)') 'wgt',jn
1218            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1219            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1220         END DO
1221         CALL iom_close (inum)
1222 
1223         ! find min and max indices in grid
1224         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1225         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1226         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1227         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1228
1229         ! and therefore dimensions of the input box
1230         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1231         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1232
1233         ! shift indexing of source grid
1234         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1235         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1236
1237         ! create input grid, give it a halo to allow gradient calculations
1238         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1239         ! a more robust solution will be given in next release
1240         ipk =  SIZE(sd%fnow, 3)
1241         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1242         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1243         !
1244         nxt_wgt = nxt_wgt + 1
1245         !
1246      ELSE
1247         CALL ctl_stop( '    fld_weight : unable to read the file ' )
1248      ENDIF
1249
1250      DEALLOCATE (ddims )
1251      !
1252   END SUBROUTINE fld_weight
1253
1254
1255   SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm,   &
1256      &                          jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
1257      !!---------------------------------------------------------------------
1258      !!                    ***  ROUTINE apply_seaoverland  ***
1259      !!
1260      !! ** Purpose :   avoid spurious fluxes in coastal or near-coastal areas
1261      !!                due to the wrong usage of "land" values from the coarse
1262      !!                atmospheric model when spatial interpolation is required
1263      !!      D. Delrosso INGV         
1264      !!----------------------------------------------------------------------
1265      INTEGER,                   INTENT(in   ) :: itmpi,itmpj,itmpz                    ! lengths
1266      INTEGER,                   INTENT(in   ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
1267      INTEGER, DIMENSION(3),     INTENT(in   ) :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
1268      REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo                              ! input/output array for seaoverland application
1269      CHARACTER (len=100),       INTENT(in   ) :: clmaskfile                           ! land/sea mask file name
1270      !
1271      INTEGER :: inum,jni,jnj,jnz,jc   ! local indices
1272      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1             ! local array for land point detection
1273      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfieldn   ! array of forcing field with undeff for land points
1274      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfield    ! array of forcing field
1275      !!---------------------------------------------------------------------
1276      !
1277      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) )
1278      !
1279      ! Retrieve the land sea mask data
1280      CALL iom_open( clmaskfile, inum )
1281      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1282      CASE(1)
1283         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1284      CASE DEFAULT
1285         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1286      END SELECT
1287      CALL iom_close( inum )
1288      !
1289      DO jnz=1,rec1_lsm(3)             !! Loop over k dimension
1290         !
1291         DO jni = 1, itmpi                               !! copy the original field into a tmp array
1292            DO jnj = 1, itmpj                            !! substituting undeff over land points
1293               zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1294               IF( zslmec1(jni,jnj,jnz) == 1. )   zfieldn(jni,jnj) = undeff_lsm
1295            END DO
1296         END DO
1297         !
1298         CALL seaoverland( zfieldn, itmpi, itmpj, zfield )
1299         DO jc = 1, nn_lsm
1300            CALL seaoverland( zfield, itmpi, itmpj, zfield )
1301         END DO
1302         !
1303         !   Check for Undeff and substitute original values
1304         IF( ANY(zfield==undeff_lsm) ) THEN
1305            DO jni = 1, itmpi
1306               DO jnj = 1, itmpj
1307                  IF( zfield(jni,jnj)==undeff_lsm )   zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1308               END DO
1309            END DO
1310         ENDIF
1311         !
1312         zfieldo(:,:,jnz) = zfield(:,:)
1313         !
1314      END DO                           !! End Loop over k dimension
1315      !
1316      DEALLOCATE ( zslmec1, zfieldn, zfield )
1317      !
1318   END SUBROUTINE apply_seaoverland 
1319
1320
1321   SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield )
1322      !!---------------------------------------------------------------------
1323      !!                    ***  ROUTINE seaoverland  ***
1324      !!
1325      !! ** Purpose :   create shifted matrices for seaoverland application 
1326      !!      D. Delrosso INGV
1327      !!----------------------------------------------------------------------
1328      INTEGER                      , INTENT(in   ) :: ileni,ilenj   ! lengths
1329      REAL, DIMENSION (ileni,ilenj), INTENT(in   ) :: zfieldn       ! array of forcing field with undeff for land points
1330      REAL, DIMENSION (ileni,ilenj), INTENT(  out) :: zfield        ! array of forcing field
1331      !
1332      REAL   , DIMENSION (ileni,ilenj)   :: zmat1, zmat2, zmat3, zmat4  ! local arrays
1333      REAL   , DIMENSION (ileni,ilenj)   :: zmat5, zmat6, zmat7, zmat8  !   -     -
1334      REAL   , DIMENSION (ileni,ilenj)   :: zlsm2d                      !   -     -
1335      REAL   , DIMENSION (ileni,ilenj,8) :: zlsm3d                      !   -     -
1336      LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d                 ! logical mask for undeff detection
1337      LOGICAL, DIMENSION (ileni,ilenj)   :: ll_msknan2d                 ! logical mask for undeff detection
1338      !!----------------------------------------------------------------------
1339      zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/)     , DIM=2 )
1340      zmat1 = eoshift( zmat8   , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/)       , DIM=1 )
1341      zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/)     , DIM=1 )
1342      zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 )
1343      zmat3 = eoshift( zmat4   , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/)       , DIM=1 )
1344      zmat5 = eoshift( zmat4   , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/)   , DIM=1 )
1345      zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 )
1346      zmat7 = eoshift( zmat8   , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/)   , DIM=1 )
1347      !
1348      zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1349      ll_msknan3d = .NOT.( zlsm3d  == undeff_lsm )
1350      ll_msknan2d = .NOT.( zfieldn == undeff_lsm )  ! FALSE where is Undeff (land)
1351      zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) )
1352      WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp )   zlsm2d = undeff_lsm
1353      zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d )
1354      !
1355   END SUBROUTINE seaoverland
1356
1357
1358   SUBROUTINE fld_interp( num, clvar, kw, kk, dta,  &
1359                          &         nrec, lsmfile)     
1360      !!---------------------------------------------------------------------
1361      !!                    ***  ROUTINE fld_interp  ***
1362      !!
1363      !! ** Purpose :   apply weights to input gridded data to create data
1364      !!                on model grid
1365      !!----------------------------------------------------------------------
1366      INTEGER                   , INTENT(in   ) ::   num     ! stream number
1367      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
1368      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
1369      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
1370      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
1371      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
1372      CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
1373      !
1374      INTEGER, DIMENSION(3) ::   rec1, recn           ! temporary arrays for start and length
1375      INTEGER, DIMENSION(3) ::   rec1_lsm, recn_lsm   ! temporary arrays for start and length in case of seaoverland
1376      INTEGER ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2    ! temporary indices
1377      INTEGER ::   jk, jn, jm, jir, jjr               ! loop counters
1378      INTEGER ::   ni, nj                             ! lengths
1379      INTEGER ::   jpimin,jpiwid                      ! temporary indices
1380      INTEGER ::   jpimin_lsm,jpiwid_lsm              ! temporary indices
1381      INTEGER ::   jpjmin,jpjwid                      ! temporary indices
1382      INTEGER ::   jpjmin_lsm,jpjwid_lsm              ! temporary indices
1383      INTEGER ::   jpi1,jpi2,jpj1,jpj2                ! temporary indices
1384      INTEGER ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
1385      INTEGER ::   itmpi,itmpj,itmpz                     ! lengths
1386      REAL(wp),DIMENSION(:,:,:), ALLOCATABLE ::   ztmp_fly_dta                 ! local array of values on input grid     
1387      !!----------------------------------------------------------------------
1388      !
1389      !! for weighted interpolation we have weights at four corners of a box surrounding
1390      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
1391      !! or by a grid value and gradients at the corner point (bicubic case)
1392      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1393
1394      !! sub grid from non-model input grid which encloses all grid points in this nemo process
1395      jpimin = ref_wgts(kw)%botleft(1)
1396      jpjmin = ref_wgts(kw)%botleft(2)
1397      jpiwid = ref_wgts(kw)%jpiwgt
1398      jpjwid = ref_wgts(kw)%jpjwgt
1399
1400      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1401      rec1(1) = MAX( jpimin-1, 1 )
1402      rec1(2) = MAX( jpjmin-1, 1 )
1403      rec1(3) = 1
1404      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1405      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1406      recn(3) = kk
1407
1408      !! where we need to put it in the non-nemo grid fly_dta
1409      !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1410      !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1411      jpi1 = 2 + rec1(1) - jpimin
1412      jpj1 = 2 + rec1(2) - jpjmin
1413      jpi2 = jpi1 + recn(1) - 1
1414      jpj2 = jpj1 + recn(2) - 1
1415
1416
1417      IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1418      !! indeces for ztmp_fly_dta
1419      ! --------------------------
1420         rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1)  ! starting index for enlarged external data, x direction
1421         rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1)  ! starting index for enlarged external data, y direction
1422         rec1_lsm(3) = 1                    ! vertical dimension
1423         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
1424         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
1425         recn_lsm(3) = kk                   ! number of vertical levels in the input file
1426
1427      !  Avoid out of bound
1428         jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1429         jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1430         jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1431         jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1432
1433         jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1434         jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1435         jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1436         jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1437
1438
1439         itmpi=jpi2_lsm-jpi1_lsm+1
1440         itmpj=jpj2_lsm-jpj1_lsm+1
1441         itmpz=kk
1442         ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1443         ztmp_fly_dta(:,:,:) = 0.0
1444         SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1445         CASE(1)
1446              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   &
1447                 &                                                                nrec, rec1_lsm, recn_lsm)
1448         CASE DEFAULT
1449              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   &
1450                 &                                                                nrec, rec1_lsm, recn_lsm)
1451         END SELECT
1452         CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),                  &
1453                 &                                      jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm,                  &
1454                 &                                      itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1455
1456
1457         ! Relative indeces for remapping
1458         ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1459         ii_lsm2 = (ii_lsm1+recn(1))-1
1460         ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1461         ij_lsm2 = (ij_lsm1+recn(2))-1
1462
1463         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1464         ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1465         DEALLOCATE(ztmp_fly_dta)
1466
1467      ELSE
1468         
1469         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1470         CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1471      ENDIF
1472     
1473
1474      !! first four weights common to both bilinear and bicubic
1475      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1476      !! note that we have to offset by 1 into fly_dta array because of halo
1477      dta(:,:,:) = 0.0
1478      DO jk = 1,4
1479        DO jn = 1, jpj
1480          DO jm = 1,jpi
1481            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1482            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1483            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1484          END DO
1485        END DO
1486      END DO
1487
1488      IF(ref_wgts(kw)%numwgt .EQ. 16) THEN
1489
1490        !! fix up halo points that we couldnt read from file
1491        IF( jpi1 == 2 ) THEN
1492           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1493        ENDIF
1494        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1495           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1496        ENDIF
1497        IF( jpj1 == 2 ) THEN
1498           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1499        ENDIF
1500        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1501           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1502        ENDIF
1503
1504        !! if data grid is cyclic we can do better on east-west edges
1505        !! but have to allow for whether first and last columns are coincident
1506        IF( ref_wgts(kw)%cyclic ) THEN
1507           rec1(2) = MAX( jpjmin-1, 1 )
1508           recn(1) = 1
1509           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1510           jpj1 = 2 + rec1(2) - jpjmin
1511           jpj2 = jpj1 + recn(2) - 1
1512           IF( jpi1 == 2 ) THEN
1513              rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1514              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1515              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1516           ENDIF
1517           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1518              rec1(1) = 1 + ref_wgts(kw)%overlap
1519              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1520              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1521           ENDIF
1522        ENDIF
1523
1524        ! gradient in the i direction
1525        DO jk = 1,4
1526          DO jn = 1, jpj
1527            DO jm = 1,jpi
1528              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1529              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1530              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1531                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1532            END DO
1533          END DO
1534        END DO
1535
1536        ! gradient in the j direction
1537        DO jk = 1,4
1538          DO jn = 1, jpj
1539            DO jm = 1,jpi
1540              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1541              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1542              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         &
1543                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1544            END DO
1545          END DO
1546        END DO
1547
1548         ! gradient in the ij direction
1549         DO jk = 1,4
1550            DO jn = 1, jpj
1551               DO jm = 1,jpi
1552                  ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1553                  nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1554                  dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1555                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1556                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
1557               END DO
1558            END DO
1559         END DO
1560         !
1561      ENDIF
1562      !
1563   END SUBROUTINE fld_interp
1564
1565
1566   FUNCTION fld_filename( sdjf, kday, kmonth, kyear )
1567      !!---------------------------------------------------------------------
1568      !!                    ***  FUNCTION fld_filename ***
1569      !!
1570      !! ** Purpose :   define the filename according to a given date
1571      !!---------------------------------------------------------------------
1572      TYPE(FLD), INTENT(in) ::   sdjf         ! input field related variables
1573      INTEGER  , INTENT(in) ::   kday, kmonth, kyear
1574      !
1575      CHARACTER(len = 256) ::   clname, fld_filename
1576      !!---------------------------------------------------------------------
1577
1578     
1579      ! build the new filename if not climatological data
1580      clname=TRIM(sdjf%clrootname)
1581      !
1582      ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
1583      IF( .NOT. sdjf%ln_clim ) THEN   
1584                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear    ! add year
1585         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a, "m",i2.2)' ) TRIM( clname          ), kmonth   ! add month
1586      ELSE
1587         ! build the new filename if climatological data
1588         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth   ! add month
1589      ENDIF
1590      IF(    sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
1591         &                               WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), kday     ! add day
1592
1593      fld_filename = clname
1594     
1595   END FUNCTION fld_filename
1596
1597
1598   FUNCTION ksec_week( cdday )
1599      !!---------------------------------------------------------------------
1600      !!                    ***  FUNCTION ksec_week ***
1601      !!
1602      !! ** Purpose :   seconds between 00h of the beginning of the week and half of the current time step
1603      !!---------------------------------------------------------------------
1604      CHARACTER(len=*), INTENT(in)   ::   cdday   ! first 3 letters of the first day of the weekly file
1605      !!
1606      INTEGER                        ::   ksec_week      ! output variable
1607      INTEGER                        ::   ijul, ishift   ! local integer
1608      CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
1609      !!----------------------------------------------------------------------
1610      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1611      DO ijul = 1, 7
1612         IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1613      END DO
1614      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1615      !
1616      ishift = ijul * NINT(rday)
1617      !
1618      ksec_week = nsec_monday + ishift
1619      ksec_week = MOD( ksec_week, 7*NINT(rday) )
1620      !
1621   END FUNCTION ksec_week
1622
1623   !!======================================================================
1624END MODULE fldread
Note: See TracBrowser for help on using the repository browser.