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_r11943_MERGE_2019/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/fldread.F90 @ 12250

Last change on this file since 12250 was 12250, checked in by smasson, 4 years ago

rev12240_dev_r11943_MERGE_2019: same as [12246], add modifications from dev_r12114_ticket_2263, results unchanged except SPITZ12 as explained in #2263

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