New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
fldread.F90 in NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/fldread.F90 @ 12351

Last change on this file since 12351 was 12351, checked in by smueller, 4 years ago

Rectification of a defect that halved the maximum length of the interval between restarts in model runs with enabled temporal interpolation of forcing fields (field-information flag ln_tint set to .true.)

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