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