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/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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