Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

extrap.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: 19.6 KB

Line
1	SUBROUTINE extrap (pfild, pmask, pwork,
2	$ kmask, kxlon, kylat,
3	$ knbor, cdextmet, cdper, kper,
4	$ krdwt, knb)
5	C****
6	C *****************************
7	C * OASIS ROUTINE - LEVEL 3 *
8	C * ------------- ------- *
9	C *****************************
10	C
11	C**** extrap - Extrapolation routine
12	C
13	C Purpose:
14	C -------
15	C Extrapolate field on land from sea values using nearest
16	C neighbors filling.
17	C
18	C** Interface:
19	C ---------
20	C CALL *extrap (pfild, pmask, pwork,
21	C $ kmask, kxlon, kylat,
22	C $ knbor, cdextmet, cdper, kper,
23	C $ krdwt, knb)
24	C
25	C Input:
26	C -----
27	C pfild : field to be extrapolated (real 2D)
28	C pmask : mask value (real)
29	C pwork : work array (real 2D)
30	C kmask : mask array (integer 2D)
31	C kxlon : number of longitudes (integer)
32	C kylat : number of latitudes (integer)
33	C knbor : nbr of neighbors for extrapolation (integer)
34	C cdextmet : extrapolation method (character)
35	C cdper : grid period.:P-period.,R-region. (character)
36	C kper : number of overlapping grid points (if any)
37	C krdwt : read/write flag for the weights and address
38	C knb : flag to identify appropriate NINENN dataset
39	C
40	C Output:
41	C ------
42	C pfild : field extrapolated
43	C
44	C Workspace:
45	C ---------
46	C zwork, zmask, ix, iy
47	C
48	C Externals:
49	C ---------
50	C None
51	C
52	C Reference:
53	C ---------
54	C See OASIS manual (1998)
55	C
56	C History:
57	C -------
58	C Version Programmer Date Description
59	C ------- ---------- ---- -----------
60	C 1.0 L. Terray 94/01/01 created
61	C 2.0 L. Terray 95/12/20 modified: new structure
62	C 2.2 L. Terray 97/12/31 modified: general cleaning
63	C 2.3 L. Terray 99/03/01 corrected: bug with
64	C latitudes fully masked
65	C 2.3 S. Valcke 99/03/30 add READ/WRITE for NINENN
66	C weights
67	C 2.3 S. Valcke 99/04/30 added: printing levels
68	C 2.3 S. Valcke 99/10/01 corrected: perodic
69	C overlapping grid
70	C
71	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72	C
73	C* ---------------- Include files and USE of modules---------------------------
74	C
75	USE mod_kinds_oasis
76	USE mod_parameter
77	USE mod_extrapol
78	USE mod_unit
79	USE mod_printing
80	C
81	C* ---------------------------- Argument declarations -------------------
82	C
83	REAL (kind=ip_realwp_p) pfild(kxlon,kylat), pwork(kxlon,kylat)
84	INTEGER (kind=ip_intwp_p) kmask(kxlon,kylat), kper
85	CHARACTER*8 cdextmet, cdper
86	C
87	C* ---------------------------- Local declarations ----------------------
88	C
89	REAL (kind=ip_realwp_p) zmask(9)
90	INTEGER (kind=ip_intwp_p) ix(9), iy(9)
91	#ifdef key_openmp
92	C TS CHECK FOR EXISTENCE OF WEIGTH FILE
93	INTEGER (kind=ip_intwp_p) :: access, status
94	#endif
95	CHARACTER*8 cladress, clweight, clincrem
96	C
97	C* ---------------------------- Poema verses ----------------------------
98	C
99	C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100	C
101	C* 1. Initialization and allocation of local arrays
102	C ---------------------------------------------
103	C
104	IF (nlogprt .GE. 2) THEN
105	WRITE (UNIT = nulou,FMT = *) ' '
106	WRITE (UNIT = nulou,FMT = *) ' '
107	WRITE (UNIT = nulou,FMT = *)
108	$ ' ROUTINE extrap - Level 3'
109	WRITE (UNIT = nulou,FMT = *)
110	$ ' ************ *****'
111	WRITE (UNIT = nulou,FMT = *) ' '
112	WRITE (UNIT = nulou,FMT = *)
113	$ ' Extrapolation on land with sea values'
114	WRITE (UNIT = nulou,FMT = *) ' '
115	WRITE (UNIT = nulou,FMT = *) ' '
116	ENDIF
117	incre = 0
118	zmask(:)=0
119	ix(:)=0
120	iy(:)=0
121	C
122	IF (.not.allocated(iaddress)) THEN
123	ALLOCATE (iaddress(ig_maxnfn,9,ig_maxgrd),stat=il_err)
124	IF (il_err.NE.0) CALL prtout
125	$ ('Error in "iaddress" allocation of extrap ',il_err,1)
126	iaddress(:,:,:)=0
127	ENDIF
128	IF (.not.allocated(iincre)) THEN
129	ALLOCATE (iincre(ig_maxnfn,ig_maxgrd),stat=il_err)
130	IF (il_err.NE.0) CALL prtout ('Error in iincre
131	$ allocation of extrap ',il_err,1)
132	iincre(:,:)=0
133	ENDIF
134	IF (.not.allocated(zweights)) THEN
135	ALLOCATE (zweights(ig_maxnfn,9,ig_maxgrd),stat=il_err)
136	IF (il_err.NE.0) CALL prtout ('Error in "zweights"
137	$ allocation of extrap ',il_err,1)
138	zweights(:,:,:)=0
139	ENDIF
140	C
141	C* To avoid problems in floating point tests
142	C
143	zwmsk = ABS(pmask) - 1.0
144	C
145	C
146	C* 2. Calculating weights and filling masked values
147	C ---------------------------------------------
148	C
149	IF (cdextmet .EQ. 'NINENN') THEN
150	IF (nlogprt .GE. 2) THEN
151	CALL prtout
152	$ ('NINENN Dataset number for this field is:',knb,1)
153	ENDIF
154	C
155	C* Define locators
156	C
157	WRITE(clweight,'(''WEIGHTS'',I1)') knb
158	WRITE(cladress,'(''ADRESSE'',I1)') knb
159	WRITE(clincrem,'(''INCREME'',I1)') knb
160	C
161	C
162	#ifdef key_openmp
163	status=0
164	status = access ( 'nweights', ' ' )
165	if ( status .eq. 0 ) write(,) "file exists in extrap"
166	if ( status .ne. 0 ) write(,) 'no such file', status
167	IF (status.ne.0 .AND. .NOT. lweight(knb)) THEN
168	#else
169	IF (krdwt .EQ. 1 .AND. .NOT. lweight(knb)) THEN
170	#endif
171	C
172	C* For krdwt=1, the weights and addresses are calculated by Oasis
173	C
174	C* Initialize iteration number, weight and address variables to zero
175	DO 201 jj = 1, kylat
176	DO 202 ji = 1, kxlon
177	iind=(jj-1)*kxlon + ji
178	iincre(knb,iind) = 0
179	DO 203 jl = 1, 9
180	iaddress(knb,jl, iind) = 0
181	zweights(knb,jl, iind) = 0.
182	203 CONTINUE
183	202 CONTINUE
184	201 CONTINUE
185	C
186	200 CONTINUE
187	incre = incre + 1
188	IF (incre .gt. 5000) THEN
189	WRITE (UNIT = nulou,FMT = *)
190	$ 'Problem in extrap.f'
191	WRITE (UNIT = nulou,FMT = *)
192	$ 'Number of iteration greater than 5000'
193	CALL HALTE ('Stop in extrap')
194	ENDIF
195	C
196	!$omp parallel do default (shared)
197	!$omp+ private (jj)
198	DO 210 jj = 1, kylat
199	C
200	C* Test if latitude 1 or kylat or n and n+1 are fully masked
201	C
202	iinf = iminim (kmask(1, jj),kxlon)
203	IF (jj .eq. 1 .OR. jj .eq. kylat) THEN
204	isup=iinf
205	ELSE
206	isup = iminim (kmask(1, jj+1),kxlon)
207	ENDIF
208	C
209	C* If yes, set number of neighbors to 3
210	C
211	IF (iinf .EQ. 1 .AND. isup .EQ. 1) THEN
212	ivoisin = MIN0(3,knbor)
213	ELSE
214	ivoisin = knbor
215	ENDIF
216	DO 220 ji = 1, kxlon
217	iind=(jj-1)*kxlon + ji
218	pwork(ji,jj) = pfild(ji,jj)
219	C
220	C* Case 1: value is not masked
221	C
222	IF (ABS(pfild(ji,jj)) .LT. zwmsk) THEN
223	GOTO 220
224	ELSE
225	C
226	C* Case 2: value at point P is masked
227	C
228	C
229	C We search over the eight closest neighbors
230	C
231	C j+1 7 8 9
232	C j 4 5 6 Current point 5 --> (i,j)
233	C j-1 1 2 3
234	C
235	C i-1 i i+1
236	C ix and iy are index arrays for the neighbors coordinates
237	C
238	inbor = 0
239	ideb = 1
240	ifin = 9
241	C
242	C* Fill up ix array
243	C
244	ix(1) = MAX0 (1,ji - 1)
245	ix(2) = ji
246	ix(3) = MIN0 (kxlon,ji + 1)
247	ix(4) = MAX0 (1,ji - 1)
248	ix(5) = ji
249	ix(6) = MIN0 (kxlon,ji + 1)
250	ix(7) = MAX0 (1,ji - 1)
251	ix(8) = ji
252	ix(9) = MIN0 (kxlon,ji + 1)
253	C
254	C* Fill up iy array
255	C
256	iy(1) = MAX0 (1,jj - 1)
257	iy(2) = MAX0 (1,jj - 1)
258	iy(3) = MAX0 (1,jj - 1)
259	iy(4) = jj
260	iy(5) = jj
261	iy(6) = jj
262	iy(7) = MIN0 (kylat,jj + 1)
263	iy(8) = MIN0 (kylat,jj + 1)
264	iy(9) = MIN0 (kylat,jj + 1)
265	C
266	C* Account for periodicity in longitude
267	C
268	IF (cdper .EQ. 'P') THEN
269	IF (ji .EQ. kxlon) THEN
270	ix(3) = 1 + kper
271	ix(6) = 1 + kper
272	ix(9) = 1 + kper
273	ELSE IF (ji .EQ. 1) THEN
274	ix(1) = kxlon - kper
275	ix(4) = kxlon - kper
276	ix(7) = kxlon - kper
277	ELSE
278	CONTINUE
279	ENDIF
280	ENDIF
281	C
282	C* Correct latitude bounds if southernmost or northernmost points
283	C
284	IF (jj .EQ. 1) ideb = 4
285	IF (jj .EQ. kylat) ifin = 6
286	C
287	C* Grid not periodic
288	C
289	IF (cdper .EQ. 'R') THEN
290	C* ji = 1
291	IF (ji .EQ. 1) THEN
292	ix(1) = ji
293	ix(2) = ji + 1
294	ix(3) = ji
295	ix(4) = ji + 1
296	ix(5) = ji
297	ix(6) = ji + 1
298	ENDIF
299	C* ji = kxlon
300	IF (ji .EQ. kxlon) THEN
301	ix(1) = ji -1
302	ix(2) = ji
303	ix(3) = ji - 1
304	ix(4) = ji
305	ix(5) = ji - 1
306	ix(6) = ji
307	ENDIF
308	C
309	C* Latitude index in both cases
310	C
311	IF (ji .EQ. 1 .OR. ji .EQ. kxlon) THEN
312	ideb = 1
313	ifin = 6
314	iy(1) = MAX0 (1,jj - 1)
315	iy(2) = MAX0 (1,jj - 1)
316	iy(3) = jj
317	iy(4) = jj
318	iy(5) = MIN0(kylat,jj + 1)
319	iy(6) = MIN0(kylat,jj + 1)
320	C
321	C* Correct latitude bounds if southernmost or northernmost points
322	C
323	IF (jj .EQ. 1) ideb = 3
324	IF (jj .EQ. kylat) ifin = 4
325	ENDIF
326	ENDIF
327	C
328	C* Find unmasked neighbors
329	C
330	DO 230 jl = ideb, ifin
331	zmask(jl) = 0.
332	ilon = ix(jl)
333	ilat = iy(jl)
334	IF (ABS(pfild(ilon,ilat)) .LT. zwmsk) THEN
335	zmask(jl) = 1.
336	inbor = inbor + 1
337	ENDIF
338	230 CONTINUE
339	C
340	C* Not enough points around point P are unmasked; interpolation on P
341	C will be done in a future call to extrap.
342	C
343	IF (inbor .LT. ivoisin) THEN
344	GOTO 220
345	ELSE
346	C
347	C* Some points around P are not masked so we use them to extrapolate
348	C* and define the iteration number, weight and address variables
349	C
350	pwork(ji,jj) = 0.
351	iincre(knb, iind) = incre
352	DO 240 jl = ideb, ifin
353	ilon = ix(jl)
354	ilat = iy(jl)
355	pwork(ji,jj) = pwork(ji,jj)
356	$ + pfild(ilon,ilat) * zmask(jl)
357	$ / FLOAT(inbor)
358	iaddress(knb,jl,iind)=(ilat-1)*kxlon+ilon
359	zweights(knb,jl,iind)=zmask(jl)/FLOAT(inbor)
360	240 CONTINUE
361	C
362	ENDIF
363	ENDIF
364	220 CONTINUE
365	210 CONTINUE
366	C
367	C* 3. Writing back unmasked field in pfild and writing weights to file
368	C ----------------------------------------------------------------
369	C
370	C* pfild then contains:
371	C - Values which were not masked
372	C - Interpolated values from the inbor neighbors
373	C - Values which are not yet interpolated
374	C
375	idoit = 0
376	DO 310 jj = 1, kylat
377	DO 320 ji = 1, kxlon
378	IF (ABS(pwork(ji,jj)) .GT. zwmsk) THEN
379	idoit = idoit + 1
380	ENDIF
381	pfild(ji,jj) = pwork(ji,jj)
382	320 CONTINUE
383	310 CONTINUE
384	C
385	IF (idoit .ne. 0) GOTO 200
386	C
387	C* Write weights, addresses and iteration numbers in file
388	C
389	C* Weights
390	#ifdef key_openmp
391	C OPEN statement removed from iniiof
392	C This is to check the existence of weights file in extrap
393
394	OPEN (UNIT = nulgn,FILE = cwninenn,
395	$ FORM ='UNFORMATTED',ERR = 340,IOSTAT = iost)
396	340 CONTINUE
397	#endif
398	C
399	CALL locwrite (clweight, zweights, ig_maxnfn9kxlon*
400	$ kylat, nulgn, iflag)
401	IF (iflag .NE. 0) THEN
402	WRITE (UNIT = nulou,FMT = *)
403	$ 'Problem in writing on UNIT = ', nulgn
404	WRITE (UNIT = nulou,FMT = *)
405	$ 'String locator is = ', clweight
406	CALL HALTE ('Stop in extrap writing weights')
407	ENDIF
408	C
409	C* Adresses
410	C
411	CALL locwrint (cladress, iaddress, ig_maxnfn9kxlon*
412	$ kylat, nulgn, iflag)
413	IF (iflag .NE. 0) THEN
414	WRITE (UNIT = nulou,FMT = *)
415	$ 'Problem in writing on UNIT = ', nulgn
416	WRITE (UNIT = nulou,FMT = *)
417	$ 'String locator is = ', cladress
418	CALL HALTE ('Stop in extrap writing addresses')
419	ENDIF
420	C
421	C* Iteration numbers
422	C
423	CALL locwrint (clincrem, iincre, ig_maxnfnkxlonkylat,
424	$ nulgn, iflag)
425	IF (iflag .NE. 0) THEN
426	WRITE (UNIT = nulou,FMT = *)
427	$ 'Problem in writing on UNIT = ', nulgn
428	WRITE (UNIT = nulou,FMT = *)
429	$ 'String locator is = ', clincrem
430	CALL HALTE ('Stop in extrap writing iteration number')
431	ENDIF
432	C
433	C* Printing
434	C
435	IF (nlogprt .GE. 2) THEN
436	CALL prtout
437	$ ('Wrote weights,addresses and iteration numbers on unit = ',
438	$ nulgn, 1)
439	ENDIF
440	ENDIF
441	C
442	#ifdef key_openmp
443	IF ( .NOT. lweight(knb)) THEN
444	#else
445	IF (krdwt .EQ. 0 .AND. .NOT. lweight(knb)) THEN
446	#endif
447	C
448	C
449	C* 4. Reading weights and filling masked values
450	C ---------------------------------------------
451	C
452	C* For krdwt=0, the weights, addresses and iteration numbers are read by Oasis
453	C
454	CALL locread (clweight, zweights, ig_maxnfn9kxlon*
455	$ kylat, nulgn, iflag)
456	IF (iflag .NE. 0) THEN
457	WRITE (UNIT = nulou,FMT = *)
458	$ 'Problem in reading on UNIT = ', nulgn
459	WRITE (UNIT = nulou,FMT = *)
460	$ 'String locator is = ', clweight
461	CALL HALTE ('Stop in extrap reading weights')
462	ENDIF
463	C
464	CALL locrint (cladress, iaddress, ig_maxnfn9kxlon*
465	$ kylat, nulgn, iflag)
466	IF (iflag .NE. 0) THEN
467	WRITE (UNIT = nulou,FMT = *)
468	$ 'Problem in reading on UNIT = ', nulgn
469	WRITE (UNIT = nulou,FMT = *)
470	$ 'String locator is = ', cladress
471	CALL HALTE ('Stop in extrap reading addresses')
472	ENDIF
473	C
474	CALL locrint (clincrem, iincre, ig_maxnfnkxlonkylat,
475	$ nulgn, iflag)
476	IF (iflag .NE. 0) THEN
477	WRITE (UNIT = nulou,FMT = *)
478	$ 'Problem in reading on UNIT = ', nulgn
479	WRITE (UNIT = nulou,FMT = *)
480	$ 'String locator is = ', clincrem
481	CALL HALTE('Stop in extrap reading iteration numbers')
482	ENDIF
483	ENDIF
484	C
485	C* Extrapolation
486	C
487	#ifdef key_openmp
488	IF ( lweight(knb)) THEN
489	#else
490	IF (krdwt .EQ. 0 .OR. lweight(knb)) THEN
491	#endif
492	cvg>>>
493	idoitold = 1000000
494	cvg<<<
495	400 CONTINUE
496
497	incre = incre + 1
498	DO 410 jj = 1, kylat
499	DO 420 ji = 1, kxlon
500	iind=(jj-1)*kxlon + ji
501	pwork(ji,jj) = pfild(ji,jj)
502	C
503	C* Case 1: value is not masked
504	C
505	IF (ABS(pfild(ji,jj)) .LT. zwmsk) THEN
506	GOTO 420
507	ELSE
508	C
509	C* Case 2: value at point P is masked
510	C
511	IF (iincre(knb,iind) .EQ. incre) THEN
512	pwork(ji,jj) = 0.
513	DO 250 jl = 1,9
514	IF (zweights(knb,jl,iind).EQ.0.) GO TO 250
515	idivi=iaddress(knb, jl, iind)/kxlon
516	imult=idivi*kxlon
517	IF (iaddress(knb,jl, iind) .EQ. imult) THEN
518	ilat=idivi
519	ilon=kxlon
520	ELSE
521	ilat=idivi+1
522	ilon=iaddress(knb,jl, iind)-imult
523	ENDIF
524	pwork(ji,jj) = pwork(ji,jj)
525	$ + pfild(ilon,ilat)
526	$ *zweights(knb,jl,iind)
527	250 CONTINUE
528	ENDIF
529	ENDIF
530	420 CONTINUE
531	410 CONTINUE
532
533	C
534	C* 5. Writing back unmasked field in pfild
535	C ------------------------------------
536	C
537	C* pfild then contains:
538	C - Values which were not masked
539	C - Interpolated values from the inbor neighbors
540	C - Values which are not yet interpolated
541	C
542	idoit = 0
543	DO 510 jj = 1, kylat
544	DO 520 ji = 1, kxlon
545	IF (ABS(pwork(ji,jj)) .GT. zwmsk) THEN
546	idoit = idoit + 1
547	ENDIF
548	pfild(ji,jj) = pwork(ji,jj)
549	520 CONTINUE
550	510 CONTINUE
551	cvg>>>
552	IF (idoit .eq. idoitold) THEN
553	WRITE (UNIT = nulou,FMT = *)
554	$ 'extrap: nweights does not fit to your masks file'
555	WRITE (UNIT = nulou,FMT = *)
556	$ 'extrap: create new file nweights ',
557	$ '($NIO = 1 in namcouple)'
558	CALL HALTE ('Stop in extrap: no progress in iteration')
559	ENDIF
560	idoitold = idoit
561	cvg<<<
562	IF (idoit .ne. 0) GOTO 400
563	ENDIF
564	C
565	C* Printing
566	C
567	IF (nlogprt .GE. 2) THEN
568	CALL prtout
569	$ ('Number of extrapolation steps incre =', incre, 1)
570	CALL FLUSH(nulou)
571	ENDIF
572	ENDIF
573	C
574	C* Put flag indicating that weights have been read or calculated to 1
575	C
576	lweight(knb)= .TRUE.
577	C
578	C* 4. End of routine
579	C --------------
580	C
581	IF (nlogprt .GE. 2) THEN
582	WRITE (UNIT = nulou,FMT = *) ' '
583	WRITE (UNIT = nulou,FMT = *)
584	$ ' --------- End of routine extrap ---------'
585	CALL FLUSH (nulou)
586	ENDIF
587	RETURN
588	END

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: CPL/oasis3/trunk/src/mod/oasis3/src/extrap.F @ 1677

Download in other formats: