source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/SBC/fldread.F90 @ 13257

Last change on this file since 13257 was 13257, checked in by orioltp, 3 months ago

Updated with trunk at r13245 and small change allocating variables in icb_oce.F90.

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