model/src/ini_cori.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cori.F,v 1.23 2006/07/13 03:00:24 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 "EESUPPORT.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_MNC
#include "MNC_PARAMS.h"
#endif
#ifdef ALLOW_MONITOR
#include "MONITOR.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
      INTEGER myThid
CEOP

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

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
        CALL READ_REC_XY_RS( 'fCoriC.bin', fCori,   1, 0, myThid )
        CALL READ_REC_XY_RS( 'fCoriG.bin', fCoriG,  1, 0, myThid )
        CALL READ_REC_XY_RS( 'fCorCs.bin', fCoriCos,1, 0, myThid )
        IF ( useCubedSphereExchange ) THEN
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
        ENDIF

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

#ifdef ALLOW_MONITOR
#ifdef ALLOW_USE_MPI
      IF ( .NOT.useSingleCPUIO .OR. mpiMyId.EQ.0 ) THEN
#endif /* ALLOW_USE_MPI */
        _BEGIN_MASTER(myThid)
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  */

        _END_MASTER(myThid)
#ifdef ALLOW_USE_MPI
      ENDIF
#endif /* ALLOW_USE_MPI */

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

#ifdef ALLOW_USE_MPI
      IF ( .NOT.useSingleCPUIO .OR. mpiMyId.EQ.0 ) THEN
#endif /* ALLOW_USE_MPI */
        _BEGIN_MASTER(myThid)

        mon_write_stdout = .FALSE.
        mon_write_mnc    = .FALSE.

        _END_MASTER(myThid)
#ifdef ALLOW_USE_MPI
      ENDIF
#endif /* ALLOW_USE_MPI */
#endif /* ALLOW_MONITOR */

      RETURN
      END
1	jmc	1.24	C $Header: /u/gcmpack/MITgcm/model/src/ini_cori.F,v 1.23 2006/07/13 03:00:24 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	jmc	1.23	#include "EESUPPORT.h"
20	cnh	1.1	#include "PARAMS.h"
21	jmc	1.23	#include "GRID.h"
22	edhill	1.20	#ifdef ALLOW_MNC
23			#include "MNC_PARAMS.h"
24			#endif
25	edhill	1.18	#ifdef ALLOW_MONITOR
26			#include "MONITOR.h"
27			#endif
28	cnh	1.1
29	cnh	1.12	C !INPUT/OUTPUT PARAMETERS:
30	cnh	1.1	INTEGER myThid
31	edhill	1.18	CEOP
32	cnh	1.1
33	cnh	1.12	C !LOCAL VARIABLES:
34	jmc	1.23	C bi,bj :: Tile Indices counters
35			C i, j :: Loop counters
36			C facGrid :: Factor for grid to meter conversion
37			INTEGER bi,bj
38			INTEGER i, j
39	cnh	1.1	_RL facGrid
40
41	edhill	1.18	C Initialise coriolis parameter
42	cnh	1.3	IF ( useConstantF ) THEN
43	edhill	1.18	C Constant F case
44			DO bj = myByLo(myThid), myByHi(myThid)
45			DO bi = myBxLo(myThid), myBxHi(myThid)
46	jmc	1.23	DO j=1-Oly,sNy+Oly
47			DO i=1-Olx,sNx+Olx
48			fCori(i,j,bi,bj) = f0
49			fCoriG(i,j,bi,bj) = f0
50			fCoriCos(i,j,bi,bj)=0. _d 0
51	edhill	1.18	ENDDO
52	jmc	1.23	ENDDO
53	cnh	1.1	ENDDO
54			ENDDO
55	cnh	1.3	ELSEIF ( useBetaPlaneF ) THEN
56	edhill	1.18	C Beta plane case
57			facGrid = 1. _d 0
58	jmc	1.23	IF ( usingSphericalPolarGrid
59			& .OR. usingCurvilinearGrid ) facGrid = deg2rad*rSphere
60	edhill	1.18	DO bj = myByLo(myThid), myByHi(myThid)
61			DO bi = myBxLo(myThid), myBxHi(myThid)
62	jmc	1.23	DO j=1-Oly,sNy+Oly
63			DO i=1-Olx,sNx+Olx
64			fCori(i,j,bi,bj) = f0+beta_yC(i,j,bi,bj)facGrid
65			fCoriG(i,j,bi,bj) = f0+beta* yG(i,j,bi,bj)*facGrid
66			fCoriCos(i,j,bi,bj)=0. _d 0
67	edhill	1.18	ENDDO
68	jmc	1.23	ENDDO
69	cnh	1.3	ENDDO
70			ENDDO
71			ELSEIF ( useSphereF ) THEN
72	edhill	1.18	C Spherical case
73			C Note in this case we assume yC is in degrees.
74			DO bj = myByLo(myThid), myByHi(myThid)
75			DO bi = myBxLo(myThid), myBxHi(myThid)
76	jmc	1.23	DO j=1-Oly,sNy+Oly
77			DO i=1-Olx,sNx+Olx
78			fCori(i,j,bi,bj) =
79	edhill	1.18	& 2. _d 0omegasin(_yC(i,j,bi,bj)*deg2rad)
80	jmc	1.23	fCoriG(i,j,bi,bj) =
81	edhill	1.18	& 2. _d 0omegasin(yG(i,j,bi,bj)*deg2rad)
82	jmc	1.23	fCoriCos(i,j,bi,bj)=
83	edhill	1.18	& 2. _d 0omegacos(_yC(i,j,bi,bj)*deg2rad)
84			ENDDO
85	jmc	1.23	ENDDO
86	cnh	1.3	ENDDO
87			ENDDO
88	jmc	1.21	c CALL WRITE_FLD_XY_RL('fCoriC',' ',fCori , 0,myThid)
89			c CALL WRITE_FLD_XY_RL('fCoriG',' ',fCoriG , 0,myThid)
90			c CALL WRITE_FLD_XY_RL('fCorCs',' ',fCoriCos,0,myThid)
91	cnh	1.3	ELSE
92	edhill	1.18	C Special custom form
93			DO bj = myByLo(myThid), myByHi(myThid)
94			DO bi = myBxLo(myThid), myBxHi(myThid)
95	jmc	1.23	DO j=1-Oly,sNy+Oly
96			DO i=1-Olx,sNx+Olx
97			fCori(i,j,bi,bj) = 0. _d 0
98			fCoriG(i,j,bi,bj) = 0. _d 0
99			fCoriCos(i,j,bi,bj)=0. _d 0
100	edhill	1.18	ENDDO
101	jmc	1.23	ENDDO
102	cnh	1.3	ENDDO
103			ENDDO
104	jmc	1.21	CALL READ_REC_XY_RS( 'fCoriC.bin', fCori, 1, 0, myThid )
105			CALL READ_REC_XY_RS( 'fCoriG.bin', fCoriG, 1, 0, myThid )
106			CALL READ_REC_XY_RS( 'fCorCs.bin', fCoriCos,1, 0, myThid )
107			IF ( useCubedSphereExchange ) THEN
108			C- deal with the 2 missing corners (for fCoriG):
109			DO bj = myByLo(myThid), myByHi(myThid)
110			DO bi = myBxLo(myThid), myBxHi(myThid)
111			C- Notes: this will only works with 6 tiles (1 per face) and
112			C with 2 polar faces + 4 equatorials:
113			IF (bi.LE.3 .OR. bi.GE.5) THEN
114			fCoriG(sNx+1,1,bi,bj) = fCoriG(1,1,bi,bj)
115			ELSE
116			fCoriG(sNx+1,1,bi,bj) = -fCoriG(1,1,bi,bj)
117			ENDIF
118			IF (bi.GE.3) THEN
119			fCoriG(1,sNy+1,bi,bj) = fCoriG(1,1,bi,bj)
120			fCoriG(sNx+1,sNy+1,bi,bj) = fCoriG(sNx+1,1,bi,bj)
121			ELSE
122			fCoriG(1,sNy+1,bi,bj) = -fCoriG(1,1,bi,bj)
123			fCoriG(sNx+1,sNy+1,bi,bj) = -fCoriG(sNx+1,1,bi,bj)
124			ENDIF
125			ENDDO
126			ENDDO
127			ENDIF
128
129	jmc	1.24	CALL EXCH_XY_RS( fCori, myThid )
130			CALL EXCH_XY_RS( fCoriCos, myThid )
131			c CALL EXCH_Z_XY_RS( fCoriG, myThid )
132			CALL EXCH_Z_3D_RS( fCoriG, 1, myThid )
133	cnh	1.3	ENDIF
134	adcroft	1.11
135	edhill	1.16	#ifdef ALLOW_MONITOR
136	jmc	1.23	#ifdef ALLOW_USE_MPI
137			IF ( .NOT.useSingleCPUIO .OR. mpiMyId.EQ.0 ) THEN
138			#endif /* ALLOW_USE_MPI */
139			_BEGIN_MASTER(myThid)
140			C-- only the master thread is allowed to switch On/Off mon_write_stdout
141			C & mon_write_mnc (since it's the only thread that uses those flags):
142
143			IF (monitor_stdio) THEN
144			mon_write_stdout = .TRUE.
145			ELSE
146			mon_write_stdout = .FALSE.
147			ENDIF
148			mon_write_mnc = .FALSE.
149	edhill	1.18	#ifdef ALLOW_MNC
150	jmc	1.23	IF (useMNC .AND. monitor_mnc) THEN
151			DO i = 1,MAX_LEN_MBUF
152			mon_fname(i:i) = ' '
153			ENDDO
154			mon_fname(1:12) = 'monitor_grid'
155			CALL MNC_CW_SET_UDIM(mon_fname, 1, myThid)
156			mon_write_mnc = .TRUE.
157			ENDIF
158			#endif /* ALLOW_MNC */
159
160			_END_MASTER(myThid)
161			#ifdef ALLOW_USE_MPI
162	edhill	1.18	ENDIF
163	jmc	1.23	#endif /* ALLOW_USE_MPI */
164
165	adcroft	1.11	CALL MON_PRINTSTATS_RS(1,fCori,'fCori',myThid)
166			CALL MON_PRINTSTATS_RS(1,fCoriG,'fCoriG',myThid)
167	adcroft	1.15	CALL MON_PRINTSTATS_RS(1,fCoriCos,'fCoriCos',myThid)
168	edhill	1.18
169	jmc	1.23	#ifdef ALLOW_USE_MPI
170			IF ( .NOT.useSingleCPUIO .OR. mpiMyId.EQ.0 ) THEN
171			#endif /* ALLOW_USE_MPI */
172			_BEGIN_MASTER(myThid)
173
174			mon_write_stdout = .FALSE.
175			mon_write_mnc = .FALSE.
176
177			_END_MASTER(myThid)
178			#ifdef ALLOW_USE_MPI
179			ENDIF
180			#endif /* ALLOW_USE_MPI */
181			#endif /* ALLOW_MONITOR */
182	cnh	1.1
183			RETURN
184			END