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 branches/devmercator2010/NEMO/OPA_SRC/SBC – NEMO

source: branches/devmercator2010/NEMO/OPA_SRC/SBC/fldread.F90 @ 2077

Last change on this file since 2077 was 2077, checked in by cbricaud, 13 years ago

commit change from DEV_r1784_3DF

  • Property svn:keywords set to Id
File size: 54.6 KB
Line 
1MODULE fldread
2   !!======================================================================
3   !!                       ***  MODULE  fldread  ***
4   !! Ocean forcing:  read input field for surface boundary condition
5   !!=====================================================================
6   !! History :  9.0  !  06-06  (G. Madec) Original code
7   !!                 !  05-08  (S. Alderson) Modified for Interpolation in memory
8   !!                 !         from input grid to model grid
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   fld_read      : read input fields used for the computation of the
13   !!                   surface boundary condition
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE phycst          ! ???
18   USE in_out_manager  ! I/O manager
19   USE iom             ! I/O manager library
20   USE geo2ocean       ! for vector rotation on to model grid
21
22   IMPLICIT NONE
23   PRIVATE   
24
25   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
26      CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file
27      INTEGER              ::   nfreqh      ! frequency of each flux file
28      CHARACTER(len = 34)  ::   clvar       ! generic name of the variable in the NetCDF flux file
29      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F)
30      LOGICAL              ::   ln_clim     ! climatology or not (T/F)
31      CHARACTER(len = 7)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly'
32      CHARACTER(len = 34)  ::   wname       ! generic name of a NetCDF weights file to be used, blank if not
33      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation
34                                           ! a string starting with "U" or "V" for each component   
35                                           ! chars 2 onwards identify which components go together 
36   END TYPE FLD_N
37
38   TYPE, PUBLIC ::   FLD        !: Input field related variables
39      CHARACTER(len = 256)            ::   clrootname   ! generic name of the NetCDF file
40      CHARACTER(len = 256)            ::   clname       ! current name of the NetCDF file
41      INTEGER                         ::   nfreqh       ! frequency of each flux file
42      CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file
43      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F)
44      LOGICAL                         ::   ln_clim      ! climatology or not (T/F)
45      CHARACTER(len = 7)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly'
46      INTEGER                         ::   num          ! iom id of the jpfld files to be read
47      INTEGER                         ::   nswap_sec    ! swapping time in second since Jan. 1st 00h of nit000 year
48      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
49      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year)
50      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:)   ::   fnow       ! input fields interpolated to now time step
51      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields
52      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
53                                                        ! into the WGTLIST structure
54      CHARACTER(len = 34)             ::   vcomp        ! symbolic name for a vector component that needs rotation
55      LOGICAL ,  DIMENSION(2)         ::   rotn         ! flag to indicate whether field has been rotated
56   END TYPE FLD
57
58!$AGRIF_DO_NOT_TREAT
59
60   !! keep list of all weights variables so they're only read in once
61   !! need to add AGRIF directives not to process this structure
62   !! also need to force wgtname to include AGRIF nest number
63   TYPE         ::   WGT        !: Input weights related variables
64      CHARACTER(len = 256)                    ::   wgtname      ! current name of the NetCDF weight file
65      INTEGER , DIMENSION(2)                  ::   ddims        ! shape of input grid
66      INTEGER , DIMENSION(2)                  ::   botleft      ! top left corner of box in input grid containing
67                                                                ! current processor grid
68      INTEGER , DIMENSION(2)                  ::   topright     ! top right corner of box
69      INTEGER                                 ::   jpiwgt       ! width of box on input grid
70      INTEGER                                 ::   jpjwgt       ! height of box on input grid
71      INTEGER                                 ::   numwgt       ! number of weights (4=bilinear, 16=bicubic)
72      INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in
73      INTEGER                                 ::   offset       ! =0 when cyclic grid has coincident first/last columns,
74                                                                ! =1 when they assumed to be one grid spacing apart     
75                                                                ! =-1 otherwise
76      LOGICAL                                 ::   cyclic       ! east-west cyclic or not
77      INTEGER, DIMENSION(:,:,:), POINTER      ::   data_jpi     ! array of source integers
78      INTEGER, DIMENSION(:,:,:), POINTER      ::   data_jpj     ! array of source integers
79      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid
80      REAL(wp), DIMENSION(:,:,:), POINTER     ::   fly_dta      ! array of values on input grid
81      REAL(wp), DIMENSION(:,:,:), POINTER     ::   col2         ! temporary array for reading in columns
82   END TYPE WGT
83
84   INTEGER,     PARAMETER             ::   tot_wgts = 10
85   TYPE( WGT ), DIMENSION(tot_wgts)   ::   ref_wgts     ! array of wgts
86   INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array
87
88!$AGRIF_END_DO_NOT_TREAT
89
90   PUBLIC   fld_read, fld_fill   ! called by sbc... modules
91
92   !!----------------------------------------------------------------------
93   !!   OPA 9.0 , LOCEAN-IPSL (2006)
94   !! $Id$
95   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
96   !!----------------------------------------------------------------------
97
98CONTAINS
99
100   SUBROUTINE fld_read( kt, kn_fsbc, sd )
101      !!---------------------------------------------------------------------
102      !!                    ***  ROUTINE fld_read  ***
103      !!                   
104      !! ** Purpose :   provide at each time step the surface ocean fluxes
105      !!                (momentum, heat, freshwater and runoff)
106      !!
107      !! ** Method  :   READ each input fields in NetCDF files using IOM
108      !!      and intepolate it to the model time-step.
109      !!         Several assumptions are made on the input file:
110      !!      blahblahblah....
111      !!----------------------------------------------------------------------
112      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
113      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)
114      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
115      !!
116      CHARACTER (LEN=34)                     ::   acomp     ! dummy weight name
117      INTEGER                                ::   kf, nf    ! dummy indices
118      INTEGER                                ::   imf       ! size of the structure sd
119      REAL(wp), DIMENSION(jpi,jpj)           ::   utmp, vtmp! temporary arrays for vector rotation
120
121      INTEGER  ::   jf         ! dummy indices
122      INTEGER  ::   jk         ! dummy indices
123      INTEGER  ::   ipk        ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
124      INTEGER  ::   kw         ! index into wgts array
125      INTEGER  ::   ireclast   ! last record to be read in the current year file
126      INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend
127      LOGICAL  ::   llnxtyr    ! open next year  file?
128      LOGICAL  ::   llnxtmth   ! open next month file?
129      LOGICAL  ::   llstop     ! stop is the file does not exist
130      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
131      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
132      CHARACTER(LEN=1000) ::   clfmt   ! write format
133      !!---------------------------------------------------------------------
134      !
135      imf = SIZE( sd )
136      !                                         ! ===================== !
137      DO jf = 1, imf                            !    LOOP OVER FIELD    !
138         !                                      ! ===================== !
139         !
140         IF( kt == nit000 )   CALL fld_init( sd(jf) )
141         !
142         ! read/update the after data?
143         IF( nsec_year + nsec1jan000 > sd(jf)%nswap_sec ) THEN
144
145            IF( sd(jf)%ln_tint ) THEN         ! time interpolation: swap before record field
146!CDIR COLLAPSE
147               sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)
148               sd(jf)%rotn(1)       = sd(jf)%rotn(2)
149            ENDIF
150
151            ! update record informations
152            CALL fld_rec( sd(jf) )
153
154            ! do we have to change the year/month/day of the forcing field??
155            IF( sd(jf)%ln_tint ) THEN
156               ! if we do time interpolation we will need to open next year/month/day file before the end of the current one
157               ! if so, we are still before the end of the year/month/day when calling fld_rec so sd(jf)%nrec_a(1) will be
158               ! larger than the record number that should be read for current year/month/day (for ex. 13 for monthly mean file)
159
160               ! last record to be read in the current file
161               IF( sd(jf)%nfreqh == -1 ) THEN
162                  IF(     sd(jf)%cltype == 'monthly'   ) THEN  ;   ireclast = 1
163                  ELSE                                         ;   ireclast = 12
164                  ENDIF
165               ELSE                             
166                  IF(     sd(jf)%cltype == 'monthly'   ) THEN  ;   ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh 
167                  ELSEIF( sd(jf)%cltype == 'daily'     ) THEN  ;   ireclast = 24                      / sd(jf)%nfreqh
168                  ELSE                                         ;   ireclast = 24 * nyear_len(     1 ) / sd(jf)%nfreqh 
169                  ENDIF
170               ENDIF
171             
172               ! do we need next file data?
173               IF( sd(jf)%nrec_a(1) > ireclast ) THEN
174
175                  sd(jf)%nrec_a(1) = 1              ! force to read the first record of the next file
176
177                  IF( .NOT. sd(jf)%ln_clim ) THEN   ! close the current file and open a new one.
178
179                     llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth)      ! open next month file?
180                     llnxtyr  = sd(jf)%cltype == 'yearly'  .OR. (nmonth == 12 .AND. llnxtmth)   ! open next year  file?
181
182                     ! if the run finishes at the end of the current year/month/day, we will allow next year/month/day file to be
183                     ! not present. If the run continue further than the current year/month/day, next year/month/day file must exist
184                     isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1))   ! second at the end of the run
185                     llstop = isecend > sd(jf)%nswap_sec                                   ! read more than 1 record of next year
186
187                     CALL fld_clopn( sd(jf), nyear  + COUNT((/llnxtyr /))                                           ,         &
188                        &                    nmonth + COUNT((/llnxtmth/)) - 12                 * COUNT((/llnxtyr /)),         &
189                        &                    nday   + 1                   - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
190
191                     IF( sd(jf)%num == 0 .AND. .NOT. llstop ) THEN    ! next year file does not exist
192                        CALL ctl_warn('next year/month/day file: '//TRIM(sd(jf)%clname)//     &
193                                &     ' not present -> back to current year/month/day')
194                        CALL fld_clopn( sd(jf), nyear, nmonth, nday )       ! back to the current year/month/day
195                        sd(jf)%nrec_a(1) = ireclast     ! force to read the last record to be read in the current year file
196                     ENDIF
197
198                  ENDIF
199               ENDIF
200       
201            ELSE
202               ! if we are not doing time interpolation, we must change the year/month/day of the file just after switching
203               ! to the NEW year/month/day. If it is the case, we are at the beginning of the year/month/day when calling
204               ! fld_rec so sd(jf)%nrec_a(1) = 1
205               IF( sd(jf)%nrec_a(1) == 1 .AND. .NOT. sd(jf)%ln_clim )   CALL fld_clopn( sd(jf), nyear, nmonth, nday )
206            ENDIF
207
208            ! read after data
209            IF( LEN(TRIM(sd(jf)%wgtname)) > 0 ) THEN
210               CALL wgt_list( sd(jf), kw )
211               ipk =  SIZE(sd(jf)%fdta,3)
212               CALL fld_interp( sd(jf)%num, sd(jf)%clvar , kw , ipk, sd(jf)%fdta(:,:,:,2) , sd(jf)%nrec_a(1) )
213            ELSE
214               SELECT CASE( SIZE(sd(jf)%fdta,3) )
215               CASE(1)
216                  CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,1,2), sd(jf)%nrec_a(1) )
217               CASE(jpk)
218                  CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,:,2), sd(jf)%nrec_a(1) )
219               END SELECT
220            ENDIF
221            sd(jf)%rotn(2) = .FALSE.
222
223         ENDIF
224         !                                      ! ===================== !
225      END DO                                    !  END LOOP OVER FIELD  !
226      !                                         ! ===================== !
227
228      IF( kt == nit000 .AND. lwp ) CALL wgt_print()
229
230      !! Vector fields may need to be rotated onto the local grid direction
231      !! This has to happen before the time interpolations
232      !! (sga: following code should be modified so that pairs arent searched for each time
233
234      DO jf = 1, imf
235         !! find vector rotations required
236         IF( LEN(TRIM(sd(jf)%vcomp)) > 0 ) THEN
237             !! east-west component has symbolic name starting with 'U'
238             IF( sd(jf)%vcomp(1:1) == 'U' ) THEN
239                !! found an east-west component, look for the north-south component
240                !! which has same symbolic name but with 'U' replaced with 'V'
241                nf = LEN_TRIM( sd(jf)%vcomp )
242                IF( nf == 1) THEN
243                   acomp = 'V'
244                ELSE
245                   acomp = 'V' // sd(jf)%vcomp(2:nf)
246                ENDIF
247                kf = -1
248                DO nf = 1, imf
249                  IF( TRIM(sd(nf)%vcomp) == TRIM(acomp) ) kf = nf
250                END DO
251                IF( kf > 0 ) THEN
252                   !! fields jf,kf are two components which need to be rotated together
253                   DO nf = 1,2
254                      !! check each time level of this pair
255                      IF( .NOT. sd(jf)%rotn(nf) .AND. .NOT. sd(kf)%rotn(nf) ) THEN
256                         utmp(:,:) = 0.0
257                         vtmp(:,:) = 0.0
258                         !
259                         ipk = SIZE( sd(kf)%fdta(:,:,:,nf) ,3 )
260                         DO jk = 1,ipk
261                            CALL rot_rep( sd(jf)%fdta(:,:,jk,nf),sd(kf)%fdta(:,:,jk,nf),'T', 'en->i', utmp(:,:) )
262                            CALL rot_rep( sd(jf)%fdta(:,:,jk,nf),sd(kf)%fdta(:,:,jk,nf),'T', 'en->j', vtmp(:,:) )
263                            sd(jf)%fdta(:,:,jk,nf) = utmp(:,:)
264                            sd(kf)%fdta(:,:,jk,nf) = vtmp(:,:)
265                         END DO
266                         !
267                         sd(jf)%rotn(nf) = .TRUE.
268                         sd(kf)%rotn(nf) = .TRUE.
269                         IF( lwp .AND. kt == nit000 ) &
270                                   WRITE(numout,*) 'fld_read: vector pair (',  &
271                                                   TRIM(sd(jf)%clvar),',',TRIM(sd(kf)%clvar), &
272                                                   ') rotated on to model grid'
273                      ENDIF
274                   END DO
275                ENDIF
276             ENDIF
277         ENDIF
278      END DO
279
280      !                                         ! ===================== !
281      DO jf = 1, imf                            !    LOOP OVER FIELD    !
282         !                                      ! ===================== !
283         !
284         ! update field at each kn_fsbc time-step
285         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN   
286            !
287            IF( sd(jf)%ln_tint ) THEN
288               IF(lwp .AND. kt - nit000 <= 100 ) THEN
289                  clfmt = "('fld_read: var ', a, ' kt = ', i8,' Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
290                     &    "' records b/a: ', i4.4, '/', i4.4, ' (', f7.2,'/', f7.2, ' days)')"
291                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, nyear, nmonth, nday,   &
292                     & 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
293               ENDIF
294               !
295               ztinta =  REAL( nsec_year + nsec1jan000 - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
296               ztintb =  1. - ztinta
297!CDIR COLLAPSE
298               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
299            ELSE
300               IF(lwp .AND. kt - nit000 <= 100 ) THEN
301                  clfmt = "('fld_read: var ', a, ' kt = ', i8,' Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
302                     &    "' record: ', i4.4, ' at ', f7.2, ' day')"
303                  WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, nyear, nmonth, nday, sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_a(2),wp)/rday
304               ENDIF
305!CDIR COLLAPSE
306               sd(jf)%fnow(:,:,:) = sd(jf)%fdta(:,:,:,2)   ! piecewise constant field
307 
308            ENDIF
309            !
310         ENDIF
311
312         IF( kt == nitend )   CALL iom_close( sd(jf)%num )   ! Close the input files
313
314         !                                      ! ===================== !
315      END DO                                    !  END LOOP OVER FIELD  !
316      !                                         ! ===================== !
317   END SUBROUTINE fld_read
318
319
320   SUBROUTINE fld_init( sdjf )
321      !!---------------------------------------------------------------------
322      !!                    ***  ROUTINE fld_init  ***
323      !!
324      !! ** Purpose :  - if time interpolation, read before data
325      !!               - open current year file
326      !!
327      !! ** Method  :   
328      !!----------------------------------------------------------------------
329      TYPE(FLD), INTENT(inout) ::   sdjf        ! input field related variables
330      !!
331      LOGICAL :: llprevyr       ! are we reading previous year  file?
332      LOGICAL :: llprevmth      ! are we reading previous month file?
333      LOGICAL :: llprevday      ! are we reading previous day   file?
334      LOGICAL :: llprev         ! llprevyr .OR. llprevmth .OR. llprevday
335      INTEGER :: idvar          ! variable id
336      INTEGER :: inrec          ! number of record existing for this variable
337      INTEGER :: kwgt
338      INTEGER :: jk             !vertical loop variable
339      INTEGER :: ipk            !number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
340      CHARACTER(LEN=1000) ::   clfmt   ! write format
341      !!---------------------------------------------------------------------
342
343      ! some default definitions...
344      sdjf%num = 0   ! default definition for non-opened file
345      IF( sdjf%ln_clim )   sdjf%clname = TRIM( sdjf%clrootname )   ! file name defaut definition, never change in this case
346      llprevyr  = .FALSE.
347      llprevmth = .FALSE.
348      llprevday = .FALSE.
349           
350      ! define record informations
351      CALL fld_rec( sdjf )
352
353      IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
354         
355         IF( sdjf%nrec_b(1) == 0  ) THEN   ! we redefine record sdjf%nrec_b(1) with the last record of previous year file
356            IF( sdjf%nfreqh == -1 ) THEN   ! monthly mean
357               IF( sdjf%cltype == 'monthly' ) THEN   ! monthly file
358                  sdjf%nrec_b(1) = 1                                                       ! force to read the unique record
359                  llprevmth = .TRUE.                                                       ! use previous month file?
360                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
361               ELSE                                  ! yearly file
362                  sdjf%nrec_b(1) = 12                                                      ! force to read december mean
363                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
364               ENDIF
365            ELSE   
366               IF( sdjf%cltype == 'monthly' ) THEN   ! monthly file
367                  sdjf%nrec_b(1) = 24 * nmonth_len(nmonth-1) / sdjf%nfreqh                 ! last record of previous month
368                  llprevmth = .NOT. sdjf%ln_clim                                           ! use previous month file?
369                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
370               ELSEIF( sdjf%cltype == 'daily' ) THEN ! daily file
371                  sdjf%nrec_b(1) = 24 / sdjf%nfreqh                                        ! last record of previous day
372                  llprevday = .NOT. sdjf%ln_clim                                           ! use previous day   file?
373                  llprevmth = llprevday .AND. nday   == 1                                  ! use previous month file?
374                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
375               ELSE                                  ! yearly file
376                  sdjf%nrec_b(1) = 24 * nyear_len(0) / sdjf%nfreqh                         ! last record of year month
377                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
378               ENDIF
379            ENDIF
380         ENDIF
381         llprev = llprevyr .OR. llprevmth .OR. llprevday
382
383         CALL fld_clopn( sdjf, nyear  - COUNT((/llprevyr /))                                              ,               &
384            &                  nmonth - COUNT((/llprevmth/)) + 12                   * COUNT((/llprevyr /)),               &
385            &                  nday   - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)), .NOT. llprev )
386
387         ! if previous year/month/day file does not exist, we switch to the current year/month/day
388         IF( llprev .AND. sdjf%num == 0 ) THEN
389            CALL ctl_warn( 'previous year/month/day file: '//TRIM(sdjf%clname)//' not present -> back to current year/month/day')
390            ! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
391            llprev = .false.
392            sdjf%nrec_b(1) = 1
393            CALL fld_clopn( sdjf, nyear, nmonth, nday )
394         ENDIF
395         
396         IF( llprev ) THEN   ! check if the last record sdjf%nrec_n(1) exists in the file
397            idvar = iom_varid( sdjf%num, sdjf%clvar )                                        ! id of the variable sdjf%clvar
398            IF( idvar <= 0 )   RETURN
399            inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar )   ! size of the last dim of idvar
400            sdjf%nrec_b(1) = MIN( sdjf%nrec_b(1), inrec )   ! make sure we select an existing record
401         ENDIF
402
403         ! read before data into sdjf%fdta(:,:,2) because we will swap data in the following part of fld_read
404         
405         IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
406            CALL wgt_list( sdjf, kwgt )
407            ipk = SIZE(sdjf%fdta,3)
408            CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, ipk, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
409         ELSE
410            SELECT CASE ( SIZE(sdjf%fdta,3) )
411            CASE(1)
412                CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_b(1) )
413            CASE(jpk)
414                CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_b(1) )
415            END SELECT
416         ENDIF
417         sdjf%rotn(2) = .FALSE.
418
419         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i4, ' at time = ', f7.2, ' days')"
420         IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_b(1), REAL(sdjf%nrec_b(2),wp)/rday
421
422         IF( llprev )   CALL iom_close( sdjf%num )   ! close previous year file (-> redefine sdjf%num to 0)
423
424      ENDIF
425
426
427      IF( sdjf%num == 0 )   CALL fld_clopn( sdjf, nyear, nmonth, nday )   ! make sure current year/month/day file is opened
428
429      sdjf%nswap_sec = nsec_year + nsec1jan000 - 1   ! force read/update the after data in the following part of fld_read
430     
431   END SUBROUTINE fld_init
432
433
434   SUBROUTINE fld_rec( sdjf )
435      !!---------------------------------------------------------------------
436      !!                    ***  ROUTINE fld_rec  ***
437      !!
438      !! ** Purpose :   compute nrec_a, nrec_b and nswap_sec
439      !!
440      !! ** Method  :   
441      !!----------------------------------------------------------------------
442      TYPE(FLD), INTENT(inout) ::   sdjf        ! input field related variables
443      !!
444      INTEGER  ::   irec        ! record number
445      INTEGER  ::   isecd       ! rday
446      REAL(wp) ::   ztmp        ! temporary variable
447      INTEGER  ::   ifreq_sec   ! frequency mean (in seconds)
448      !!----------------------------------------------------------------------
449      !
450      IF( sdjf%nfreqh == -1 ) THEN      ! monthly mean
451         !
452         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
453            !
454            !                  INT( ztmp )
455            !                     /|\
456            !                    1 |    *----
457            !                    0 |----(             
458            !                      |----+----|--> time
459            !                      0   /|\   1   (nday/nmonth_len(nmonth))
460            !                           |   
461            !                           |   
462            !       forcing record :  nmonth
463            !                           
464            ztmp  = 0.e0
465            IF(  REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) .GT. 0.5 ) ztmp  = 1.0
466         ELSE
467            ztmp  = 0.e0
468         ENDIF
469         irec = nmonth + INT( ztmp )
470
471         IF( sdjf%ln_tint ) THEN   ;   sdjf%nswap_sec = nmonth_half(irec) + nsec1jan000   ! swap at the middle of the month
472         ELSE                      ;   sdjf%nswap_sec = nmonth_end (irec) + nsec1jan000   ! swap at the end    of the month
473         ENDIF
474
475         IF( sdjf%cltype == 'monthly' ) THEN
476
477            sdjf%nrec_b(:) = (/ 0, nmonth_half(irec - 1 ) + nsec1jan000 /)
478            sdjf%nrec_a(:) = (/ 1, nmonth_half(irec     ) + nsec1jan000 /)
479
480            IF( ztmp  == 1. ) THEN
481              sdjf%nrec_b(1) = 1
482              sdjf%nrec_a(1) = 2
483            ENDIF
484
485         ELSE
486
487            sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /)   ! define after  record number and time
488            irec = irec - 1                                                ! move back to previous record
489            sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /)   ! define before record number and time
490
491         ENDIF
492         !
493      ELSE                              ! higher frequency mean (in hours)
494         !
495         ifreq_sec = sdjf%nfreqh * 3600   ! frequency mean (in seconds)
496         ! number of second since the beginning of the file
497         IF(     sdjf%cltype == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month,wp)   ! since 00h on the 1st day of the current month
498         ELSEIF( sdjf%cltype == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day  ,wp)   ! since 00h of the current day
499         ELSE                                      ;   ztmp = REAL(nsec_year ,wp)   ! since 00h on Jan 1 of the current year
500         ENDIF
501         IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record
502            !
503            !                  INT( ztmp )
504            !                     /|\
505            !                    2 |        *-----(
506            !                    1 |  *-----(
507            !                    0 |--(             
508            !                      |--+--|--+--|--+--|--> time
509            !                      0 /|\ 1 /|\ 2 /|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
510            !                         |     |     |
511            !                         |     |     |
512            !       forcing record :  1     2     3
513            !                   
514            ztmp= ztmp / ifreq_sec + 0.5
515         ELSE                 
516            !
517            !                  INT( ztmp )
518            !                     /|\
519            !                    2 |           *-----(
520            !                    1 |     *-----(
521            !                    0 |-----(             
522            !                      |--+--|--+--|--+--|--> time
523            !                      0 /|\ 1 /|\ 2 /|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
524            !                         |     |     |
525            !                         |     |     |
526            !       forcing record :  1     2     3
527            !                           
528            ztmp= ztmp / ifreq_sec
529         ENDIF
530         irec = 1 + INT( ztmp )
531
532         isecd = NINT(rday)
533         ! after record index and second since Jan. 1st 00h of nit000 year
534         sdjf%nrec_a(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /)
535         IF( sdjf%cltype == 'monthly' )   &   ! add the number of seconds between 00h Jan 1 and the end of previous month
536            sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * SUM(nmonth_len(1:nmonth -1))   ! ok if nmonth=1
537         IF( sdjf%cltype == 'daily'   )   &   ! add the number of seconds between 00h Jan 1 and the end of previous day
538            sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * ( nday_year - 1 )
539
540         ! before record index and second since Jan. 1st 00h of nit000 year
541         irec = irec - 1.                           ! move back to previous record
542         sdjf%nrec_b(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /)
543         IF( sdjf%cltype == 'monthly' )   &   ! add the number of seconds between 00h Jan 1 and the end of previous month
544            sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * SUM(nmonth_len(1:nmonth -1))   ! ok if nmonth=1
545         IF( sdjf%cltype == 'daily'   )   &   ! add the number of seconds between 00h Jan 1 and the end of previous day
546            sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * ( nday_year - 1 )
547
548         ! swapping time in second since Jan. 1st 00h of nit000 year
549         IF( sdjf%ln_tint ) THEN   ;   sdjf%nswap_sec =  sdjf%nrec_a(2)                     ! swap at the middle of the record
550         ELSE                      ;   sdjf%nswap_sec =  sdjf%nrec_a(2) + ifreq_sec / 2     ! swap at the end    of the record
551         ENDIF       
552         !
553      ENDIF
554      !
555   END SUBROUTINE fld_rec
556
557
558   SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
559      !!---------------------------------------------------------------------
560      !!                    ***  ROUTINE fld_clopn  ***
561      !!
562      !! ** Purpose :   update the file name and open the file
563      !!
564      !! ** Method  :   
565      !!----------------------------------------------------------------------
566      TYPE(FLD), INTENT(inout)           ::   sdjf     ! input field related variables
567      INTEGER  , INTENT(in   )           ::   kyear    ! year value
568      INTEGER  , INTENT(in   )           ::   kmonth   ! month value
569      INTEGER  , INTENT(in   )           ::   kday     ! day value
570      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
571
572      IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
573      ! build the new filename if not climatological data
574      IF( .NOT. sdjf%ln_clim ) THEN   ;   WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear    ! add year
575         IF( sdjf%cltype /= 'yearly' )    WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname     ), kmonth   ! add month
576         IF( sdjf%cltype == 'daily'  )    WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname     ), kday     ! add day
577      ELSE
578         ! build the new filename if climatological data
579         IF( sdjf%cltype == 'monthly' )   WRITE(sdjf%clname, '(a,"_m" ,i2.2)' ) TRIM( sdjf%clrootname ), kmonth   ! add month
580      ENDIF
581      CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
582      !
583   END SUBROUTINE fld_clopn
584
585
586   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
587      !!---------------------------------------------------------------------
588      !!                    ***  ROUTINE fld_fill  ***
589      !!
590      !! ** Purpose :   fill sdf with sdf_n and control print
591      !!
592      !! ** Method  :   
593      !!----------------------------------------------------------------------
594      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
595      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
596      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
597      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   !
598      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    !
599      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      !
600      !
601      INTEGER  ::   jf       ! dummy indices
602      !!---------------------------------------------------------------------
603
604      DO jf = 1, SIZE(sdf)
605         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
606         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
607         sdf(jf)%clvar      = sdf_n(jf)%clvar
608         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
609         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
610         sdf(jf)%cltype     = sdf_n(jf)%cltype
611         sdf(jf)%wgtname = " "
612         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
613         sdf(jf)%vcomp   = sdf_n(jf)%vcomp
614      END DO
615
616      IF(lwp) THEN      ! control print
617         WRITE(numout,*)
618         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
619         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
620         WRITE(numout,*) '          '//TRIM( cdnam )//' Namelist'
621         WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
622         DO jf = 1, SIZE(sdf)
623            WRITE(numout,*) '               root filename: '  , TRIM( sdf(jf)%clrootname ),   &
624               &                          ' variable name: '  , TRIM( sdf(jf)%clvar      )
625            WRITE(numout,*) '               frequency: '      ,       sdf(jf)%nfreqh      ,   &
626               &                          ' time interp: '    ,       sdf(jf)%ln_tint     ,   &
627               &                          ' climatology: '    ,       sdf(jf)%ln_clim     ,   &
628               &                          ' weights    : '    , TRIM( sdf(jf)%wgtname    ),   &
629               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   &
630               &                          ' data type: '      ,       sdf(jf)%cltype
631         END DO
632      ENDIF
633     
634   END SUBROUTINE fld_fill
635
636
637   SUBROUTINE wgt_list( sd, kwgt )
638      !!---------------------------------------------------------------------
639      !!                    ***  ROUTINE wgt_list  ***
640      !!
641      !! ** Purpose :   search array of WGTs and find a weights file
642      !!                entry, or return a new one adding it to the end
643      !!                if it is a new entry, the weights data is read in and
644      !!                restructured (fld_weight)
645      !!
646      !! ** Method  :   
647      !!----------------------------------------------------------------------
648      TYPE( FLD ),      INTENT(in)    ::   sd        ! field with name of weights file
649      INTEGER,      INTENT(inout)     ::   kwgt      ! index of weights
650      !!
651      INTEGER                         ::   kw
652      INTEGER                         ::   nestid
653      LOGICAL                         ::   found
654      !!----------------------------------------------------------------------
655      !
656      !! search down linked list
657      !! weights filename is either present or we hit the end of the list
658      found = .FALSE.
659
660      !! because agrif nest part of filenames are now added in iom_open
661      !! to distinguish between weights files on the different grids, need to track
662      !! nest number explicitly
663      nestid = 0
664#if defined key_agrif
665      nestid = Agrif_Fixed()
666#endif
667      DO kw = 1, nxt_wgt-1
668         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
669             ref_wgts(kw)%nestid == nestid) THEN
670            kwgt = kw
671            found = .TRUE.
672            EXIT
673         ENDIF
674      END DO
675      IF( .NOT.found ) THEN
676         kwgt = nxt_wgt
677         CALL fld_weight( sd )
678      ENDIF
679
680   END SUBROUTINE wgt_list
681
682   SUBROUTINE wgt_print( )
683      !!---------------------------------------------------------------------
684      !!                    ***  ROUTINE wgt_print  ***
685      !!
686      !! ** Purpose :   print the list of known weights
687      !!
688      !! ** Method  :   
689      !!----------------------------------------------------------------------
690      !!
691      INTEGER                         ::   kw
692      !!----------------------------------------------------------------------
693      !
694
695      DO kw = 1, nxt_wgt-1
696         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
697         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
698         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
699         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
700         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
701         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
702         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
703         IF( ref_wgts(kw)%cyclic ) THEN
704            WRITE(numout,*) '       cyclical'
705            IF( ref_wgts(kw)%offset > 0 ) WRITE(numout,*) '                 with offset'
706         ELSE
707            WRITE(numout,*) '       not cyclical'
708         ENDIF
709         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
710      END DO
711
712   END SUBROUTINE wgt_print
713
714   SUBROUTINE fld_weight( sd )
715      !!---------------------------------------------------------------------
716      !!                    ***  ROUTINE fld_weight  ***
717      !!
718      !! ** Purpose :   create a new WGT structure and fill in data from 
719      !!                file, restructuring as required
720      !!
721      !! ** Method  :   
722      !!----------------------------------------------------------------------
723      TYPE( FLD ),      INTENT(in)            ::   sd            ! field with name of weights file
724      !!
725      INTEGER                                 ::   jn            ! dummy loop indices
726      INTEGER                                 ::   inum          ! temporary logical unit
727      INTEGER                                 ::   id            ! temporary variable id
728      INTEGER                                 ::   ipk           ! temporary vertical dimension
729      CHARACTER (len=5)                       ::   aname
730      INTEGER , DIMENSION(3)                  ::   ddims
731      INTEGER , DIMENSION(jpi, jpj)           ::   data_src
732      REAL(wp), DIMENSION(jpi, jpj)           ::   data_tmp
733      REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::   line2, lines  ! temporary array to read 2 lineumns
734      CHARACTER (len=34)                      ::   lonvar
735      LOGICAL                                 ::   cyclical
736      REAL(wp)                                ::   resid, dlon   ! temporary array to read 2 lineumns
737      INTEGER                                 ::   offset        ! temporary integer
738      !!----------------------------------------------------------------------
739      !
740      IF( nxt_wgt > tot_wgts ) THEN
741        CALL ctl_stop("fld_weights: weights array size exceeded, increase tot_wgts")
742      ENDIF
743      !
744      !! new weights file entry, add in extra information
745      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
746      !! input data file is representative of all other files to be opened and processed with the
747      !! current weights file
748
749      !! open input data file (non-model grid)
750      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
751
752      !! get dimensions
753      id = iom_varid( inum, sd%clvar, ddims )
754
755      !! check for an east-west cyclic grid
756      !! try to guess name of longitude variable
757
758      lonvar = 'nav_lon'
759      id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
760      IF( id <= 0 ) THEN
761         lonvar = 'lon'
762         id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
763      ENDIF
764
765      offset = -1
766      cyclical = .FALSE.
767      IF( id > 0 ) THEN
768         !! found a longitude variable
769         !! now going to assume that grid is regular so we can read a single row
770
771         !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
772         !! worse, we cant pass line2(:,1) to iom_get since this is treated as a 1d array which doesnt match input file
773         ALLOCATE( lines(ddims(1),2) )
774         CALL iom_get(inum, jpdom_unknown, lonvar, lines(:,:), 1, kstart=(/1,1/), kcount=(/ddims(1),2/) )
775
776         !! find largest grid spacing
777         lines(1:ddims(1)-1,2) = lines(2:ddims(1),1) - lines(1:ddims(1)-1,1)
778         dlon = MAXVAL( lines(1:ddims(1)-1,2) )
779
780         resid = ABS(ABS(lines(ddims(1),1)-lines(1,1))-360.0)
781         IF( resid < rsmall ) THEN
782            !! end rows overlap in longitude
783            offset = 0
784            cyclical = .TRUE.
785         ELSEIF( resid < 2.0*dlon ) THEN
786            !! also call it cyclic if difference between end points is less than twice dlon from 360
787            offset = 1
788            cyclical = .TRUE.
789         ENDIF
790
791         DEALLOCATE( lines )
792
793      ELSE
794         !! guessing failed
795         !! read in first and last columns of data variable
796         !! since we dont know the name of the longitude variable (or even if there is one)
797         !! we assume that if these two columns are equal, file is cyclic east-west
798
799         !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
800         !! worse, we cant pass line2(1,:) to iom_get since this is treated as a 1d array which doesnt match input file
801         ALLOCATE( lines(2,ddims(2)), line2(2,ddims(2)) )
802         CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/1,1/), kcount=(/2,ddims(2)/) )
803         lines(2,:) = line2(1,:)
804
805         CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/ddims(1)-1,1/), kcount=(/2,ddims(2)/) )
806         lines(1,:) = line2(2,:)
807
808         resid = SUM( ABS(lines(1,:) - lines(2,:)) )
809         IF( resid < ddims(2)*rsmall ) THEN
810            offset = 0
811            cyclical = .TRUE.
812         ENDIF
813
814         DEALLOCATE( lines, line2 )
815      ENDIF
816
817      !! close it
818      CALL iom_close( inum )
819
820      !! now open the weights file
821
822      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
823      IF ( inum > 0 ) THEN
824
825         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
826         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
827         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
828         ref_wgts(nxt_wgt)%offset = -1
829         ref_wgts(nxt_wgt)%cyclic = .FALSE.
830         IF( cyclical ) THEN
831            ref_wgts(nxt_wgt)%offset = offset
832            ref_wgts(nxt_wgt)%cyclic = .TRUE.
833         ENDIF
834         ref_wgts(nxt_wgt)%nestid = 0
835#if defined key_agrif
836         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
837#endif
838         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
839         !! for each weight wgtNN there is an integer array srcNN which gives the point in
840         !! the input data grid which is to be multiplied by the weight
841         !! they are both arrays on the model grid so the result of the multiplication is
842         !! added into an output array on the model grid as a running sum
843
844         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
845         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
846         IF( id <= 0) THEN
847            ref_wgts(nxt_wgt)%numwgt = 4
848         ELSE
849            ref_wgts(nxt_wgt)%numwgt = 16
850         ENDIF
851
852         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
853         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
854         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
855
856         DO jn = 1,4
857            aname = ' '
858            WRITE(aname,'(a3,i2.2)') 'src',jn
859            data_tmp(:,:) = 0
860            CALL iom_get ( inum, jpdom_unknown, aname, data_tmp(1:nlci,1:nlcj), &
861                           kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) )
862            data_src(:,:) = INT(data_tmp(:,:))
863            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
864            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
865         END DO
866
867         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
868            aname = ' '
869            WRITE(aname,'(a3,i2.2)') 'wgt',jn
870            ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn) = 0.0
871            CALL iom_get ( inum, jpdom_unknown, aname, ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn), &
872                           kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) )
873         END DO
874         CALL iom_close (inum)
875 
876         ! find min and max indices in grid
877         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:))
878         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:))
879         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:))
880         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:))
881
882         ! and therefore dimensions of the input box
883         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
884         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
885
886         ! shift indexing of source grid
887         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
888         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
889
890         ! create input grid, give it a halo to allow gradient calculations
891         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
892         ! a more robust solution will be given in next release
893         ipk =  SIZE(sd%fdta,3)
894         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
895         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
896
897         nxt_wgt = nxt_wgt + 1
898
899      ELSE
900         CALL ctl_stop( '    fld_weight : unable to read the file ' )
901      ENDIF
902
903   END SUBROUTINE fld_weight
904
905   SUBROUTINE fld_interp(num, clvar, kw, kk, dta, nrec)
906      !!---------------------------------------------------------------------
907      !!                    ***  ROUTINE fld_interp  ***
908      !!
909      !! ** Purpose :   apply weights to input gridded data to create data
910      !!                on model grid
911      !!
912      !! ** Method  :   
913      !!----------------------------------------------------------------------
914      INTEGER,          INTENT(in)                           ::   num                 ! stream number
915      CHARACTER(LEN=*), INTENT(in)                           ::   clvar               ! variable name
916      INTEGER,          INTENT(in)                           ::   kw                  ! weights number
917      INTEGER,          INTENT(in)                           ::   kk                  ! vertical dimension of kk
918      REAL(wp),         INTENT(inout), DIMENSION(jpi,jpj,kk) ::   dta                 ! output field on model grid
919      INTEGER,          INTENT(in)                           ::   nrec                ! record number to read (ie time slice)
920      !!
921      INTEGER, DIMENSION(3)                                  ::   rec1,recn           ! temporary arrays for start and length
922      INTEGER                                                ::  jk, jn, jm           ! loop counters
923      INTEGER                                                ::  ni, nj               ! lengths
924      INTEGER                                                ::  jpimin,jpiwid        ! temporary indices
925      INTEGER                                                ::  jpjmin,jpjwid        ! temporary indices
926      INTEGER                                                ::  jpi1,jpi2,jpj1,jpj2  ! temporary indices
927      !!----------------------------------------------------------------------
928      !
929
930      !! for weighted interpolation we have weights at four corners of a box surrounding
931      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
932      !! or by a grid value and gradients at the corner point (bicubic case)
933      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
934
935      !! sub grid where we already have weights
936      jpimin = ref_wgts(kw)%botleft(1)
937      jpjmin = ref_wgts(kw)%botleft(2)
938      jpiwid = ref_wgts(kw)%jpiwgt
939      jpjwid = ref_wgts(kw)%jpjwgt
940
941      !! what we need to read into sub grid in order to calculate gradients
942      rec1(1) = MAX( jpimin-1, 1 )
943      rec1(2) = MAX( jpjmin-1, 1 )
944      rec1(3) = 1
945      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
946      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
947      recn(3) = kk
948
949      !! where we need to read it to
950      jpi1 = 2 + rec1(1) - jpimin
951      jpj1 = 2 + rec1(2) - jpjmin
952      jpi2 = jpi1 + recn(1) - 1
953      jpj2 = jpj1 + recn(2) - 1
954
955      ref_wgts(kw)%fly_dta(:,:,:) = 0.0
956      SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
957      CASE(1)
958           CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
959      CASE(jpk) 
960           CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
961      END SELECT 
962
963      !! first four weights common to both bilinear and bicubic
964      !! note that we have to offset by 1 into fly_dta array because of halo
965      dta(:,:,:) = 0.0
966      DO jk = 1,4
967        DO jn = 1, nlcj
968          DO jm = 1,nlci
969            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
970            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
971            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,jk)
972          END DO
973        END DO
974      END DO
975
976      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
977
978        !! fix up halo points that we couldnt read from file
979        IF( jpi1 == 2 ) THEN
980           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
981        ENDIF
982        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
983           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
984        ENDIF
985        IF( jpj1 == 2 ) THEN
986           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
987        ENDIF
988        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
989           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
990        ENDIF
991
992        !! if data grid is cyclic we can do better on east-west edges
993        !! but have to allow for whether first and last columns are coincident
994        IF( ref_wgts(kw)%cyclic ) THEN
995           rec1(2) = MAX( jpjmin-1, 1 )
996           recn(1) = 2
997           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
998           jpj1 = 2 + rec1(2) - jpjmin
999           jpj2 = jpj1 + recn(2) - 1
1000           IF( jpi1 == 2 ) THEN
1001              rec1(1) = ref_wgts(kw)%ddims(1) - 1
1002              SELECT CASE( SIZE( ref_wgts(kw)%col2(:,jpj1:jpj2,:),3) )
1003              CASE(1)
1004                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn)
1005              CASE(jpk)         
1006                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn)
1007              END SELECT     
1008              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2,:)
1009           ENDIF
1010           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1011              rec1(1) = 1
1012              SELECT CASE( SIZE( ref_wgts(kw)%col2(:,jpj1:jpj2,:),3) )
1013              CASE(1)
1014                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn)
1015              CASE(jpk)
1016                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn)
1017              END SELECT
1018              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2,:)
1019           ENDIF
1020        ENDIF
1021
1022        ! gradient in the i direction
1023        DO jk = 1,4
1024          DO jn = 1, nlcj
1025            DO jm = 1,nlci
1026              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1027              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1028              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1029                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1030            END DO
1031          END DO
1032        END DO
1033
1034        ! gradient in the j direction
1035        DO jk = 1,4
1036          DO jn = 1, nlcj
1037            DO jm = 1,nlci
1038              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1039              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1040              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         &
1041                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1042            END DO
1043          END DO
1044        END DO
1045
1046        ! gradient in the ij direction
1047        DO jk = 1,4
1048          DO jn = 1, jpj
1049            DO jm = 1,jpi
1050              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1051              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1052              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1053                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1054                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
1055            END DO
1056          END DO
1057        END DO
1058
1059      END IF
1060
1061   END SUBROUTINE fld_interp
1062 
1063END MODULE fldread
Note: See TracBrowser for help on using the repository browser.