/[MITgcm]/MITgcm/pkg/therm_seaice/ice_forcing.F
ViewVC logotype

Contents of /MITgcm/pkg/therm_seaice/ice_forcing.F

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

Revision 1.1 - (show annotations) (download)
Thu Nov 21 19:11:42 2002 UTC (21 years, 5 months ago) by cheisey
Branch: MAIN
CVS Tags: checkpoint47b_post, checkpoint51f_post, checkpoint48i_post, checkpoint51l_post, checkpoint47a_post, checkpoint51k_post, checkpoint48d_post, checkpoint50b_post, checkpoint47i_post, checkpoint51o_post, checkpoint48g_post, branchpoint-genmake2, checkpoint50c_pre, checkpoint50, checkpoint51j_post, branch-exfmods-tag, checkpoint51q_post, checkpoint47e_post, checkpoint50f_post, checkpoint50a_post, checkpoint48e_post, checkpoint47c_post, checkpoint50f_pre, checkpoint52a_pre, checkpoint48b_post, checkpoint47j_post, checkpoint47d_pre, checkpoint50d_pre, checkpoint47h_post, checkpoint51d_post, checkpoint51, checkpoint51r_post, checkpoint48c_pre, checkpoint52, checkpoint52b_pre, checkpoint48c_post, checkpoint50d_post, checkpoint51o_pre, checkpoint47f_post, checkpoint51t_post, checkpoint51b_pre, checkpoint52a_post, checkpoint51i_post, checkpoint50e_post, checkpoint47g_post, checkpoint50h_post, checkpoint50c_post, checkpoint51a_post, checkpoint51n_pre, checkpoint47d_post, checkpoint50e_pre, checkpoint50b_pre, checkpoint48d_pre, checkpoint50i_post, checkpoint51p_post, checkpoint51n_post, checkpoint51e_post, checkpoint51b_post, checkpoint48a_post, checkpoint51h_pre, checkpoint48f_post, checkpoint51i_pre, checkpoint51l_pre, checkpoint50g_post, checkpoint51u_post, checkpoint51c_post, checkpoint51g_post, checkpoint51m_post, ecco_c52_e35, checkpoint48, checkpoint49, checkpoint51f_pre, checkpoint51s_post, checkpoint48h_post
Branch point for: checkpoint51n_branch, branch-nonh, tg2-branch, branch-genmake2, branch-exfmods-curt 
Two packages:  bulk_force (Bulk forcing)
and therm_seaice (thermodynamic_seaice) - adopted from LANL CICE.v2.0.2
Earlier integration from Stephaine Dutkiewicz
and Patrick Heimbach.

Two ifdef statements for compile time,
ALLOW_THERM_SEAICE and ALLOW_BULK_FORCE

Two switches in data.pkg to turn on at run-time:

cat data.pkg
# Packages
 &PACKAGES
 useBulkForce=.TRUE.,
 useThermSeaIce=.TRUE.,
 &

WARNING:  useSEAICE and useThermSEAICE are mutually exclusive.

The bulk package requires an additional parameter file
with two namelists, data.ice and data.blk.

c ADAPTED FROM:
c LANL CICE.v2.0.2
c-----------------------------------------------------------------------
c.. thermodynamics (vertical physics) based on M. Winton 3-layer model
c.. See Bitz, C. M. and W. H. Lipscomb, 1999:  "An energy-conserving
c..       thermodynamic sea ice model for climate study."  J. Geophys.
c..       Res., 104, 15669 - 15677.
c..     Winton, M., 1999:  "A reformulated three-layer sea ice model."
c..       Submitted to J. Atmos. Ocean. Technol.

c.. authors Elizabeth C. Hunke and William Lipscomb
c..         Fluid Dynamics Group, Los Alamos National Laboratory
c-----------------------------------------------------------------------
1 cswdice ---- takes place of external_forcing_surf when
2 c the ice model is working
3
4 #include "CPP_OPTIONS.h"
5
6 C !ROUTINE: ICE_FORCING
7 C !INTERFACE:
8 SUBROUTINE ICE_FORCING(
9 I bi, bj, iMin, iMax, jMin, jMax,
10 I myThid )
11 C *==========================================================*
12 C | SUBROUTINE ICE_FORCING
13 C | o Determines forcing terms based on external fields
14 C | relaxation terms etc, including effects of ice.
15 C | cswd - Apr 02 -- add influence of partial ice cover
16 C *==========================================================*
17
18 C !USES:
19 IMPLICIT NONE
20 C === Global variables ===
21 #include "SIZE.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "FFIELDS.h"
25 #include "DYNVARS.h"
26 #include "GRID.h"
27 #ifdef ALLOW_BULK_FORCE
28 #include "BULKF_ICE_CONSTANTS.h"
29 #include "BULKF.h"
30 #include "BULKF_DIAG.h"
31 #ifdef CONSERV_BULKF
32 #include "BULKF_CONSERV.h"
33 #endif
34 #ifdef ALLOW_THERM_SEAICE
35 #include "ICE.h"
36 #include "ICE_DIAGS.h"
37 #endif
38 #endif
39
40
41 C !INPUT/OUTPUT PARAMETERS:
42 C === Routine arguments ===
43 C myThid :: Thread no. that called this routine.
44 INTEGER myThid
45 INTEGER bi,bj
46 INTEGER iMin, iMax
47 INTEGER jMin, jMax
48
49 C !LOCAL VARIABLES:
50 C === Local variables ===
51 INTEGER i,j
52 _RL qleft
53 _RL fsalt
54 _RL ffresh
55 _RL Tf, cphm, frzmlt, compact
56
57 #ifdef ALLOW_THERM_SEAICE
58
59
60 DO j = jMin, jMax
61 DO i = iMin, iMax
62
63 Tf=-mu_Tf*salt(i,j,1,bi,bj)
64 cQQ from v3.0
65 cphm = cpwater*rhosw*drF(1)*hFacC(i,j,1,bi,bj)
66 frzmlt = (Tf-theta(i,j,1,bi,bj))*cphm/deltaTtracer
67 cQQ from v2
68 cQQ frzmlt = (Tf-theta(i,j,1,bi,bj))*cpwater*30.d0/dt
69 frzmlt = min(max(frzmlt,-1.d4),1.d4)
70
71 IF (iceMask(i,j,bi,bj).eq.0.and.frzmlt.lt.1e-6) THEN
72 c no freezing
73 c Zonal wind stress fu:
74 surfaceTendencyU(i,j,bi,bj) = fu(i,j,bi,bj)
75 & *horiVertRatio*recip_rhoNil*recip_dRf(1)
76 c Meridional wind stress fv:
77 surfaceTendencyV(i,j,bi,bj) = fV(i,j,bi,bj)
78 & *horiVertRatio*recip_rhoNil*recip_dRf(1)
79 c Net heat flux Qnet:
80 surfaceTendencyT(i,j,bi,bj) = -Qnet(i,j,bi,bj)
81 & *recip_Cp*recip_rhoNil*recip_dRf(1)
82 c relax to sst equatorward of relaxlat
83 if (abs(yC(i,j,bi,bj)).le.relaxlat) then
84 surfaceTendencyT(i,j,bi,bj) =
85 & surfaceTendencyT(i,j,bi,bj)
86 & - lambdaThetaClimRelax*
87 & (theta(i,j,1,bi,bj)-SST(i,j,bi,bj))
88 endif
89
90 c Net salt flux
91 #ifdef USE_NATURAL_BCS
92 c Freshwater flux EmPmR:
93 surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)
94 & *recip_dRf(1)*salt(i,j,1,bi,bj)
95 #else
96 c Freshwater flux EmPmR:
97 surfaceTendencyS(i,j,bi,bj) = EmPmR(i,j,bi,bj)
98 & *recip_dRf(1)*35.
99 #endif
100 if (abs(yC(i,j,bi,bj)).le.relaxlat) then
101 surfaceTendencyS(i,j,bi,bj) =
102 & surfaceTendencyS(i,j,bi,bj)
103 & - lambdaSaltClimRelax*
104 & (salt(i,j,1,bi,bj)-SSS(i,j,bi,bj))
105 endif
106 iceMask(i,j,bi,bj)=0.0
107 iceHeight(i,j,bi,bj)=0.0
108 snowHeight(i,j,bi,bj)=0.0
109 snow(i,j,bi,bj)=0.0
110 Tsrf(i,j,bi,bj)=theta(i,j,1,bi,bj)
111 Tice1(i,j,bi,bj)=0.0
112 Tice2(i,j,bi,bj)=0.0
113 Qice1(i,j,bi,bj)=0.0
114 Qice2(i,j,bi,bj)=0.0
115
116 #ifdef ALLOW_TIMEAVE
117 ICE_qleft_AVE(i,j,bi,bj)=ICE_qleft_AVE(i,j,bi,bj)
118 & +Qnet(i,j,bi,bj)*deltaTclock
119 ICE_fresh_AVE(i,j,bi,bj)=ICE_fresh_AVE(i,j,bi,bj)
120 & +(EmPmR(i,j,bi,bj))*deltaTclock
121 #endif /*ALLOW_TIMEAVE*/
122 c
123 ELSE
124 qleft=0.D0
125 fsalt=0.D0
126 ffresh=0.D0
127 compact=0.d0
128 if (hfacC(i,j,1,bi,bj).gt.0.d0) then
129 c if new ice should have formed last timestep
130 if (icemask(i,j,bi,bj).eq.0.d0) then
131 call ice_start(i,j,bi, bj, myThid,
132 & frzmlt,ffresh,fsalt,Tf, compact)
133 qleft=0.d0
134 theta(i,j,1,bi,bj)=Tf
135 salt(i,j,1,bi,bj)=salt(i,j,1,bi,bj)+
136 & (ffresh-fsalt)/RHONIL
137 & *recip_dRf(1)*35.*deltaTtracer
138 endif
139 CALL ICE_THERM(i,j,bi,bj,qleft,fsalt,ffresh,
140 & compact, myThid)
141 c create more ice
142 cQQ if (compact.gt.0.d0.and.compact.lt.1.d0.and.
143 cQQ & frzmlt.ge.1e-6) then
144 cQQ call ice_extrastart(i,j,bi, bj, myThid,qleft,
145 cQQ & frzmlt,ffresh,fsalt,Tf,Tair(i,j,bi,bj), compact)
146 cQQ endif
147 endif
148 c redistribute ice if too thin
149 if (compact.gt.0.d0.and.iceHeight(i,j,bi,bj).lt.hicemin) then
150 compact=compact*iceHeight(i,j,bi,bj)/hicemin
151 iceHeight(i,j,bi,bj)=hicemin
152 endif
153 c
154 c Zonal wind stress fu:
155 surfaceTendencyU(i,j,bi,bj) = 0.D0
156 c Meridional wind stress fv:
157 surfaceTendencyV(i,j,bi,bj) = 0.D0
158 c Net heat flux Qnet:
159 surfaceTendencyT(i,j,bi,bj) =
160 & -( compact*qleft + (1.d0-compact)*Qnet(i,j,bi,bj) )
161 & *recip_Cp*recip_rhoNil*recip_dRf(1)
162 cQQQQ relaxing below ice
163 c relax to sst equatorward of relaxlat
164 c if (abs(yC(i,j,bi,bj)).le.relaxlat) then
165 c surfaceTendencyT(i,j,bi,bj) =
166 c & surfaceTendencyT(i,j,bi,bj)
167 c & - lambdaThetaClimRelax*
168 c & (theta(i,j,1,bi,bj)-SST(i,j,bi,bj))
169 c endif
170 cQQQQQ
171
172 c Freshwater flux:
173 ffresh=-(ffresh-fsalt)/RHONIL
174 surfaceTendencyS(i,j,bi,bj) =
175 & ( compact*(ffresh+rain(i,j,bi,bj))
176 & + (1.d0-compact)*EmPmR(i,j,bi,bj))
177 & *recip_dRf(1)*35.
178 cQQQQQ relax below ice
179 c relax to SSS equatorward of relaxlat
180 c if (abs(yC(i,j,bi,bj)).le.relaxlat) then
181 c surfaceTendencyS(i,j,bi,bj) =
182 c & surfaceTendencyS(i,j,bi,bj)
183 c & - lambdaSaltClimRelax*
184 c & (salt(i,j,1,bi,bj)-SSS(i,j,bi,bj))
185 c endif
186 cQQQQQQ
187
188 #ifdef ALLOW_TIMEAVE
189 ICE_qleft_AVE(i,j,bi,bj)=ICE_qleft_AVE(i,j,bi,bj)
190 & +( compact*qleft +
191 & (1.d0-compact)*Qnet(i,j,bi,bj))*deltaTclock
192 ICE_fresh_AVE(i,j,bi,bj)=ICE_fresh_AVE(i,j,bi,bj)
193 & +( compact*(ffresh+rain(i,j,bi,bj)) +
194 & (1.d0-compact)*EmPmR(i,j,bi,bj))*deltaTclock
195 #endif /*ALLOW_TIMEAVE*/
196
197 iceMask(i,j,bi,bj)=compact
198
199 ENDIF !iceMask
200 cswdice --- end add ---
201
202 ENDDO
203 ENDDO
204
205 #endif /*ALLOW_THERM_SEAICE*/
206
207 RETURN
208 END
  ViewVC Help
Powered by ViewVC 1.1.22