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/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5232

Last change on this file since 5232 was 5232, checked in by davestorkey, 9 years ago

Svn keywords deactivated using "svn propdel" in
branch 2015/dev_r5021_UKMO1_CICE_coupling.

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