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

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2777

Last change on this file since 2777 was 2777, checked in by smasson, 13 years ago

LIM3 compiling and (partly?) running in v3_3_1, see ticket#817

  • Property svn:keywords set to Id
File size: 61.5 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
11   !!----------------------------------------------------------------------
12   !!   fld_read      : read input fields used for the computation of the
13   !!                   surface boundary condition
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE phycst          ! ???
18   USE in_out_manager  ! I/O manager
19   USE iom             ! I/O manager library
20   USE geo2ocean       ! for vector rotation on to model grid
21   USE lib_mpp         ! MPP library
22   USE ioipsl, ONLY :   ymds2ju, ju2ymds   ! for calendar
23
24   IMPLICIT NONE
25   PRIVATE   
26
27   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
28      CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file
29      INTEGER              ::   nfreqh      ! frequency of each flux file
30      CHARACTER(len = 34)  ::   clvar       ! generic name of the variable in the NetCDF flux file
31      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F)
32      LOGICAL              ::   ln_clim     ! climatology or not (T/F)
33      CHARACTER(len = 8)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly'
34      CHARACTER(len = 34)  ::   wname       ! generic name of a NetCDF weights file to be used, blank if not
35      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation
36      !                                     ! a string starting with "U" or "V" for each component   
37      !                                     ! chars 2 onwards identify which components go together 
38   END TYPE FLD_N
39
40   TYPE, PUBLIC ::   FLD        !: Input field related variables
41      CHARACTER(len = 256)            ::   clrootname   ! generic name of the NetCDF file
42      CHARACTER(len = 256)            ::   clname       ! current name of the NetCDF file
43      INTEGER                         ::   nfreqh       ! frequency of each flux file
44      CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file
45      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F)
46      LOGICAL                         ::   ln_clim      ! climatology or not (T/F)
47      CHARACTER(len = 8)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly'
48      INTEGER                         ::   num          ! iom id of the jpfld files to be read
49      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
50      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year)
51      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow       ! input fields interpolated to now time step
52      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields
53      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
54      !                                                 ! into the WGTLIST structure
55      CHARACTER(len = 34)             ::   vcomp        ! symbolic name for a vector component that needs rotation
56      LOGICAL                         ::   rotn         ! flag to indicate whether field has been rotated
57   END TYPE FLD
58
59!$AGRIF_DO_NOT_TREAT
60
61   !! keep list of all weights variables so they're only read in once
62   !! need to add AGRIF directives not to process this structure
63   !! also need to force wgtname to include AGRIF nest number
64   TYPE         ::   WGT        !: Input weights related variables
65      CHARACTER(len = 256)                    ::   wgtname      ! current name of the NetCDF weight file
66      INTEGER , DIMENSION(2)                  ::   ddims        ! shape of input grid
67      INTEGER , DIMENSION(2)                  ::   botleft      ! top left corner of box in input grid containing
68      !                                                         ! current processor grid
69      INTEGER , DIMENSION(2)                  ::   topright     ! top right corner of box
70      INTEGER                                 ::   jpiwgt       ! width of box on input grid
71      INTEGER                                 ::   jpjwgt       ! height of box on input grid
72      INTEGER                                 ::   numwgt       ! number of weights (4=bilinear, 16=bicubic)
73      INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in
74      INTEGER                                 ::   overlap      ! =0 when cyclic grid has no overlapping EW columns
75      !                                                         ! =>1 when they have one or more overlapping columns     
76      !                                                         ! =-1 not cyclic
77      LOGICAL                                 ::   cyclic       ! east-west cyclic or not
78      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpi     ! array of source integers
79      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpj     ! array of source integers
80      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid
81      REAL(wp), DIMENSION(:,:,:), POINTER     ::   fly_dta      ! array of values on input grid
82      REAL(wp), DIMENSION(:,:,:), POINTER     ::   col          ! temporary array for reading in columns
83   END TYPE WGT
84
85   INTEGER,     PARAMETER             ::   tot_wgts = 10
86   TYPE( WGT ), DIMENSION(tot_wgts)   ::   ref_wgts     ! array of wgts
87   INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array
88
89!$AGRIF_END_DO_NOT_TREAT
90
91   PUBLIC   fld_read, fld_fill   ! called by sbc... modules
92
93   !!----------------------------------------------------------------------
94   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
95   !! $Id$
96   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
97   !!----------------------------------------------------------------------
98CONTAINS
99
100   SUBROUTINE fld_read( kt, kn_fsbc, sd )
101      !!---------------------------------------------------------------------
102      !!                    ***  ROUTINE fld_read  ***
103      !!                   
104      !! ** Purpose :   provide at each time step the surface ocean fluxes
105      !!                (momentum, heat, freshwater and runoff)
106      !!
107      !! ** Method  :   READ each input fields in NetCDF files using IOM
108      !!      and intepolate it to the model time-step.
109      !!         Several assumptions are made on the input file:
110      !!      blahblahblah....
111      !!----------------------------------------------------------------------
112      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
113      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)
114      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
115      !!
116      INTEGER  ::   imf        ! size of the structure sd
117      INTEGER  ::   jf         ! dummy indices
118      INTEGER  ::   ireclast   ! last record to be read in the current year file
119      INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend
120      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
121      LOGICAL  ::   llnxtyr    ! open next year  file?
122      LOGICAL  ::   llnxtmth   ! open next month file?
123      LOGICAL  ::   llstop     ! stop is the file does not exist
124      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
125      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
126      CHARACTER(LEN=1000) ::   clfmt   ! write format
127      !!---------------------------------------------------------------------
128      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
129      isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1))   ! middle of sbc time step
130      imf = SIZE( sd )
131      !
132      IF( kt == nit000 ) THEN                      ! initialization
133         DO jf = 1, imf 
134            CALL fld_init( kn_fsbc, sd(jf) )       ! read each before field (put them in after as they will be swapped)
135         END DO
136         IF( lwp ) CALL wgt_print()                ! control print
137         CALL fld_rot( kt, sd )                    ! rotate vector fiels if needed
138      ENDIF
139      !                                            ! ====================================== !
140      IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! update field at each kn_fsbc time-step !
141         !                                         ! ====================================== !
142         !
143         DO jf = 1, imf                            ! ---   loop over field   --- !
144           
145            IF( isecsbc > sd(jf)%nrec_a(2) .OR. kt == nit000 ) THEN  ! read/update the after data?
146
147               IF( sd(jf)%ln_tint ) THEN                             ! swap before record field and informations
148                  sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)
149!CDIR COLLAPSE
150                  sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)
151               ENDIF
152
153               CALL fld_rec( kn_fsbc, sd(jf) )                       ! update record informations
154
155               ! do we have to change the year/month/week/day of the forcing field??
156               IF( sd(jf)%ln_tint ) THEN
157                  ! if we do time interpolation we will need to open next year/month/week/day file before the end of the current
158                  ! 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)
159                  ! will be larger than the record number that should be read for current year/month/week/day
160
161                  ! last record to be read in the current file
162                  IF    ( sd(jf)%nfreqh == -12 ) THEN                 ;   ireclast = 1    !  yearly mean
163                  ELSEIF( sd(jf)%nfreqh ==  -1 ) THEN                                     ! monthly mean
164                     IF(     sd(jf)%cltype      == 'monthly' ) THEN   ;   ireclast = 1
165                     ELSE                                             ;   ireclast = 12
166                     ENDIF
167                  ELSE                                                                    ! higher frequency mean (in hours)
168                     IF(     sd(jf)%cltype      == 'monthly' ) THEN   ;   ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh 
169                     ELSEIF( sd(jf)%cltype(1:4) == 'week'    ) THEN   ;   ireclast = 24 * 7                  / sd(jf)%nfreqh
170                     ELSEIF( sd(jf)%cltype      == 'daily'   ) THEN   ;   ireclast = 24                      / sd(jf)%nfreqh
171                     ELSE                                             ;   ireclast = 24 * nyear_len(     1 ) / sd(jf)%nfreqh 
172                     ENDIF
173                  ENDIF
174
175                  ! do we need next file data?
176                  IF( sd(jf)%nrec_a(1) > ireclast ) THEN
177
178                     sd(jf)%nrec_a(1) = 1              ! force to read the first record of the next file
179
180                     IF( .NOT. sd(jf)%ln_clim ) THEN   ! close the current file and open a new one.
181
182                        llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth)      ! open next month file?
183                        llnxtyr  = sd(jf)%cltype == 'yearly'  .OR. (nmonth == 12 .AND. llnxtmth)   ! open next year  file?
184
185                        ! if the run finishes at the end of the current year/month/week/day, we will allow next
186                        ! year/month/week/day file to be not present. If the run continue further than the current
187                        ! year/month/week/day, next year/month/week/day file must exist
188                        isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1))   ! second at the end of the run
189                        llstop = isecend > sd(jf)%nrec_a(2)                                   ! read more than 1 record of next year
190
191                        CALL fld_clopn( sd(jf), nyear  + COUNT((/llnxtyr /))                                           ,         &
192                           &                    nmonth + COUNT((/llnxtmth/)) - 12                 * COUNT((/llnxtyr /)),         &
193                           &                    nday   + 1                   - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
194
195                        IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN    ! next year file does not exist
196                           CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)//     &
197                              &     ' not present -> back to current year/month/day')
198                           CALL fld_clopn( sd(jf), nyear, nmonth, nday )       ! back to the current year/month/day
199                           sd(jf)%nrec_a(1) = ireclast     ! force to read the last record to be read in the current year file
200                        ENDIF
201
202                     ENDIF
203                  ENDIF
204
205               ELSE
206                  ! if we are not doing time interpolation, we must change the year/month/week/day of the file just after
207                  ! switching to the NEW year/month/week/day. If it is the case, we are at the beginning of the
208                  ! year/month/week/day when calling fld_rec so sd(jf)%nrec_a(1) = 1
209                  IF( sd(jf)%nrec_a(1) == 1 .AND. .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) )   &
210                     &   CALL fld_clopn( sd(jf), nyear, nmonth, nday )
211               ENDIF
212
213               ! read after data
214               CALL fld_get( sd(jf) )
215
216            ENDIF
217         END DO                                    ! --- end loop over field --- !
218
219         CALL fld_rot( kt, sd )                    ! rotate vector fiels if needed
220
221         DO jf = 1, imf                            ! ---   loop over field   --- !
222            !
223            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
224               IF(lwp .AND. kt - nit000 <= 100 ) THEN
225                  clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
226                     &    "', records b/a: ', i4.4, '/', i4.4, ' (days ', f7.2,'/', f7.2, ')')"
227                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &
228                     & 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
229               ENDIF
230               ! temporal interpolation weights
231               ztinta =  REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
232               ztintb =  1. - ztinta
233!CDIR COLLAPSE
234               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
235            ELSE   ! nothing to do...
236               IF(lwp .AND. kt - nit000 <= 100 ) THEN
237                  clfmt = "('fld_read: var ', a, ' kt = ', i8,' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
238                     &    "', record: ', i4.4, ' (days ', f7.2, ' <-> ', f7.2, ')')"
239                  WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
240                     &                 sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
241               ENDIF
242            ENDIF
243            !
244            IF( kt == nitend - kn_fsbc + 1 )   CALL iom_close( sd(jf)%num )   ! Close the input files
245
246         END DO                                    ! --- end loop over field --- !
247         !
248         !                                         ! ====================================== !
249      ENDIF                                        ! update field at each kn_fsbc time-step !
250      !                                            ! ====================================== !
251      !
252   END SUBROUTINE fld_read
253
254
255   SUBROUTINE fld_init( kn_fsbc, sdjf )
256      !!---------------------------------------------------------------------
257      !!                    ***  ROUTINE fld_init  ***
258      !!
259      !! ** Purpose :  - if time interpolation, read before data
260      !!               - open current year file
261      !!----------------------------------------------------------------------
262      INTEGER  , INTENT(in   ) ::   kn_fsbc   ! sbc computation period (in time step)
263      TYPE(FLD), INTENT(inout) ::   sdjf      ! input field related variables
264      !!
265      LOGICAL :: llprevyr              ! are we reading previous year  file?
266      LOGICAL :: llprevmth             ! are we reading previous month file?
267      LOGICAL :: llprevweek            ! are we reading previous week  file?
268      LOGICAL :: llprevday             ! are we reading previous day   file?
269      LOGICAL :: llprev                ! llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
270      INTEGER :: idvar                 ! variable id
271      INTEGER :: inrec                 ! number of record existing for this variable
272      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
273      INTEGER :: isec_week             ! number of seconds since start of the weekly file
274      CHARACTER(LEN=1000) ::   clfmt   ! write format
275      !!---------------------------------------------------------------------
276     
277      ! some default definitions...
278      sdjf%num = 0   ! default definition for non-opened file
279      IF( sdjf%ln_clim )   sdjf%clname = TRIM( sdjf%clrootname )   ! file name defaut definition, never change in this case
280      llprevyr   = .FALSE.
281      llprevmth  = .FALSE.
282      llprevweek = .FALSE.
283      llprevday  = .FALSE.
284      isec_week  = 0
285           
286      IF( sdjf%cltype(1:4) == 'week' .AND. nn_leapy == 0 )   &
287         &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdjf%clrootname)//') needs nn_leapy = 1')
288      IF( sdjf%cltype(1:4) == 'week' .AND. sdjf%ln_clim  )   &
289         &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdjf%clrootname)//') needs ln_clim = .FALSE.')
290
291      ! define record informations
292      CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. )  ! return before values in sdjf%nrec_a (as we will swap it later)
293
294      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
295
296      IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
297
298         IF( sdjf%nrec_a(1) == 0  ) THEN   ! we redefine record sdjf%nrec_a(1) with the last record of previous year file
299            IF    ( sdjf%nfreqh == -12 ) THEN   ! yearly mean
300               IF( sdjf%cltype == 'yearly' ) THEN             ! yearly file
301                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
302                  llprevyr  = .NOT. sdjf%ln_clim                                           ! use previous year  file?
303               ELSE
304                  CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clname) )
305               ENDIF
306            ELSEIF( sdjf%nfreqh ==  -1 ) THEN   ! monthly mean
307               IF( sdjf%cltype == 'monthly' ) THEN            ! monthly file
308                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
309                  llprevmth = .TRUE.                                                       ! use previous month file?
310                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
311               ELSE                                           ! yearly file
312                  sdjf%nrec_a(1) = 12                                                      ! force to read december mean
313                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
314               ENDIF
315            ELSE                                ! higher frequency mean (in hours)
316               IF    ( sdjf%cltype      == 'monthly' ) THEN   ! monthly file
317                  sdjf%nrec_a(1) = 24 * nmonth_len(nmonth-1) / sdjf%nfreqh                 ! last record of previous month
318                  llprevmth = .TRUE.                                                       ! use previous month file?
319                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
320               ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ! weekly file
321                  llprevweek = .TRUE.                                                      ! use previous week  file?
322                  sdjf%nrec_a(1) = 24 * 7 / sdjf%nfreqh                                    ! last record of previous week
323                  isec_week = NINT(rday) * 7                                               ! add a shift toward previous week
324               ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ! daily file
325                  sdjf%nrec_a(1) = 24 / sdjf%nfreqh                                        ! last record of previous day
326                  llprevday = .TRUE.                                                       ! use previous day   file?
327                  llprevmth = llprevday .AND. nday   == 1                                  ! use previous month file?
328                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
329               ELSE                                           ! yearly file
330                  sdjf%nrec_a(1) = 24 * nyear_len(0) / sdjf%nfreqh                         ! last record of previous year
331                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
332               ENDIF
333            ENDIF
334         ENDIF
335         IF ( sdjf%cltype(1:4) == 'week' ) THEN
336            isec_week = isec_week + ksec_week( sdjf%cltype(6:8) )   ! second since the beginning of the week
337            llprevmth = isec_week > nsec_month                      ! longer time since the beginning of the week than the month
338            llprevyr  = llprevmth .AND. nmonth == 1
339         ENDIF
340         llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
341         !
342         iyear  = nyear  - COUNT((/llprevyr /))
343         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
344         iday   = nday   - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
345         !
346         CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev )
347
348         ! if previous year/month/day file does not exist, we switch to the current year/month/day
349         IF( llprev .AND. sdjf%num <= 0 ) THEN
350            CALL ctl_warn( 'previous year/month/week/day file: '//TRIM(sdjf%clname)//   &
351               &           ' not present -> back to current year/month/week/day' )
352            ! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
353            llprev = .FALSE.
354            sdjf%nrec_a(1) = 1
355            CALL fld_clopn( sdjf, nyear, nmonth, nday )
356         ENDIF
357         
358         IF( llprev ) THEN   ! check if the last record sdjf%nrec_n(1) exists in the file
359            idvar = iom_varid( sdjf%num, sdjf%clvar )                                        ! id of the variable sdjf%clvar
360            IF( idvar <= 0 )   RETURN
361            inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar )   ! size of the last dim of idvar
362            sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec )   ! make sure we select an existing record
363         ENDIF
364
365         ! read before data
366         CALL fld_get( sdjf )  ! read before values in after arrays(as we will swap it later)
367
368         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i4, ' at time = ', f7.2, ' days')"
369         IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
370
371         IF( llprev )   CALL iom_close( sdjf%num )          ! force to close previous year file (-> redefine sdjf%num to 0)
372
373      ENDIF
374
375      ! make sure current year/month/day file is opened
376      IF( sdjf%num <= 0 ) THEN
377         !
378         IF ( sdjf%cltype(1:4) == 'week' ) THEN
379            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
380            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
381            llprevyr   = llprevmth .AND. nmonth == 1
382         ELSE
383            isec_week  = 0
384            llprevmth  = .FALSE.
385            llprevyr   = .FALSE.
386         ENDIF
387         !
388         iyear  = nyear  - COUNT((/llprevyr /))
389         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
390         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
391         !
392         CALL fld_clopn( sdjf, iyear, imonth, iday )
393      ENDIF 
394      !
395   END SUBROUTINE fld_init
396
397
398   SUBROUTINE fld_rec( kn_fsbc, sdjf, ldbefore )
399      !!---------------------------------------------------------------------
400      !!                    ***  ROUTINE fld_rec  ***
401      !!
402      !! ** Purpose : Compute
403      !!              if sdjf%ln_tint = .TRUE.
404      !!                  nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping)
405      !!              if sdjf%ln_tint = .FALSE.
406      !!                  nrec_a(1): record number
407      !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record (for print only)
408      !!----------------------------------------------------------------------
409      INTEGER  , INTENT(in   )           ::   kn_fsbc   ! sbc computation period (in time step)
410      TYPE(FLD), INTENT(inout)           ::   sdjf      ! input field related variables
411      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldbefore  ! sent back before record values (default = .FALSE.)
412                                                        ! used only if sdjf%ln_tint = .TRUE.
413      !!
414      LOGICAL  ::   llbefore    ! local definition of ldbefore
415      INTEGER  ::   iendrec     ! end of this record (in seconds)
416      INTEGER  ::   imth        ! month number
417      INTEGER  ::   ifreq_sec   ! frequency mean (in seconds)
418      INTEGER  ::   isec_week   ! number of seconds since the start of the weekly file
419      REAL(wp) ::   ztmp        ! temporary variable
420      !!----------------------------------------------------------------------
421      !
422      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
423      !
424      IF( PRESENT(ldbefore) ) THEN   ;   llbefore = ldbefore .AND. sdjf%ln_tint   ! needed only if sdjf%ln_tint = .TRUE.
425      ELSE                           ;   llbefore = .FALSE.
426      ENDIF
427      !
428      !                                      ! =========== !
429      IF    ( sdjf%nfreqh == -12 ) THEN      ! yearly mean
430         !                                   ! =========== !
431         !
432         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
433            !
434            !                  INT( ztmp )
435            !                     /|\
436            !                    1 |    *----
437            !                    0 |----(             
438            !                      |----+----|--> time
439            !                      0   /|\   1   (nday/nyear_len(1))
440            !                           |   
441            !                           |   
442            !       forcing record :    1
443            !                           
444            ztmp = REAL( nday, wp ) / REAL( nyear_len(1), wp ) + 0.5
445            sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
446            ! swap at the middle of the year
447            IF( llbefore ) THEN   ;   sdjf%nrec_a(2) = nsec1jan000 - NINT(0.5 * rday) * nyear_len(0)
448            ELSE                  ;   sdjf%nrec_a(2) = nsec1jan000 + NINT(0.5 * rday) * nyear_len(1)   
449            ENDIF
450         ELSE                                    ! no time interpolation
451            sdjf%nrec_a(1) = 1
452            sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000   ! swap at the end    of the year
453            sdjf%nrec_b(2) = nsec1jan000                               ! beginning of the year (only for print)
454         ENDIF
455         !
456         !                                   ! ============ !
457      ELSEIF( sdjf%nfreqh ==  -1 ) THEN      ! monthly mean !
458         !                                   ! ============ !
459         !
460         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
461            !
462            !                  INT( ztmp )
463            !                     /|\
464            !                    1 |    *----
465            !                    0 |----(             
466            !                      |----+----|--> time
467            !                      0   /|\   1   (nday/nmonth_len(nmonth))
468            !                           |   
469            !                           |   
470            !       forcing record :  nmonth
471            !                           
472            ztmp = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5
473            imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
474            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
475            ELSE                                  ;   sdjf%nrec_a(1) = imth
476            ENDIF
477            sdjf%nrec_a(2) = nmonth_half(   imth ) + nsec1jan000   ! swap at the middle of the month
478         ELSE                                    ! no time interpolation
479            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1
480            ELSE                                  ;   sdjf%nrec_a(1) = nmonth
481            ENDIF
482            sdjf%nrec_a(2) =  nmonth_end(nmonth  ) + nsec1jan000   ! swap at the end    of the month
483            sdjf%nrec_b(2) =  nmonth_end(nmonth-1) + nsec1jan000   ! beginning of the month (only for print)
484         ENDIF
485         !
486         !                                   ! ================================ !
487      ELSE                                   ! higher frequency mean (in hours)
488         !                                   ! ================================ !
489         !
490         ifreq_sec = sdjf%nfreqh * 3600                                                 ! frequency mean (in seconds)
491         IF( sdjf%cltype(1:4) == 'week' )   isec_week = ksec_week( sdjf%cltype(6:8) )   ! since the first day of the current week
492         ! number of second since the beginning of the file
493         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month,wp)  ! since the first day of the current month
494         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ztmp = REAL(isec_week ,wp)  ! since the first day of the current week
495         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day  ,wp)  ! since 00h of the current day
496         ELSE                                           ;   ztmp = REAL(nsec_year ,wp)  ! since 00h on Jan 1 of the current year
497         ENDIF
498         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1)   ! shift time to be centrered in the middle of sbc time step
499         ztmp = ztmp + 0.01 * rdttra(1)                          ! add 0.01 time step to avoid truncation error
500         IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record
501            !
502            !                  INT( ztmp )
503            !                     /|\
504            !                    2 |        *-----(
505            !                    1 |  *-----(
506            !                    0 |--(             
507            !                      |--+--|--+--|--+--|--> time
508            !                      0 /|\ 1 /|\ 2 /|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
509            !                         |     |     |
510            !                         |     |     |
511            !       forcing record :  1     2     3
512            !                   
513            ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
514         ELSE                                   ! no time interpolation
515            !
516            !                  INT( ztmp )
517            !                     /|\
518            !                    2 |           *-----(
519            !                    1 |     *-----(
520            !                    0 |-----(             
521            !                      |--+--|--+--|--+--|--> time
522            !                      0 /|\ 1 /|\ 2 /|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
523            !                         |     |     |
524            !                         |     |     |
525            !       forcing record :  1     2     3
526            !                           
527            ztmp= ztmp / REAL(ifreq_sec, wp)
528         ENDIF
529         sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))   ! record nomber to be read
530
531         iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000       ! end of this record (in second)
532         ! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
533         IF( sdjf%cltype      == 'monthly' )   iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
534         IF( sdjf%cltype(1:4) == 'week'    )   iendrec = iendrec + ( nsec_year - isec_week )
535         IF( sdjf%cltype      == 'daily'   )   iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
536         IF( sdjf%ln_tint ) THEN
537             sdjf%nrec_a(2) = iendrec - ifreq_sec / 2        ! swap at the middle of the record
538         ELSE
539             sdjf%nrec_a(2) = iendrec                        ! swap at the end    of the record
540             sdjf%nrec_b(2) = iendrec - ifreq_sec            ! beginning of the record (only for print)
541         ENDIF
542         !
543      ENDIF
544      !
545   END SUBROUTINE fld_rec
546
547
548   SUBROUTINE fld_get( sdjf )
549      !!---------------------------------------------------------------------
550      !!                    ***  ROUTINE fld_clopn  ***
551      !!
552      !! ** Purpose :   read the data
553      !!----------------------------------------------------------------------
554      TYPE(FLD), INTENT(inout) ::   sdjf   ! input field related variables
555      !!
556      INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
557      INTEGER                  ::   iw     ! index into wgts array
558      !!---------------------------------------------------------------------
559           
560      ipk = SIZE( sdjf%fnow, 3 )
561      IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
562         CALL wgt_list( sdjf, iw )
563         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
564         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) )
565         ENDIF
566      ELSE
567         SELECT CASE( ipk )
568         CASE(1)   
569            IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
570            ELSE                      ;   CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,1  ), sdjf%nrec_a(1) )
571            ENDIF
572         CASE DEFAULT
573            IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
574            ELSE                      ;   CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) )
575            ENDIF
576         END SELECT
577      ENDIF
578      !
579      sdjf%rotn = .false.   ! vector not yet rotated
580
581   END SUBROUTINE fld_get
582
583
584   SUBROUTINE fld_rot( kt, sd )
585      !!---------------------------------------------------------------------
586      !!                    ***  ROUTINE fld_clopn  ***
587      !!
588      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
589      !!----------------------------------------------------------------------
590      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
591      USE wrk_nemo, ONLY: utmp => wrk_2d_4, vtmp => wrk_2d_5      ! 2D workspace
592      !!
593      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
594      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
595      !!
596      INTEGER                      ::   ju, jv, jk   ! loop indices
597      INTEGER                      ::   imf          ! size of the structure sd
598      INTEGER                      ::   ill          ! character length
599      INTEGER                      ::   iv           ! indice of V component
600      CHARACTER (LEN=100)          ::   clcomp       ! dummy weight name
601      !!---------------------------------------------------------------------
602
603      IF(wrk_in_use(2, 4,5) ) THEN
604         CALL ctl_stop('fld_rot: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
605      END IF
606
607      !! (sga: following code should be modified so that pairs arent searched for each time
608      !
609      imf = SIZE( sd )
610      DO ju = 1, imf
611         ill = LEN_TRIM( sd(ju)%vcomp )
612         IF( ill > 0 .AND. .NOT. sd(ju)%rotn ) THEN   ! find vector rotations required             
613             IF( sd(ju)%vcomp(1:1) == 'U' ) THEN      ! east-west component has symbolic name starting with 'U'
614                ! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
615                clcomp = 'V' // sd(ju)%vcomp(2:ill)   ! works even if ill == 1
616                iv = -1
617                DO jv = 1, imf
618                  IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
619                END DO
620                IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together
621                   DO jk = 1, SIZE( sd(ju)%fnow, 3 )
622                      IF( sd(ju)%ln_tint )THEN
623                         CALL rot_rep( sd(ju)%fdta(:,:,jk,2), sd(iv)%fdta(:,:,jk,2), 'T', 'en->i', utmp(:,:) )
624                         CALL rot_rep( sd(ju)%fdta(:,:,jk,2), sd(iv)%fdta(:,:,jk,2), 'T', 'en->j', vtmp(:,:) )
625                         sd(ju)%fdta(:,:,jk,2) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,2) = vtmp(:,:)
626                      ELSE
627                         CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->i', utmp(:,:) )
628                         CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->j', vtmp(:,:) )
629                         sd(ju)%fnow(:,:,jk  ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk  ) = vtmp(:,:)
630                      ENDIF
631                   END DO
632                   sd(ju)%rotn = .TRUE.               ! vector was rotated
633                   IF( lwp .AND. kt == nit000 )   WRITE(numout,*)   &
634                      &   'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
635                ENDIF
636             ENDIF
637          ENDIF
638       END DO
639      !
640      IF(wrk_not_released(2, 4,5) )    CALL ctl_stop('fld_rot: ERROR: failed to release workspace arrays.')
641      !
642   END SUBROUTINE fld_rot
643
644
645   SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
646      !!---------------------------------------------------------------------
647      !!                    ***  ROUTINE fld_clopn  ***
648      !!
649      !! ** Purpose :   update the file name and open the file
650      !!----------------------------------------------------------------------
651      TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
652      INTEGER          , INTENT(in   ) ::   kyear    ! year value
653      INTEGER          , INTENT(in   ) ::   kmonth   ! month value
654      INTEGER          , INTENT(in   ) ::   kday     ! day value
655      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
656      !!----------------------------------------------------------------------
657
658      IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
659      ! build the new filename if not climatological data
660      sdjf%clname=TRIM(sdjf%clrootname)
661      !
662      ! note that sdjf%ln_clim is is only acting on presence of the year in the file
663      IF( .NOT. sdjf%ln_clim ) THEN   
664                                         WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear    ! add year
665         IF( sdjf%cltype /= 'yearly' )   WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname     ), kmonth   ! add month
666      ELSE
667         ! build the new filename if climatological data
668         IF( sdjf%cltype /= 'yearly' )   WRITE(sdjf%clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth   ! add month
669      ENDIF
670      IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
671            &                            WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname     ), kday     ! add day
672      !
673      CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
674      !
675   END SUBROUTINE fld_clopn
676
677
678   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
679      !!---------------------------------------------------------------------
680      !!                    ***  ROUTINE fld_fill  ***
681      !!
682      !! ** Purpose :   fill sdf with sdf_n and control print
683      !!----------------------------------------------------------------------
684      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
685      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
686      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
687      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   !
688      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    !
689      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      !
690      !
691      INTEGER  ::   jf       ! dummy indices
692      !!---------------------------------------------------------------------
693
694      DO jf = 1, SIZE(sdf)
695         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
696         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
697         sdf(jf)%clvar      = sdf_n(jf)%clvar
698         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
699         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
700         sdf(jf)%cltype     = sdf_n(jf)%cltype
701         sdf(jf)%wgtname = " "
702         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
703         sdf(jf)%vcomp   = sdf_n(jf)%vcomp
704      END DO
705
706      IF(lwp) THEN      ! control print
707         WRITE(numout,*)
708         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
709         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
710         WRITE(numout,*) '          '//TRIM( cdnam )//' Namelist'
711         WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
712         DO jf = 1, SIZE(sdf)
713            WRITE(numout,*) '               root filename: '  , TRIM( sdf(jf)%clrootname ),   &
714               &                          ' variable name: '  , TRIM( sdf(jf)%clvar      )
715            WRITE(numout,*) '               frequency: '      ,       sdf(jf)%nfreqh      ,   &
716               &                          ' time interp: '    ,       sdf(jf)%ln_tint     ,   &
717               &                          ' climatology: '    ,       sdf(jf)%ln_clim     ,   &
718               &                          ' weights    : '    , TRIM( sdf(jf)%wgtname    ),   &
719               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   &
720               &                          ' data type: '      ,       sdf(jf)%cltype
721            call flush(numout)
722         END DO
723      ENDIF
724     
725   END SUBROUTINE fld_fill
726
727
728   SUBROUTINE wgt_list( sd, kwgt )
729      !!---------------------------------------------------------------------
730      !!                    ***  ROUTINE wgt_list  ***
731      !!
732      !! ** Purpose :   search array of WGTs and find a weights file
733      !!                entry, or return a new one adding it to the end
734      !!                if it is a new entry, the weights data is read in and
735      !!                restructured (fld_weight)
736      !!----------------------------------------------------------------------
737      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
738      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
739      !!
740      INTEGER ::   kw, nestid   ! local integer
741      LOGICAL ::   found        ! local logical
742      !!----------------------------------------------------------------------
743      !
744      !! search down linked list
745      !! weights filename is either present or we hit the end of the list
746      found = .FALSE.
747
748      !! because agrif nest part of filenames are now added in iom_open
749      !! to distinguish between weights files on the different grids, need to track
750      !! nest number explicitly
751      nestid = 0
752#if defined key_agrif
753      nestid = Agrif_Fixed()
754#endif
755      DO kw = 1, nxt_wgt-1
756         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
757             ref_wgts(kw)%nestid == nestid) THEN
758            kwgt = kw
759            found = .TRUE.
760            EXIT
761         ENDIF
762      END DO
763      IF( .NOT.found ) THEN
764         kwgt = nxt_wgt
765         CALL fld_weight( sd )
766      ENDIF
767      !
768   END SUBROUTINE wgt_list
769
770
771   SUBROUTINE wgt_print( )
772      !!---------------------------------------------------------------------
773      !!                    ***  ROUTINE wgt_print  ***
774      !!
775      !! ** Purpose :   print the list of known weights
776      !!----------------------------------------------------------------------
777      INTEGER ::   kw   !
778      !!----------------------------------------------------------------------
779      !
780      DO kw = 1, nxt_wgt-1
781         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
782         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
783         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
784         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
785         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
786         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
787         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
788         IF( ref_wgts(kw)%cyclic ) THEN
789            WRITE(numout,*) '       cyclical'
790            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap
791         ELSE
792            WRITE(numout,*) '       not cyclical'
793         ENDIF
794         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
795      END DO
796      !
797   END SUBROUTINE wgt_print
798
799
800   SUBROUTINE fld_weight( sd )
801      !!---------------------------------------------------------------------
802      !!                    ***  ROUTINE fld_weight  ***
803      !!
804      !! ** Purpose :   create a new WGT structure and fill in data from 
805      !!                file, restructuring as required
806      !!----------------------------------------------------------------------
807      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released
808      USE wrk_nemo, ONLY:   data_tmp =>  wrk_2d_1     ! 2D real    workspace
809      USE wrk_nemo, ONLY:   data_src => iwrk_2d_1     ! 2D integer workspace
810      !!
811      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
812      !!
813      INTEGER                ::   jn            ! dummy loop indices
814      INTEGER                ::   inum          ! temporary logical unit
815      INTEGER                ::   id            ! temporary variable id
816      INTEGER                ::   ipk           ! temporary vertical dimension
817      CHARACTER (len=5)      ::   aname
818      INTEGER , DIMENSION(3) ::   ddims
819      LOGICAL                ::   cyclical
820      INTEGER                ::   zwrap      ! local integer
821      !!----------------------------------------------------------------------
822      !
823      IF(  wrk_in_use(2, 1)  .OR.  iwrk_in_use(2,1) ) THEN
824         CALL ctl_stop('fld_weight: requested workspace arrays are unavailable')   ;   RETURN
825      ENDIF
826      !
827      IF( nxt_wgt > tot_wgts ) THEN
828        CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
829      ENDIF
830      !
831      !! new weights file entry, add in extra information
832      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
833      !! input data file is representative of all other files to be opened and processed with the
834      !! current weights file
835
836      !! open input data file (non-model grid)
837      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
838
839      !! get dimensions
840      id = iom_varid( inum, sd%clvar, ddims )
841
842      !! close it
843      CALL iom_close( inum )
844
845      !! now open the weights file
846
847      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
848      IF ( inum > 0 ) THEN
849
850         !! determine whether we have an east-west cyclic grid
851         !! from global attribute called "ew_wrap" in the weights file
852         !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
853         !! since this is the most common forcing configuration
854
855         CALL iom_getatt(inum, 'ew_wrap', zwrap)
856         IF( zwrap >= 0 ) THEN
857            cyclical = .TRUE.
858         ELSE IF( zwrap == -999 ) THEN
859            cyclical = .TRUE.
860            zwrap = 0
861         ELSE
862            cyclical = .FALSE.
863         ENDIF
864
865         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
866         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
867         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
868         ref_wgts(nxt_wgt)%overlap = zwrap
869         ref_wgts(nxt_wgt)%cyclic = cyclical
870         ref_wgts(nxt_wgt)%nestid = 0
871#if defined key_agrif
872         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
873#endif
874         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
875         !! for each weight wgtNN there is an integer array srcNN which gives the point in
876         !! the input data grid which is to be multiplied by the weight
877         !! they are both arrays on the model grid so the result of the multiplication is
878         !! added into an output array on the model grid as a running sum
879
880         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
881         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
882         IF( id <= 0) THEN
883            ref_wgts(nxt_wgt)%numwgt = 4
884         ELSE
885            ref_wgts(nxt_wgt)%numwgt = 16
886         ENDIF
887
888         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
889         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
890         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
891
892         DO jn = 1,4
893            aname = ' '
894            WRITE(aname,'(a3,i2.2)') 'src',jn
895            data_tmp(:,:) = 0
896            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
897            data_src(:,:) = INT(data_tmp(:,:))
898            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
899            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
900         END DO
901
902         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
903            aname = ' '
904            WRITE(aname,'(a3,i2.2)') 'wgt',jn
905            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
906            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
907         END DO
908         CALL iom_close (inum)
909 
910         ! find min and max indices in grid
911         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
912         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
913         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
914         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
915
916         ! and therefore dimensions of the input box
917         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
918         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
919
920         ! shift indexing of source grid
921         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
922         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
923
924         ! create input grid, give it a halo to allow gradient calculations
925         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
926         ! a more robust solution will be given in next release
927         ipk =  SIZE(sd%fnow, 3)
928         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
929         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
930
931         nxt_wgt = nxt_wgt + 1
932
933      ELSE
934         CALL ctl_stop( '    fld_weight : unable to read the file ' )
935      ENDIF
936
937      IF(  wrk_not_released(2, 1) .OR.    &
938          iwrk_not_released(2, 1)  )   CALL ctl_stop('fld_weight: failed to release workspace arrays')
939      !
940   END SUBROUTINE fld_weight
941
942
943   SUBROUTINE fld_interp( num, clvar, kw, kk, dta, nrec )
944      !!---------------------------------------------------------------------
945      !!                    ***  ROUTINE fld_interp  ***
946      !!
947      !! ** Purpose :   apply weights to input gridded data to create data
948      !!                on model grid
949      !!----------------------------------------------------------------------
950      INTEGER                   , INTENT(in   ) ::   num     ! stream number
951      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
952      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
953      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
954      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
955      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
956      !!
957      INTEGER, DIMENSION(3) ::   rec1,recn   ! temporary arrays for start and length
958      INTEGER ::  jk, jn, jm           ! loop counters
959      INTEGER ::  ni, nj               ! lengths
960      INTEGER ::  jpimin,jpiwid        ! temporary indices
961      INTEGER ::  jpjmin,jpjwid        ! temporary indices
962      INTEGER ::  jpi1,jpi2,jpj1,jpj2  ! temporary indices
963      !!----------------------------------------------------------------------
964      !
965      !! for weighted interpolation we have weights at four corners of a box surrounding
966      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
967      !! or by a grid value and gradients at the corner point (bicubic case)
968      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
969
970      !! sub grid from non-model input grid which encloses all grid points in this nemo process
971      jpimin = ref_wgts(kw)%botleft(1)
972      jpjmin = ref_wgts(kw)%botleft(2)
973      jpiwid = ref_wgts(kw)%jpiwgt
974      jpjwid = ref_wgts(kw)%jpjwgt
975
976      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
977      rec1(1) = MAX( jpimin-1, 1 )
978      rec1(2) = MAX( jpjmin-1, 1 )
979      rec1(3) = 1
980      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
981      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
982      recn(3) = kk
983
984      !! where we need to put it in the non-nemo grid fly_dta
985      !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
986      !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
987      jpi1 = 2 + rec1(1) - jpimin
988      jpj1 = 2 + rec1(2) - jpjmin
989      jpi2 = jpi1 + recn(1) - 1
990      jpj2 = jpj1 + recn(2) - 1
991
992      ref_wgts(kw)%fly_dta(:,:,:) = 0.0
993      SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
994      CASE(1)
995           CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
996      CASE DEFAULT
997           CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
998      END SELECT 
999
1000      !! first four weights common to both bilinear and bicubic
1001      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1002      !! note that we have to offset by 1 into fly_dta array because of halo
1003      dta(:,:,:) = 0.0
1004      DO jk = 1,4
1005        DO jn = 1, jpj
1006          DO jm = 1,jpi
1007            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1008            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1009            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1010          END DO
1011        END DO
1012      END DO
1013
1014      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1015
1016        !! fix up halo points that we couldnt read from file
1017        IF( jpi1 == 2 ) THEN
1018           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1019        ENDIF
1020        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1021           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1022        ENDIF
1023        IF( jpj1 == 2 ) THEN
1024           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1025        ENDIF
1026        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1027           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1028        ENDIF
1029
1030        !! if data grid is cyclic we can do better on east-west edges
1031        !! but have to allow for whether first and last columns are coincident
1032        IF( ref_wgts(kw)%cyclic ) THEN
1033           rec1(2) = MAX( jpjmin-1, 1 )
1034           recn(1) = 1
1035           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1036           jpj1 = 2 + rec1(2) - jpjmin
1037           jpj2 = jpj1 + recn(2) - 1
1038           IF( jpi1 == 2 ) THEN
1039              rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1040              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1041              CASE(1)
1042                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1043              CASE DEFAULT
1044                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1045              END SELECT     
1046              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1047           ENDIF
1048           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1049              rec1(1) = 1 + ref_wgts(kw)%overlap
1050              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1051              CASE(1)
1052                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1053              CASE DEFAULT
1054                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1055              END SELECT
1056              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1057           ENDIF
1058        ENDIF
1059
1060        ! gradient in the i direction
1061        DO jk = 1,4
1062          DO jn = 1, jpj
1063            DO jm = 1,jpi
1064              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1065              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1066              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1067                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1068            END DO
1069          END DO
1070        END DO
1071
1072        ! gradient in the j direction
1073        DO jk = 1,4
1074          DO jn = 1, jpj
1075            DO jm = 1,jpi
1076              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1077              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1078              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         &
1079                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1080            END DO
1081          END DO
1082        END DO
1083
1084         ! gradient in the ij direction
1085         DO jk = 1,4
1086            DO jn = 1, jpj
1087               DO jm = 1,jpi
1088                  ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1089                  nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1090                  dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1091                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1092                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
1093               END DO
1094            END DO
1095         END DO
1096         !
1097      END IF
1098      !
1099   END SUBROUTINE fld_interp
1100
1101
1102   FUNCTION ksec_week( cdday )
1103      !!---------------------------------------------------------------------
1104      !!                    ***  FUNCTION kshift_week ***
1105      !!
1106      !! ** Purpose : 
1107      !!---------------------------------------------------------------------
1108      CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file
1109      !!
1110      INTEGER                        ::   ksec_week  ! output variable
1111      INTEGER                        ::   ijul       !temp variable
1112      INTEGER                        ::   ishift     !temp variable
1113      CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
1114      !!----------------------------------------------------------------------
1115      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1116      DO ijul = 1, 7
1117         IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1118      END DO
1119      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1120      !
1121      ishift = ijul * NINT(rday)
1122      !
1123      ksec_week = nsec_week + ishift
1124      ksec_week = MOD( ksec_week, 7*NINT(rday) )
1125      !
1126   END FUNCTION ksec_week
1127
1128   !!======================================================================
1129END MODULE fldread
Note: See TracBrowser for help on using the repository browser.