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

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

dev_r12114_ticket_2263: new version of fldread calendar, see #2263

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