model/src/ini_cori.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cori.F,v 1.27 2008/05/28 03:06:57 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_CORI

C     !INTERFACE:
      SUBROUTINE INI_CORI( myThid )
C     !DESCRIPTION:
C     Initialise coriolis term.

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_EXCH2
# include "W2_EXCH2_SIZE.h"
# include "W2_EXCH2_TOPOLOGY.h"
#endif
#ifdef ALLOW_MNC
# include "MNC_PARAMS.h"
#endif
#ifdef ALLOW_MONITOR
# include "MONITOR.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
      INTEGER myThid
CEOP

C     === Functions ====
      LOGICAL  MASTER_CPU_IO
      EXTERNAL MASTER_CPU_IO

C     !LOCAL VARIABLES:
C     bi,bj   :: Tile Indices counters
C     i, j    :: Loop counters
C     facGrid :: Factor for grid to meter conversion
      INTEGER bi,bj
      INTEGER i, j
      _RL facGrid
#ifndef OLD_GRID_IO
      INTEGER myTile, iG, iLen
      CHARACTER*(MAX_LEN_FNAM) fName
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER  ILNBLNK
      EXTERNAL ILNBLNK
#endif


C     Initialise coriolis parameter
      IF     ( useConstantF ) THEN
C       Constant F case
        DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j=1-Oly,sNy+Oly
              DO i=1-Olx,sNx+Olx
                fCori(i,j,bi,bj)  = f0
                fCoriG(i,j,bi,bj) = f0
                fCoriCos(i,j,bi,bj)=0. _d 0
              ENDDO
            ENDDO
          ENDDO
        ENDDO
      ELSEIF ( useBetaPlaneF ) THEN
C       Beta plane case
        facGrid = 1. _d 0
        IF ( usingSphericalPolarGrid
     &     .OR. usingCurvilinearGrid ) facGrid = deg2rad*rSphere
        DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j=1-Oly,sNy+Oly
              DO i=1-Olx,sNx+Olx
                fCori(i,j,bi,bj)  = f0+beta*_yC(i,j,bi,bj)*facGrid
                fCoriG(i,j,bi,bj) = f0+beta* yG(i,j,bi,bj)*facGrid
                fCoriCos(i,j,bi,bj)=0. _d 0
              ENDDO
            ENDDO
          ENDDO
        ENDDO
      ELSEIF ( useSphereF ) THEN
C       Spherical case
C       Note in this case we assume yC is in degrees.
        DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j=1-Oly,sNy+Oly
              DO i=1-Olx,sNx+Olx
                fCori(i,j,bi,bj)  =
     &                 2. _d 0*omega*sin(_yC(i,j,bi,bj)*deg2rad)
                fCoriG(i,j,bi,bj) =
     &                 2. _d 0*omega*sin(yG(i,j,bi,bj)*deg2rad)
                fCoriCos(i,j,bi,bj)=
     &                 2. _d 0*omega*cos(_yC(i,j,bi,bj)*deg2rad)
              ENDDO
            ENDDO
          ENDDO
        ENDDO
c       CALL WRITE_FLD_XY_RL('fCoriC',' ',fCori ,  0,myThid)
c       CALL WRITE_FLD_XY_RL('fCoriG',' ',fCoriG , 0,myThid)
c       CALL WRITE_FLD_XY_RL('fCorCs',' ',fCoriCos,0,myThid)
      ELSE
C       Special custom form
        DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j=1-Oly,sNy+Oly
              DO i=1-Olx,sNx+Olx
                fCori(i,j,bi,bj)  = 0. _d 0
                fCoriG(i,j,bi,bj) = 0. _d 0
                fCoriCos(i,j,bi,bj)=0. _d 0
              ENDDO
            ENDDO
          ENDDO
        ENDDO
        _BARRIER
        CALL READ_REC_XY_RS( 'fCoriC.bin', fCori,   1, 0, myThid )
        CALL READ_REC_XY_RS( 'fCorCs.bin', fCoriCos,1, 0, myThid )
        IF ( .NOT.useCubedSphereExchange ) THEN
         CALL READ_REC_XY_RS('fCoriG.bin', fCoriG,  1, 0, myThid )
        ELSE
#ifdef OLD_GRID_IO
         CALL READ_REC_XY_RS('fCoriG.bin', fCoriG,  1, 0, myThid )
         _BARRIER
C-       deal with the 2 missing corners (for fCoriG):
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
C-  Notes: this will only works with 6 tiles (1 per face) and
C    with 2 polar faces + 4 equatorials:
           IF (bi.LE.3 .OR. bi.GE.5) THEN
             fCoriG(sNx+1,1,bi,bj) = fCoriG(1,1,bi,bj)
           ELSE
             fCoriG(sNx+1,1,bi,bj) = -fCoriG(1,1,bi,bj)
           ENDIF
           IF (bi.GE.3) THEN
             fCoriG(1,sNy+1,bi,bj) = fCoriG(1,1,bi,bj)
             fCoriG(sNx+1,sNy+1,bi,bj) = fCoriG(sNx+1,1,bi,bj)
           ELSE
             fCoriG(1,sNy+1,bi,bj) = -fCoriG(1,1,bi,bj)
             fCoriG(sNx+1,sNy+1,bi,bj) = -fCoriG(sNx+1,1,bi,bj)
           ENDIF
          ENDDO
         ENDDO
#else  /* OLD_GRID_IO */
         _BEGIN_MASTER(myThid)
         DO bj = 1,nSy
          DO bi = 1,nSx
           iG = bi+(myXGlobalLo-1)/sNx
           myTile = iG
#ifdef ALLOW_EXCH2
           myTile = W2_myTileList(bi)
           iG = exch2_myface(myTile)
#endif
           WRITE(fName,'(2A,I3.3,A)') 'fCoriG','.face',iG,'.bin'
           iLen = ILNBLNK(fName)
           WRITE(msgBuf,'(A,I6,2A)')
     &       ' Reading tile:', myTile, ' from file ', fName(1:iLen)
           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                         SQUEEZE_RIGHT , myThid )
#ifdef ALLOW_MDSIO
           CALL MDS_FACEF_READ_RS( fName, readBinaryPrec, 1,
     &                             fCoriG, bi, bj, myThid )
#else /* ALLOW_MDSIO */
           WRITE(msgBuf,'(A)') 'INI_CORI: Needs to compile MDSIO pkg'
           CALL PRINT_ERROR( msgBuf, myThid )
           STOP 'ABNORMAL END: S/R INI_CORI'
#endif /* ALLOW_MDSIO */

          ENDDO
         ENDDO
         _END_MASTER(myThid)
#endif /* OLD_GRID_IO */
        ENDIF

        CALL EXCH_XY_RS( fCori, myThid )
        CALL EXCH_XY_RS( fCoriCos, myThid )
        CALL EXCH_Z_3D_RS( fCoriG, 1, myThid )
      ENDIF

#ifdef ALLOW_MONITOR
      IF ( MASTER_CPU_IO(myThid) ) THEN
C--   only the master thread is allowed to switch On/Off mon_write_stdout
C     & mon_write_mnc (since it's the only thread that uses those flags):

        IF (monitor_stdio) THEN
          mon_write_stdout = .TRUE.
        ELSE
          mon_write_stdout = .FALSE.
        ENDIF
        mon_write_mnc = .FALSE.
#ifdef ALLOW_MNC
        IF (useMNC .AND. monitor_mnc) THEN
          DO i = 1,MAX_LEN_MBUF
            mon_fname(i:i) = ' '
          ENDDO
          mon_fname(1:12) = 'monitor_grid'
          CALL MNC_CW_SET_UDIM(mon_fname, 1, myThid)
          mon_write_mnc = .TRUE.
        ENDIF
#endif /*  ALLOW_MNC  */

      ENDIF

      CALL MON_PRINTSTATS_RS(1,fCori,'fCori',myThid)
      CALL MON_PRINTSTATS_RS(1,fCoriG,'fCoriG',myThid)
      CALL MON_PRINTSTATS_RS(1,fCoriCos,'fCoriCos',myThid)

      IF ( MASTER_CPU_IO(myThid) ) THEN
        mon_write_stdout = .FALSE.
        mon_write_mnc    = .FALSE.
      ENDIF
#endif /* ALLOW_MONITOR */

      RETURN
      END
1	jmc	1.28	C $Header: /u/gcmpack/MITgcm/model/src/ini_cori.F,v 1.27 2008/05/28 03:06:57 jmc Exp $
2	adcroft	1.10	C $Name: $
3	cnh	1.1
4	jmc	1.17	#include "PACKAGES_CONFIG.h"
5	cnh	1.6	#include "CPP_OPTIONS.h"
6	cnh	1.1
7	cnh	1.12	CBOP
8			C !ROUTINE: INI_CORI
9	edhill	1.18
10	cnh	1.12	C !INTERFACE:
11	cnh	1.1	SUBROUTINE INI_CORI( myThid )
12	edhill	1.18	C !DESCRIPTION:
13	jmc	1.23	C Initialise coriolis term.
14	cnh	1.12
15			C !USES:
16	adcroft	1.7	IMPLICIT NONE
17	cnh	1.1	#include "SIZE.h"
18			#include "EEPARAMS.h"
19			#include "PARAMS.h"
20	jmc	1.23	#include "GRID.h"
21	jmc	1.27	#ifdef ALLOW_EXCH2
22	jmc	1.28	# include "W2_EXCH2_SIZE.h"
23	jmc	1.27	# include "W2_EXCH2_TOPOLOGY.h"
24			#endif
25	edhill	1.20	#ifdef ALLOW_MNC
26	jmc	1.27	# include "MNC_PARAMS.h"
27	edhill	1.20	#endif
28	edhill	1.18	#ifdef ALLOW_MONITOR
29	jmc	1.27	# include "MONITOR.h"
30	edhill	1.18	#endif
31	cnh	1.1
32	cnh	1.12	C !INPUT/OUTPUT PARAMETERS:
33	cnh	1.1	INTEGER myThid
34	edhill	1.18	CEOP
35	cnh	1.1
36	jmc	1.26	C === Functions ====
37			LOGICAL MASTER_CPU_IO
38			EXTERNAL MASTER_CPU_IO
39
40	cnh	1.12	C !LOCAL VARIABLES:
41	jmc	1.23	C bi,bj :: Tile Indices counters
42			C i, j :: Loop counters
43			C facGrid :: Factor for grid to meter conversion
44			INTEGER bi,bj
45			INTEGER i, j
46	cnh	1.1	_RL facGrid
47	jmc	1.27	#ifndef OLD_GRID_IO
48			INTEGER myTile, iG, iLen
49			CHARACTER*(MAX_LEN_FNAM) fName
50			CHARACTER*(MAX_LEN_MBUF) msgBuf
51			INTEGER ILNBLNK
52			EXTERNAL ILNBLNK
53			#endif
54
55	cnh	1.1
56	edhill	1.18	C Initialise coriolis parameter
57	cnh	1.3	IF ( useConstantF ) THEN
58	edhill	1.18	C Constant F case
59			DO bj = myByLo(myThid), myByHi(myThid)
60			DO bi = myBxLo(myThid), myBxHi(myThid)
61	jmc	1.23	DO j=1-Oly,sNy+Oly
62			DO i=1-Olx,sNx+Olx
63			fCori(i,j,bi,bj) = f0
64			fCoriG(i,j,bi,bj) = f0
65			fCoriCos(i,j,bi,bj)=0. _d 0
66	edhill	1.18	ENDDO
67	jmc	1.23	ENDDO
68	cnh	1.1	ENDDO
69			ENDDO
70	cnh	1.3	ELSEIF ( useBetaPlaneF ) THEN
71	edhill	1.18	C Beta plane case
72			facGrid = 1. _d 0
73	jmc	1.23	IF ( usingSphericalPolarGrid
74			& .OR. usingCurvilinearGrid ) facGrid = deg2rad*rSphere
75	edhill	1.18	DO bj = myByLo(myThid), myByHi(myThid)
76			DO bi = myBxLo(myThid), myBxHi(myThid)
77	jmc	1.23	DO j=1-Oly,sNy+Oly
78			DO i=1-Olx,sNx+Olx
79			fCori(i,j,bi,bj) = f0+beta_yC(i,j,bi,bj)facGrid
80			fCoriG(i,j,bi,bj) = f0+beta* yG(i,j,bi,bj)*facGrid
81			fCoriCos(i,j,bi,bj)=0. _d 0
82	edhill	1.18	ENDDO
83	jmc	1.23	ENDDO
84	cnh	1.3	ENDDO
85			ENDDO
86			ELSEIF ( useSphereF ) THEN
87	edhill	1.18	C Spherical case
88			C Note in this case we assume yC is in degrees.
89			DO bj = myByLo(myThid), myByHi(myThid)
90			DO bi = myBxLo(myThid), myBxHi(myThid)
91	jmc	1.23	DO j=1-Oly,sNy+Oly
92			DO i=1-Olx,sNx+Olx
93			fCori(i,j,bi,bj) =
94	edhill	1.18	& 2. _d 0omegasin(_yC(i,j,bi,bj)*deg2rad)
95	jmc	1.23	fCoriG(i,j,bi,bj) =
96	edhill	1.18	& 2. _d 0omegasin(yG(i,j,bi,bj)*deg2rad)
97	jmc	1.23	fCoriCos(i,j,bi,bj)=
98	edhill	1.18	& 2. _d 0omegacos(_yC(i,j,bi,bj)*deg2rad)
99			ENDDO
100	jmc	1.23	ENDDO
101	cnh	1.3	ENDDO
102			ENDDO
103	jmc	1.21	c CALL WRITE_FLD_XY_RL('fCoriC',' ',fCori , 0,myThid)
104			c CALL WRITE_FLD_XY_RL('fCoriG',' ',fCoriG , 0,myThid)
105			c CALL WRITE_FLD_XY_RL('fCorCs',' ',fCoriCos,0,myThid)
106	cnh	1.3	ELSE
107	edhill	1.18	C Special custom form
108			DO bj = myByLo(myThid), myByHi(myThid)
109			DO bi = myBxLo(myThid), myBxHi(myThid)
110	jmc	1.23	DO j=1-Oly,sNy+Oly
111			DO i=1-Olx,sNx+Olx
112			fCori(i,j,bi,bj) = 0. _d 0
113			fCoriG(i,j,bi,bj) = 0. _d 0
114			fCoriCos(i,j,bi,bj)=0. _d 0
115	edhill	1.18	ENDDO
116	jmc	1.23	ENDDO
117	cnh	1.3	ENDDO
118			ENDDO
119	jmc	1.26	_BARRIER
120	jmc	1.21	CALL READ_REC_XY_RS( 'fCoriC.bin', fCori, 1, 0, myThid )
121			CALL READ_REC_XY_RS( 'fCorCs.bin', fCoriCos,1, 0, myThid )
122	jmc	1.27	IF ( .NOT.useCubedSphereExchange ) THEN
123			CALL READ_REC_XY_RS('fCoriG.bin', fCoriG, 1, 0, myThid )
124			ELSE
125			#ifdef OLD_GRID_IO
126			CALL READ_REC_XY_RS('fCoriG.bin', fCoriG, 1, 0, myThid )
127	jmc	1.26	_BARRIER
128	jmc	1.21	C- deal with the 2 missing corners (for fCoriG):
129			DO bj = myByLo(myThid), myByHi(myThid)
130			DO bi = myBxLo(myThid), myBxHi(myThid)
131			C- Notes: this will only works with 6 tiles (1 per face) and
132			C with 2 polar faces + 4 equatorials:
133			IF (bi.LE.3 .OR. bi.GE.5) THEN
134			fCoriG(sNx+1,1,bi,bj) = fCoriG(1,1,bi,bj)
135			ELSE
136			fCoriG(sNx+1,1,bi,bj) = -fCoriG(1,1,bi,bj)
137			ENDIF
138			IF (bi.GE.3) THEN
139			fCoriG(1,sNy+1,bi,bj) = fCoriG(1,1,bi,bj)
140			fCoriG(sNx+1,sNy+1,bi,bj) = fCoriG(sNx+1,1,bi,bj)
141			ELSE
142			fCoriG(1,sNy+1,bi,bj) = -fCoriG(1,1,bi,bj)
143			fCoriG(sNx+1,sNy+1,bi,bj) = -fCoriG(sNx+1,1,bi,bj)
144			ENDIF
145			ENDDO
146			ENDDO
147	jmc	1.27	#else /* OLD_GRID_IO */
148			_BEGIN_MASTER(myThid)
149			DO bj = 1,nSy
150			DO bi = 1,nSx
151			iG = bi+(myXGlobalLo-1)/sNx
152			myTile = iG
153			#ifdef ALLOW_EXCH2
154			myTile = W2_myTileList(bi)
155			iG = exch2_myface(myTile)
156			#endif
157			WRITE(fName,'(2A,I3.3,A)') 'fCoriG','.face',iG,'.bin'
158			iLen = ILNBLNK(fName)
159			WRITE(msgBuf,'(A,I6,2A)')
160			& ' Reading tile:', myTile, ' from file ', fName(1:iLen)
161			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
162			& SQUEEZE_RIGHT , myThid )
163			#ifdef ALLOW_MDSIO
164			CALL MDS_FACEF_READ_RS( fName, readBinaryPrec, 1,
165			& fCoriG, bi, bj, myThid )
166			#else /* ALLOW_MDSIO */
167			WRITE(msgBuf,'(A)') 'INI_CORI: Needs to compile MDSIO pkg'
168			CALL PRINT_ERROR( msgBuf, myThid )
169			STOP 'ABNORMAL END: S/R INI_CORI'
170			#endif /* ALLOW_MDSIO */
171
172			ENDDO
173			ENDDO
174			_END_MASTER(myThid)
175			#endif /* OLD_GRID_IO */
176	jmc	1.21	ENDIF
177
178	jmc	1.24	CALL EXCH_XY_RS( fCori, myThid )
179			CALL EXCH_XY_RS( fCoriCos, myThid )
180			CALL EXCH_Z_3D_RS( fCoriG, 1, myThid )
181	cnh	1.3	ENDIF
182	adcroft	1.11
183	edhill	1.16	#ifdef ALLOW_MONITOR
184	jmc	1.26	IF ( MASTER_CPU_IO(myThid) ) THEN
185	jmc	1.23	C-- only the master thread is allowed to switch On/Off mon_write_stdout
186			C & mon_write_mnc (since it's the only thread that uses those flags):
187
188			IF (monitor_stdio) THEN
189			mon_write_stdout = .TRUE.
190			ELSE
191			mon_write_stdout = .FALSE.
192			ENDIF
193			mon_write_mnc = .FALSE.
194	edhill	1.18	#ifdef ALLOW_MNC
195	jmc	1.23	IF (useMNC .AND. monitor_mnc) THEN
196			DO i = 1,MAX_LEN_MBUF
197			mon_fname(i:i) = ' '
198			ENDDO
199			mon_fname(1:12) = 'monitor_grid'
200			CALL MNC_CW_SET_UDIM(mon_fname, 1, myThid)
201			mon_write_mnc = .TRUE.
202			ENDIF
203			#endif /* ALLOW_MNC */
204
205	edhill	1.18	ENDIF
206	jmc	1.23
207	adcroft	1.11	CALL MON_PRINTSTATS_RS(1,fCori,'fCori',myThid)
208			CALL MON_PRINTSTATS_RS(1,fCoriG,'fCoriG',myThid)
209	adcroft	1.15	CALL MON_PRINTSTATS_RS(1,fCoriCos,'fCoriCos',myThid)
210	edhill	1.18
211	jmc	1.26	IF ( MASTER_CPU_IO(myThid) ) THEN
212	jmc	1.23	mon_write_stdout = .FALSE.
213			mon_write_mnc = .FALSE.
214			ENDIF
215			#endif /* ALLOW_MONITOR */
216	cnh	1.1
217			RETURN
218			END