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

Last change on this file was 13546, checked in by smasson, 2 months ago

trunk: supress the use of undefined values, see #2535

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