source: NEMO/branches/2019/dev_r12114_ticket_2263/src/OCE/SBC/fldread.F90 @ 12127

Last change on this file since 12127 was 12127, checked in by smasson, 16 months ago

dev_r12114_ticket_2263: replace integer kt_offset by real pt_offset, see #2263

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