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

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

Last change on this file was 15023, checked in by gsamson, 3 years ago

merge ticket2680_C1D_PAPA branch back into the trunk; see ticket #2680 for details

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