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

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

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

Last change on this file since 12340 was 12340, checked in by acc, 2 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:


Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

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