model/src/ini_cori.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cori.F,v 1.20 2004/10/13 04:37:37 edhill 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"
#ifdef ALLOW_MNC
#include "MNC_PARAMS.h"
#endif
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_MONITOR
#include "MONITOR.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
C     bi,bj  - Loop counters
C     I,J,K
C     facGrid - Factor for grid to meter conversion
      INTEGER bi, bj
      INTEGER  I,  J, K
      _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)
c           DO K=1,Nr
              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.
                ENDDO
              ENDDO
c           ENDDO
          ENDDO
        ENDDO
      ELSEIF ( useBetaPlaneF ) THEN
C       Beta plane case
        facGrid = 1. _d 0
        IF ( usingSphericalPolarGrid ) facGrid = deg2rad*rSphere
        DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
c           DO K=1,Nr
              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.
                ENDDO
              ENDDO
c           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)
c           DO K=1,Nr
              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
c           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)
c           DO K=1,Nr
              DO J=1-Oly,sNy+Oly
                DO I=1-Olx,sNx+Olx
                  fCori(i,j,bi,bj)=0.
                  fCoriG(i,j,bi,bj)=0.
                  fCoriCos(i,j,bi,bj)=0.
                ENDDO
              ENDDO
c           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

        _EXCH_XY_R4(fCori,myThid)
        CALL EXCH_Z_XY_RS(fCoriG,myThid)
        _EXCH_XY_R4(fCoriCos,myThid)
      ENDIF

#ifdef ALLOW_MONITOR
      mon_write_stdout = .FALSE.
      mon_write_mnc    = .FALSE.
      IF (monitor_stdio) THEN
        mon_write_stdout = .TRUE.
      ENDIF
      
#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  */
      
      CALL MON_PRINTSTATS_RS(1,fCori,'fCori',myThid)
      CALL MON_PRINTSTATS_RS(1,fCoriG,'fCoriG',myThid)
      CALL MON_PRINTSTATS_RS(1,fCoriCos,'fCoriCos',myThid)

      mon_write_stdout = .FALSE.
      mon_write_mnc    = .FALSE.
#endif

      RETURN
      END
1	jmc	1.21	C $Header: /u/gcmpack/MITgcm/model/src/ini_cori.F,v 1.20 2004/10/13 04:37:37 edhill 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			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	edhill	1.20	#ifdef ALLOW_MNC
21			#include "MNC_PARAMS.h"
22			#endif
23	cnh	1.1	#include "GRID.h"
24			#include "DYNVARS.h"
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	cnh	1.1	C bi,bj - Loop counters
35			C I,J,K
36			C facGrid - Factor for grid to meter conversion
37			INTEGER bi, bj
38			INTEGER I, J, K
39			_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.21	c DO K=1,Nr
47	edhill	1.18	DO J=1-Oly,sNy+Oly
48			DO I=1-Olx,sNx+Olx
49			fCori(i,j,bi,bj)=f0
50			fCoriG(i,j,bi,bj)=f0
51			fCoriCos(i,j,bi,bj)=0.
52			ENDDO
53			ENDDO
54	jmc	1.21	c ENDDO
55	cnh	1.1	ENDDO
56			ENDDO
57	cnh	1.3	ELSEIF ( useBetaPlaneF ) THEN
58	edhill	1.18	C Beta plane case
59			facGrid = 1. _d 0
60			IF ( usingSphericalPolarGrid ) facGrid = deg2rad*rSphere
61			DO bj = myByLo(myThid), myByHi(myThid)
62			DO bi = myBxLo(myThid), myBxHi(myThid)
63	jmc	1.21	c DO K=1,Nr
64	edhill	1.18	DO J=1-Oly,sNy+Oly
65			DO I=1-Olx,sNx+Olx
66			fCori(i,j,bi,bj)=f0+beta_yC(i,j,bi,bj)facGrid
67			fCoriG(i,j,bi,bj)=f0+betayG(i,j,bi,bj)facGrid
68			fCoriCos(i,j,bi,bj)=0.
69			ENDDO
70			ENDDO
71	jmc	1.21	c ENDDO
72	cnh	1.3	ENDDO
73			ENDDO
74			ELSEIF ( useSphereF ) THEN
75	edhill	1.18	C Spherical case
76			C Note in this case we assume yC is in degrees.
77			DO bj = myByLo(myThid), myByHi(myThid)
78			DO bi = myBxLo(myThid), myBxHi(myThid)
79	jmc	1.21	c DO K=1,Nr
80	edhill	1.18	DO J=1-Oly,sNy+Oly
81			DO I=1-Olx,sNx+Olx
82			fCori(i,j,bi,bj)=
83			& 2. _d 0omegasin(_yC(i,j,bi,bj)*deg2rad)
84			fCoriG(i,j,bi,bj)=
85			& 2. _d 0omegasin(yG(i,j,bi,bj)*deg2rad)
86			fCoriCos(i,j,bi,bj)=
87			& 2. _d 0omegacos(_yC(i,j,bi,bj)*deg2rad)
88			ENDDO
89			ENDDO
90	jmc	1.21	c ENDDO
91	cnh	1.3	ENDDO
92			ENDDO
93	jmc	1.21	c CALL WRITE_FLD_XY_RL('fCoriC',' ',fCori , 0,myThid)
94			c CALL WRITE_FLD_XY_RL('fCoriG',' ',fCoriG , 0,myThid)
95			c CALL WRITE_FLD_XY_RL('fCorCs',' ',fCoriCos,0,myThid)
96	cnh	1.3	ELSE
97	edhill	1.18	C Special custom form
98			DO bj = myByLo(myThid), myByHi(myThid)
99			DO bi = myBxLo(myThid), myBxHi(myThid)
100	jmc	1.21	c DO K=1,Nr
101	edhill	1.18	DO J=1-Oly,sNy+Oly
102			DO I=1-Olx,sNx+Olx
103			fCori(i,j,bi,bj)=0.
104			fCoriG(i,j,bi,bj)=0.
105			fCoriCos(i,j,bi,bj)=0.
106			ENDDO
107			ENDDO
108	jmc	1.21	c ENDDO
109	cnh	1.3	ENDDO
110			ENDDO
111	jmc	1.21	CALL READ_REC_XY_RS( 'fCoriC.bin', fCori, 1, 0, myThid )
112			CALL READ_REC_XY_RS( 'fCoriG.bin', fCoriG, 1, 0, myThid )
113			CALL READ_REC_XY_RS( 'fCorCs.bin', fCoriCos,1, 0, myThid )
114			IF ( useCubedSphereExchange ) THEN
115			C- deal with the 2 missing corners (for fCoriG):
116			DO bj = myByLo(myThid), myByHi(myThid)
117			DO bi = myBxLo(myThid), myBxHi(myThid)
118			C- Notes: this will only works with 6 tiles (1 per face) and
119			C with 2 polar faces + 4 equatorials:
120			IF (bi.LE.3 .OR. bi.GE.5) THEN
121			fCoriG(sNx+1,1,bi,bj) = fCoriG(1,1,bi,bj)
122			ELSE
123			fCoriG(sNx+1,1,bi,bj) = -fCoriG(1,1,bi,bj)
124			ENDIF
125			IF (bi.GE.3) THEN
126			fCoriG(1,sNy+1,bi,bj) = fCoriG(1,1,bi,bj)
127			fCoriG(sNx+1,sNy+1,bi,bj) = fCoriG(sNx+1,1,bi,bj)
128			ELSE
129			fCoriG(1,sNy+1,bi,bj) = -fCoriG(1,1,bi,bj)
130			fCoriG(sNx+1,sNy+1,bi,bj) = -fCoriG(sNx+1,1,bi,bj)
131			ENDIF
132			ENDDO
133			ENDDO
134			ENDIF
135
136			_EXCH_XY_R4(fCori,myThid)
137			CALL EXCH_Z_XY_RS(fCoriG,myThid)
138			_EXCH_XY_R4(fCoriCos,myThid)
139	cnh	1.3	ENDIF
140	adcroft	1.11
141	edhill	1.16	#ifdef ALLOW_MONITOR
142	edhill	1.18	mon_write_stdout = .FALSE.
143			mon_write_mnc = .FALSE.
144			IF (monitor_stdio) THEN
145			mon_write_stdout = .TRUE.
146			ENDIF
147
148			#ifdef ALLOW_MNC
149			IF (useMNC .AND. monitor_mnc) THEN
150			DO i = 1,MAX_LEN_MBUF
151			mon_fname(i:i) = ' '
152			ENDDO
153			mon_fname(1:12) = 'monitor_grid'
154			CALL MNC_CW_SET_UDIM(mon_fname, 1, myThid)
155			mon_write_mnc = .TRUE.
156			ENDIF
157			#endif /* ALLOW_MNC */
158
159	adcroft	1.11	CALL MON_PRINTSTATS_RS(1,fCori,'fCori',myThid)
160			CALL MON_PRINTSTATS_RS(1,fCoriG,'fCoriG',myThid)
161	adcroft	1.15	CALL MON_PRINTSTATS_RS(1,fCoriCos,'fCoriCos',myThid)
162	edhill	1.18
163			mon_write_stdout = .FALSE.
164			mon_write_mnc = .FALSE.
165	heimbach	1.13	#endif
166	cnh	1.1
167			RETURN
168			END