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

source: branches/devmercator2010_1/NEMO/OPA_SRC/SBC/fldread.F90 @ 2134

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

remove debug print and comment

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