pkg/aim_v23/aim_aim2sioce.F

C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/aim_aim2sioce.F,v 1.8 2007/10/01 13:34:43 jmc Exp $
C $Name:  $

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

CBOP
C     !ROUTINE: AIM_AIM2SIOCE
C     !INTERFACE:
      SUBROUTINE AIM_AIM2SIOCE(
     I               land_frc, siceFrac,
     O               prcAtm, snowPrc,
     I               bi, bj, myTime, myIter, myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R AIM_AIM2SIOCE
C     | o Interface between AIM and thSIce pkg or (coupled) ocean
C     *==========================================================*
C     | o compute surface fluxes over ocean (ice-free + ice covered)
C     |   for diagnostics, thsice package and (slab, coupled) ocean
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     == Global variables ===
C-- size for MITgcm & Physics package :
#include "AIM_SIZE.h"

#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"

C-- Physics package
#include "AIM_PARAMS.h"
#include "com_physcon.h"
#include "com_physvar.h"

#ifdef ALLOW_THSICE
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"
#endif

C updated fields (in commom blocks):
C  if using thSIce:
C      Qsw(inp)   :: SW radiation through the sea-ice down to the ocean (+=up)
C      Qsw(out)   :: SW radiation down to the ocean (ice-free + ice-covered)(+=up)
C      Qnet(out)  :: Net heat flux out of the ocean (ice-free ocean only)(+=up)
C             and the Ice-Covered contribution will be added in S/R THSICE_STEP_FWD
C      EmPmR(out) :: Net fresh water flux out off the ocean (ice-free ocean only)
C             and the Ice-Covered contribution will be added in S/R THSICE_STEP_FWD
C      sHeating(in/out) :: air - seaice surface heat flux left to melt the ice
C      icFrwAtm   :: Evaporation over sea-ice [kg/m2/s] (>0 if evaporate)
C      icFlxSW    :: net SW heat flux through the ice to the ocean [W/m2] (+=dw)
C  if not using thSIce:
C      Qsw(out)   :: SW radiation down to the ocean (ice-free + ice-covered)(+=up)
C      Qnet(out)  :: Net heat flux out of the ocean (ice-free + ice-covered)(+=up)
C      EmPmR(out) :: Net fresh water flux out off the ocean (ice-free + ice-covered)

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     land_frc :: land fraction [0-1]
C     siceFrac :: sea-ice fraction (relative to full grid-cell) [0-1]
C     prcAtm   :: total precip from the atmosphere [kg/m2/s]
C     snowPrc  :: snow precip over sea-ice [kg/m2/s]
C     bi,bj    :: Tile indices
C     myTime   :: Current time of simulation ( s )
C     myIter   :: Current iteration number in simulation
C     myThid   :: My Thread Id number
      _RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL siceFrac(sNx,sNy)
      _RL prcAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL snowPrc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER bi, bj, myIter, myThid
      _RL myTime
CEOP

#ifdef ALLOW_AIM
C     == Local variables ==
C     i,j,I2      :: loop counters
C     convPrcEvp  :: units conversion factor for Precip & Evap:
C                 :: from AIM units (g/m2/s) to model EmPmR units ( kg/m2/s )
      _RL convPrcEvp
      _RL icFrac, opFrac
      INTEGER i,j,I2

C--   Initialisation :

C--   Atmospheric Physics Fluxes

C     from g/m2/s to kg/m2/s :
      convPrcEvp = 1. _d -3

      DO j=1,sNy
       DO i=1,sNx
        IF ( land_frc(i,j,bi,bj).GE.1. _d 0 ) THEN
C-    Full Land grid-cell: set all fluxes to zero (this has no effect on the
C        model integration and just put this to get meaningfull diagnostics)
         prcAtm(i,j)     = 0. _d 0
         Qnet(i,j,bi,bj) = 0. _d 0
         EmPmR(i,j,bi,bj)= 0. _d 0
         Qsw(i,j,bi,bj)  = 0. _d 0
        ELSE
         I2 = i+(j-1)*sNx

C-    Total Precip (no distinction between ice-covered / ice-free fraction):
         prcAtm(i,j) = ( PRECNV(I2,myThid)
     &                 + PRECLS(I2,myThid) )

C-    Net surface heat flux over ice-free ocean (+=down)
C     note: with aim_splitSIOsFx=F, ice-free & ice covered contribution are
C     already merged together and Qnet is the mean heat flux over the grid box.
         Qnet(i,j,bi,bj) =
     &                         SSR(I2,2,myThid)
     &                       - SLR(I2,2,myThid)
     &                       - SHF(I2,2,myThid)
     &                       - EVAP(I2,2,myThid)*ALHC

C-    E-P over ice-free ocean [m/s]: (same as above is aim_splitSIOsFx=F)
         EmPmR(i,j,bi,bj) = ( EVAP(I2,2,myThid)
     &                      - prcAtm(i,j) ) * convPrcEvp

C-    Net short wave (ice-free ocean) into the ocean (+=down)
         Qsw(i,j,bi,bj) = SSR(I2,2,myThid)

        ENDIF
       ENDDO
      ENDDO

#ifdef ALLOW_THSICE
      IF ( useThSIce ) THEN
       DO j=1,sNy
        DO i=1,sNx
         I2 = i+(j-1)*sNx
C-    Mixed-Layer Ocean: (for thsice slab_ocean and coupler)
C     NOTE: masking is now applied much earlier, during initialisation
c        IF (land_frc(i,j,bi,bj).EQ.1. _d 0) hOceMxL(i,j,bi,bj) = 0.

C-    Evaporation over sea-ice: (for thsice)
         icFrwAtm(i,j,bi,bj) = EVAP(I2,3,myThid)*convPrcEvp

C-    short-wave downward heat flux (ice-free ocean + ice-covered):
C     note: at this point we already called THSICE_IMPL_TEMP to solve for
C     seaice temp and SW flux through the ice. SW is not modified after, and
C     can therefore combine the open-ocean & ice-covered ocean SW fluxes.
         icFrac = iceMask(i,j,bi,bj)
         opFrac = 1. _d 0 - icFrac
         Qsw(i,j,bi,bj) = icFrac*icFlxSW(i,j,bi,bj)
     &                  + opFrac*Qsw(i,j,bi,bj)

        ENDDO
       ENDDO

       IF ( aim_energPrecip ) THEN
C--   Add energy flux related to Precip. (snow, T_rain) over sea-ice
         DO j=1,sNy
          DO i=1,sNx
           IF ( iceMask(i,j,bi,bj).GT.0. _d 0 ) THEN
            I2 = i+(j-1)*sNx
            IF ( EnPrec(I2,myThid).GE.0. _d 0 ) THEN
C-    positive => add to surface heating
              sHeating(i,j,bi,bj) = sHeating(i,j,bi,bj)
     &                            + EnPrec(I2,myThid)*prcAtm(i,j)
              snowPrc(i,j) = 0. _d 0
            ELSE
C-    negative => make snow
              snowPrc(i,j) = prcAtm(i,j)*convPrcEvp
            ENDIF
           ELSE
              snowPrc(i,j) = 0. _d 0
           ENDIF
          ENDDO
         ENDDO
       ENDIF

      ELSEIF ( aim_splitSIOsFx ) THEN
#else /* ALLOW_THSICE */
      IF ( aim_splitSIOsFx ) THEN
#endif /* ALLOW_THSICE */
C-    aim_splitSIOsFx=T: fluxes over sea-ice (3) & ice-free ocean (2) were
C     computed separately and here we merge the 2 fractions
       DO j=1,sNy
        DO i=1,sNx
         I2 = i+(j-1)*sNx
         IF ( siceFrac(i,j) .GT. 0. ) THEN
          icFrac = siceFrac(i,j)/(1. _d 0 - land_frc(i,j,bi,bj))
          opFrac = 1. _d 0 - icFrac

C-    Net surface heat flux over sea-ice + ice-free ocean (+=down)
          Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*opFrac
     &                    + (  SSR(I2,3,myThid)
     &                       - SLR(I2,3,myThid)
     &                       - SHF(I2,3,myThid)
     &                       - EVAP(I2,3,myThid)*ALHC
     &                      )*icFrac
C-    E-P over sea-ice + ice-free ocean [m/s]:
          EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*opFrac
     &                     + ( EVAP(I2,3,myThid)
     &                       - prcAtm(i,j) ) * convPrcEvp * icFrac

C-    Net short wave (ice-free ocean) into the ocean (+=down)
          Qsw(i,j,bi,bj) = opFrac*Qsw(i,j,bi,bj)
     &                   + icFrac*SSR(I2,3,myThid)

         ENDIF
        ENDDO
       ENDDO

C--   end of If useThSIce / elseif aim_splitSIOsFx blocks
      ENDIF

      IF ( aim_energPrecip ) THEN
C--   Ice free fraction: Add energy flux related to Precip. (snow, T_rain):
        DO j=1,sNy
         DO i=1,sNx
          I2 = i+(j-1)*sNx
          Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
     &                    + EnPrec(I2,myThid)*prcAtm(i,j)
         ENDDO
        ENDDO
      ENDIF

      DO j=1,sNy
        DO i=1,sNx
C-    Total Precip : convert units
          prcAtm(i,j) = prcAtm(i,j) * convPrcEvp
C-    Oceanic convention: Heat flux are > 0 upward ; reverse sign.
          Qsw(i,j,bi,bj) = -Qsw(i,j,bi,bj)
          Qnet(i,j,bi,bj)= -Qnet(i,j,bi,bj)
        ENDDO
      ENDDO

#endif /* ALLOW_AIM */

      RETURN
      END
1	jmc	1.9	C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/aim_aim2sioce.F,v 1.8 2007/10/01 13:34:43 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "AIM_OPTIONS.h"
5			#ifdef ALLOW_THSICE
6			#include "THSICE_OPTIONS.h"
7			#endif
8
9			CBOP
10			C !ROUTINE: AIM_AIM2SIOCE
11			C !INTERFACE:
12			SUBROUTINE AIM_AIM2SIOCE(
13	jmc	1.2	I land_frc, siceFrac,
14	jmc	1.9	O prcAtm, snowPrc,
15	jmc	1.1	I bi, bj, myTime, myIter, myThid)
16
17			C !DESCRIPTION: \bv
18			C ==========================================================
19			C \| S/R AIM_AIM2SIOCE
20			C \| o Interface between AIM and thSIce pkg or (coupled) ocean
21			C ==========================================================
22			C \| o compute surface fluxes over ocean (ice-free + ice covered)
23			C \| for diagnostics, thsice package and (slab, coupled) ocean
24			C ==========================================================
25			C \ev
26
27			C !USES:
28			IMPLICIT NONE
29
30			C == Global variables ===
31			C-- size for MITgcm & Physics package :
32	jmc	1.3	#include "AIM_SIZE.h"
33	jmc	1.1
34			#include "EEPARAMS.h"
35			#include "PARAMS.h"
36			#include "FFIELDS.h"
37
38			C-- Physics package
39			#include "AIM_PARAMS.h"
40			#include "com_physcon.h"
41			#include "com_physvar.h"
42
43			#ifdef ALLOW_THSICE
44	jmc	1.3	#include "THSICE_SIZE.h"
45	jmc	1.1	#include "THSICE_PARAMS.h"
46			#include "THSICE_VARS.h"
47			#endif
48
49	jmc	1.3	C updated fields (in commom blocks):
50			C if using thSIce:
51			C Qsw(inp) :: SW radiation through the sea-ice down to the ocean (+=up)
52			C Qsw(out) :: SW radiation down to the ocean (ice-free + ice-covered)(+=up)
53			C Qnet(out) :: Net heat flux out of the ocean (ice-free ocean only)(+=up)
54			C and the Ice-Covered contribution will be added in S/R THSICE_STEP_FWD
55			C EmPmR(out) :: Net fresh water flux out off the ocean (ice-free ocean only)
56			C and the Ice-Covered contribution will be added in S/R THSICE_STEP_FWD
57			C sHeating(in/out) :: air - seaice surface heat flux left to melt the ice
58	jmc	1.5	C icFrwAtm :: Evaporation over sea-ice [kg/m2/s] (>0 if evaporate)
59			C icFlxSW :: net SW heat flux through the ice to the ocean [W/m2] (+=dw)
60	jmc	1.3	C if not using thSIce:
61			C Qsw(out) :: SW radiation down to the ocean (ice-free + ice-covered)(+=up)
62			C Qnet(out) :: Net heat flux out of the ocean (ice-free + ice-covered)(+=up)
63			C EmPmR(out) :: Net fresh water flux out off the ocean (ice-free + ice-covered)
64
65	jmc	1.1	C !INPUT/OUTPUT PARAMETERS:
66			C == Routine arguments ==
67			C land_frc :: land fraction [0-1]
68	jmc	1.2	C siceFrac :: sea-ice fraction (relative to full grid-cell) [0-1]
69	jmc	1.1	C prcAtm :: total precip from the atmosphere [kg/m2/s]
70	jmc	1.9	C snowPrc :: snow precip over sea-ice [kg/m2/s]
71	jmc	1.7	C bi,bj :: Tile indices
72	jmc	1.1	C myTime :: Current time of simulation ( s )
73			C myIter :: Current iteration number in simulation
74	jmc	1.3	C myThid :: My Thread Id number
75	jmc	1.1	_RS land_frc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76	jmc	1.2	_RL siceFrac(sNx,sNy)
77	jmc	1.9	_RL prcAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78			_RL snowPrc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79	jmc	1.1	INTEGER bi, bj, myIter, myThid
80			_RL myTime
81			CEOP
82
83			#ifdef ALLOW_AIM
84			C == Local variables ==
85	jmc	1.8	C i,j,I2 :: loop counters
86	jmc	1.9	C convPrcEvp :: units conversion factor for Precip & Evap:
87	jmc	1.8	C :: from AIM units (g/m2/s) to model EmPmR units ( kg/m2/s )
88			_RL convPrcEvp
89	jmc	1.1	_RL icFrac, opFrac
90			INTEGER i,j,I2
91
92			C-- Initialisation :
93
94			C-- Atmospheric Physics Fluxes
95
96	jmc	1.2	C from g/m2/s to kg/m2/s :
97	jmc	1.8	convPrcEvp = 1. _d -3
98	jmc	1.1
99			DO j=1,sNy
100	jmc	1.4	DO i=1,sNx
101			IF ( land_frc(i,j,bi,bj).GE.1. _d 0 ) THEN
102			C- Full Land grid-cell: set all fluxes to zero (this has no effect on the
103			C model integration and just put this to get meaningfull diagnostics)
104			prcAtm(i,j) = 0. _d 0
105			Qnet(i,j,bi,bj) = 0. _d 0
106			EmPmR(i,j,bi,bj)= 0. _d 0
107			Qsw(i,j,bi,bj) = 0. _d 0
108			ELSE
109	jmc	1.1	I2 = i+(j-1)*sNx
110
111	jmc	1.3	C- Total Precip (no distinction between ice-covered / ice-free fraction):
112	jmc	1.1	prcAtm(i,j) = ( PRECNV(I2,myThid)
113			& + PRECLS(I2,myThid) )
114
115			C- Net surface heat flux over ice-free ocean (+=down)
116	jmc	1.3	C note: with aim_splitSIOsFx=F, ice-free & ice covered contribution are
117			C already merged together and Qnet is the mean heat flux over the grid box.
118			Qnet(i,j,bi,bj) =
119	jmc	1.1	& SSR(I2,2,myThid)
120			& - SLR(I2,2,myThid)
121			& - SHF(I2,2,myThid)
122			& - EVAP(I2,2,myThid)*ALHC
123
124	jmc	1.3	C- E-P over ice-free ocean [m/s]: (same as above is aim_splitSIOsFx=F)
125	jmc	1.1	EmPmR(i,j,bi,bj) = ( EVAP(I2,2,myThid)
126	jmc	1.8	& - prcAtm(i,j) ) * convPrcEvp
127	jmc	1.1
128			C- Net short wave (ice-free ocean) into the ocean (+=down)
129			Qsw(i,j,bi,bj) = SSR(I2,2,myThid)
130
131	jmc	1.4	ENDIF
132			ENDDO
133	jmc	1.1	ENDDO
134
135			#ifdef ALLOW_THSICE
136			IF ( useThSIce ) THEN
137			DO j=1,sNy
138			DO i=1,sNx
139			I2 = i+(j-1)*sNx
140	jmc	1.3	C- Mixed-Layer Ocean: (for thsice slab_ocean and coupler)
141	jmc	1.7	C NOTE: masking is now applied much earlier, during initialisation
142			c IF (land_frc(i,j,bi,bj).EQ.1. _d 0) hOceMxL(i,j,bi,bj) = 0.
143	jmc	1.1
144	jmc	1.3	C- Evaporation over sea-ice: (for thsice)
145	jmc	1.8	icFrwAtm(i,j,bi,bj) = EVAP(I2,3,myThid)*convPrcEvp
146	jmc	1.1
147	jmc	1.3	C- short-wave downward heat flux (ice-free ocean + ice-covered):
148	jmc	1.5	C note: at this point we already called THSICE_IMPL_TEMP to solve for
149	jmc	1.3	C seaice temp and SW flux through the ice. SW is not modified after, and
150			C can therefore combine the open-ocean & ice-covered ocean SW fluxes.
151	jmc	1.1	icFrac = iceMask(i,j,bi,bj)
152			opFrac = 1. _d 0 - icFrac
153	jmc	1.5	Qsw(i,j,bi,bj) = icFrac*icFlxSW(i,j,bi,bj)
154			& + opFrac*Qsw(i,j,bi,bj)
155	jmc	1.1
156			ENDDO
157			ENDDO
158
159	jmc	1.3	IF ( aim_energPrecip ) THEN
160	jmc	1.1	C-- Add energy flux related to Precip. (snow, T_rain) over sea-ice
161	jmc	1.3	DO j=1,sNy
162			DO i=1,sNx
163			IF ( iceMask(i,j,bi,bj).GT.0. _d 0 ) THEN
164			I2 = i+(j-1)*sNx
165			IF ( EnPrec(I2,myThid).GE.0. _d 0 ) THEN
166	jmc	1.1	C- positive => add to surface heating
167	jmc	1.3	sHeating(i,j,bi,bj) = sHeating(i,j,bi,bj)
168			& + EnPrec(I2,myThid)*prcAtm(i,j)
169	jmc	1.9	snowPrc(i,j) = 0. _d 0
170	jmc	1.3	ELSE
171	jmc	1.1	C- negative => make snow
172	jmc	1.9	snowPrc(i,j) = prcAtm(i,j)*convPrcEvp
173	jmc	1.3	ENDIF
174			ELSE
175	jmc	1.9	snowPrc(i,j) = 0. _d 0
176	jmc	1.3	ENDIF
177			ENDDO
178			ENDDO
179			ENDIF
180
181	jmc	1.2	ELSEIF ( aim_splitSIOsFx ) THEN
182			#else /* ALLOW_THSICE */
183			IF ( aim_splitSIOsFx ) THEN
184			#endif /* ALLOW_THSICE */
185	jmc	1.3	C- aim_splitSIOsFx=T: fluxes over sea-ice (3) & ice-free ocean (2) were
186			C computed separately and here we merge the 2 fractions
187	jmc	1.2	DO j=1,sNy
188			DO i=1,sNx
189			I2 = i+(j-1)*sNx
190			IF ( siceFrac(i,j) .GT. 0. ) THEN
191			icFrac = siceFrac(i,j)/(1. _d 0 - land_frc(i,j,bi,bj))
192			opFrac = 1. _d 0 - icFrac
193
194			C- Net surface heat flux over sea-ice + ice-free ocean (+=down)
195			Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*opFrac
196			& + ( SSR(I2,3,myThid)
197			& - SLR(I2,3,myThid)
198			& - SHF(I2,3,myThid)
199			& - EVAP(I2,3,myThid)*ALHC
200			& )*icFrac
201			C- E-P over sea-ice + ice-free ocean [m/s]:
202			EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*opFrac
203			& + ( EVAP(I2,3,myThid)
204	jmc	1.8	& - prcAtm(i,j) ) * convPrcEvp * icFrac
205	jmc	1.2
206			C- Net short wave (ice-free ocean) into the ocean (+=down)
207			Qsw(i,j,bi,bj) = opFrac*Qsw(i,j,bi,bj)
208			& + icFrac*SSR(I2,3,myThid)
209
210	jmc	1.1	ENDIF
211			ENDDO
212			ENDDO
213	jmc	1.3
214			C-- end of If useThSIce / elseif aim_splitSIOsFx blocks
215	jmc	1.2	ENDIF
216
217			IF ( aim_energPrecip ) THEN
218	jmc	1.3	C-- Ice free fraction: Add energy flux related to Precip. (snow, T_rain):
219	jmc	1.2	DO j=1,sNy
220			DO i=1,sNx
221			I2 = i+(j-1)*sNx
222			Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
223			& + EnPrec(I2,myThid)*prcAtm(i,j)
224			ENDDO
225			ENDDO
226	jmc	1.1	ENDIF
227
228			DO j=1,sNy
229			DO i=1,sNx
230			C- Total Precip : convert units
231	jmc	1.8	prcAtm(i,j) = prcAtm(i,j) * convPrcEvp
232	jmc	1.1	C- Oceanic convention: Heat flux are > 0 upward ; reverse sign.
233			Qsw(i,j,bi,bj) = -Qsw(i,j,bi,bj)
234			Qnet(i,j,bi,bj)= -Qnet(i,j,bi,bj)
235			ENDDO
236			ENDDO
237
238			#endif /* ALLOW_AIM */
239
240			RETURN
241			END