source: TOOLS/MOZAIC/src/MOZAIC/interpol.f90 @ 3326

Last change on this file since 3326 was 3326, checked in by omamce, 7 years ago

O.M. : Utility to generate interpolatio weights for OASIS-MCT
File size: 14.5 KB
Line 
1! -*- Mode: f90 -*-
2!!!
3!!! lecture d'un fichier NetCDF et interpolation type Mosaic
4!!! Suppose un calcul préalable des poids
5!!!
6!!! Modifiez le module 'modele' selon votre cas de figure
7!!!
8MODULE declare_1
9   !!
10   USE defprec
11   IMPLICIT NONE
12   !! Data types
13   INTEGER, PARAMETER :: r4 = r_4, r8 = r_8
14   INTEGER, PARAMETER :: i1=1, i2 = 2, i4 = i_4, i8 = i_8
15   INTEGER, PARAMETER :: rl = r_std, il = i_std
16   INTEGER (kind=il), PARAMETER :: nout = 6 ! Fichier de sortie des messages
17   !!
18END MODULE declare_1
19!!
20MODULE modele_1
21  USE declare_1
22  IMPLICIT NONE
23  CHARACTER (len=80) :: cfile_in = 'routing.nc' ! Fichier entree
24  CHARACTER (len=80) :: cfile_ou = 'runoff_oce.nc' ! Fichier sortie
25  ! OPA -> LMD
26!$$$  INTEGER (kind=il), PARAMETER :: jpi1 = 92, jpj1 = 76, jpk1 = 1      ! Dimension modele entree
27!$$$  INTEGER (kind=il), PARAMETER :: jpi2 = 72, jpj2 = 46, jpk2 = 1      ! Dimension modele sortie
28!$$$  CHARACTER (len=8) :: clweight =  "WEIGHTS1" , cladress = "ADRESSE1" ! Reference des poids et adresses à utiliser
29!$$$  INTEGER (kind=il), PARAMETER :: jpn   = 15_il ! Nombre maxi de voisins
30  ! LMD -> OPA
31  INTEGER (kind=il), PARAMETER :: jpi1 = 72_il, jpj1 = 46_il, jpk1 = 1_il ! Dimension modele entree
32  INTEGER (kind=il), PARAMETER :: jpi2 = 92_il, jpj2 = 76_il, jpk2 = 1_il ! Dimension modele sortie
33  CHARACTER (len=8) :: clweight =  "WEIGHTS4" , cladress = "ADRESSE4" ! Reference des poids et adresses à utiliser
34  INTEGER (kind=il), PARAMETER :: jpn   = 6_il ! Nombre maxi de voisins
35  CHARACTER (len=80) :: cfile_ma = '/home/users/marti/GRAF/DATA/orca4.nc' ! Fichier geo sortie
36  !
37  INTEGER (kind=il), PARAMETER :: jpdim =  4_il ! Nombre maxi de dimensions
38  INTEGER (kind=il), PARAMETER :: jpvar = 30_il ! Nombre maxi de variables
39  CHARACTER (len=8) :: ctaux = 'FRTU', ctauy = 'FRTV' ! Nom des composantes du stress
40  REAL (kind=r8) :: zfreeze = -1.77_r8 ! Point de congélation
41  !
42  REAL    (kind=rl), DIMENSION ( jpn,jpi2*jpj2) :: weight
43  INTEGER (kind=il), DIMENSION ( jpn,jpi2*jpj2) :: kadd
44  REAL    (kind=rl), DIMENSION ( jpi2*jpj2) :: glon, glat
45END MODULE modele_1
46!!
47PROGRAM interpol
48   USE declare_1
49   USE modele_1
50   USE fliocom
51   USE errioipsl
52   IMPLICIT NONE
53   INTEGER (kind=i4),PARAMETER :: in = i4
54   !!
55   REAL (kind=r8), DIMENSION (jpi1*jpj1) :: zvar_in8
56   REAL (kind=r8), DIMENSION (jpi2*jpj2) :: zvar_ou8
57   INTEGER (kind=i4), DIMENSION (flio_max_dims) :: dim_len1, dim_len2
58   CHARACTER (len=20), DIMENSION (flio_max_dims) :: cldim
59   INTEGER (kind=i4) :: jdim, jvar
60   INTEGER (kind=i4), DIMENSION (jpvar) :: xtype, ndims, natts
61   REAL (kind=r8), DIMENSION (jpvar) :: spval = 0.0_r8
62   CHARACTER (len=20), DIMENSION (jpvar) :: clvar
63   INTEGER (kind=i4), DIMENSION (jpvar, flio_max_dims) :: dimsid
64   INTEGER (kind=in) :: jd, jv, jn, jo 
65   INTEGER (kind=il) :: ji, jj, jt, jnt, tid, xid, yid, tid_v, kerr
66   INTEGER (kind=i4) :: ncid_in, ncid_ou, ncid_ma, ja
67   INTEGER (kind=il) :: len1, len2, nwei = 13
68   CHARACTER (len=8) :: c8
69   INTEGER (kind=i4), DIMENSION ( jpdim) :: istart1, istart2, icount1, icount2
70   CHARACTER (len=80) :: cname
71   REAL (kind=r8), DIMENSION ( :), ALLOCATABLE :: zt
72   INTEGER (kind=il)       , DIMENSION ( :), ALLOCATABLE :: it
73   CHARACTER (len=80) :: cfname
74   ! For flio
75   CHARACTER (len=80), DIMENSION(1) :: catt
76   LOGICAL :: l_ex
77   !!
78   !! Lecture poids et adresses
79   WRITE ( unit = cfname, fmt = '("/home/climdata/ocean/CPL/ORCA4xLMD7245/mozaic.w.runoff.i", 1I1, ".r", 1I1)' ) il, rl
80   OPEN ( unit = nwei, file = TRIM(cfname), form = 'unformatted', status = 'old', action = 'read')
81   ALLOCATE ( it ( jpi2*jpj2*jpn) )
82   CALL locrint (cladress, it, jpi2*jpj2*jpn, nwei, kerr, nout) ;  IF ( kerr /= 0 ) STOP 055
83   DO jo = 1, jpi2*jpj2
84      DO jn = 1, jpn
85         kadd  (jn, jo) = it ( jn + (jo-1)*jpn)
86      END DO
87   ENDDO
88   DEALLOCATE ( it)
89   !
90   ALLOCATE ( zt ( jpi2*jpj2*jpn) )
91   CALL locread (clweight, zt, jpi2*jpj2*jpn, nwei, kerr, nout) ;  IF ( kerr /= 0 ) STOP 056
92   DO jo = 1, jpi2*jpj2
93      DO jn = 1, jpn
94         weight (jn, jo) = zt (jn + (jo-1)*jpn)
95      END DO
96   ENDDO
97   DEALLOCATE ( zt)
98   !! Lecture geo
99   CALL flioopfd (TRIM(cfile_ma), ncid_ma)
100   CALL fliogetv (ncid_ma, 'nav_lon', glon)
101   CALL fliogetv (ncid_ma, 'nav_lat', glat)
102   CALL flioclo  (ncid_ma)
103   !! Ouverture fichier entree
104   CALL flioopfd (TRIM(cfile_in), ncid_in, nb_dim=jdim, nb_gat=natts(1))
105   !! Dimensions
106   CALL flioinqf (ncid_in, ln_dim=dim_len1(1:jdim), nb_var=jvar )
107   CALL flioinqn (ncid_in, cn_dim=cldim(1:jdim), cn_var=clvar(1:jvar))
108   SeekDim: DO jd = 1, jdim
109      WRITE ( unit = c8, fmt = '(1I8)' ) jd
110      IF      ( dim_len1 (jd) == jpi1 ) THEN
111         dim_len2 (jd) = jpi2
112      ELSE IF ( dim_len1 (jd) == jpj1 ) THEN
113         dim_len2 (jd) = jpj2
114      ELSE IF ( dim_len1 (jd) == jpk1 ) THEN
115         dim_len2 (jd) = jpk2
116      ELSE
117         dim_len2 (jd) = dim_len1 (jd)
118      END IF
119      ! Determination du nombre de pas de temps
120      jnt = 0
121      WRITE(unit=c8,fmt='(2I4)') jd
122      IF ( INDEX ( cldim(jd), 'time') /= 0 .OR. INDEX ( cldim(jd), 'TIME') /= 0 .OR. INDEX ( cldim(jd), 'Time') /= 0) THEN
123         jnt = dim_len1(jd) ; tid = jd ; dim_len2 (jd) = -1 ! Unlimited dimension
124      END IF
125      IF ( INDEX ( cldim(jd), 'X') /= 0 .OR. INDEX ( cldim(jd), 'x') /= 0 ) xid = jd 
126      IF ( INDEX ( cldim(jd), 'Y') /= 0 .OR. INDEX ( cldim(jd), 'y') /= 0 ) yid = jd 
127   END DO SeekDim
128   WRITE ( nout, *) 'Nombre de dimensions : ', jdim
129   WRITE ( nout, *) 'Longueur temps       : ', jnt
130   WRITE ( nout, *) 'Dimension temps : ', tid, cldim(tid)
131   !! Ouverture fichier sortie
132   CALL fliocrfd (TRIM(cfile_ou), cldim(1:jdim), dim_len2(1:jdim), ncid_in )
133   CALL flioputa (ncid_in, '?', 'Comment', TRIM(c_comment) )
134   !! Attributs globaux
135   DO ja = 1, natts (1) 
136      WRITE ( unit = c8, fmt = '(2I4)' ) ja
137      CALL flioinqn (ncid_in, ia_start=ja, ia_count=1, cn_gat=catt) ; cname=catt(1)
138      CALL fliocpya (ncid_in, '?', cname, ncid_ou, '?' )
139   END DO
140   !! Definition variables
141   SeekVar: DO jv = 1, jvar
142      CALL flioinqv (ncid_in, TRIM(clvar(jv)), l_ex, nb_atts=natts(jv), v_t=xtype(jv), nb_dims=ndims(jv))
143      CALL flioinqv (ncid_in, TRIM(clvar(jv)), l_ex, id_dims=dimsid(jv,1:ndims(jv)))
144      WRITE ( nout, FMT='(1A20,8I5)' ) TRIM(clvar(jv)), xtype(jv), dimsid(jv,1:ndims(jv)), natts(jv), jv
145      IF ( INDEX ( clvar(jv), 'time') /= 0 .OR. INDEX ( clvar(jv), 'TIME') /= 0 .OR. &
146         & INDEX ( clvar(jv), 'Time') /= 0) tid_v = jv
147      CALL fliodefv (ncid_ou, TRIM(clvar(jv)), dimsid(jv,1:ndims(jv)), v_t=xtype(jv))
148      !! Traitements des attributs
149      DO ja = 1, natts(jv)
150         WRITE ( unit = c8, fmt = '(2I4)' ) jv, ja
151         CALL flioinqv (ncid_in, TRIM(clvar(jv)), l_ex, ia_start=ja, ia_count=1, cn_atts=catt) ; cname=catt(1)
152         CALL fliocpya (ncid_in, TRIM(clvar(jv)), cname, ncid_ou, TRIM(clvar(jv)) )
153      END DO
154      IF ( TRIM(cname) == 'missing_value' .OR. TRIM(cname) == 'MISSING_VALUE' ) THEN
155         CALL fliogeta ( ncid_in, TRIM(clvar(jv)), cname, spval(jv))
156         WRITE ( unit = nout, fmt = *) 'Missing value ' // TRIM (clvar(jv)) // " : ", spval(jv)
157      ENDIF
158   END DO SeekVar
159   !!
160   Temps: DO jt = 1_il, 50 !jnt
161      Var: DO jv = 1_il, jvar
162         WRITE(unit=c8,fmt='(2I4)') jt, jv
163         istart1 = 1_il ; icount1 = 1_il
164         istart2 = 1_il ; icount2 = 1_il
165         SELECT CASE ( ndims(jv) )
166         CASE (1)
167            IF      ( jv == tid_v ) THEN
168               istart1 (1_il)= jt ; icount1 (1_il) = 1_il
169               istart2 (1_il)= jt ; icount2 (1_il) = 1_il
170            ELSE IF ( jv == xid ) THEN
171               istart1 (1_il) = 1_il ; icount1 (1_il) = jpi1
172               istart2 (1_il) = 1_il ; icount2 (1_il) = jpi2
173               CALL flioputv (ncid_ou, TRIM(clvar(jv)),  (/ (REAL(ji,kind=r8),ji=1,jpi2) /) ) ! Ecrit Y
174               CYCLE var
175            ELSE IF ( jv == yid ) THEN
176               istart1 (1_il) = 1_il ; icount1 (1) = jpj1
177               istart2 (1_il) = 1_il ; icount2 (1) = jpj2
178               CALL flioputv (ncid_ou, TRIM(clvar(jv)),  (/ (REAL(jj,kind=r8),jj=1,jpj2) /) ) ! Ecrit Y
179               CYCLE var
180            ELSE
181               IF ( jt > 1_il) CYCLE var
182            END IF
183         CASE ( 2) 
184            IF ( jt > 1 ) CYCLE var
185            istart1 (1:2) = (/  1_il, 1_il /) ; icount1 (1:2) = (/jpi1, jpj1 /)
186            istart2 (1:2) = (/  1_il, 1_il /) ; icount2 (1:2) = (/jpi2, jpj2 /)
187         CASE ( 3) 
188            istart1 (1:3) = (/  1_il, 1_il, jt /)
189            icount1 (1:3) = (/jpi1, jpj1, 1_il /)
190            istart2 (1:3) = (/  1_il, 1_il, jt /)
191            icount2 (1:3) = (/jpi2, jpj2, 1_il /)
192         CASE ( 4) 
193            istart1 (1:4) = (/  1_il, 1_il, 1_il, jt /) ; icount1 (1:4)  = (/jpi1, jpj1, 1_il, 1_il /)
194            istart2 (1:4) = (/  1_il, 1_il, 1_il, jt /) ; icount2 (1:4)  = (/jpi2, jpj2, 1_il, 1_il /)
195         END SELECT
196         !!
197         len1 = product ( icount1 ) ; len2 = product ( icount2 )
198         WRITE ( nout, FMT = '("No ", 1I3, ", Nom : ", 1A20, ", Time : ", 1I3, ", Lens : ", 2I6, ", Type : ", 1I4)' ) &
199            &  jv, TRIM(clvar(jv)), jt, len1, len2, xtype(jv)
200         !
201         ! Lecture
202         CALL fliogetv ( ncid_in, TRIM(clvar(jv)), zvar_in8(1:len1), istart1, icount1)
203         zvar_in8(1:len1) = zvar_in8(1:len1) + 0.0_r8
204         WHERE ( zvar_in8(1:len1) == spval ( jv) ) zvar_in8(1:len1) = 0.0_r8
205!$$$      WRITE (nout,*) 'extrema entree : ', MINVAL ( zvar_in8(1:len1)) , MAXVAL ( zvar_in8(1:len1))
206         ! Interpolation eventuelle
207         IF ( ndims(jv) == 1 ) THEN
208            zvar_ou8(1:len1) = zvar_in8(1:len1)
209         ELSE
210!$$$          IF ( jv == 6 ) THEN
211!$$$              WRITE(nout,*) 'Attention : traitement glace '
212!$$$              DO ji = 1, len1
213!$$$                IF ( zvar_in8(ji) <= zfreeze ) THEN
214!$$$                     zvar_in8(ji) = 1.0_rl
215!$$$                 ELSE
216!$$$                     zvar_in8(ji) = 0.0_rl
217!$$$                 END IF
218!$$$               END DO
219!$$$          ENDIF
220            IF      ( TRIM(clvar(jv)) == 'nav_lon' .OR. TRIM(clvar(jv)) == 'longitude' .OR. & 
221               &      TRIM(clvar(jv)) == 'NAV_LON' .OR. TRIM(clvar(jv)) == 'LONGITUDE' ) THEN
222               zvar_ou8(1:len2) = glon
223            ELSE IF ( TRIM(clvar(jv)) == 'nav_lat' .OR. TRIM(clvar(jv)) == 'latitude' .OR. &
224               &      TRIM(clvar(jv)) == 'NAV_LAT' .OR. TRIM(clvar(jv)) == 'LATITUDE' ) THEN
225               zvar_ou8(1:len2) = glat
226            ELSE
227               CALL inter ( zvar_in8(1:len1), zvar_ou8(1:len2))
228            ENDIF
229         END IF
230!$$$      WRITE (nout,*) 'extrema sortie : ', MINVAL ( zvar_ou8(1:len2)) , MAXVAL ( zvar_ou8(1:len2))
231         WRITE (nout,FMT='("extrema : ",4E12.3)' ) MINVAL ( zvar_in8(1:len1)), MINVAL ( zvar_ou8(1:len2)) &
232            &  , MAXVAL ( zvar_in8(1:len1)), MAXVAL ( zvar_ou8(1:len2))
233         ! Ecriture
234         CALL flioputv ( ncid_ou, TRIM(clvar(jv)), zvar_ou8(1:len2))
235      END DO Var
236   END DO Temps
237
238   !!
239   CALL flioclo ( ncid_in)
240   CALL flioclo ( ncid_ou)
241
242   STOP
243CONTAINS
244   SUBROUTINE inter ( pin, pout)
245      USE declare_1
246      USE modele_1
247      IMPLICIT NONE
248      REAL (kind=r8), INTENT ( in), DIMENSION ( jpi1*jpj1) :: pin
249      REAL (kind=r8), INTENT (out), DIMENSION ( jpi2*jpj2) :: pout
250      INTEGER (kind=il) :: ji, jn, in
251      !!
252      pout = 0.0_rl
253      !
254      DO ji = 1, jpi2*jpj2
255         DO jn = 1, jpn
256            in  = kadd ( jn, ji)
257            IF ( in > 0 ) pout ( ji) = pout ( ji) + weight ( jn, ji) * pin ( in )
258         END DO
259      END DO
260      !!
261   END SUBROUTINE inter
262   !!
263   SUBROUTINE locread ( cdfldn, pfield, kdimax, knulre, kflgre, kout)
264      USE declare_1
265      !!
266      IMPLICIT NONE
267      !
268      CHARACTER (LEN = 8), INTENT ( in) :: cdfldn ! Name of field to search for
269      INTEGER (kind=il), INTENT (  in) :: kdimax, knulre    ! Field dimension, unit to read
270      INTEGER (kind=il), INTENT (  in) :: kout              ! Standard output
271      INTEGER (kind=il), INTENT ( out) :: kflgre            ! Error status code
272      REAL (kind=r8)  , INTENT ( out), DIMENSION ( :) :: pfield ! Field
273      !
274      CHARACTER (len = 8) clecfl
275      INTEGER (kind=il) :: ji, kerr
276      !
277      !*    1. Initialization
278      !        --------------
279      !
280      WRITE (UNIT = kout,FMT = '( "Locread :  Read field ", A8, " in unit = ", 1I6)' ) cdfldn, knulre
281      !
282      !     2. Find field in file
283      !        ------------------
284      !
285      REWIND ( UNIT = knulre)
286200   CONTINUE
287      !* Find string
288      READ (UNIT = knulre, ERR = 200, END = 210, IOSTAT = kerr) clecfl
289      IF ( clecfl /= cdfldn ) GO TO  200
290      !* Read associated field
291      READ (UNIT = knulre, IOSTAT = kerr) ( pfield (ji),ji = 1, kdimax)
292      IF ( kerr /= 0 ) WRITE ( unit = nout, fmt = *) 'Problem reading field'
293      kflgre = 0
294      GO TO 220
295      !* Problem in reading
296210   kflgre = 1
297      WRITE ( unit = nout, fmt = *) 'Problem in reading'
298220   CONTINUE
299      !
300      !*    3. End of routine
301      !        --------------
302      !
303      RETURN
304   END SUBROUTINE locread
305   !
306   SUBROUTINE locrint ( cdfldn, kfield, kdimax, knulre, kflgre, kout)
307      !
308      USE declare_1
309      !!
310      IMPLICIT none
311      !
312      CHARACTER (LEN = 8), INTENT ( in) :: cdfldn ! Name of field to search for
313      INTEGER (kind=il), INTENT (  in) :: kdimax, knulre    ! Field dimension, unit to read
314      INTEGER (kind=il), INTENT (  in) :: kout              ! Standard output
315      INTEGER (kind=il), INTENT ( out) :: kflgre            ! Error status code
316      INTEGER (kind=il), INTENT ( inout), DIMENSION (:) :: kfield ! Field
317      !
318      CHARACTER (len = 8) :: clecfl
319      INTEGER (kind=il) :: ji, kerr
320      !
321      !*    1. Initialization
322      !        --------------
323      !
324      WRITE (UNIT = kout, FMT = '( "Locrint :  Read field ", A8, " in unit = ", 1I6, " Dim : ", 1I9)' ) cdfldn, knulre, kdimax
325      !
326      !     2. Find field in file
327      !        ------------------
328      !
329      REWIND ( UNIT = knulre )
330200   CONTINUE
331      !
332      !* Find string
333      READ (UNIT = knulre, ERR = 200, END = 210) clecfl
334!$$$    WRITE ( unit = kout, fmt = *) clecfl
335      !
336      IF ( clecfl /= cdfldn ) GO TO  200
337      !* Read associated field
338      READ (UNIT = knulre, IOSTAT = kerr) ( kfield (ji),ji = 1, kdimax)
339      IF ( kerr /= 0 ) WRITE ( unit = nout, fmt = *) 'Problem reading field'
340      !* Reading done and ok
341      kflgre = 0
342      GO TO 220
343      !* Problem in reading
344210   kflgre = 1
345      WRITE ( unit = nout, fmt = *) 'Problem in reading. ji : ', ji, kfield(1), kfield(kdimax)
346      !
347220   CONTINUE
348      !
349      !*    3. End of routine
350      !        --------------
351      !
352      RETURN
353   END SUBROUTINE locrint
354   !
355END PROGRAM interpol
Note: See TracBrowser for help on using the repository browser.