source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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