New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
fldread.F90 in branches/devmercator2010/NEMO/OPA_SRC/SBC – NEMO

source: branches/devmercator2010/NEMO/OPA_SRC/SBC/fldread.F90 @ 2076

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

add change dev_1784_EVP

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