/[MITgcm]/MITgcm/pkg/thsice/thsice_step_fwd.F
ViewVC logotype

Annotation of /MITgcm/pkg/thsice/thsice_step_fwd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.2 - (hide annotations) (download)
Wed Dec 31 17:44:32 2003 UTC (20 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint52l_pre, checkpoint52e_pre, hrcube4, checkpoint52j_post, checkpoint52e_post, hrcube_1, checkpoint52l_post, checkpoint52k_post, checkpoint52f_post, hrcube5, checkpoint52i_post, checkpoint52j_pre, checkpoint52i_pre, checkpoint52h_pre, checkpoint52f_pre, hrcube_2, hrcube_3
Changes since 1.1: +26 -20 lines 
change how Albedo is computed over sea-ice:
 - separate albedo S/R (thsice_albedo.F) from thsice_therm.F file.
 - use GISS formulation for snow-covered area ; turn on snow-aging ;
 - bare-ice albedo as in LANL code ;
 - addd diagnostic of surface albedo.
1 jmc 1.2 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.1 2003/11/23 01:20:13 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "THSICE_OPTIONS.h"
5    
6     C !ROUTINE: THSICE_STEP_FWD
7     C !INTERFACE:
8     SUBROUTINE THSICE_STEP_FWD(
9     I bi, bj, iMin, iMax, jMin, jMax,
10     I myTime, myIter, myThid )
11     C *==========================================================*
12     C | SUBROUTINE THSICE_STEP_FWD
13     C | o Step Forward Therm-SeaIce model.
14     C *==========================================================*
15    
16     C !USES:
17     IMPLICIT NONE
18     C === Global variables ===
19     #include "SIZE.h"
20     #include "EEPARAMS.h"
21     #include "PARAMS.h"
22     #include "FFIELDS.h"
23     #include "DYNVARS.h"
24     #include "GRID.h"
25     #ifdef ALLOW_BULK_FORCE
26     #include "BULKF.h"
27     #endif
28     #include "THSICE_SIZE.h"
29     #include "THSICE_PARAMS.h"
30     #include "THSICE.h"
31     #include "THSICE_DIAGS.h"
32    
33     C !INPUT/OUTPUT PARAMETERS:
34     C === Routine arguments ===
35     C myIter :: iteration counter for this thread
36     C myTime :: time counter for this thread
37     C myThid :: thread number for this instance of the routine.
38     INTEGER bi,bj
39     INTEGER iMin, iMax
40     INTEGER jMin, jMax
41     _RL myTime
42     INTEGER myIter
43     INTEGER myThid
44    
45     #ifdef ALLOW_THSICE
46     C !LOCAL VARIABLES:
47     C === Local variables ===
48 jmc 1.2 C snowAge :: snow age (s)
49     C albedo :: surface albedo [0-1]
50     C fSWabs :: net Short-Wave (+ = down) at surface (W m-2)
51 jmc 1.1 C Fbot :: the oceanic heat flux already incorporated (ice_therm)
52     C flxAtm :: net heat flux from the atmosphere ( >0 downward)
53     C evpAtm :: evaporation to the atmosphere
54     C frwAtm :: net fresh-water flux (E-P-R) to the atmosphere (m/s)
55     C qleft :: net heat flux from the ice to the ocean
56     C ffresh :: fresh-water flux from the ice to the ocean
57     C fsalt :: mass salt flux to the ocean
58     INTEGER i,j
59 jmc 1.2 _RL snowAge
60     _RL albedo
61     _RL fSWabs
62     _RL qleft, qNewIce
63 jmc 1.1 _RL fsalt
64     _RL ffresh
65     _RL Tf, cphm, frzmlt
66     _RL Fbot, esurp
67     _RL flxAtm, evpAtm, frwAtm
68     _RL openFrac, iceFrac, qicAv
69     _RL oceHs, oceV2s, oceSs, oceTs
70     _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr)
71    
72     LOGICAL dBug
73    
74     dBug = .FALSE.
75     1010 FORMAT(A,1P4E11.3)
76    
77     DO j = jMin, jMax
78     DO i = iMin, iMax
79     c dBug = ( bi.EQ.3 .AND. i.EQ.13 .AND. j.EQ.13 )
80    
81     Tf = -mu_Tf*salt(i,j,1,bi,bj)
82     cphm = cpwater*rhosw*drF(1)*hFacC(i,j,1,bi,bj)
83 jmc 1.2 oceTs = theta(i,j,1,bi,bj)
84     frzmlt = (Tf-oceTs)*cphm/thSIce_deltaT
85     compact= iceMask(i,j,bi,bj)
86     hIce = iceHeight(i,j,bi,bj)
87     hSnow = snowHeight(i,j,bi,bj)
88 jmc 1.1 Fbot = 0. _d 0
89     snow(i,j,bi,bj) = 0. _d 0
90     saltFlux(i,j,bi,bj) = 0. _d 0
91    
92 jmc 1.2 IF (dBug .AND. (frzmlt.GT.0. .OR. compact.GT.0.) ) THEN
93 jmc 1.1 WRITE(6,1010) 'ThSI_FWD:-0- iceMask,hIc,hSn,Qnet=',
94 jmc 1.2 & compact, hIce, hSnow, Qnet(i,j,bi,bj)
95 jmc 1.1 WRITE(6,1010) 'ThSI_FWD: ocTs,Tf,frzmlt=',
96 jmc 1.2 & oceTs,Tf,frzmlt
97 jmc 1.1 ENDIF
98    
99     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
100     C part.1 : ice-covered fraction ;
101     C can only reduce the ice-fraction but not increase it.
102     C-------
103     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
104     oceHs = hfacC(i,j,1,bi,bj)*drF(1)
105     oceSs = salt (i,j,1,bi,bj)
106     oceV2s = ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
107     & + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
108     & + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
109     & + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj) )*0.5 _d 0
110 jmc 1.2 snowAge = sage(i,j,bi,bj)
111     c snowAge = thSIce_deltaT
112 jmc 1.1 Tsf = Tsrf(i,j,bi,bj)
113     Tice(1) = Tice1(i,j,bi,bj)
114     Tice(2) = Tice2(i,j,bi,bj)
115     qicen(1)= Qice1(i,j,bi,bj)
116     qicen(2)= Qice2(i,j,bi,bj)
117 jmc 1.2 CALL THSICE_ALBEDO(hIce,hSnow,Tsf,snowAge,albedo)
118     fSWabs = solar(i,j,bi,bj)*(1. _d 0 - albedo)
119 jmc 1.1 CALL THSICE_THERM(
120 jmc 1.2 I fSWabs, oceHs, oceV2s, oceSs, oceTs,
121 jmc 1.1 U compact, hIce, hSnow, Tsf, Tice, qicen,
122     O qleft, ffresh, fsalt, Fbot,
123     O flxAtm, evpAtm,
124     I i,j, bi,bj, myThid)
125    
126     C-- Diagnostic of Atmospheric Fluxes over sea-ice :
127     frwAtm = evpAtm - snow(i,j,bi,bj)*rhos/rhofw
128     C note: Any flux of mass (here fresh water) that enter or leave the system
129     C with a non zero energy HAS TO be counted: add snow precip.
130     flxAtm = flxAtm - Lfresh*snow(i,j,bi,bj)*rhos
131    
132     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
133     IF (dBug) WRITE(6,1010) 'ThSI_FWD: iceFrac,flxAtm,evpAtm,flxSnw=',
134     & iceMask(i,j,bi,bj),flxAtm,evpAtm,-Lfresh*snow(i,j,bi,bj)*rhos
135     IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qleft,fsalt,ffresh=',
136     & compact,qleft,fsalt,ffresh
137     #ifdef CHECK_ENERGY_CONSERV
138     iceFrac = iceMask(i,j,bi,bj)
139     CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,
140     I iceFrac, compact, hIce, hSnow, qicen,
141     I qleft, ffresh, fsalt, flxAtm, frwAtm,
142     I myTime, myIter, myThid )
143     #endif /* CHECK_ENERGY_CONSERV */
144     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
145    
146     C-- Update Sea-Ice state :
147     c iceMask(i,j,bi,bj)=compact
148     iceheight(i,j,bi,bj) = hIce
149     snowheight(i,j,bi,bj)= hSnow
150     Tsrf(i,j,bi,bj) =Tsf
151     Tice1(i,j,bi,bj)=Tice(1)
152     Tice2(i,j,bi,bj)=Tice(2)
153     Qice1(i,j,bi,bj)=qicen(1)
154     Qice2(i,j,bi,bj)=qicen(2)
155    
156     C-- Net fluxes :
157     ffresh = ffresh/rhofw
158     ffresh = -ffresh-rain(i,j,bi,bj)-runoff(i,j,bi,bj)
159     frwAtm = frwAtm-rain(i,j,bi,bj)-runoff(i,j,bi,bj)
160     iceFrac = iceMask(i,j,bi,bj)
161     openFrac= 1. _d 0-iceFrac
162     #ifdef ALLOW_TIMEAVE
163     ICE_Qnet_AVE(i,j,bi,bj) = ICE_Qnet_AVE(i,j,bi,bj)
164     & + ( -iceFrac*flxAtm + openFrac*Qnet(i,j,bi,bj)
165     & )*thSIce_deltaT
166     ICE_FWfx_AVE(i,j,bi,bj) = ICE_FWfx_AVE(i,j,bi,bj)
167     & + ( iceFrac*frwAtm + openFrac*EmPmR(i,j,bi,bj)
168     & )*thSIce_deltaT
169 jmc 1.2 ICE_albedo_AVE(i,j,bi,bj) = ICE_albedo_AVE(i,j,bi,bj)
170     & + iceFrac*albedo*thSIce_deltaT
171 jmc 1.1 #endif /*ALLOW_TIMEAVE*/
172     Qnet(i,j,bi,bj)=-iceFrac*qleft + openFrac*Qnet(i,j,bi,bj)
173     EmPmR(i,j,bi,bj)=iceFrac*ffresh+openFrac*EmPmR(i,j,bi,bj)
174     saltFlux(i,j,bi,bj)=-iceFrac*fsalt
175    
176     IF (dBug) WRITE(6,1010)'ThSI_FWD:-1- compact,hIc,hSn,Qnet=',
177     & compact,hIce,hSnow,Qnet(i,j,bi,bj)
178    
179     ELSEIF (hFacC(i,j,1,bi,bj).gt.0. _d 0) THEN
180    
181     #ifdef ALLOW_TIMEAVE
182     ICE_Qnet_AVE(i,j,bi,bj) = ICE_Qnet_AVE(i,j,bi,bj)
183     & +Qnet(i,j,bi,bj)*thSIce_deltaT
184     ICE_FWfx_AVE(i,j,bi,bj) = ICE_FWfx_AVE(i,j,bi,bj)
185     & +EmPmR(i,j,bi,bj)*thSIce_deltaT
186     #endif /*ALLOW_TIMEAVE*/
187    
188     ENDIF
189    
190     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
191     C part.2 : freezing of sea-water
192     C over ice-free fraction and what is left from ice-covered fraction
193     C-------
194     esurp = frzmlt - Fbot*iceMask(i,j,bi,bj)
195     IF (esurp.GT.0. _d 0) THEN
196     iceFrac = compact
197     IF ( compact.GT.0. _d 0 ) THEN
198     qicen(1)= Qice1(i,j,bi,bj)
199     qicen(2)= Qice2(i,j,bi,bj)
200     ELSE
201     qicen(1)= -cpwater*Tmlt1
202     & + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)
203     qicen(2)= -cpice *Tf + Lfresh
204     ENDIF
205     qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0
206     CALL THSICE_START( myThid,
207     I esurp, qicAv, Tf,
208     O qNewIce, ffresh, fsalt,
209     U oceTs, compact, hIce, hSnow )
210     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
211     IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,qNewIce,fsalt,ffresh='
212     & ,compact,qNewIce,fsalt,ffresh
213     #ifdef CHECK_ENERGY_CONSERV
214     flxAtm = 0.
215     frwAtm = 0.
216     CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,
217     I iceFrac, compact, hIce, hSnow, qicen,
218     I qNewIce, ffresh, fsalt, flxAtm, frwAtm,
219     I myTime, myIter, myThid )
220     #endif /* CHECK_ENERGY_CONSERV */
221     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
222     C-- Update Sea-Ice state :
223     IF ( compact.GT.0. _d 0 .AND. iceFrac.EQ.0. _d 0) THEN
224     Tsrf(i,j,bi,bj) = Tf
225     Tice1(i,j,bi,bj) = Tf
226     Tice2(i,j,bi,bj) = Tf
227     Qice1(i,j,bi,bj) = qicen(1)
228     Qice2(i,j,bi,bj) = qicen(2)
229     ENDIF
230     iceheight(i,j,bi,bj) = hIce
231     snowheight(i,j,bi,bj)= hSnow
232     C-- Net fluxes :
233     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - qNewIce
234     EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- ffresh/rhofw
235     saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt
236    
237     IF (dBug) WRITE(6,1010)'ThSI_FWD:-2- compact,hIc,hSn,Qnet=',
238     & compact,hIce,hSnow,Qnet(i,j,bi,bj)
239     C-- - if esurp > 0 : end
240     ENDIF
241    
242     IF ( compact .GT. 0. _d 0 ) THEN
243     iceMask(i,j,bi,bj)=compact
244 jmc 1.2 IF ( hSnow .EQ. 0. _d 0 ) sage(i,j,bi,bj) = 0. _d 0
245 jmc 1.1 ELSE
246     iceMask(i,j,bi,bj) = 0. _d 0
247     iceHeight(i,j,bi,bj)= 0. _d 0
248     snowHeight(i,j,bi,bj)=0. _d 0
249 jmc 1.2 sage(i,j,bi,bj) = 0. _d 0
250     Tsrf(i,j,bi,bj) = oceTs
251 jmc 1.1 Tice1(i,j,bi,bj) = 0. _d 0
252     Tice2(i,j,bi,bj) = 0. _d 0
253     Qice1(i,j,bi,bj) = 0. _d 0
254     Qice2(i,j,bi,bj) = 0. _d 0
255     ENDIF
256    
257     #ifndef CHECK_ENERGY_CONSERV
258     #ifdef ALLOW_TIMEAVE
259     ICE_qleft_AVE(i,j,bi,bj)=ICE_qleft_AVE(i,j,bi,bj)
260     & + ( Qnet(i,j,bi,bj)
261     & )*thSIce_deltaT
262     ICE_fresh_AVE(i,j,bi,bj)=ICE_fresh_AVE(i,j,bi,bj)
263     & + ( EmPmR(i,j,bi,bj)
264     & )*thSIce_deltaT
265     ICE_salFx_AVE(i,j,bi,bj)=ICE_salFx_AVE(i,j,bi,bj)
266     & +saltFlux(i,j,bi,bj)*thSIce_deltaT
267     #endif /* ALLOW_TIMEAVE */
268     #endif /* CHECK_ENERGY_CONSERV */
269    
270     ENDDO
271     ENDDO
272    
273     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
274     #endif /* ALLOW_THSICE */
275    
276     RETURN
277     END
  ViewVC Help
Powered by ViewVC 1.1.22