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.

solisl_fdir.h90 in tags/start/NEMO/OPA_SRC/SOL – NEMO

Context Navigation

source: tags/start/NEMO/OPA_SRC/SOL/solisl_fdir.h90 @ 328

Last change on this file since 328 was 3, checked in by opalod, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.6 KB

Line
1	!!----------------------------------------------------------------------
2	!! * solisl_fdir.h90 *
3	!!----------------------------------------------------------------------
4
5	!!----------------------------------------------------------------------
6	!! islmat : Compute the matrix associated with islands
7	!! and save it in a direct access file
8	!! islbsf : Compute the barotropic streamfunction of each island
9	!! and save it in a direct access file
10	!!----------------------------------------------------------------------
11
12	SUBROUTINE islmat
13	!!----------------------------------------------------------------------
14	!! * ROUTINE islmat *
15	!!
16	!! ** Purpose : Compute and invert the island matrix
17	!!
18	!! ** Method : aisl(jni,jnj) matrix constituted by bsfisl circulation
19	!!
20	!! ** Action : - aisl : island matrix
21	!! - aislm1 : invert of the island matrix
22	!!
23	!! ** file : islands: contain bsfisl, aisl and aislm1
24	!!
25	!! History :
26	!! 1.0 ! 88-10 (G. Madec) Original code
27	!! 7.0 ! 96-01 (G. Madec) suppression of common work arrays
28	!! 8.5 ! 02-08 (G. Madec) Free form, F90
29	!!----------------------------------------------------------------------
30	!! * Local declarations
31	INTEGER :: ji, jj, jni, jnj, ios, jl
32	REAL(wp) :: zwx(jpi,jpj,2)
33	!!----------------------------------------------------------------------
34
35
36	! I. Island matrix lecture in numisl (if it exists)
37	! ==================================
38	! file already open in islbsf routine
39	! Lecture
40	READ( numisl, REC=jpisl+2, IOSTAT=ios ) aisl, aislm1
41	IF( ios /= 0 ) GOTO 110
42
43	! Control print
44	IF(lwp) THEN
45	WRITE(numout,*)
46	WRITE(numout,*) ' islmat: lecture aisl/aislm1 in numisl done'
47	WRITE(numout,*) ' ~~~~~~'
48	WRITE(numout,*)
49	WRITE(numout,*) ' island matrix'
50	WRITE(numout,*) ' -------------'
51	WRITE(numout,*)
52
53	DO jnj = 1, jpisl
54	WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl )
55	END DO
56
57	WRITE(numout,*)
58	WRITE(numout,*) ' inverse of the island matrix'
59	WRITE(numout,*) ' ----------------------------'
60	WRITE(numout,*)
61
62	DO jnj = 1, jpisl
63	WRITE(numout,'(12e11.3)') ( aislm1(jni,jnj), jni=1, jpisl )
64	END DO
65	ENDIF
66
67	CLOSE(numisl)
68
69	RETURN
70
71	110 CONTINUE
72
73
74	! II. Island matrix computation
75	! =============================
76
77	DO jnj = 1, jpisl
78
79	! 1.1 Circulation of bsf(jnj) around island jni
80
81	DO jj = 2, jpj
82	zwx(:,jj,1) = -( bsfisl(:,jj,jnj) - bsfisl(:,jj-1,jnj) )
83	END DO
84	DO ji = 1, jpi
85	zwx(ji,1,1) = 0.e0
86	END DO
87
88	DO jj = 1, jpj
89	DO ji = 2, jpi
90	zwx(ji,jj,2) = ( bsfisl(ji,jj,jnj) - bsfisl(ji-1,jj,jnj) )
91	END DO
92	END DO
93	DO jj = 1, jpj
94	zwx(1,jj,2) = 0.e0
95	END DO
96
97	! 1.2 Island matrix
98
99	DO jni = 1, jpisl
100	aisl(jni,jnj) = 0.
101	DO jl = 1, 2
102	DO jj = 1, jpj
103	DO ji = 1, jpi
104	aisl(jni,jnj) = aisl(jni,jnj) + acisl2(ji,jj,jl,jni) * zwx(ji,jj,jl)
105	END DO
106	END DO
107	END DO
108	END DO
109
110	END DO
111	# if defined key_mpp
112	CALL mpp_sum( aisl, jpisl*jpisl )
113	# endif
114
115	! 1.3 Control print
116
117	IF(lwp) THEN
118	WRITE(numout,*)
119	WRITE(numout,*) ' islmat: island matrix'
120	WRITE(numout,*) ' ~~~~~~'
121	WRITE(numout,*)
122
123	DO jnj = 1, jpisl
124	WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl )
125	END DO
126	ENDIF
127
128
129	! 2. Invertion of the island matrix
130	! ---------------------------------
131
132	! 2.1 Call of an imsl routine for the matrix invertion
133
134	CALL linrg( jpisl, aisl, jpisl, aislm1, jpisl )
135
136	! 2.2 Control print
137
138	IF(lwp) THEN
139	WRITE(numout,*)
140	WRITE(numout,*) ' islmat: inverse of the island matrix'
141	WRITE(numout,*) ' ~~~~~~'
142	WRITE(numout,*)
143
144	DO jnj = 1, jpisl
145	WRITE(numout, '(12e11.3)') ( aislm1(jni,jnj), jni=1, jpisl )
146	END DO
147	ENDIF
148
149
150	! 3. Output of aisl and aislm1 in numisl
151	! --------------------------------------
152
153	IF(lwp) WRITE(numisl,REC=jpisl+2) aisl, aislm1
154
155	CLOSE( numisl )
156
157	END SUBROUTINE islmat
158
159
160	SUBROUTINE islbsf
161	!!----------------------------------------------------------------------
162	!! * ROUTINE islbsf *
163	!!
164	!! ** Purpose : Compute the barotropic stream function associated with
165	!! each island using FETI or preconditioned conjugate gradient
166	!! method or read them in numisl.
167	!!
168	!! ** Method : Direct acces file
169	!!
170	!! ** file : numisl: barotropic stream function associated
171	!! with each island of the domain
172	!!
173	!! ** Action : - bsfisl, the streamfunction which takes the value 1
174	!! over island ni, and 0 over the others islands
175	!! - update 'numisl' file
176	!!
177	!! History :
178	!! ! 87-10 (G. Madec) Original code
179	!! ! 91-11 (G. Madec)
180	!! ! 93-03 (G. Madec) release 7.1
181	!! ! 93-04 (M. Guyon) loops and suppress pointers
182	!! ! 96-11 (A. Weaver) correction to preconditioning
183	!! ! 98-02 (M. Guyon) FETI method
184	!! ! 99-11 (M. Imbard) NetCDF FORMAT with IOIPSL
185	!! 8.5 ! 02-08 (G. Madec) Free form, F90
186	!!----------------------------------------------------------------------
187	!! * Local declarations
188	INTEGER :: ji, jj, jni, jii, jnp
189	INTEGER :: iimlu, ijmlu, inmlu, ios, iju, iloc
190	INTEGER :: ii, ij, icile, icut, inmax, indic
191	REAL(wp) :: zepsr, zeplu, zgwgt
192	REAL(wp) :: zep(jpisl)
193	!!----------------------------------------------------------------------
194
195
196	! 0. Initializations
197	! ==================
198
199	! set to zero the convergence indicator
200	icile = 0
201
202	! set to zero of bsfisl
203	bsfisl(:,:,1:jpisl) = 0.e0
204
205
206	! I. Lecture of bsfisl in numisl (if it exists)
207	! =============================================
208
209	! open islands file
210	ibloc = 4096
211	ilglo = ibloc((jpiglojpjglo*jpbyt-1 )/ibloc+1)
212	clname = 'islands'
213	CALL ctlopn( numisl, 'islands', 'UNKNOWN', 'UNFORMATTED', 'DIRECT', &
214	ilglo,numout,lwp,0)
215
216	icut = 0
217	iimlu = 0
218	ijmlu = 0
219	inmlu = 0
220	zeplu = 0.
221	READ (numisl, REC=1, IOSTAT=ios) iimlu, ijmlu, inmlu, zeplu
222	IF( ios /= 0 ) GOTO 110
223	DO jni = 1, jpisl
224	CALL read2( numisl, bsfisl(1,1,jni), 1, jni+1 )
225	END DO
226
227	110 CONTINUE
228
229	! the read precision is not the required one : icut=888
230	IF( zeplu > epsisl ) icut = 888
231
232	! the read domain does not correspond to the model one : icut=999
233	IF( iimlu /= jpiglo .OR. ijmlu /= jpjglo .OR. inmlu /= jpisl ) icut = 999
234
235	! Control print
236	IF( icut == 999 ) THEN
237	IF(lwp) THEN
238	WRITE(numout,*)
239	WRITE(numout,*) ' islbsf : lecture bsfisl in numisl failed'
240	WRITE(numout,*) ' ~~~~~~'
241	WRITE(numout,*) ' icut= ', icut
242	WRITE(numout,*) ' imlu= ', iimlu, ' jmlu= ', ijmlu, &
243	' nilu= ', inmlu, ' epsisl lu= ', zeplu
244	WRITE(numout,*) ' the bsfisl are computed from zero'
245	ENDIF
246	ELSEIF( icut == 888 ) THEN
247	IF(lwp) THEN
248	WRITE(numout,*)
249	WRITE(numout,*) ' islbsf : lecture bsfisl in numisl done'
250	WRITE(numout,*) ' ~~~~~~'
251	WRITE(numout,*) ' the required accuracy is not reached'
252	WRITE(numout,*) ' epsisl lu= ', zeplu,' epsisl required ', epsisl
253	WRITE(numout,*) ' the bsfisl are computed from the read values'
254	ENDIF
255	ELSE
256	IF(lwp) THEN
257	WRITE(numout,*)
258	WRITE(numout,*) ' islbsf : lecture bsfisl in numisl done'
259	WRITE(numout,*) ' ~~~~~~'
260	ENDIF
261	RETURN
262	ENDIF
263
264
265	! II. Compute the bsfisl (if icut=888 or 999)
266	! ======================
267
268	! save nmax
269	inmax = nmax
270
271	! set the number of iteration of island computation
272	nmax = nmisl
273
274	! Loop over islands
275	! -----------------
276
277	DO jni = 1, jpisl
278
279
280	! 1. Initalizations of island computation
281	! ---------------------------------------
282
283	! 1.0 Set the pcg solution gcb to zero
284	gcb(:,:) = 0.e0
285
286	! 1.1 Set first guess gcx either to zero or to the read bsfisl
287
288	IF( icut == 999 ) THEN
289	gcx(:,:) = 0.e0
290	ELSEIF( icut == 888 ) THEN
291	IF(lwp) WRITE(numout,*) ' islbsf: bsfisl read are used as first guess'
292	! C a u t i o n : bsfisl masked because it contains 1 along island seaside
293	gcx(:,:) = bsfisl(:,:,jni) * bmask(:,:)
294	ENDIF
295
296	! 1.2 Right hand side of the stream FUNCTION equation
297
298	# if defined key_mpp
299
300	! north fold treatment
301	IF( npolj == 3 ) iloc = jpiglo -(nimpp-1+nimppt(nono+1)-1)
302	IF( npolj == 4 ) iloc = jpiglo - 2*(nimpp-1)
303	t2p1(:,1,1) = 0.
304	! north and south grid-points
305	DO jii = 1, 2
306	DO jnp = 1, mnisl(jii,jni)
307	ii = miisl(jnp,jii,jni)
308	ij = mjisl(jnp,jii,jni)
309	IF( ( npolj == 3 .OR. npolj == 4 ) .AND. ( ij == nlcj-1 .AND. jii == 1) ) THEN
310	iju=iloc-ii+1
311	t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
312	ELSE
313	gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
314	ENDIF
315	END DO
316	END DO
317
318	! east and west grid-points
319
320	DO jii = 3, 4
321	DO jnp = 1, mnisl(jii,jni)
322	ii = miisl(jnp,jii,jni)
323	ij = mjisl(jnp,jii,jni)
324	gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
325	END DO
326	END DO
327	CALL mpplnks( gcb )
328
329	# else
330
331	! north and south grid-points
332	DO jii = 1, 2
333	DO jnp = 1, mnisl(jii,jni)
334	ii = miisl(jnp,jii,jni)
335	ij = mjisl(jnp,jii,jni)
336	IF( ( nperio == 3 .OR. nperio == 4 ) .AND. ( ij == jpj-1 .AND. jii == 1) ) THEN
337	gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
338	ELSE
339	gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
340	ENDIF
341	END DO
342	END DO
343
344	! east and west grid-points
345	DO jii = 3, 4
346	DO jnp = 1, mnisl(jii,jni)
347	ii = miisl(jnp,jii,jni)
348	ij = mjisl(jnp,jii,jni)
349	IF( bmask(ii-jii+3,ij) /= 0. ) THEN
350	gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
351	ELSE
352	! east-west cyclic boundary conditions
353	IF( ii-jii+3 == 1 ) THEN
354	gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
355	ENDIF
356	ENDIF
357	END DO
358	END DO
359
360	# endif
361
362	! 1.4 Preconditioned right hand side and absolute precision
363
364	IF( nsolv == 3 ) THEN
365	! FETI method
366	ncut = 0
367	rnorme=0.
368	DO jj = 1, jpj
369	DO ji = 1, jpi
370	gcb (ji,jj) = bmask(ji,jj) * gcb(ji,jj)
371	rnorme = rnorme + gcb(ji,jj) * gcb(ji,jj)
372	END DO
373	END DO
374
375	CALL mpp_sum( rnorme )
376
377	IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
378	epsr = epsisl * epsisl * rnorme
379	indic = 0
380	ELSE
381	ncut = 0
382	rnorme = 0.
383	DO jj = 1, jpj
384	DO ji = 1, jpi
385	gcb (ji,jj) = gcdprc(ji,jj) * gcb(ji,jj)
386	zgwgt = gcdmat(ji,jj) * gcb(ji,jj)
387	rnorme = rnorme + gcb(ji,jj) * zgwgt
388	END DO
389	END DO
390	#if defined key_mpp
391	CALL mpp_sum( rnorme )
392	#endif
393	IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
394	epsr = epsisl * epsisl * rnorme
395	indic = 0
396	ENDIF
397
398
399	! 3. PCG solver for gcp.gcx=gcb in monotask
400	! -------------------------------------------
401
402	IF(nsolv == 3) THEN
403	! FETI method
404	! precision to compute Islands matrix A
405
406	epsilo = epsisl
407	CALL solfet( indic )
408	! --->
409	! precision to compute grad PS
410
411	epsilo = eps
412	ELSE
413
414	CALL solpcg( indic )
415
416	ENDIF
417
418
419	! 4. Save the solution in bsfisl
420	! ------------------------------
421
422	DO jj = 1, jpj
423	DO ji = 1, jpi
424	bsfisl(ji,jj,jni) = gcx(ji,jj)
425	END DO
426	END DO
427
428
429	! 5. Boundary conditions
430	! ----------------------
431
432	! set to 1. coastal gridpoints of the island
433	DO jnp = 1, mnisl(0,jni)
434	ii = miisl(jnp,0,jni)
435	ij = mjisl(jnp,0,jni)
436	bsfisl(ii,ij,jni) = 1.
437	END DO
438
439	! cyclic boundary conditions
440	IF( nperio == 1 .OR. nperio == 4 ) THEN
441	DO jj = 1, jpj
442	bsfisl( 1 ,jj,jni) = bsfisl(jpim1,jj,jni)
443	bsfisl(jpi,jj,jni) = bsfisl( 2 ,jj,jni)
444	END DO
445	ENDIF
446	IF( nperio == 3 .OR. nperio == 4 ) THEN
447	DO ji = 1, jpim1
448	iju = jpi-ji+1
449	bsfisl(ji,jpj-1,jni) = bsfisl(iju,jpj-2,jni)
450	bsfisl(ji, jpj ,jni) = bsfisl(iju,jpj-3,jni)
451	END DO
452	ENDIF
453	#if defined key_mpp
454	CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. )
455	#endif
456
457
458	! 6. Control print
459	! ----------------
460
461	IF(lwp)WRITE(numout,*)
462	IF(lwp)WRITE(numout,*) ' islbsf: island number: ', jni
463	IF(lwp)WRITE (numout,9290) niter, res, SQRT(epsr)/epsisl
464	zep(jni) = MAX(epsisl, res/(SQRT(epsr)/epsisl))
465	IF(indic < 0) THEN
466	icile = icile-1
467	IF(lwp)WRITE(numout,*) ' pcg do not converge for island: ', jni
468	IF(lwp)WRITE(numout,*) ' Precision reached: ',zep(jni)
469	ENDIF
470
471	9290 FORMAT(' niter :',i4,' , res :',e20.10,' , gcb :',e20.10)
472
473
474	! End loop for islands
475	! --------------------
476
477	END DO
478
479
480	! 7. Reset PCG
481	! ------------
482
483	! reset the number of iteration for pcg
484	nmax = inmax
485
486	! reset to zero pcg arrays
487	gcx (:,:) = 0.e0
488	gcxb (:,:) = 0.e0
489	gcb (:,:) = 0.e0
490	gcr (:,:) = 0.e0
491	gcdes(:,:) = 0.e0
492	gccd (:,:) = 0.e0
493
494
495	! III. Output of bsfisl in numisl
496	! ===============================
497
498	IF( icile == 0 .AND. icut /= 0 ) THEN
499	IF(lwp) THEN
500	WRITE(numout,*)
501	WRITE(numout,*)' islbsf : write bsfisl in numisl ', numisl
502	WRITE(numout,*)' ~~~~~~'
503	ENDIF
504
505	IF(lwp) WRITE(numisl,REC=1) jpiglo, jpjglo, jpisl, epsisl
506	DO jni = 1, jpisl
507	CALL write2( numisl, bsfisl(1,1,jni), 1, jni+1 )
508	END DO
509	ENDIF
510
511	IF( icile < 0 ) THEN
512	IF(lwp) THEN
513	WRITE(numout,*)
514	WRITE(numout,*) ' islbsf : number of island without convergence: ',ABS(icile)
515	WRITE(numout,*) ' ~~~~~~'
516	ENDIF
517	zepsr = epsisl
518	DO jni = 1, jpisl
519	IF(lwp)WRITE(numout,*) ' isl ',jni,' precision reached ', zep(jni)
520	zepsr = MAX( zep(jni), zepsr )
521	END DO
522	IF( zepsr == 0. ) zepsr = epsisl
523	IF(lwp) THEN
524	WRITE(numout,*) ' save value of precision reached: ',zepsr
525	WRITE(numout,*)
526	WRITE(numout,*)' islbsf : save bsfisl in numisl ',numisl
527	WRITE(numout,*)' ~~~~~~'
528	ENDIF
529
530	IF(lwp) WRITE(numisl,REC=1) jpiglo, jpjglo, jpisl, zepsr
531	DO jni = 1, jpisl
532	CALL write2( numisl, bsfisl(1,1,jni), 1, jni+1 )
533	END DO
534	nstop = nstop + 1
535	ENDIF
536
537	END SUBROUTINE islbsf

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

Download in other formats:

Original Format