pkg/thsice/thsice_ini_vars.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_ini_vars.F,v 1.12 2008/05/22 01:59:51 dimitri Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"
#ifdef ALLOW_AIM
# include "AIM_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: THSICE_INI_VARS
C     !INTERFACE:
      SUBROUTINE THSICE_INI_VARS( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R THSICE_INI_VARS
C     | o initialize THermo_SeaICE variables
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "DYNVARS.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"
#ifdef ALLOW_AIM
# include "AIM_FFIELDS.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid :: My Thread Id. number
      INTEGER myThid
CEOP

#ifdef ALLOW_THSICE
C     == Local variables ==
C     bi,bj  :: Loop counters
C     i,j    :: Loop counters
      INTEGER bi, bj
      INTEGER i, j
      INTEGER prec
c     CHARACTER*(MAX_LEN_FNAM) fn
      _RL v2Loc


c     set up ice arrays to zero if starting ice
      DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
C-        state variables :
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            iceMask(i,j,bi,bj)  = 0. _d 0
            iceHeight(i,j,bi,bj)= 0. _d 0
            snowHeight(i,j,bi,bj)=0. _d 0
            Tsrf(i,j,bi,bj)     = 0. _d 0
            Tice1(i,j,bi,bj)    = 0. _d 0
            Tice2(i,j,bi,bj)    = 0. _d 0
            Qice1(i,j,bi,bj)    = 0. _d 0
            Qice2(i,j,bi,bj)    = 0. _d 0
            snowAge(i,j,bi,bj)  = 0. _d 0
           ENDDO
          ENDDO
C-        fluxes :
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            sHeating(i,j,bi,bj) = 0. _d 0
            flxCndBt(i,j,bi,bj) = 0. _d 0
            snowPrc(i,j,bi,bj)  = 0. _d 0
            siceAlb(i,j,bi,bj)  = 0. _d 0
            icFlxSW (i,j,bi,bj) = 0. _d 0
            icFlxAtm(i,j,bi,bj) = 0. _d 0
            icFrwAtm(i,j,bi,bj) = 0. _d 0
C-        needed when using advection/diffusion:
            oceFWfx(i,j,bi,bj)  = 0. _d 0
            oceSflx(i,j,bi,bj)  = 0. _d 0
            oceQnet(i,j,bi,bj)  = 0. _d 0
           ENDDO
          ENDDO
C-        oceanic mixed layer state :
          v2Loc = vMxL_default*vMxL_default
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             hOceMxL(i,j,bi,bj) = hMxL_default
             tOceMxL(i,j,bi,bj) = 0. _d 0
             sOceMxL(i,j,bi,bj) = sMxL_default
             v2ocMxL(i,j,bi,bj) = v2Loc
           ENDDO
          ENDDO
#ifdef ALLOW_AIM
          IF ( useAIM ) THEN
C-        Mask mixed layer depth : depth is used in thsice slab_ocean
C         and this mask is used in thsice_advdiff and if coupled
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             IF ( aim_landFr(i,j,bi,bj).EQ.1. _d 0 )
     &       hOceMxL(i,j,bi,bj) = 0.
            ENDDO
           ENDDO
          ENDIF
#endif /* ALLOW_AIM */
        ENDDO
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF ( startIceModel.LE.0 .AND. nIter0.NE.0 ) THEN
C--     Read ice pickup fields
       _BARRIER
        prec=precFloat64
        CALL THSICE_READ_PICKUP( prec, nIter0, myThid )

      ELSE
C--     Read initial conditions:
        _BEGIN_MASTER( myThid )
        IF ( thSIceThick_InitFile .NE. ' ' ) THEN
         CALL READ_REC_XY_RL(thSIceThick_InitFile,iceHeight,1,0,myThid)
        ENDIF
        IF ( thSIceFract_InitFile .NE. ' ' ) THEN
         CALL READ_REC_XY_RL(thSIceFract_InitFile,iceMask,1,0,myThid)
        ELSE
C     default to 100% ice concentration when thickness is non-zero
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             IF (iceHeight(i,j,bi,bj) .NE. 0. _d 0)
     &            iceMask(i,j,bi,bj) = 1. _d 0
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDIF
        IF ( thSIceSnowH_InitFile .NE. ' ' ) THEN
         CALL READ_REC_XY_RL(thSIceSnowH_InitFile,snowHeight,1,0,myThid)
        ELSE
C     default to <= 20 cm snow thickness where there is ice
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             snowHeight(i,j,bi,bj) = 0.2 _d 0
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDIF
        IF ( thSIceSnowA_InitFile .NE. ' ' ) THEN
         CALL READ_REC_XY_RL(thSIceSnowA_InitFile,snowAge,1,0,myThid)
        ELSE
C     default to 3 days for snow age
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             IF (snowHeight(i,j,bi,bj) .NE. 0. _d 0)
     &            snowAge(i,j,bi,bj) = 259200. _d 0
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDIF
        IF ( thSIceEnthp_InitFile .NE. ' ' ) THEN
         CALL READ_REC_XY_RL(thSIceEnthp_InitFile,Qice1,1,0,myThid)
         CALL READ_REC_XY_RL(thSIceEnthp_InitFile,Qice2,2,0,myThid)
        ELSE
C     default to 3.4e5 J/kg for enthalpy where there is sea ice
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             IF (iceMask(i,j,bi,bj) .NE. 0. _d 0) THEN
              Qice1(i,j,bi,bj) = 3.4 _d 5
              Qice2(i,j,bi,bj) = 3.4 _d 5
             ENDIF
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDIF
        IF ( thSIceTsurf_InitFile .NE. ' ' ) THEN
         CALL READ_REC_XY_RS(thSIceTsurf_InitFile,Tsrf,1,0,myThid)
        ENDIF
        _END_MASTER(myThid)
      ENDIF

        _EXCH_XY_R8(iceMask,myThid)
        _EXCH_XY_R8(iceHeight, myThid)
        _EXCH_XY_R8(snowHeight,myThid)
        _EXCH_XY_R8(Tsrf,   myThid)
        _EXCH_XY_R8(Tice1,  myThid)
        _EXCH_XY_R8(Tice2,  myThid)
        _EXCH_XY_R8(Qice1,  myThid)
        _EXCH_XY_R8(Qice2,  myThid)
        _EXCH_XY_R8(snowAge,myThid)

C--     Initialise Sea-Ice Loading for SeaIce-Dynamics :
        IF ( useSEAICE ) THEN
          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
               sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
     &                               + iceHeight(i,j,bi,bj)*rhoi
     &                               )*iceMask(i,j,bi,bj)
             ENDDO
            ENDDO
           ENDDO
          ENDDO
        ENDIF

#endif /* ALLOW_THSICE */

      RETURN
      END
1	dimitri	1.13	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_ini_vars.F,v 1.12 2008/05/22 01:59:51 dimitri Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5	jmc	1.8	#ifdef ALLOW_AIM
6			# include "AIM_OPTIONS.h"
7			#endif
8	jmc	1.1
9	jmc	1.2	CBOP
10			C !ROUTINE: THSICE_INI_VARS
11			C !INTERFACE:
12	jmc	1.1	SUBROUTINE THSICE_INI_VARS( myThid )
13	jmc	1.2
14			C !DESCRIPTION: \bv
15	jmc	1.1	C ==========================================================
16	jmc	1.2	C \| S/R THSICE_INI_VARS
17	jmc	1.1	C \| o initialize THermo_SeaICE variables
18			C ==========================================================
19	jmc	1.2	C \ev
20
21			C !USES:
22	jmc	1.1	IMPLICIT NONE
23
24			C === Global variables ===
25			#include "SIZE.h"
26			#include "EEPARAMS.h"
27			#include "PARAMS.h"
28	jmc	1.7	#include "FFIELDS.h"
29	dimitri	1.12	#include "DYNVARS.h"
30	jmc	1.1	#include "THSICE_PARAMS.h"
31	jmc	1.2	#include "THSICE_VARS.h"
32	jmc	1.8	#ifdef ALLOW_AIM
33			# include "AIM_FFIELDS.h"
34			#endif
35	jmc	1.1
36	jmc	1.2	C !INPUT/OUTPUT PARAMETERS:
37	jmc	1.1	C == Routine arguments ==
38	jmc	1.8	C myThid :: My Thread Id. number
39	jmc	1.1	INTEGER myThid
40	jmc	1.2	CEOP
41	jmc	1.1
42			#ifdef ALLOW_THSICE
43			C == Local variables ==
44	jmc	1.2	C bi,bj :: Loop counters
45			C i,j :: Loop counters
46	jmc	1.1	INTEGER bi, bj
47	jmc	1.8	INTEGER i, j
48	jmc	1.1	INTEGER prec
49	jmc	1.5	c CHARACTER*(MAX_LEN_FNAM) fn
50	jmc	1.2	_RL v2Loc
51	jmc	1.1
52
53			c set up ice arrays to zero if starting ice
54	jmc	1.8	DO bj = myByLo(myThid), myByHi(myThid)
55			DO bi = myBxLo(myThid), myBxHi(myThid)
56	jmc	1.2	C- state variables :
57	jmc	1.6	DO j=1-OLy,sNy+OLy
58			DO i=1-OLx,sNx+OLx
59	jmc	1.2	iceMask(i,j,bi,bj) = 0. _d 0
60			iceHeight(i,j,bi,bj)= 0. _d 0
61	jmc	1.1	snowHeight(i,j,bi,bj)=0. _d 0
62	jmc	1.2	Tsrf(i,j,bi,bj) = 0. _d 0
63			Tice1(i,j,bi,bj) = 0. _d 0
64			Tice2(i,j,bi,bj) = 0. _d 0
65			Qice1(i,j,bi,bj) = 0. _d 0
66			Qice2(i,j,bi,bj) = 0. _d 0
67			snowAge(i,j,bi,bj) = 0. _d 0
68			ENDDO
69			ENDDO
70			C- fluxes :
71	jmc	1.6	DO j=1-OLy,sNy+OLy
72			DO i=1-OLx,sNx+OLx
73	jmc	1.2	sHeating(i,j,bi,bj) = 0. _d 0
74			flxCndBt(i,j,bi,bj) = 0. _d 0
75			snowPrc(i,j,bi,bj) = 0. _d 0
76	jmc	1.4	siceAlb(i,j,bi,bj) = 0. _d 0
77	jmc	1.6	icFlxSW (i,j,bi,bj) = 0. _d 0
78			icFlxAtm(i,j,bi,bj) = 0. _d 0
79			icFrwAtm(i,j,bi,bj) = 0. _d 0
80	jmc	1.11	C- needed when using advection/diffusion:
81	heimbach	1.9	oceFWfx(i,j,bi,bj) = 0. _d 0
82			oceSflx(i,j,bi,bj) = 0. _d 0
83			oceQnet(i,j,bi,bj) = 0. _d 0
84	jmc	1.2	ENDDO
85			ENDDO
86			C- oceanic mixed layer state :
87			v2Loc = vMxL_default*vMxL_default
88	jmc	1.6	DO j=1-OLy,sNy+OLy
89			DO i=1-OLx,sNx+OLx
90	jmc	1.2	hOceMxL(i,j,bi,bj) = hMxL_default
91			tOceMxL(i,j,bi,bj) = 0. _d 0
92			sOceMxL(i,j,bi,bj) = sMxL_default
93			v2ocMxL(i,j,bi,bj) = v2Loc
94	jmc	1.1	ENDDO
95			ENDDO
96	jmc	1.8	#ifdef ALLOW_AIM
97			IF ( useAIM ) THEN
98			C- Mask mixed layer depth : depth is used in thsice slab_ocean
99			C and this mask is used in thsice_advdiff and if coupled
100			DO j=1-OLy,sNy+OLy
101			DO i=1-OLx,sNx+OLx
102			IF ( aim_landFr(i,j,bi,bj).EQ.1. _d 0 )
103			& hOceMxL(i,j,bi,bj) = 0.
104			ENDDO
105			ENDDO
106			ENDIF
107			#endif /* ALLOW_AIM */
108	jmc	1.1	ENDDO
109	jmc	1.8	ENDDO
110	jmc	1.1
111	jmc	1.3	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
112
113	jmc	1.2	IF ( startIceModel.LE.0 .AND. nIter0.NE.0 ) THEN
114	jmc	1.1	C-- Read ice pickup fields
115			_BARRIER
116			prec=precFloat64
117	jmc	1.2	CALL THSICE_READ_PICKUP( prec, nIter0, myThid )
118	jmc	1.1
119	jmc	1.3	ELSE
120			C-- Read initial conditions:
121			_BEGIN_MASTER( myThid )
122	dimitri	1.12	IF ( thSIceThick_InitFile .NE. ' ' ) THEN
123			CALL READ_REC_XY_RL(thSIceThick_InitFile,iceHeight,1,0,myThid)
124			ENDIF
125	jmc	1.3	IF ( thSIceFract_InitFile .NE. ' ' ) THEN
126			CALL READ_REC_XY_RL(thSIceFract_InitFile,iceMask,1,0,myThid)
127	dimitri	1.12	ELSE
128	dimitri	1.13	C default to 100% ice concentration when thickness is non-zero
129	dimitri	1.12	DO bj = myByLo(myThid), myByHi(myThid)
130			DO bi = myBxLo(myThid), myBxHi(myThid)
131			DO j=1-OLy,sNy+OLy
132			DO i=1-OLx,sNx+OLx
133			IF (iceHeight(i,j,bi,bj) .NE. 0. _d 0)
134			& iceMask(i,j,bi,bj) = 1. _d 0
135			ENDDO
136			ENDDO
137			ENDDO
138			ENDDO
139	jmc	1.3	ENDIF
140			IF ( thSIceSnowH_InitFile .NE. ' ' ) THEN
141			CALL READ_REC_XY_RL(thSIceSnowH_InitFile,snowHeight,1,0,myThid)
142	dimitri	1.12	ELSE
143	dimitri	1.13	C default to <= 20 cm snow thickness where there is ice
144	dimitri	1.12	DO bj = myByLo(myThid), myByHi(myThid)
145			DO bi = myBxLo(myThid), myBxHi(myThid)
146			DO j=1-OLy,sNy+OLy
147			DO i=1-OLx,sNx+OLx
148	dimitri	1.13	snowHeight(i,j,bi,bj) = 0.2 _d 0
149	dimitri	1.12	ENDDO
150			ENDDO
151			ENDDO
152			ENDDO
153	jmc	1.3	ENDIF
154			IF ( thSIceSnowA_InitFile .NE. ' ' ) THEN
155			CALL READ_REC_XY_RL(thSIceSnowA_InitFile,snowAge,1,0,myThid)
156	dimitri	1.12	ELSE
157	dimitri	1.13	C default to 3 days for snow age
158	dimitri	1.12	DO bj = myByLo(myThid), myByHi(myThid)
159			DO bi = myBxLo(myThid), myBxHi(myThid)
160			DO j=1-OLy,sNy+OLy
161			DO i=1-OLx,sNx+OLx
162			IF (snowHeight(i,j,bi,bj) .NE. 0. _d 0)
163			& snowAge(i,j,bi,bj) = 259200. _d 0
164			ENDDO
165			ENDDO
166			ENDDO
167			ENDDO
168	jmc	1.3	ENDIF
169			IF ( thSIceEnthp_InitFile .NE. ' ' ) THEN
170			CALL READ_REC_XY_RL(thSIceEnthp_InitFile,Qice1,1,0,myThid)
171			CALL READ_REC_XY_RL(thSIceEnthp_InitFile,Qice2,2,0,myThid)
172	dimitri	1.12	ELSE
173	dimitri	1.13	C default to 3.4e5 J/kg for enthalpy where there is sea ice
174	dimitri	1.12	DO bj = myByLo(myThid), myByHi(myThid)
175			DO bi = myBxLo(myThid), myBxHi(myThid)
176			DO j=1-OLy,sNy+OLy
177			DO i=1-OLx,sNx+OLx
178			IF (iceMask(i,j,bi,bj) .NE. 0. _d 0) THEN
179			Qice1(i,j,bi,bj) = 3.4 _d 5
180			Qice2(i,j,bi,bj) = 3.4 _d 5
181			ENDIF
182			ENDDO
183			ENDDO
184			ENDDO
185			ENDDO
186	jmc	1.3	ENDIF
187			IF ( thSIceTsurf_InitFile .NE. ' ' ) THEN
188			CALL READ_REC_XY_RS(thSIceTsurf_InitFile,Tsrf,1,0,myThid)
189			ENDIF
190			_END_MASTER(myThid)
191			ENDIF
192
193	jmc	1.2	_EXCH_XY_R8(iceMask,myThid)
194	jmc	1.1	_EXCH_XY_R8(iceHeight, myThid)
195	jmc	1.2	_EXCH_XY_R8(snowHeight,myThid)
196			_EXCH_XY_R8(Tsrf, myThid)
197			_EXCH_XY_R8(Tice1, myThid)
198			_EXCH_XY_R8(Tice2, myThid)
199			_EXCH_XY_R8(Qice1, myThid)
200			_EXCH_XY_R8(Qice2, myThid)
201			_EXCH_XY_R8(snowAge,myThid)
202	jmc	1.1
203	jmc	1.7	C-- Initialise Sea-Ice Loading for SeaIce-Dynamics :
204			IF ( useSEAICE ) THEN
205			DO bj = myByLo(myThid), myByHi(myThid)
206			DO bi = myBxLo(myThid), myBxHi(myThid)
207			DO j=1-OLy,sNy+OLy
208			DO i=1-OLx,sNx+OLx
209			sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
210			& + iceHeight(i,j,bi,bj)*rhoi
211			& )*iceMask(i,j,bi,bj)
212			ENDDO
213			ENDDO
214			ENDDO
215			ENDDO
216			ENDIF
217
218	jmc	1.1	#endif /* ALLOW_THSICE */
219
220			RETURN
221			END