pkg/seaice/seaice_init_fixed.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_fixed.F,v 1.5 2009/06/25 14:36:15 heimbach Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_INIT_FIXED( myThid )
C     *==========================================================*
C     | SUBROUTINE SEAICE_INIT_FIXED
C     | o Initialization of sea ice model.
C     *==========================================================*
C     *==========================================================*
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SEAICE.h"
CML#include "SEAICE_GRID.h"
#include "SEAICE_DIAGS.h"
#include "SEAICE_PARAMS.h"
#include "FFIELDS.h"

C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid
CEndOfInterface

C     === Local variables ===
C     i,j,k,bi,bj - Loop counters

      INTEGER i, j, k, bi, bj
      INTEGER kSurface
#ifndef SEAICE_CGRID
      _RS  mask_uice
#endif
cif(
cif   Helper variable for determining the fraction of sw radiation
cif   penetrating the model's shallowest layer
      INTEGER dummyIter
      _RL dummyTime
      _RL swfracba(2)
      _RL FACTORM
      INTEGER IMAX
cif)

      IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
       kSurface        = Nr
      ELSE
       kSurface        = 1
      ENDIF

C     Initialize MNC variable information for SEAICE
      IF ( useMNC .AND.
     &    (seaice_tave_mnc.OR.seaice_dump_mnc.OR.SEAICE_mon_mnc)
     &   ) THEN
        CALL SEAICE_MNC_INIT( myThid )
      ENDIF

cif(
#ifdef SHORTWAVE_HEATING
      IMAX      = 2
      FACTORM   = -1.0
      dummyTime = 1.0
      dummyIter = 0
      swfracba(1) = abs(rF(1))
      swfracba(2) = abs(rF(2))
      CALL SWFRAC(
     I       IMAX,FACTORM,
     U       swfracba,
     I       dummyTime,dummyIter,myThid)
      SWFRACB = swfracba(2)
#endif
cif)

C--   Initialize grid info
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          HEFFM(i,j,bi,bj)       = 0. _d 0
         ENDDO
        ENDDO
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          HEFFM(i,j,bi,bj)= 1. _d 0
          IF (_hFacC(i,j,kSurface,bi,bj).eq.0.)
     &         HEFFM(i,j,bi,bj)= 0. _d 0
         ENDDO
        ENDDO
        DO j=1-OLy+1,sNy+OLy
         DO i=1-OLx+1,sNx+OLx
#ifndef SEAICE_CGRID
          UVM(i,j,bi,bj)=0. _d 0
          mask_uice=HEFFM(i,j,  bi,bj)+HEFFM(i-1,j-1,bi,bj)
     &             +HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j,  bi,bj)
          IF(mask_uice.GT.3.5 _d 0) UVM(i,j,bi,bj)=1. _d 0
#endif /* SEAICE_CGRID */
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef SEAICE_CGRID
C     coefficients for metric terms
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          k1AtC(I,J,bi,bj) = 0.0 _d 0
          k1AtZ(I,J,bi,bj) = 0.0 _d 0
          k2AtC(I,J,bi,bj) = 0.0 _d 0
          k2AtZ(I,J,bi,bj) = 0.0 _d 0
         ENDDO
        ENDDO
        IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
C     This is the only case where tan(phi) is not zero. In this case
C     C and U points, and Z and V points have the same phi, so that we
C     only need a copy here. Do not use tan(YC) and tan(YG), because these
C     can be the geographical coordinates and not the correct grid
C     coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
           k2AtZ(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
          ENDDO
         ENDDO
        ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
C     compute metric term coefficients from finite difference approximation
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx-1
           k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
     &          * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
     &          * _recip_dxF(I,J,bi,bj)
          ENDDO
         ENDDO
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx+1,sNx+OLx
           k1AtZ(I,J,bi,bj) = _recip_dyU(I,J,bi,bj)
     &          * ( _dyC(I,J,bi,bj) - _dyC(I-1,J,bi,bj) )
     &          * _recip_dxV(I,J,bi,bj)
          ENDDO
         ENDDO
         DO j=1-OLy,sNy+OLy-1
          DO i=1-OLx,sNx+OLx
           k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
     &          * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
     &          * _recip_dyF(I,J,bi,bj)
          ENDDO
         ENDDO
         DO j=1-OLy+1,sNy+OLy
          DO i=1-OLx,sNx+OLx
           k2AtC(I,J,bi,bj) = _recip_dxV(I,J,bi,bj)
     &          * ( _dxC(I,J,bi,bj) - _dxC(I,J-1,bi,bj) )
     &          * _recip_dyU(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDIF
       ENDDO
      ENDDO
#endif /* SEAICE_CGRID */

#ifndef SEAICE_CGRID
C--   Choose a proxy level for geostrophic velocity,
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          KGEO(i,j,bi,bj)   = 0
         ENDDO
        ENDDO
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
#ifdef SEAICE_BICE_STRESS
          KGEO(i,j,bi,bj) = 1
#else /* SEAICE_BICE_STRESS */
          IF (klowc(i,j,bi,bj) .LT. 2) THEN
           KGEO(i,j,bi,bj) = 1
          ELSE
           KGEO(i,j,bi,bj) = 2
           DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND.
     &          KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
              KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
           ENDDO
          ENDIF
#endif /* SEAICE_BICE_STRESS */
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#endif /* SEAICE_CGRID */

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
        CALL SEAICE_DIAGNOSTICS_INIT( myThid )
      ENDIF
#endif

#ifdef ALLOW_TIMEAVE
C     Initialize averages to zero
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        CALL TIMEAVE_RESET(FUtave   ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(FVtave   ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(EmPmRtave,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(QNETtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(QSWtave  ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(UICEtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(VICEtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(HEFFtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(AREAtave ,1,bi,bj,myThid)
        DO k=1,Nr
         SEAICE_TimeAve(k,bi,bj)=ZERO
        ENDDO
       ENDDO
      ENDDO
#endif /* ALLOW_TIMEAVE */

      RETURN
      END
1	jmc	1.6	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_fixed.F,v 1.5 2009/06/25 14:36:15 heimbach Exp $
2	heimbach	1.1	C $Name: $
3
4			#include "SEAICE_OPTIONS.h"
5
6			CStartOfInterface
7			SUBROUTINE SEAICE_INIT_FIXED( myThid )
8	jmc	1.6	C ==========================================================
9			C \| SUBROUTINE SEAICE_INIT_FIXED
10			C \| o Initialization of sea ice model.
11			C ==========================================================
12			C ==========================================================
13	heimbach	1.1	IMPLICIT NONE
14	jmc	1.6
15	heimbach	1.1	C === Global variables ===
16			#include "SIZE.h"
17			#include "EEPARAMS.h"
18			#include "PARAMS.h"
19			#include "GRID.h"
20			#include "SEAICE.h"
21			CML#include "SEAICE_GRID.h"
22			#include "SEAICE_DIAGS.h"
23			#include "SEAICE_PARAMS.h"
24			#include "FFIELDS.h"
25
26			C === Routine arguments ===
27			C myThid - Thread no. that called this routine.
28			INTEGER myThid
29			CEndOfInterface
30	jmc	1.6
31	heimbach	1.1	C === Local variables ===
32			C i,j,k,bi,bj - Loop counters
33
34	mlosch	1.4	INTEGER i, j, k, bi, bj
35			INTEGER kSurface
36	jmc	1.6	#ifndef SEAICE_CGRID
37	mlosch	1.4	_RS mask_uice
38	jmc	1.6	#endif
39	heimbach	1.3	cif(
40			cif Helper variable for determining the fraction of sw radiation
41			cif penetrating the model's shallowest layer
42			INTEGER dummyIter
43			_RL dummyTime
44			_RL swfracba(2)
45			_RL FACTORM
46			INTEGER IMAX
47			cif)
48	heimbach	1.1
49	mlosch	1.4	IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
50			kSurface = Nr
51			ELSE
52			kSurface = 1
53	heimbach	1.1	ENDIF
54
55	jmc	1.6	C Initialize MNC variable information for SEAICE
56			IF ( useMNC .AND.
57			& (seaice_tave_mnc.OR.seaice_dump_mnc.OR.SEAICE_mon_mnc)
58			& ) THEN
59			CALL SEAICE_MNC_INIT( myThid )
60			ENDIF
61
62	heimbach	1.3	cif(
63			#ifdef SHORTWAVE_HEATING
64			IMAX = 2
65			FACTORM = -1.0
66			dummyTime = 1.0
67			dummyIter = 0
68	jmc	1.6	swfracba(1) = abs(rF(1))
69	heimbach	1.3	swfracba(2) = abs(rF(2))
70			CALL SWFRAC(
71			I IMAX,FACTORM,
72			U swfracba,
73			I dummyTime,dummyIter,myThid)
74			SWFRACB = swfracba(2)
75			#endif
76			cif)
77
78	mlosch	1.4	C-- Initialize grid info
79			DO bj=myByLo(myThid),myByHi(myThid)
80			DO bi=myBxLo(myThid),myBxHi(myThid)
81			DO j=1-OLy,sNy+OLy
82			DO i=1-OLx,sNx+OLx
83			HEFFM(i,j,bi,bj) = 0. _d 0
84			ENDDO
85			ENDDO
86			DO j=1-OLy,sNy+OLy
87			DO i=1-OLx,sNx+OLx
88			HEFFM(i,j,bi,bj)= 1. _d 0
89			IF (_hFacC(i,j,kSurface,bi,bj).eq.0.)
90			& HEFFM(i,j,bi,bj)= 0. _d 0
91			ENDDO
92			ENDDO
93			DO j=1-OLy+1,sNy+OLy
94			DO i=1-OLx+1,sNx+OLx
95	heimbach	1.5	#ifndef SEAICE_CGRID
96	mlosch	1.4	UVM(i,j,bi,bj)=0. _d 0
97			mask_uice=HEFFM(i,j, bi,bj)+HEFFM(i-1,j-1,bi,bj)
98			& +HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j, bi,bj)
99			IF(mask_uice.GT.3.5 _d 0) UVM(i,j,bi,bj)=1. _d 0
100			#endif /* SEAICE_CGRID */
101			ENDDO
102			ENDDO
103			ENDDO
104			ENDDO
105
106			#ifdef SEAICE_CGRID
107	jmc	1.6	C coefficients for metric terms
108	mlosch	1.4	DO bj=myByLo(myThid),myByHi(myThid)
109			DO bi=myBxLo(myThid),myBxHi(myThid)
110			DO j=1-OLy,sNy+OLy
111			DO i=1-OLx,sNx+OLx
112			k1AtC(I,J,bi,bj) = 0.0 _d 0
113			k1AtZ(I,J,bi,bj) = 0.0 _d 0
114			k2AtC(I,J,bi,bj) = 0.0 _d 0
115			k2AtZ(I,J,bi,bj) = 0.0 _d 0
116			ENDDO
117			ENDDO
118			IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
119	jmc	1.6	C This is the only case where tan(phi) is not zero. In this case
120			C C and U points, and Z and V points have the same phi, so that we
121	mlosch	1.4	C only need a copy here. Do not use tan(YC) and tan(YG), because these
122	jmc	1.6	C can be the geographical coordinates and not the correct grid
123	mlosch	1.4	C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
124			DO j=1-OLy,sNy+OLy
125			DO i=1-OLx,sNx+OLx
126			k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
127			k2AtZ(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
128			ENDDO
129			ENDDO
130			ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
131			C compute metric term coefficients from finite difference approximation
132			DO j=1-OLy,sNy+OLy
133			DO i=1-OLx,sNx+OLx-1
134			k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
135			& * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
136			& * _recip_dxF(I,J,bi,bj)
137			ENDDO
138			ENDDO
139			DO j=1-OLy,sNy+OLy
140			DO i=1-OLx+1,sNx+OLx
141	jmc	1.6	k1AtZ(I,J,bi,bj) = _recip_dyU(I,J,bi,bj)
142	mlosch	1.4	& * ( _dyC(I,J,bi,bj) - _dyC(I-1,J,bi,bj) )
143			& * _recip_dxV(I,J,bi,bj)
144			ENDDO
145			ENDDO
146			DO j=1-OLy,sNy+OLy-1
147			DO i=1-OLx,sNx+OLx
148	jmc	1.6	k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
149	mlosch	1.4	& * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
150			& * _recip_dyF(I,J,bi,bj)
151			ENDDO
152			ENDDO
153			DO j=1-OLy+1,sNy+OLy
154			DO i=1-OLx,sNx+OLx
155	jmc	1.6	k2AtC(I,J,bi,bj) = _recip_dxV(I,J,bi,bj)
156	mlosch	1.4	& * ( _dxC(I,J,bi,bj) - _dxC(I,J-1,bi,bj) )
157			& * _recip_dyU(I,J,bi,bj)
158			ENDDO
159			ENDDO
160			ENDIF
161			ENDDO
162			ENDDO
163			#endif /* SEAICE_CGRID */
164
165			#ifndef SEAICE_CGRID
166			C-- Choose a proxy level for geostrophic velocity,
167			DO bj=myByLo(myThid),myByHi(myThid)
168			DO bi=myBxLo(myThid),myBxHi(myThid)
169			DO j=1-OLy,sNy+OLy
170			DO i=1-OLx,sNx+OLx
171			KGEO(i,j,bi,bj) = 0
172			ENDDO
173			ENDDO
174			DO j=1-OLy,sNy+OLy
175			DO i=1-OLx,sNx+OLx
176			#ifdef SEAICE_BICE_STRESS
177			KGEO(i,j,bi,bj) = 1
178			#else /* SEAICE_BICE_STRESS */
179			IF (klowc(i,j,bi,bj) .LT. 2) THEN
180			KGEO(i,j,bi,bj) = 1
181			ELSE
182			KGEO(i,j,bi,bj) = 2
183			DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND.
184			& KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
185			KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
186			ENDDO
187			ENDIF
188			#endif /* SEAICE_BICE_STRESS */
189			ENDDO
190			ENDDO
191			ENDDO
192			ENDDO
193			#endif /* SEAICE_CGRID */
194
195			#ifdef ALLOW_DIAGNOSTICS
196			IF ( useDiagnostics ) THEN
197			CALL SEAICE_DIAGNOSTICS_INIT( myThid )
198			ENDIF
199			#endif
200
201	heimbach	1.1	#ifdef ALLOW_TIMEAVE
202			C Initialize averages to zero
203			DO bj = myByLo(myThid), myByHi(myThid)
204			DO bi = myBxLo(myThid), myBxHi(myThid)
205			CALL TIMEAVE_RESET(FUtave ,1,bi,bj,myThid)
206			CALL TIMEAVE_RESET(FVtave ,1,bi,bj,myThid)
207			CALL TIMEAVE_RESET(EmPmRtave,1,bi,bj,myThid)
208			CALL TIMEAVE_RESET(QNETtave ,1,bi,bj,myThid)
209			CALL TIMEAVE_RESET(QSWtave ,1,bi,bj,myThid)
210			CALL TIMEAVE_RESET(UICEtave ,1,bi,bj,myThid)
211			CALL TIMEAVE_RESET(VICEtave ,1,bi,bj,myThid)
212			CALL TIMEAVE_RESET(HEFFtave ,1,bi,bj,myThid)
213			CALL TIMEAVE_RESET(AREAtave ,1,bi,bj,myThid)
214			DO k=1,Nr
215			SEAICE_TimeAve(k,bi,bj)=ZERO
216			ENDDO
217			ENDDO
218			ENDDO
219			#endif /* ALLOW_TIMEAVE */
220
221			RETURN
222			END