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 @ 2131

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

add changes from dev_1784_WEEK

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