itd/code/seaice_budget_ocean.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_budget_ocean.F,v 1.19 2012/11/09 22:15:18 heimbach Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_BUDGET_OCEAN(
     I     UG,
     I     TSURF,
     O     netHeatFlux, SWHeatFlux,
     I     bi, bj, myTime, myIter, myThid )
C     *================================================================*
C     | SUBROUTINE seaice_budget_ocean
C     | o Calculate surface heat fluxes over open ocean
C     |   see Hibler, MWR, 108, 1943-1973, 1980
C     |   If SEAICE_EXTERNAL_FLUXES is defined this routine simply
C     |   copies the global fields to the seaice-local fields.
C     *================================================================*
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "FFIELDS.h"
#ifndef SEAICE_EXTERNAL_FLUXES
# include "PARAMS.h"
# include "GRID.h"
# include "SEAICE_SIZE.h"
# include "SEAICE_PARAMS.h"
# ifdef ALLOW_EXF
#  include "EXF_OPTIONS.h"
#  include "EXF_FIELDS.h"
# endif
#endif

C     === Routine arguments ===
C     INPUT:
C     UG      :: thermal wind of atmosphere
C     TSURF   :: ocean surface temperature in Kelvin
C     bi,bj   :: loop indices
C     myTime  :: Simulation time
C     myIter  :: Simulation timestep number
C     myThid  :: Thread no. that called this routine.
C     OUTPUT:
C     netHeatFlux :: net surface heat flux over open water or under ice
C     SWHeatFlux  :: short wave heat flux over open water or under ice
      _RL UG         (1:sNx,1:sNy)
      _RL TSURF      (1:sNx,1:sNy)
      _RL netHeatFlux(1:sNx,1:sNy)
      _RL SWHeatFlux (1:sNx,1:sNy)
      _RL myTime
      INTEGER bi, bj, myIter, myThid
CEndOfInterface

C     === Local variables ===
C     i,j - Loop counters
      INTEGER i, j
#ifndef SEAICE_EXTERNAL_FLUXES
      _RL  QS1, D1, D1W, D3, TMELT

C     local copies of global variables
      _RL tsurfLoc   (1:sNx,1:sNy)
      _RL atempLoc   (1:sNx,1:sNy)
      _RL lwdownLoc  (1:sNx,1:sNy)

C     auxiliary variable
      _RL ssq, sstdegC
      _RL recip_rhoConstFresh, recip_lhEvap

C NOW DEFINE ASSORTED CONSTANTS
C SATURATION VAPOR PRESSURE CONSTANT
      QS1=0.622 _d 0/1013.0 _d 0
C SENSIBLE HEAT CONSTANT
      D1=SEAICE_dalton*SEAICE_cpAir*SEAICE_rhoAir
C WATER LATENT HEAT CONSTANT
      D1W=SEAICE_dalton*SEAICE_lhEvap*SEAICE_rhoAir
C STEFAN BOLTZMAN CONSTANT TIMES EMISSIVITY
      D3=SEAICE_emissivity*SEAICE_boltzmann
C MELTING TEMPERATURE OF ICE
      TMELT = celsius2K
C
      recip_lhEvap = 1./SEAICE_lhEvap
      recip_rhoConstFresh = 1./rhoConstFresh

      DO J=1,sNy
       DO I=1,sNx
        netHeatFlux(I,J) = 0. _d 0
        SWHeatFlux (I,J) = 0. _d 0
C
C     MAX_TICE does not exist anly longer, lets see if it works without
C       tsurfLoc (I,J) = MIN(celsius2K+MAX_TICE,TSURF(I,J))
        tsurfLoc (I,J) = TSURF(I,J)
# ifdef ALLOW_ATM_TEMP
C     Is this necessary?
        atempLoc (I,J) = MAX(celsius2K+MIN_ATEMP,ATEMP(I,J,bi,bj))
# endif
# ifdef ALLOW_DOWNWARD_RADIATION
        lwdownLoc(I,J) = MAX(MIN_LWDOWN,LWDOWN(I,J,bi,bj))
# endif
       ENDDO
      ENDDO
#endif /* SEAICE_EXTERNAL_FLUXES */

C NOW DETERMINE OPEN WATER HEAT BUD. ASSUMING TSURF=WATER TEMP.
C WATER ALBEDO IS ASSUMED TO BE THE CONSTANT SEAICE_waterAlbedo
      DO J=1,sNy
       DO I=1,sNx
#ifdef SEAICE_EXTERNAL_FLUXES
        netHeatFlux(I,J) = Qnet(I,J,bi,bj)
        SWHeatFlux (I,J) =  Qsw(I,J,bi,bj)
#else /* SEAICE_EXTERNAL_FLUXES undefined */
C     This is an example of how one could implement surface fluxes
C     over the ocean (if one dislikes the fluxes computed in pkg/exf).
C     In this example, the exf-fields are re-used so that they no longer
C     have the same values as at the time when they are saved for
C     diagnostics (e.g., hl, hs, lwflux, sflux). To properly
C     diagnose them, one has to save them again as different (SI-)fields.
# ifdef ALLOW_DOWNWARD_RADIATION
C     net upward short wave heat flux
        SWHeatFlux(I,J) = (SEAICE_waterAlbedo - 1. _d 0)
     &       *swdown(I,J,bi,bj)
C     lwup = emissivity*stefanBoltzmann*Tsrf^4 + (1-emissivity)*lwdown
C     the second terms is the reflected incoming long wave radiation
C     so that the net upward long wave heat flux is:
        lwflux(I,J,bi,bj) = - lwdownLoc(I,J)*SEAICE_emissivity
     &       + D3*tsurfLoc(I,J)**4
        sstdegC = tsurfLoc(I,J) - TMELT
C     downward sensible heat
        hs(I,J,bi,bj) = D1*UG(I,J)*(atempLoc(I,J)-tsurfLoc(I,J))
C     saturation humidity
        ssq = QS1*6.11 _d 0 *EXP( 17.2694 _d 0
     &                           *sstdegC/(sstdegC+237.3 _d 0) )
C     downward latent heat
        hl(I,J,bi,bj) = D1W*UG(I,J)*(AQH(I,J,bi,bj)-ssq)
C     net heat is positive upward
        netHeatFlux(I,J)=SWHeatFlux(I,J)
     &       + lwflux(I,J,bi,bj)
     &       - hs(I,J,bi,bj) - hl(I,J,bi,bj)
C     compute evaporation here again because latent heat is different
C     from its previous value
        evap(i,j,bi,bj) = -hl(I,J,bi,bj)
     &       *recip_lhEvap*recip_rhoConstFresh
C     Salt flux from Precipitation and Evaporation.
        sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)
#  ifdef ALLOW_RUNOFF
        sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
#  endif
        sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskC(i,j,1,bi,bj)
        empmr(i,j,bi,bj) = sflux(i,j,bi,bj)*rhoConstFresh
# endif /* ALLOW_DOWNWARD_RADIATION */
#endif /* SEAICE_EXTERNAL_FLUXES */
       ENDDO
      ENDDO

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_budget_ocean.F,v 1.19 2012/11/09 22:15:18 heimbach Exp $
2	C $Name: $
3
4	#include "SEAICE_OPTIONS.h"
5
6	CStartOfInterface
7	SUBROUTINE SEAICE_BUDGET_OCEAN(
8	I UG,
9	I TSURF,
10	O netHeatFlux, SWHeatFlux,
11	I bi, bj, myTime, myIter, myThid )
12	C ================================================================
13	C \| SUBROUTINE seaice_budget_ocean
14	C \| o Calculate surface heat fluxes over open ocean
15	C \| see Hibler, MWR, 108, 1943-1973, 1980
16	C \| If SEAICE_EXTERNAL_FLUXES is defined this routine simply
17	C \| copies the global fields to the seaice-local fields.
18	C ================================================================
19	IMPLICIT NONE
20
21	C === Global variables ===
22	#include "SIZE.h"
23	#include "EEPARAMS.h"
24	#include "FFIELDS.h"
25	#ifndef SEAICE_EXTERNAL_FLUXES
26	# include "PARAMS.h"
27	# include "GRID.h"
28	# include "SEAICE_SIZE.h"
29	# include "SEAICE_PARAMS.h"
30	# ifdef ALLOW_EXF
31	# include "EXF_OPTIONS.h"
32	# include "EXF_FIELDS.h"
33	# endif
34	#endif
35
36	C === Routine arguments ===
37	C INPUT:
38	C UG :: thermal wind of atmosphere
39	C TSURF :: ocean surface temperature in Kelvin
40	C bi,bj :: loop indices
41	C myTime :: Simulation time
42	C myIter :: Simulation timestep number
43	C myThid :: Thread no. that called this routine.
44	C OUTPUT:
45	C netHeatFlux :: net surface heat flux over open water or under ice
46	C SWHeatFlux :: short wave heat flux over open water or under ice
47	_RL UG (1:sNx,1:sNy)
48	_RL TSURF (1:sNx,1:sNy)
49	_RL netHeatFlux(1:sNx,1:sNy)
50	_RL SWHeatFlux (1:sNx,1:sNy)
51	_RL myTime
52	INTEGER bi, bj, myIter, myThid
53	CEndOfInterface
54
55	C === Local variables ===
56	C i,j - Loop counters
57	INTEGER i, j
58	#ifndef SEAICE_EXTERNAL_FLUXES
59	_RL QS1, D1, D1W, D3, TMELT
60
61	C local copies of global variables
62	_RL tsurfLoc (1:sNx,1:sNy)
63	_RL atempLoc (1:sNx,1:sNy)
64	_RL lwdownLoc (1:sNx,1:sNy)
65
66	C auxiliary variable
67	_RL ssq, sstdegC
68	_RL recip_rhoConstFresh, recip_lhEvap
69
70	C NOW DEFINE ASSORTED CONSTANTS
71	C SATURATION VAPOR PRESSURE CONSTANT
72	QS1=0.622 _d 0/1013.0 _d 0
73	C SENSIBLE HEAT CONSTANT
74	D1=SEAICE_daltonSEAICE_cpAirSEAICE_rhoAir
75	C WATER LATENT HEAT CONSTANT
76	D1W=SEAICE_daltonSEAICE_lhEvapSEAICE_rhoAir
77	C STEFAN BOLTZMAN CONSTANT TIMES EMISSIVITY
78	D3=SEAICE_emissivity*SEAICE_boltzmann
79	C MELTING TEMPERATURE OF ICE
80	TMELT = celsius2K
81	C
82	recip_lhEvap = 1./SEAICE_lhEvap
83	recip_rhoConstFresh = 1./rhoConstFresh
84
85	DO J=1,sNy
86	DO I=1,sNx
87	netHeatFlux(I,J) = 0. _d 0
88	SWHeatFlux (I,J) = 0. _d 0
89	C
90	C MAX_TICE does not exist anly longer, lets see if it works without
91	C tsurfLoc (I,J) = MIN(celsius2K+MAX_TICE,TSURF(I,J))
92	tsurfLoc (I,J) = TSURF(I,J)
93	# ifdef ALLOW_ATM_TEMP
94	C Is this necessary?
95	atempLoc (I,J) = MAX(celsius2K+MIN_ATEMP,ATEMP(I,J,bi,bj))
96	# endif
97	# ifdef ALLOW_DOWNWARD_RADIATION
98	lwdownLoc(I,J) = MAX(MIN_LWDOWN,LWDOWN(I,J,bi,bj))
99	# endif
100	ENDDO
101	ENDDO
102	#endif /* SEAICE_EXTERNAL_FLUXES */
103
104	C NOW DETERMINE OPEN WATER HEAT BUD. ASSUMING TSURF=WATER TEMP.
105	C WATER ALBEDO IS ASSUMED TO BE THE CONSTANT SEAICE_waterAlbedo
106	DO J=1,sNy
107	DO I=1,sNx
108	#ifdef SEAICE_EXTERNAL_FLUXES
109	netHeatFlux(I,J) = Qnet(I,J,bi,bj)
110	SWHeatFlux (I,J) = Qsw(I,J,bi,bj)
111	#else /* SEAICE_EXTERNAL_FLUXES undefined */
112	C This is an example of how one could implement surface fluxes
113	C over the ocean (if one dislikes the fluxes computed in pkg/exf).
114	C In this example, the exf-fields are re-used so that they no longer
115	C have the same values as at the time when they are saved for
116	C diagnostics (e.g., hl, hs, lwflux, sflux). To properly
117	C diagnose them, one has to save them again as different (SI-)fields.
118	# ifdef ALLOW_DOWNWARD_RADIATION
119	C net upward short wave heat flux
120	SWHeatFlux(I,J) = (SEAICE_waterAlbedo - 1. _d 0)
121	& *swdown(I,J,bi,bj)
122	C lwup = emissivitystefanBoltzmannTsrf^4 + (1-emissivity)*lwdown
123	C the second terms is the reflected incoming long wave radiation
124	C so that the net upward long wave heat flux is:
125	lwflux(I,J,bi,bj) = - lwdownLoc(I,J)*SEAICE_emissivity
126	& + D3tsurfLoc(I,J)*4
127	sstdegC = tsurfLoc(I,J) - TMELT
128	C downward sensible heat
129	hs(I,J,bi,bj) = D1UG(I,J)(atempLoc(I,J)-tsurfLoc(I,J))
130	C saturation humidity
131	ssq = QS16.11 _d 0 EXP( 17.2694 _d 0
132	& *sstdegC/(sstdegC+237.3 _d 0) )
133	C downward latent heat
134	hl(I,J,bi,bj) = D1WUG(I,J)(AQH(I,J,bi,bj)-ssq)
135	C net heat is positive upward
136	netHeatFlux(I,J)=SWHeatFlux(I,J)
137	& + lwflux(I,J,bi,bj)
138	& - hs(I,J,bi,bj) - hl(I,J,bi,bj)
139	C compute evaporation here again because latent heat is different
140	C from its previous value
141	evap(i,j,bi,bj) = -hl(I,J,bi,bj)
142	& recip_lhEvaprecip_rhoConstFresh
143	C Salt flux from Precipitation and Evaporation.
144	sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)
145	# ifdef ALLOW_RUNOFF
146	sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
147	# endif
148	sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskC(i,j,1,bi,bj)
149	empmr(i,j,bi,bj) = sflux(i,j,bi,bj)*rhoConstFresh
150	# endif /* ALLOW_DOWNWARD_RADIATION */
151	#endif /* SEAICE_EXTERNAL_FLUXES */
152	ENDDO
153	ENDDO
154
155	RETURN
156	END