/[MITgcm]/MITgcm_contrib/cg2d_bench/ini_cg2d.F
ViewVC logotype

Contents of /MITgcm_contrib/cg2d_bench/ini_cg2d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (show annotations) (download)
Fri May 12 22:23:42 2006 UTC (17 years, 11 months ago) by ce107
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +5 -4 lines
Fixed for single/double precision.

1 C $Id$
2 CStartOfInterface
3 SUBROUTINE INI_CG2D
4 C /==========================================================\
5 C | SUBROUTINE INI_CG2D |
6 C | o Initialise 2d conjugate gradient solver operators. |
7 C |==========================================================|
8 C | These arrays are purely a function of the basin geom. |
9 C | We set then here once and them use then repeatedly. |
10 C | In this example we hard code a periodic channel with a |
11 C | grid spacing on 1. |
12 C \==========================================================/
13 IMPLICIT NONE
14
15 C === Global variables ===
16 #include "SIZE.h"
17 #include "EEPARAMS.h"
18 #include "PARAMS.h"
19 #include "CG2D.h"
20
21 C === Routine arguments ===
22 C myThid - Thread no. that called this routine.
23 INTEGER myThid
24 CEndOfInterface
25
26 C === Local variables ===
27 C xG, yG - Global coordinate location.
28 C zG
29 C iG, jG - Global coordinate index
30 C bi,bj - Loop counters
31 C faceArea - Temporary used to hold cell face areas.
32 C I,J,K
33 INTEGER I, J, K
34 INTEGER iG, jG
35
36 C-- Initialise laplace operator
37 C aW2d: integral Ax/dX
38 C aS2d: integral Ay/dY
39 DO J=1-OLy,sNy+OLy
40 DO I=1-OLx,sNx+OLx
41 iG = myXGlobalLo+I-1
42 jG = myYGlobalLo+J-1
43 aW2d(I,J) = 1. _d 0
44 aS2d(I,J) = 1. _d 0
45 IF ( jG .EQ. 1 .OR. jG .EQ. nY+1 ) THEN
46 aS2d(I,J) = 0. _d 0
47 ENDIF
48 ENDDO
49 ENDDO
50 C CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D')
51 C CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D')
52
53 C-- Initialise preconditioner
54 DO J=1,sNy
55 DO I=1,sNx
56 pC(I,J) = 1. _d 0
57 IF (
58 & aW2d(I,J) + aW2d(I+1,J)
59 & +aS2d(I,J) + aS2D(I,J+1)
60 & .EQ. 0. _d 0
61 & ) pC(I,J) = 0. _d 0
62 pW(I,J) = 0. _d 0
63 pS(I,J) = 0. _d 0
64 ENDDO
65 ENDDO
66 C-- Update overlap regions
67 CALL EXCH_XY_R8(pC)
68 CALL EXCH_XY_R8(pW)
69 CALL EXCH_XY_R8(pS)
70
71 C-- Set default values for initial guess
72 DO J=1,sNy
73 DO I=1,sNx
74 cg2d_x(I,J) = 0. _d 0
75 ENDDO
76 ENDDO
77 C-- Update overlap regions
78 CALL EXCH_XY_R8(cg2d_x)
79
80 RETURN
81 END

  ViewVC Help
Powered by ViewVC 1.1.22