source: CPL/oasis3/trunk/src/lib/anaism/src/namset.f @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago

Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).

File size: 16.2 KB
Line 
1      SUBROUTINE namset (px1, py1, kmsk1, kvmsk1, kngx1, kngy1, 
2     $                   cdper1, kper1,
3     $                   px2, py2, kmsk2, kvmsk2, kngx2, kngy2, 
4     $                   cdper2, kper2,
5     $                   cdtyp,
6     $                   pr1to2, k1to2, kw1to2, psgr1, psgr2,
7     $                   krdwt, kwlun, kvmsz2,
8     $                   kmskz2, knumber)
9C****
10C               *****************************
11C               * OASIS ROUTINE  -  LEVEL 3 *
12C               * -------------     ------- *
13C               *****************************
14C
15C**** *namset* - Initializations for the mesh naive method
16C
17C     Purpose:
18C     -------
19C     Read or calculate the interpolation weights used in Anaism
20C
21C**   Interface:
22C     ---------
23C       *CALL*  *namset(px1, py1, kmsk1, kvmsk1, kngx1, kngy1, cdper1, kper1,
24C                       px2, py2, kmsk2, kvmsk2, kngx2, kngy2, cdper2, kper2,
25C                       pr1to2, k1to2, kw1to2, psgr1, psgr2,
26C                       krdwt, kwlun, kvmsz2,
27C                       kmskz2, knumber)*
28C     Input:
29C     -----
30C                kngx1   : number of longitudes for source grid
31C                kngy1   : number of latitudes for source grid
32C                px1     : longitudes for source grid (real 2D)
33C                py1     : latitudes for source grid (real 2D)
34C                kmsk1   : the mask for source grid (integer 2D)
35C                kvmsk1  : the value of the mask for source grid
36C                cdper1  : source grid periodicity
37C                kper1   : number of overlapped points for source grid 
38C                kngx2   : number of longitudes for target grid
39C                kngy2   : number of latitudes for target grid
40C                px2     : longitudes for target grid (real 2D)
41C                py2     : latitudes for target grid (real 2D)
42C                kmsk2   : the mask of target grid (integer 2D)
43C                kvmsk2  : the value of the mask for target grid 
44C                cdper2  : target grid periodicity
45C                kper2   : number of overlapped points for target grid
46C                cdtyp   : type of source grid 
47C                kw1to2  : maximum number of overlapped neighbors
48C                krdwt   : read/write flag for the weights
49C                kwlun   : logical unit for the weights
50C                kvmsz2  : mask value for array kmskz2
51C                knumber : flag to identify appropriate Anaism dataset
52C
53C     Output:
54C     ------
55C                pr1to2  : weights for Anaism interpolation (real 2D)
56C                k1to2   : source grid neighbors adresses (integer 2D)
57C                kmskz2  : number of source grid neighbors (integer 2D)
58C                psgr1   : source grid square surfaces (real 2D)
59C                psgr2   : target grid square surfaces (real 2D)
60C                 
61C
62C     Workspace:
63C     ---------
64C     None
65C
66C     External:
67C     --------
68C     pcssph, ssumr, grstat, pmrhal, locwrint, locwrite, locread, locrint,
69C     icoor, jcoor
70C
71C     References:
72C     ----------
73C     O. Thual, Simple ocean-atmosphere interpolation. 
74C               Part A: The method, EPICOA 0629 (1992)
75C               Part B: Software implementation, EPICOA 0630 (1992)
76C     See also OASIS manual (1995)
77C
78C     History:
79C     -------
80C       Version   Programmer     Date      Description
81C       -------   ----------     ----      ----------- 
82C       1.1       O. Thual       93/04/15  created 
83C       2.0       L. Terray      95/10/01  modified: new structure
84C       2.3       S. Valcke      99/04/30  added: printing levels
85C       2.3       L. Terray      99/09/15  changed periodicity variables
86C
87C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88C
89C* ---------------------------- Include files ---------------------------
90C
91      USE mod_kinds_oasis
92      USE mod_unit
93      USE mod_printing
94C
95C* ---------------------------- Argument declarations -------------------
96C
97      REAL (kind=ip_realwp_p) px1(kngx1,kngy1), py1(kngx1,kngy1), 
98     $    psgr1(kngx1,kngy1)
99      REAL (kind=ip_realwp_p) px2(kngx2,kngy2), py2(kngx2,kngy2), 
100     $    psgr2(kngx2,kngy2)
101      REAL (kind=ip_realwp_p) pr1to2(kw1to2,kngx2*kngy2)
102      INTEGER (kind=ip_intwp_p) kmsk1(kngx1,kngy1), kmsk2(kngx2,kngy2)
103      INTEGER (kind=ip_intwp_p) k1to2(kw1to2,kngx2*kngy2), 
104     $    kmskz2(kngx2,kngy2)
105      CHARACTER*8 cdper1, cdper2
106      CHARACTER*1 cdtyp
107C
108C* ---------------------------- Local declarations ----------------------
109C
110      CHARACTER*8 cladress, clweight, cloverlp
111      CHARACTER*1 cltyp
112C
113C* ---------------------------- Poema verses ----------------------------
114C
115C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
116C
117C*    1. Initializations
118C        ---------------
119C
120      IF (nlogprt .GE. 2) THEN
121          WRITE (UNIT = nulou,FMT = *) ' '
122          WRITE (UNIT = nulou,FMT = *) ' '
123          WRITE (UNIT = nulou,FMT = *) 
124     $    '           Routine namset  -  Level 3'
125          WRITE (UNIT = nulou,FMT = *) 
126     $    '           **************     *******'
127          WRITE (UNIT = nulou,FMT = *) ' '
128          WRITE (UNIT = nulou,FMT = *) 
129     $    ' Set up ANAIS-MESH interpolation'
130          WRITE (UNIT = nulou,FMT = *) ' '
131          WRITE (UNIT = nulou,FMT = *) ' '
132      ENDIF
133C
134C* Define global dimensions + other local variables
135C
136      ing1 = kngx1 * kngy1
137      ing2 = kngx2 * kngy2
138      ivmsk1 = kvmsk1
139      ivmsk2 = kvmsk2
140      iflag = 0
141C
142C* Define character strings to locate needed data
143C
144      WRITE(clweight,'(''WEIGHTS'',I1)') knumber
145      WRITE(cladress,'(''ADRESSE'',I1)') knumber
146      WRITE(cloverlp,'(''OVERLAP'',I1)') knumber
147C
148C
149C*    2. Statistics of source grid
150C        -------------------------
151C
152      cltyp = cdtyp
153C* The following routines calculate some interesting info about the grids
154C
155C* Calculate surface elements (assume spherical and periodic grid)
156C
157      CALL pcssph (px1, py1, psgr1, kngx1, kngy1)
158      zsum = ssumr (psgr1, ing1)
159C
160C* Printing on ANAIS output file
161C
162      IF (nlogprt .GE. 2) THEN
163          WRITE (UNIT = nulan,FMT = *) '           ANAIS output FILE '
164          WRITE (UNIT = nulan,FMT = *) '           ***** ****** **** '
165          WRITE (UNIT = nulan,FMT = *) ' '
166          WRITE (UNIT = nulan,FMT = *) '            Routine namset '
167          WRITE (UNIT = nulan,FMT = *) '            -------------- '
168          WRITE (UNIT = nulan,FMT = *) ' '
169          WRITE (UNIT = nulan,FMT = *) 
170     $    '      Some statistics about the source grid '
171          WRITE (UNIT = nulan,FMT = *) 
172     $    '      ************************************* '
173          WRITE (UNIT = nulan,FMT = *) ' '
174          WRITE (UNIT = nulan,FMT = *) 
175     $    ' Sum of grid square surfaces = ', zsum
176          WRITE (UNIT = nulan,FMT = *) 
177     $    ' It must be equal to 4 x PI  = ', 16. * atan(1.) 
178          WRITE (UNIT = nulan,FMT = *) ' '
179      ENDIF
180C
181C* Calculate for ocean points (unmasked points) the following data :
182C                ************
183C           - total surface
184C           - average inter-point distance
185C           - inverse of the sum of surface elements squared
186C
187      CALL grstat (px1, py1, psgr1, kmsk1, kngx1, kngy1, cltyp,
188     $             zamsh1, zstot1, zhhi1, ivmsk1) 
189C 
190C* Printing
191C 
192      IF (nlogprt .GE. 2) THEN
193          WRITE (UNIT = nulan,FMT = *) 
194     $    ' Some info for ocean points only '
195          WRITE (UNIT = nulan,FMT = *) 
196     $    ' ******************************* '
197          WRITE (UNIT = nulan,FMT = *) ' '
198          WRITE (UNIT = nulan,FMT = *) 
199     $    ' Average mesh distance = ', sqrt(zamsh1)
200          WRITE (UNIT = nulan,FMT = *) ' Total surface = ', zstot1
201          WRITE (UNIT = nulan,FMT = *) 
202     $    ' Inverse of sum surf. elts. squared = ', zhhi1
203      ENDIF
204C
205C
206C*    3. Statistics of target grid
207C        -------------------------
208C
209C* The following routines calculate some interesting info about the grids
210C
211C* Calculate surface elements (assume spherical and periodic grid)
212C
213      CALL pcssph (px2, py2, psgr2, kngx2, kngy2)
214      zsum = ssumr (psgr2, ing2)
215C
216C* Printing
217C
218      IF (nlogprt .GE. 2) THEN 
219          WRITE (UNIT = nulan,FMT = *) ' '
220          WRITE (UNIT = nulan,FMT = *) 
221     $    '      Some statistics about the target grid '
222          WRITE (UNIT = nulan,FMT = *) 
223     $    '      ************************************* '
224          WRITE (UNIT = nulan,FMT = *) ' '
225          WRITE (UNIT = nulan,FMT = *) 
226     $    ' Sum of grid square surfaces = ', zsum
227          WRITE (UNIT = nulan,FMT = *) 
228     $    ' It must be equal to 4 x PI  = ', 
229     $                        16. * atan(1.) 
230          WRITE (UNIT = nulan,FMT = *) ' '
231      ENDIF
232C
233C* Calculate for ocean points (unmasked points) the following data :
234C                ************
235C           - total surface
236C           - average inter-point distance
237C           - inverse of the sum of surface elements squared
238C
239      CALL grstat (px2, py2, psgr2, kmsk2, kngx2, kngy2, cltyp,
240     $             zamsh2, zstot2, zhhi2, ivmsk2) 
241C 
242C* Printing
243C
244      IF (nlogprt .GE. 2) THEN   
245          WRITE (UNIT = nulan,FMT = *) 
246     $    ' Some info for ocean points only '
247          WRITE (UNIT = nulan,FMT = *) 
248     $    ' ******************************* '
249          WRITE (UNIT = nulan,FMT = *) ' '
250          WRITE (UNIT = nulan,FMT = *) 
251     $    ' Average mesh distance = ', sqrt(zamsh2)
252          WRITE (UNIT = nulan,FMT = *) ' Total surface = ', zstot2
253          WRITE (UNIT = nulan,FMT = *) 
254     $    ' Inverse of sum surf. elts. squared = ', zhhi2
255      ENDIF
256C 
257C
258C*    4. Preparation of the neighbors search
259C        -----------------------------------
260C
261C* Consistency checking
262C
263      IF (kw1to2 .GT. ing2) THEN
264          WRITE (UNIT = nulou,FMT = *) ' WARNING: kw1to2 .gt. ing2 '
265          WRITE (UNIT = nulou,FMT = *) ' *******  ------      ---- '
266          WRITE (UNIT = nulou,FMT = *) 
267     $        ' kw1to2 = ',kw1to2,' ing2 = ',ing2
268      ENDIF
269C
270C
271C*    5. Weights determination
272C        ---------------------
273C
274C* Initialization for read write options
275C 
276      IF (krdwt .EQ. 1) THEN
277C
278C* Writing
279C  -------
280C
281C* Calculate and write weights after checking numbers 
282C 
283C* Calculates weights 
284C
285          CALL pmrhal (pr1to2, k1to2, kw1to2,
286     $                 px1, py1, kmsk1, kngx1, kngy1, cdper1, kper1,
287     $                 px2, py2, kmsk2, kngx2, kngy2, cdper2, kper2,
288     $                 ivmsk1, ivmsk2, kmskz2, kvmsz2)
289C
290C* Write  weights + other useful stuff
291C
292C* - Write main locator string + weights
293C
294          CALL locwrite (clweight, pr1to2, kw1to2*kngx2*kngy2, 
295     $                   kwlun, iflag)
296          IF (iflag .NE. 0) THEN
297              WRITE (UNIT = nulou,FMT = *) 
298     $            'Problem in writing on UNIT = ', kwlun
299              WRITE (UNIT = nulou,FMT = *)
300     $            'String locator is = ', clweight
301              CALL HALTE ('Stop in namset')
302          ENDIF
303C 
304C* - Write overlapped neighbors adresses
305C
306          CALL locwrint (cladress, k1to2, kw1to2*kngx2*kngy2,
307     $                   kwlun, iflag)
308          IF (iflag .NE. 0) THEN
309              WRITE (UNIT = nulou,FMT = *) 
310     $            'Problem in writing on UNIT = ', kwlun
311              WRITE (UNIT = nulou,FMT = *)
312     $            'String locator is = ', cladress
313              CALL HALTE ('Stop in namset')
314          ENDIF
315C
316C* - Write overlapped neighbors numbers
317C
318          CALL locwrint (cloverlp, kmskz2, kngx2*kngy2,
319     $                   kwlun, iflag)
320          IF (iflag .NE. 0) THEN
321              WRITE (UNIT = nulou,FMT = *) 
322     $            'Problem in writing on UNIT = ', kwlun
323              WRITE (UNIT = nulou,FMT = *)
324     $            'String locator is = ', cloverlp
325              CALL HALTE ('Stop in namset')
326          ENDIF
327C
328C* - Write  checkings numbers (grid parameters)
329C
330          WRITE(kwlun) kngx1, kngy1, kngx2, kngy2
331C
332C* Job is done
333C
334          IF (nlogprt .GE. 2) THEN 
335              CALL prtout('Wrote weights on unit = ', kwlun, 1) 
336          ENDIF
337C
338C* Reading
339C  -------
340C
341        ELSE
342C
343C
344C* Reads weights and checks 
345C
346C* - Read main locator string + weights
347          CALL locread (clweight, pr1to2, kw1to2*kngx2*kngy2, 
348     $                   kwlun, iflag)
349          IF (iflag .NE. 0) THEN
350              WRITE (UNIT = nulou,FMT = *) 
351     $            'Problem in reading on UNIT = ', kwlun
352              WRITE (UNIT = nulou,FMT = *)
353     $            'String locator is = ', clweight
354              CALL HALTE ('Stop in namset')
355          ENDIF
356C
357C* - Read overlapped neighbors adresses
358C
359          CALL locrint (cladress, k1to2, kw1to2*kngx2*kngy2,
360     $                   kwlun, iflag)
361          IF (iflag .NE. 0) THEN
362              WRITE (UNIT = nulou,FMT = *) 
363     $            'Problem in reading on UNIT = ', kwlun
364              WRITE (UNIT = nulou,FMT = *)
365     $            'String locator is = ', cladress
366              CALL HALTE ('Stop in namset')
367          ENDIF
368C
369C* - Read overlapped neighbors numbers
370C
371          CALL locrint (cloverlp, kmskz2, kngx2*kngy2,
372     $                   kwlun, iflag)
373          IF (iflag .NE. 0) THEN
374              WRITE (UNIT = nulou,FMT = *) 
375     $            'Problem in reading on UNIT = ', kwlun
376              WRITE (UNIT = nulou,FMT = *)
377     $            'String locator is = ', cloverlp
378              CALL HALTE ('Stop in namset')
379          ENDIF
380C
381C* - Read  checkings numbers (grid parameters)
382C
383          READ(kwlun) ingx1, ingy1, ingx2, ingy2
384C
385C* Checks
386C
387          IF (ingx1 .NE. kngx1 .OR. ingy1 .NE. kngy1 .OR. 
388     $        ingx2 .NE. kngx2 .OR. ingy2 .NE. kngy2) THEN
389              WRITE (UNIT = nulou,FMT = *) 
390     $            ' Inconsistency in mweights file'
391              WRITE (UNIT = nulou,FMT = *) 
392     $            'ingx1 = ',ingx1,'kngx1 = ',kngx1
393              WRITE (UNIT = nulou,FMT = *) 
394     $            'ingy1 = ',ingy1,'kngy1 = ',kngy1
395              WRITE (UNIT = nulou,FMT = *) 
396     $            'ingx2 = ',ingx2,'kngx2 = ',kngx2
397              WRITE (UNIT = nulou,FMT = *) 
398     $            'ingy2 = ',ingy2,'kngy2 = ',kngy2
399              CALL HALTE ('Stop in namset')
400          ENDIF
401C
402C* Reading of weights done
403C
404          IF (nlogprt .GE. 2) THEN
405              CALL prtout
406     $        ('Reading of weights done on unit = ', kwlun, 1) 
407          ENDIF
408      ENDIF
409C
410C
411C*    6. Printing weights and adresses on ANAIS output file
412C        --------------------------------------------------
413C
414      IF (nlogprt .GE. 2) THEN
415          WRITE (UNIT = nulan,FMT = *) ' '
416          WRITE (UNIT = nulan,FMT = *) 
417     $    ' Print weights and adresses of overlapped neighbors '
418          WRITE (UNIT = nulan,FMT = *)
419     $    ' ************************************************** '
420          WRITE (UNIT = nulan,FMT = *) ' '
421          DO 610 jj = 1, kngy2 
422            DO 620 ji = 1, kngx2
423              WRITE (UNIT = nulan,FMT = 6010) ji, jj
424              WRITE (UNIT = nulan,FMT = 6020) px2(ji,jj), py2(ji,jj)
425              WRITE (UNIT = nulan,FMT = 6030) kmskz2(ji,jj)
426              IF (kmskz2(ji,jj) .NE. 0) THEN
427                  IF (kmsk2(ji,jj) .EQ. 0) THEN
428                      DO 630 jn = 1, kmskz2(ji,jj)
429                        iind = ji + kngx2 * (jj-1)
430                        iadr = k1to2(jn,iind)
431                        ix = icoor(iadr,kngx1)
432                        iy = jcoor(iadr,kngx1)
433                        WRITE(UNIT = nulan,FMT = 6040) 
434     $                  jn, pr1to2(jn,iind)
435                        WRITE(UNIT = nulan,FMT = 6050) 
436     $                  ix, iy, px1(ix,iy), py1(ix,iy)
437 630                  CONTINUE
438                  ELSE
439                      WRITE(UNIT = nulan,FMT = 6060)
440                  ENDIF
441              ENDIF
442 620        CONTINUE
443 610      CONTINUE
444          CALL FLUSH(nulan)
445      ENDIF
446C
447C* Formats
448C
449 6010 FORMAT(/,' Target grid point --  i = ',I3,' j = ',I3)
450 6020 FORMAT(' Point coordinates -- long = ',F9.4,' lat = ',F9.4)
451 6030 FORMAT(3X,'Number of overlapped source grid squares = ',I3)
452 6040 FORMAT(5X,' neighbor number = ',I3,' weight is = ',F6.4)
453 6050 FORMAT(5X,' i = ',I3,' j = ',I3,' lon = ',F9.4,' lat = ',F9.4)
454 6060 FORMAT(5X,' This is a continental point on the target grid ')
455C
456C
457C*    7. End of routine
458C        --------------
459C
460      IF (nlogprt .GE. 2) THEN
461          WRITE(UNIT = nulou,FMT = *) ' '
462          WRITE(UNIT = nulou,FMT = *) 
463     $    '          --------- End of routine namset ---------'
464          CALL FLUSH(nulou)
465      ENDIF
466      RETURN
467      END
468
469
470
471
Note: See TracBrowser for help on using the repository browser.