source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 10 years ago

update licence of all NEMO files…

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