1 |
|
|
!Completed |
2 |
|
|
MODULE ocean_slab_mod |
3 |
|
|
! |
4 |
|
|
! This module is used for both surface ocean and sea-ice when using the slab ocean, |
5 |
|
|
! "ocean=slab". |
6 |
|
|
! |
7 |
|
|
|
8 |
|
|
USE dimphy |
9 |
|
|
USE indice_sol_mod |
10 |
|
|
USE surface_data |
11 |
|
|
USE mod_grid_phy_lmdz, ONLY: klon_glo |
12 |
|
|
USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root |
13 |
|
|
|
14 |
|
|
IMPLICIT NONE |
15 |
|
|
PRIVATE |
16 |
|
|
PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice |
17 |
|
|
|
18 |
|
|
!*********************************************************************************** |
19 |
|
|
! Global saved variables |
20 |
|
|
!*********************************************************************************** |
21 |
|
|
! number of slab vertical layers |
22 |
|
|
INTEGER, PUBLIC, SAVE :: nslay |
23 |
|
|
!$OMP THREADPRIVATE(nslay) |
24 |
|
|
! timestep for coupling (update slab temperature) in timesteps |
25 |
|
|
INTEGER, PRIVATE, SAVE :: cpl_pas |
26 |
|
|
!$OMP THREADPRIVATE(cpl_pas) |
27 |
|
|
! cyang = 1/heat capacity of top layer (rho.c.H) |
28 |
|
|
REAL, PRIVATE, SAVE :: cyang |
29 |
|
|
!$OMP THREADPRIVATE(cyang) |
30 |
|
|
! depth of slab layers (1 or 2) |
31 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh |
32 |
|
|
!$OMP THREADPRIVATE(slabh) |
33 |
|
|
! slab temperature |
34 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab |
35 |
|
|
!$OMP THREADPRIVATE(tslab) |
36 |
|
|
! heat flux convergence due to Ekman |
37 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_ekman |
38 |
|
|
!$OMP THREADPRIVATE(dt_ekman) |
39 |
|
|
! heat flux convergence due to horiz diffusion |
40 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_hdiff |
41 |
|
|
!$OMP THREADPRIVATE(dt_hdiff) |
42 |
|
|
! heat flux convergence due to GM eddy advection |
43 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_gm |
44 |
|
|
!$OMP THREADPRIVATE(dt_gm) |
45 |
|
|
! Heat Flux correction |
46 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_qflux |
47 |
|
|
!$OMP THREADPRIVATE(dt_qflux) |
48 |
|
|
! fraction of ocean covered by sea ice (sic / (oce+sic)) |
49 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic |
50 |
|
|
!$OMP THREADPRIVATE(fsic) |
51 |
|
|
! temperature of the sea ice |
52 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice |
53 |
|
|
!$OMP THREADPRIVATE(tice) |
54 |
|
|
! sea ice thickness, in kg/m2 |
55 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice |
56 |
|
|
!$OMP THREADPRIVATE(seaice) |
57 |
|
|
! net surface heat flux, weighted by open ocean fraction |
58 |
|
|
! slab_bils accumulated over cpl_pas timesteps |
59 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum |
60 |
|
|
!$OMP THREADPRIVATE(bils_cum) |
61 |
|
|
! net heat flux into the ocean below the ice : conduction + solar radiation |
62 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg |
63 |
|
|
!$OMP THREADPRIVATE(slab_bilg) |
64 |
|
|
! slab_bilg over cpl_pas timesteps |
65 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum |
66 |
|
|
!$OMP THREADPRIVATE(bilg_cum) |
67 |
|
|
! wind stress saved over cpl_pas timesteps |
68 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: taux_cum |
69 |
|
|
!$OMP THREADPRIVATE(taux_cum) |
70 |
|
|
REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tauy_cum |
71 |
|
|
!$OMP THREADPRIVATE(tauy_cum) |
72 |
|
|
|
73 |
|
|
!*********************************************************************************** |
74 |
|
|
! Parameters (could be read in def file: move to slab_init) |
75 |
|
|
!*********************************************************************************** |
76 |
|
|
! snow and ice physical characteristics: |
77 |
|
|
REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp |
78 |
|
|
REAL, PARAMETER :: t_melt=273.15 ! melting ice temp |
79 |
|
|
REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3 |
80 |
|
|
REAL, PARAMETER :: ice_den=917. ! ice density |
81 |
|
|
REAL, PARAMETER :: sea_den=1025. ! sea water density |
82 |
|
|
REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice |
83 |
|
|
REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow |
84 |
|
|
REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice |
85 |
|
|
REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, snow and ice |
86 |
|
|
REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice |
87 |
|
|
|
88 |
|
|
! control of snow and ice cover & freeze / melt (heights converted to kg/m2) |
89 |
|
|
REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm |
90 |
|
|
REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean |
91 |
|
|
REAL, PARAMETER :: ice_frac_min=0.001 |
92 |
|
|
REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction |
93 |
|
|
REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness |
94 |
|
|
REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness |
95 |
|
|
! below ice_thin, priority is melt lateral / grow height |
96 |
|
|
! ice_thin is also height of new ice |
97 |
|
|
REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness |
98 |
|
|
! above ice_thick, priority is melt height / grow lateral |
99 |
|
|
REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice |
100 |
|
|
REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height |
101 |
|
|
|
102 |
|
|
! albedo and radiation parameters |
103 |
|
|
REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo |
104 |
|
|
REAL, PARAMETER :: alb_sno_del=0.3 !max snow albedo = min + del |
105 |
|
|
REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice |
106 |
|
|
REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice |
107 |
|
|
REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the |
108 |
|
|
! ice (no snow) |
109 |
|
|
REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1) |
110 |
|
|
|
111 |
|
|
! horizontal transport |
112 |
|
|
LOGICAL, PUBLIC, SAVE :: slab_hdiff |
113 |
|
|
!$OMP THREADPRIVATE(slab_hdiff) |
114 |
|
|
LOGICAL, PUBLIC, SAVE :: slab_gm |
115 |
|
|
!$OMP THREADPRIVATE(slab_gm) |
116 |
|
|
REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion |
117 |
|
|
!$OMP THREADPRIVATE(coef_hdiff) |
118 |
|
|
INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment |
119 |
|
|
!$OMP THREADPRIVATE(slab_ekman) |
120 |
|
|
!$OMP THREADPRIVATE(slab_cadj) |
121 |
|
|
|
122 |
|
|
!*********************************************************************************** |
123 |
|
|
|
124 |
|
|
CONTAINS |
125 |
|
|
! |
126 |
|
|
!*********************************************************************************** |
127 |
|
|
! |
128 |
|
|
SUBROUTINE ocean_slab_init(dtime, pctsrf_rst) |
129 |
|
|
!, seaice_rst etc |
130 |
|
|
|
131 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
132 |
|
|
USE mod_phys_lmdz_transfert_para, ONLY : gather |
133 |
|
|
USE slab_heat_transp_mod, ONLY : ini_slab_transp |
134 |
|
|
|
135 |
|
|
! Input variables |
136 |
|
|
!*********************************************************************************** |
137 |
|
|
REAL, INTENT(IN) :: dtime |
138 |
|
|
! Variables read from restart file |
139 |
|
|
REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst |
140 |
|
|
! surface fractions from start file |
141 |
|
|
|
142 |
|
|
! Local variables |
143 |
|
|
!************************************************************************************ |
144 |
|
|
INTEGER :: error |
145 |
|
|
REAL, DIMENSION(klon_glo) :: zmasq_glo |
146 |
|
|
CHARACTER (len = 80) :: abort_message |
147 |
|
|
CHARACTER (len = 20) :: modname = 'ocean_slab_intit' |
148 |
|
|
|
149 |
|
|
!*********************************************************************************** |
150 |
|
|
! Define some parameters |
151 |
|
|
!*********************************************************************************** |
152 |
|
|
! Number of slab layers |
153 |
|
|
nslay=2 |
154 |
|
|
CALL getin_p('slab_layers',nslay) |
155 |
|
|
print *,'number of slab layers : ',nslay |
156 |
|
|
! Layer thickness |
157 |
|
|
ALLOCATE(slabh(nslay), stat = error) |
158 |
|
|
IF (error /= 0) THEN |
159 |
|
|
abort_message='Pb allocation slabh' |
160 |
|
|
CALL abort_physic(modname,abort_message,1) |
161 |
|
|
ENDIF |
162 |
|
|
slabh(1)=50. |
163 |
|
|
CALL getin_p('slab_depth',slabh(1)) |
164 |
|
|
IF (nslay.GT.1) THEN |
165 |
|
|
slabh(2)=150. |
166 |
|
|
END IF |
167 |
|
|
|
168 |
|
|
! cyang = 1/heat capacity of top layer (rho.c.H) |
169 |
|
|
cyang=1/(slabh(1)*sea_den*sea_cap) |
170 |
|
|
|
171 |
|
|
! cpl_pas coupling period (update of tslab and ice fraction) |
172 |
|
|
! pour un calcul a chaque pas de temps, cpl_pas=1 |
173 |
|
|
cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour |
174 |
|
|
CALL getin_p('cpl_pas',cpl_pas) |
175 |
|
|
print *,'cpl_pas',cpl_pas |
176 |
|
|
|
177 |
|
|
! Horizontal diffusion |
178 |
|
|
slab_hdiff=.FALSE. |
179 |
|
|
CALL getin_p('slab_hdiff',slab_hdiff) |
180 |
|
|
coef_hdiff=25000. |
181 |
|
|
CALL getin_p('coef_hdiff',coef_hdiff) |
182 |
|
|
! Ekman transport |
183 |
|
|
slab_ekman=0 |
184 |
|
|
CALL getin_p('slab_ekman',slab_ekman) |
185 |
|
|
! GM eddy advection (2-layers only) |
186 |
|
|
slab_gm=.FALSE. |
187 |
|
|
CALL getin_p('slab_gm',slab_gm) |
188 |
|
|
IF (slab_ekman.LT.2) THEN |
189 |
|
|
slab_gm=.FALSE. |
190 |
|
|
ENDIF |
191 |
|
|
! Convective adjustment |
192 |
|
|
IF (nslay.EQ.1) THEN |
193 |
|
|
slab_cadj=0 |
194 |
|
|
ELSE |
195 |
|
|
slab_cadj=1 |
196 |
|
|
END IF |
197 |
|
|
CALL getin_p('slab_cadj',slab_cadj) |
198 |
|
|
|
199 |
|
|
!************************************************************************************ |
200 |
|
|
! Allocate surface fraction read from restart file |
201 |
|
|
!************************************************************************************ |
202 |
|
|
ALLOCATE(fsic(klon), stat = error) |
203 |
|
|
IF (error /= 0) THEN |
204 |
|
|
abort_message='Pb allocation tmp_pctsrf_slab' |
205 |
|
|
CALL abort_physic(modname,abort_message,1) |
206 |
|
|
ENDIF |
207 |
|
|
fsic(:)=0. |
208 |
|
|
!zmasq = continent fraction |
209 |
|
|
WHERE (1.-zmasq(:)>EPSFRA) |
210 |
|
|
fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:)) |
211 |
|
|
END WHERE |
212 |
|
|
|
213 |
|
|
!************************************************************************************ |
214 |
|
|
! Allocate saved fields |
215 |
|
|
!************************************************************************************ |
216 |
|
|
ALLOCATE(tslab(klon,nslay), stat=error) |
217 |
|
|
IF (error /= 0) CALL abort_physic & |
218 |
|
|
(modname,'pb allocation tslab', 1) |
219 |
|
|
|
220 |
|
|
ALLOCATE(bils_cum(klon), stat = error) |
221 |
|
|
IF (error /= 0) THEN |
222 |
|
|
abort_message='Pb allocation slab_bils_cum' |
223 |
|
|
CALL abort_physic(modname,abort_message,1) |
224 |
|
|
ENDIF |
225 |
|
|
bils_cum(:) = 0.0 |
226 |
|
|
|
227 |
|
|
IF (version_ocean=='sicINT') THEN ! interactive sea ice |
228 |
|
|
ALLOCATE(slab_bilg(klon), stat = error) |
229 |
|
|
IF (error /= 0) THEN |
230 |
|
|
abort_message='Pb allocation slab_bilg' |
231 |
|
|
CALL abort_physic(modname,abort_message,1) |
232 |
|
|
ENDIF |
233 |
|
|
slab_bilg(:) = 0.0 |
234 |
|
|
ALLOCATE(bilg_cum(klon), stat = error) |
235 |
|
|
IF (error /= 0) THEN |
236 |
|
|
abort_message='Pb allocation slab_bilg_cum' |
237 |
|
|
CALL abort_physic(modname,abort_message,1) |
238 |
|
|
ENDIF |
239 |
|
|
bilg_cum(:) = 0.0 |
240 |
|
|
ALLOCATE(tice(klon), stat = error) |
241 |
|
|
IF (error /= 0) THEN |
242 |
|
|
abort_message='Pb allocation slab_tice' |
243 |
|
|
CALL abort_physic(modname,abort_message,1) |
244 |
|
|
ENDIF |
245 |
|
|
ALLOCATE(seaice(klon), stat = error) |
246 |
|
|
IF (error /= 0) THEN |
247 |
|
|
abort_message='Pb allocation slab_seaice' |
248 |
|
|
CALL abort_physic(modname,abort_message,1) |
249 |
|
|
ENDIF |
250 |
|
|
END IF |
251 |
|
|
|
252 |
|
|
IF (slab_hdiff) THEN !horizontal diffusion |
253 |
|
|
ALLOCATE(dt_hdiff(klon,nslay), stat = error) |
254 |
|
|
IF (error /= 0) THEN |
255 |
|
|
abort_message='Pb allocation dt_hdiff' |
256 |
|
|
CALL abort_physic(modname,abort_message,1) |
257 |
|
|
ENDIF |
258 |
|
|
dt_hdiff(:,:) = 0.0 |
259 |
|
|
ENDIF |
260 |
|
|
|
261 |
|
|
ALLOCATE(dt_qflux(klon,nslay), stat = error) |
262 |
|
|
IF (error /= 0) THEN |
263 |
|
|
abort_message='Pb allocation dt_qflux' |
264 |
|
|
CALL abort_physic(modname,abort_message,1) |
265 |
|
|
ENDIF |
266 |
|
|
dt_qflux(:,:) = 0.0 |
267 |
|
|
|
268 |
|
|
IF (slab_gm) THEN !GM advection |
269 |
|
|
ALLOCATE(dt_gm(klon,nslay), stat = error) |
270 |
|
|
IF (error /= 0) THEN |
271 |
|
|
abort_message='Pb allocation dt_gm' |
272 |
|
|
CALL abort_physic(modname,abort_message,1) |
273 |
|
|
ENDIF |
274 |
|
|
dt_gm(:,:) = 0.0 |
275 |
|
|
ENDIF |
276 |
|
|
|
277 |
|
|
IF (slab_ekman.GT.0) THEN ! ekman transport |
278 |
|
|
ALLOCATE(dt_ekman(klon,nslay), stat = error) |
279 |
|
|
IF (error /= 0) THEN |
280 |
|
|
abort_message='Pb allocation dt_ekman' |
281 |
|
|
CALL abort_physic(modname,abort_message,1) |
282 |
|
|
ENDIF |
283 |
|
|
dt_ekman(:,:) = 0.0 |
284 |
|
|
ALLOCATE(taux_cum(klon), stat = error) |
285 |
|
|
IF (error /= 0) THEN |
286 |
|
|
abort_message='Pb allocation taux_cum' |
287 |
|
|
CALL abort_physic(modname,abort_message,1) |
288 |
|
|
ENDIF |
289 |
|
|
taux_cum(:) = 0.0 |
290 |
|
|
ALLOCATE(tauy_cum(klon), stat = error) |
291 |
|
|
IF (error /= 0) THEN |
292 |
|
|
abort_message='Pb allocation tauy_cum' |
293 |
|
|
CALL abort_physic(modname,abort_message,1) |
294 |
|
|
ENDIF |
295 |
|
|
tauy_cum(:) = 0.0 |
296 |
|
|
ENDIF |
297 |
|
|
|
298 |
|
|
! Initialize transport |
299 |
|
|
IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN |
300 |
|
|
CALL gather(zmasq,zmasq_glo) |
301 |
|
|
! Master thread/process only |
302 |
|
|
!$OMP MASTER |
303 |
|
|
IF (is_mpi_root) THEN |
304 |
|
|
CALL ini_slab_transp(zmasq_glo) |
305 |
|
|
END IF |
306 |
|
|
!$OMP END MASTER |
307 |
|
|
END IF |
308 |
|
|
|
309 |
|
|
END SUBROUTINE ocean_slab_init |
310 |
|
|
! |
311 |
|
|
!*********************************************************************************** |
312 |
|
|
! |
313 |
|
|
SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified) |
314 |
|
|
|
315 |
|
|
! this routine sends back the sea ice and ocean fraction to the main physics |
316 |
|
|
! routine. Called only with interactive sea ice |
317 |
|
|
|
318 |
|
|
! Arguments |
319 |
|
|
!************************************************************************************ |
320 |
|
|
INTEGER, INTENT(IN) :: itime ! current timestep |
321 |
|
|
INTEGER, INTENT(IN) :: jour ! day in year (not |
322 |
|
|
REAL , INTENT(IN) :: dtime ! physics timestep (s) |
323 |
|
|
REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction |
324 |
|
|
LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is |
325 |
|
|
! modified at this time step |
326 |
|
|
|
327 |
|
|
pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:)) |
328 |
|
|
pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:)) |
329 |
|
|
is_modified=.TRUE. |
330 |
|
|
|
331 |
|
|
END SUBROUTINE ocean_slab_frac |
332 |
|
|
! |
333 |
|
|
!************************************************************************************ |
334 |
|
|
! |
335 |
|
|
SUBROUTINE ocean_slab_noice( & |
336 |
|
|
itime, dtime, jour, knon, knindex, & |
337 |
|
|
p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, & |
338 |
|
|
AcoefH, AcoefQ, BcoefH, BcoefQ, & |
339 |
|
|
AcoefU, AcoefV, BcoefU, BcoefV, & |
340 |
|
|
ps, u1, v1, gustiness, tsurf_in, & |
341 |
|
|
radsol, snow, & |
342 |
|
|
qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & |
343 |
|
|
tsurf_new, dflux_s, dflux_l, slab_bils) |
344 |
|
|
|
345 |
|
|
USE calcul_fluxs_mod |
346 |
|
|
USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman1,slab_ekman2,slab_gmdiff |
347 |
|
|
USE mod_phys_lmdz_para |
348 |
|
|
|
349 |
|
|
INCLUDE "clesphys.h" |
350 |
|
|
|
351 |
|
|
! This routine |
352 |
|
|
! (1) computes surface turbulent fluxes over points with some open ocean |
353 |
|
|
! (2) reads additional Q-flux (everywhere) |
354 |
|
|
! (3) computes horizontal transport (diffusion & Ekman) |
355 |
|
|
! (4) updates slab temperature every cpl_pas ; creates new ice if needed. |
356 |
|
|
|
357 |
|
|
! Note : |
358 |
|
|
! klon total number of points |
359 |
|
|
! knon number of points with open ocean (varies with time) |
360 |
|
|
! knindex gives position of the knon points within klon. |
361 |
|
|
! In general, local saved variables have klon values |
362 |
|
|
! variables exchanged with PBL module have knon. |
363 |
|
|
|
364 |
|
|
! Input arguments |
365 |
|
|
!*********************************************************************************** |
366 |
|
|
INTEGER, INTENT(IN) :: itime ! current timestep INTEGER, |
367 |
|
|
INTEGER, INTENT(IN) :: jour ! day in year (for Q-Flux) |
368 |
|
|
INTEGER, INTENT(IN) :: knon ! number of points |
369 |
|
|
INTEGER, DIMENSION(klon), INTENT(IN) :: knindex |
370 |
|
|
REAL, INTENT(IN) :: dtime ! timestep (s) |
371 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: p1lay |
372 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm |
373 |
|
|
! drag coefficients |
374 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow |
375 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum ! near surface T, q |
376 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ |
377 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV |
378 |
|
|
! exchange coefficients for boundary layer scheme |
379 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: ps ! surface pressure |
380 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! surface wind |
381 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! surface temperature |
382 |
|
|
REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux |
383 |
|
|
|
384 |
|
|
! In/Output arguments |
385 |
|
|
!************************************************************************************ |
386 |
|
|
REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2 |
387 |
|
|
|
388 |
|
|
! Output arguments |
389 |
|
|
!************************************************************************************ |
390 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: qsurf |
391 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat |
392 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 |
393 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture |
394 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l |
395 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: slab_bils |
396 |
|
|
|
397 |
|
|
! Local variables |
398 |
|
|
!************************************************************************************ |
399 |
|
|
INTEGER :: i,ki,k |
400 |
|
|
REAL :: t_cadj |
401 |
|
|
! for surface heat fluxes |
402 |
|
|
REAL, DIMENSION(klon) :: cal, beta, dif_grnd |
403 |
|
|
! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes |
404 |
|
|
REAL, DIMENSION(klon) :: diff_sst, diff_siv |
405 |
|
|
REAL, DIMENSION(klon,nslay) :: lmt_bils |
406 |
|
|
! for surface wind stress |
407 |
|
|
REAL, DIMENSION(klon) :: u0, v0 |
408 |
|
|
REAL, DIMENSION(klon) :: u1_lay, v1_lay |
409 |
|
|
! for new ice creation |
410 |
|
|
REAL :: e_freeze, h_new, dfsic |
411 |
|
|
! horizontal diffusion and Ekman local vars |
412 |
|
|
! dimension = global domain (klon_glo) instead of // subdomains |
413 |
|
|
REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo |
414 |
|
|
! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop |
415 |
|
|
REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp |
416 |
|
|
REAL, DIMENSION(klon_glo,nslay) :: tslab_glo |
417 |
|
|
REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo |
418 |
|
|
|
419 |
|
|
!**************************************************************************************** |
420 |
|
|
! 1) Surface fluxes calculation |
421 |
|
|
! |
422 |
|
|
!**************************************************************************************** |
423 |
|
|
!cal(:) = 0. ! infinite thermal inertia |
424 |
|
|
!beta(:) = 1. ! wet surface |
425 |
|
|
!dif_grnd(:) = 0. ! no diffusion into ground |
426 |
|
|
! EV: use calbeta |
427 |
|
|
CALL calbeta(dtime, is_oce, knon, snow,qsurf, beta, cal, dif_grnd) |
428 |
|
|
|
429 |
|
|
|
430 |
|
|
|
431 |
|
|
! Suppose zero surface speed |
432 |
|
|
u0(:)=0.0 |
433 |
|
|
v0(:)=0.0 |
434 |
|
|
u1_lay(:) = u1(:) - u0(:) |
435 |
|
|
v1_lay(:) = v1(:) - v0(:) |
436 |
|
|
|
437 |
|
|
! Compute latent & sensible fluxes |
438 |
|
|
CALL calcul_fluxs(knon, is_oce, dtime, & |
439 |
|
|
tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, & |
440 |
|
|
precip_rain, precip_snow, snow, qsurf, & |
441 |
|
|
radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & |
442 |
|
|
f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & |
443 |
|
|
tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) |
444 |
|
|
|
445 |
|
|
! save total cumulated heat fluxes locally |
446 |
|
|
! radiative + turbulent + melt of falling snow |
447 |
|
|
slab_bils(:)=0. |
448 |
|
|
DO i=1,knon |
449 |
|
|
ki=knindex(i) |
450 |
|
|
slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) & |
451 |
|
|
-precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki))) |
452 |
|
|
bils_cum(ki)=bils_cum(ki)+slab_bils(ki) |
453 |
|
|
END DO |
454 |
|
|
|
455 |
|
|
! Compute surface wind stress |
456 |
|
|
CALL calcul_flux_wind(knon, dtime, & |
457 |
|
|
u0, v0, u1, v1, gustiness, cdragm, & |
458 |
|
|
AcoefU, AcoefV, BcoefU, BcoefV, & |
459 |
|
|
p1lay, temp_air, & |
460 |
|
|
flux_u1, flux_v1) |
461 |
|
|
|
462 |
|
|
! save cumulated wind stress |
463 |
|
|
IF (slab_ekman.GT.0) THEN |
464 |
|
|
DO i=1,knon |
465 |
|
|
ki=knindex(i) |
466 |
|
|
taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas |
467 |
|
|
tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas |
468 |
|
|
END DO |
469 |
|
|
ENDIF |
470 |
|
|
|
471 |
|
|
!**************************************************************************************** |
472 |
|
|
! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc |
473 |
|
|
! |
474 |
|
|
!**************************************************************************************** |
475 |
|
|
CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) |
476 |
|
|
! lmt_bils and diff_sst,siv saved by limit_slab |
477 |
|
|
! qflux = total QFlux correction (in W/m2) |
478 |
|
|
dt_qflux(:,1)=lmt_bils(:,1)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400. |
479 |
|
|
IF (nslay.GT.1) THEN |
480 |
|
|
dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay) |
481 |
|
|
END IF |
482 |
|
|
|
483 |
|
|
!**************************************************************************************** |
484 |
|
|
! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport) |
485 |
|
|
! Bring to freezing temp and make sea ice if necessary |
486 |
|
|
! |
487 |
|
|
!***********************************************o***************************************** |
488 |
|
|
tsurf_new=tsurf_in |
489 |
|
|
IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction |
490 |
|
|
! *********************************** |
491 |
|
|
! Horizontal transport |
492 |
|
|
! *********************************** |
493 |
|
|
IF (slab_ekman.GT.0) THEN |
494 |
|
|
! copy wind stress to global var |
495 |
|
|
CALL gather(taux_cum,taux_glo) |
496 |
|
|
CALL gather(tauy_cum,tauy_glo) |
497 |
|
|
END IF |
498 |
|
|
|
499 |
|
|
IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN |
500 |
|
|
CALL gather(tslab,tslab_glo) |
501 |
|
|
! Compute horiz transport on one process only |
502 |
|
|
IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus |
503 |
|
|
IF (slab_hdiff) THEN |
504 |
|
|
dt_hdiff_glo(:,:)=0. |
505 |
|
|
END IF |
506 |
|
|
IF (slab_ekman.GT.0) THEN |
507 |
|
|
dt_ekman_glo(:,:)=0. |
508 |
|
|
END IF |
509 |
|
|
IF (slab_gm) THEN |
510 |
|
|
dt_gm_glo(:,:)=0. |
511 |
|
|
END IF |
512 |
|
|
DO i=1,cpl_pas ! time splitting for numerical stability |
513 |
|
|
IF (slab_ekman.GT.0) THEN |
514 |
|
|
SELECT CASE (slab_ekman) |
515 |
|
|
CASE (1) |
516 |
|
|
CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp) |
517 |
|
|
CASE (2) |
518 |
|
|
CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm) |
519 |
|
|
CASE DEFAULT |
520 |
|
|
dt_ekman_tmp(:,:)=0. |
521 |
|
|
END SELECT |
522 |
|
|
dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:) |
523 |
|
|
! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1 |
524 |
|
|
DO k=1,nslay |
525 |
|
|
dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den) |
526 |
|
|
ENDDO |
527 |
|
|
tslab_glo=tslab_glo+dt_ekman_tmp*dtime |
528 |
|
|
IF (slab_gm) THEN ! Gent-McWilliams eddy advection |
529 |
|
|
dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) |
530 |
|
|
! convert dt from K.s-1.(kg.m-2) to K.s-1 |
531 |
|
|
DO k=1,nslay |
532 |
|
|
dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den) |
533 |
|
|
END DO |
534 |
|
|
tslab_glo=tslab_glo+dt_hdiff_tmp*dtime |
535 |
|
|
END IF |
536 |
|
|
ENDIF |
537 |
|
|
! GM included in Ekman_2 |
538 |
|
|
! IF (slab_gm) THEN ! Gent-McWilliams eddy advection |
539 |
|
|
! CALL slab_gmdiff(tslab_glo,dt_hdiff_tmp) |
540 |
|
|
! ! convert dt_gm from K.m.s-1 to K.s-1 |
541 |
|
|
! DO k=1,nslay |
542 |
|
|
! dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/slabh(k) |
543 |
|
|
! END DO |
544 |
|
|
! tslab_glo=tslab_glo+dt_hdiff_tmp*dtime |
545 |
|
|
! dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) |
546 |
|
|
! END IF |
547 |
|
|
IF (slab_hdiff) THEN ! horizontal diffusion |
548 |
|
|
! laplacian of slab T |
549 |
|
|
CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp) |
550 |
|
|
! multiply by diff coef and normalize to 50m slab equivalent |
551 |
|
|
dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh) |
552 |
|
|
dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:) |
553 |
|
|
tslab_glo=tslab_glo+dt_hdiff_tmp*dtime |
554 |
|
|
END IF |
555 |
|
|
END DO ! time splitting |
556 |
|
|
IF (slab_hdiff) THEN |
557 |
|
|
!dt_hdiff_glo saved in W/m2 |
558 |
|
|
DO k=1,nslay |
559 |
|
|
dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas |
560 |
|
|
END DO |
561 |
|
|
END IF |
562 |
|
|
IF (slab_gm) THEN |
563 |
|
|
!dt_hdiff_glo saved in W/m2 |
564 |
|
|
dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas |
565 |
|
|
END IF |
566 |
|
|
IF (slab_ekman.GT.0) THEN |
567 |
|
|
! dt_ekman_glo saved in W/m2 |
568 |
|
|
dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas |
569 |
|
|
END IF |
570 |
|
|
END IF ! master process |
571 |
|
|
!$OMP BARRIER |
572 |
|
|
! Send new fields back to all processes |
573 |
|
|
CALL Scatter(tslab_glo,tslab) |
574 |
|
|
IF (slab_hdiff) THEN |
575 |
|
|
CALL Scatter(dt_hdiff_glo,dt_hdiff) |
576 |
|
|
END IF |
577 |
|
|
IF (slab_gm) THEN |
578 |
|
|
CALL Scatter(dt_gm_glo,dt_gm) |
579 |
|
|
END IF |
580 |
|
|
IF (slab_ekman.GT.0) THEN |
581 |
|
|
CALL Scatter(dt_ekman_glo,dt_ekman) |
582 |
|
|
! clear wind stress |
583 |
|
|
taux_cum(:)=0. |
584 |
|
|
tauy_cum(:)=0. |
585 |
|
|
END IF |
586 |
|
|
ENDIF ! transport |
587 |
|
|
|
588 |
|
|
! *********************************** |
589 |
|
|
! Other heat fluxes |
590 |
|
|
! *********************************** |
591 |
|
|
! Add read QFlux |
592 |
|
|
DO k=1,nslay |
593 |
|
|
tslab(:,k)=tslab(:,k)+dt_qflux(:,k)*cyang*dtime*cpl_pas & |
594 |
|
|
*slabh(1)/slabh(k) |
595 |
|
|
END DO |
596 |
|
|
! Add cumulated surface fluxes |
597 |
|
|
tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime |
598 |
|
|
! Convective adjustment if 2 layers |
599 |
|
|
IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN |
600 |
|
|
DO i=1,klon |
601 |
|
|
IF (tslab(i,2).GT.tslab(i,1)) THEN |
602 |
|
|
! mean (mass-weighted) temperature |
603 |
|
|
t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:)) |
604 |
|
|
tslab(i,1)=t_cadj |
605 |
|
|
tslab(i,2)=t_cadj |
606 |
|
|
END IF |
607 |
|
|
END DO |
608 |
|
|
END IF |
609 |
|
|
! *********************************** |
610 |
|
|
! Update surface temperature and ice |
611 |
|
|
! *********************************** |
612 |
|
|
SELECT CASE(version_ocean) |
613 |
|
|
CASE('sicNO') ! no sea ice even below freezing ! |
614 |
|
|
DO i=1,knon |
615 |
|
|
ki=knindex(i) |
616 |
|
|
tsurf_new(i)=tslab(ki,1) |
617 |
|
|
END DO |
618 |
|
|
CASE('sicOBS') ! "realistic" case, for prescribed sea ice |
619 |
|
|
! tslab cannot be below freezing, or above it if there is sea ice |
620 |
|
|
DO i=1,knon |
621 |
|
|
ki=knindex(i) |
622 |
|
|
IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN |
623 |
|
|
tslab(ki,1)=t_freeze |
624 |
|
|
END IF |
625 |
|
|
tsurf_new(i)=tslab(ki,1) |
626 |
|
|
END DO |
627 |
|
|
CASE('sicINT') ! interactive sea ice |
628 |
|
|
DO i=1,knon |
629 |
|
|
ki=knindex(i) |
630 |
|
|
IF (fsic(ki).LT.epsfra) THEN ! Free of ice |
631 |
|
|
IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice |
632 |
|
|
! quantity of new ice formed |
633 |
|
|
e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat |
634 |
|
|
! new ice |
635 |
|
|
tice(ki)=t_freeze |
636 |
|
|
fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin) |
637 |
|
|
IF (fsic(ki).GT.ice_frac_min) THEN |
638 |
|
|
seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max) |
639 |
|
|
tslab(ki,1)=t_freeze |
640 |
|
|
ELSE |
641 |
|
|
fsic(ki)=0. |
642 |
|
|
END IF |
643 |
|
|
tsurf_new(i)=t_freeze |
644 |
|
|
ELSE |
645 |
|
|
tsurf_new(i)=tslab(ki,1) |
646 |
|
|
END IF |
647 |
|
|
ELSE ! ice present |
648 |
|
|
tsurf_new(i)=t_freeze |
649 |
|
|
IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice |
650 |
|
|
! quantity of new ice formed over open ocean |
651 |
|
|
e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) & |
652 |
|
|
/(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) |
653 |
|
|
! new ice height and fraction |
654 |
|
|
h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new |
655 |
|
|
dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new) |
656 |
|
|
h_new=MIN(e_freeze/dfsic,h_ice_max) |
657 |
|
|
! update tslab to freezing over open ocean only |
658 |
|
|
tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki)) |
659 |
|
|
! update sea ice |
660 |
|
|
seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) & |
661 |
|
|
/(dfsic+fsic(ki)) |
662 |
|
|
fsic(ki)=fsic(ki)+dfsic |
663 |
|
|
! update snow? |
664 |
|
|
END IF ! tslab below freezing |
665 |
|
|
END IF ! sea ice present |
666 |
|
|
END DO |
667 |
|
|
END SELECT |
668 |
|
|
bils_cum(:)=0.0! clear cumulated fluxes |
669 |
|
|
END IF ! coupling time |
670 |
|
|
END SUBROUTINE ocean_slab_noice |
671 |
|
|
! |
672 |
|
|
!**************************************************************************************** |
673 |
|
|
|
674 |
|
|
SUBROUTINE ocean_slab_ice( & |
675 |
|
|
itime, dtime, jour, knon, knindex, & |
676 |
|
|
tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, & |
677 |
|
|
AcoefH, AcoefQ, BcoefH, BcoefQ, & |
678 |
|
|
AcoefU, AcoefV, BcoefU, BcoefV, & |
679 |
|
|
ps, u1, v1, gustiness, & |
680 |
|
|
radsol, snow, qsurf, qsol, agesno, & |
681 |
|
|
alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & |
682 |
|
|
tsurf_new, dflux_s, dflux_l, swnet) |
683 |
|
|
|
684 |
|
|
USE calcul_fluxs_mod |
685 |
|
|
|
686 |
|
|
INCLUDE "YOMCST.h" |
687 |
|
|
INCLUDE "clesphys.h" |
688 |
|
|
|
689 |
|
|
! Input arguments |
690 |
|
|
!**************************************************************************************** |
691 |
|
|
INTEGER, INTENT(IN) :: itime, jour, knon |
692 |
|
|
INTEGER, DIMENSION(klon), INTENT(IN) :: knindex |
693 |
|
|
REAL, INTENT(IN) :: dtime |
694 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in |
695 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: p1lay |
696 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm |
697 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow |
698 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum |
699 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ |
700 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV |
701 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: ps |
702 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness |
703 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: swnet |
704 |
|
|
|
705 |
|
|
! In/Output arguments |
706 |
|
|
!**************************************************************************************** |
707 |
|
|
REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol |
708 |
|
|
REAL, DIMENSION(klon), INTENT(INOUT) :: agesno |
709 |
|
|
REAL, DIMENSION(klon), INTENT(INOUT) :: radsol |
710 |
|
|
|
711 |
|
|
! Output arguments |
712 |
|
|
!**************************************************************************************** |
713 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: qsurf |
714 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval |
715 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval |
716 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat |
717 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 |
718 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new |
719 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l |
720 |
|
|
|
721 |
|
|
! Local variables |
722 |
|
|
!**************************************************************************************** |
723 |
|
|
INTEGER :: i,ki |
724 |
|
|
REAL, DIMENSION(klon) :: cal, beta, dif_grnd |
725 |
|
|
REAL, DIMENSION(klon) :: u0, v0 |
726 |
|
|
REAL, DIMENSION(klon) :: u1_lay, v1_lay |
727 |
|
|
! intermediate heat fluxes: |
728 |
|
|
REAL :: f_cond, f_swpen |
729 |
|
|
! for snow/ice albedo: |
730 |
|
|
REAL :: alb_snow, alb_ice, alb_pond |
731 |
|
|
REAL :: frac_snow, frac_ice, frac_pond |
732 |
|
|
! for ice melt / freeze |
733 |
|
|
REAL :: e_melt, snow_evap, h_test |
734 |
|
|
! dhsic, dfsic change in ice mass, fraction. |
735 |
|
|
REAL :: dhsic, dfsic, frac_mf |
736 |
|
|
|
737 |
|
|
!**************************************************************************************** |
738 |
|
|
! 1) Flux calculation |
739 |
|
|
!**************************************************************************************** |
740 |
|
|
! Suppose zero surface speed |
741 |
|
|
u0(:)=0.0 |
742 |
|
|
v0(:)=0.0 |
743 |
|
|
u1_lay(:) = u1(:) - u0(:) |
744 |
|
|
v1_lay(:) = v1(:) - v0(:) |
745 |
|
|
|
746 |
|
|
! set beta, cal, compute conduction fluxes inside ice/snow |
747 |
|
|
slab_bilg(:)=0. |
748 |
|
|
!dif_grnd(:)=0. |
749 |
|
|
!beta(:) = 1. |
750 |
|
|
! EV: use calbeta to calculate beta and then recalculate properly cal |
751 |
|
|
CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, cal, dif_grnd) |
752 |
|
|
|
753 |
|
|
|
754 |
|
|
DO i=1,knon |
755 |
|
|
ki=knindex(i) |
756 |
|
|
IF (snow(i).GT.snow_min) THEN |
757 |
|
|
! snow-layer heat capacity |
758 |
|
|
cal(i)=2.*RCPD/(snow(i)*ice_cap) |
759 |
|
|
! snow conductive flux |
760 |
|
|
f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i) |
761 |
|
|
! all shortwave flux absorbed |
762 |
|
|
f_swpen=0. |
763 |
|
|
! bottom flux (ice conduction) |
764 |
|
|
slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki) |
765 |
|
|
! update ice temperature |
766 |
|
|
tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) & |
767 |
|
|
*(slab_bilg(ki)+f_cond)*dtime |
768 |
|
|
ELSE ! bare ice |
769 |
|
|
! ice-layer heat capacity |
770 |
|
|
cal(i)=2.*RCPD/(seaice(ki)*ice_cap) |
771 |
|
|
! conductive flux |
772 |
|
|
f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki) |
773 |
|
|
! penetrative shortwave flux... |
774 |
|
|
f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den) |
775 |
|
|
slab_bilg(ki)=f_swpen-f_cond |
776 |
|
|
END IF |
777 |
|
|
radsol(i)=radsol(i)+f_cond-f_swpen |
778 |
|
|
END DO |
779 |
|
|
! weight fluxes to ocean by sea ice fraction |
780 |
|
|
slab_bilg(:)=slab_bilg(:)*fsic(:) |
781 |
|
|
|
782 |
|
|
! calcul_fluxs (sens, lat etc) |
783 |
|
|
CALL calcul_fluxs(knon, is_sic, dtime, & |
784 |
|
|
tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, & |
785 |
|
|
precip_rain, precip_snow, snow, qsurf, & |
786 |
|
|
radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & |
787 |
|
|
f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & |
788 |
|
|
tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) |
789 |
|
|
DO i=1,knon |
790 |
|
|
IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i) |
791 |
|
|
END DO |
792 |
|
|
|
793 |
|
|
! calcul_flux_wind |
794 |
|
|
CALL calcul_flux_wind(knon, dtime, & |
795 |
|
|
u0, v0, u1, v1, gustiness, cdragm, & |
796 |
|
|
AcoefU, AcoefV, BcoefU, BcoefV, & |
797 |
|
|
p1lay, temp_air, & |
798 |
|
|
flux_u1, flux_v1) |
799 |
|
|
|
800 |
|
|
!**************************************************************************************** |
801 |
|
|
! 2) Update snow and ice surface |
802 |
|
|
!**************************************************************************************** |
803 |
|
|
! snow precip |
804 |
|
|
DO i=1,knon |
805 |
|
|
ki=knindex(i) |
806 |
|
|
IF (precip_snow(i) > 0.) THEN |
807 |
|
|
snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki))) |
808 |
|
|
END IF |
809 |
|
|
! snow and ice sublimation |
810 |
|
|
IF (evap(i) > 0.) THEN |
811 |
|
|
snow_evap = MIN (snow(i) / dtime, evap(i)) |
812 |
|
|
snow(i) = snow(i) - snow_evap * dtime |
813 |
|
|
snow(i) = MAX(0.0, snow(i)) |
814 |
|
|
seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime) |
815 |
|
|
ENDIF |
816 |
|
|
! Melt / Freeze snow from above if Tsurf>0 |
817 |
|
|
IF (tsurf_new(i).GT.t_melt) THEN |
818 |
|
|
! energy available for melting snow (in kg of melted snow /m2) |
819 |
|
|
e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. & |
820 |
|
|
/(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i)) |
821 |
|
|
! remove snow |
822 |
|
|
IF (snow(i).GT.e_melt) THEN |
823 |
|
|
snow(i)=snow(i)-e_melt |
824 |
|
|
tsurf_new(i)=t_melt |
825 |
|
|
ELSE ! all snow is melted |
826 |
|
|
! add remaining heat flux to ice |
827 |
|
|
e_melt=e_melt-snow(i) |
828 |
|
|
tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki)) |
829 |
|
|
tsurf_new(i)=tice(ki) |
830 |
|
|
END IF |
831 |
|
|
END IF |
832 |
|
|
! melt ice from above if Tice>0 |
833 |
|
|
IF (tice(ki).GT.t_melt) THEN |
834 |
|
|
! quantity of ice melted (kg/m2) |
835 |
|
|
e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. & |
836 |
|
|
/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0) |
837 |
|
|
! melt from above, height only |
838 |
|
|
dhsic=MIN(seaice(ki)-h_ice_min,e_melt) |
839 |
|
|
e_melt=e_melt-dhsic |
840 |
|
|
IF (e_melt.GT.0) THEN |
841 |
|
|
! lateral melt if ice too thin |
842 |
|
|
dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki)) |
843 |
|
|
! if all melted add remaining heat to ocean |
844 |
|
|
e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min) |
845 |
|
|
slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime |
846 |
|
|
! update height and fraction |
847 |
|
|
fsic(ki)=fsic(ki)-dfsic |
848 |
|
|
END IF |
849 |
|
|
seaice(ki)=seaice(ki)-dhsic |
850 |
|
|
! surface temperature at melting point |
851 |
|
|
tice(ki)=t_melt |
852 |
|
|
tsurf_new(i)=t_melt |
853 |
|
|
END IF |
854 |
|
|
! convert snow to ice if below floating line |
855 |
|
|
h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den |
856 |
|
|
IF (h_test.GT.0.) THEN !snow under water |
857 |
|
|
! extra snow converted to ice (with added frozen sea water) |
858 |
|
|
dhsic=h_test/(sea_den-ice_den+sno_den) |
859 |
|
|
seaice(ki)=seaice(ki)+dhsic |
860 |
|
|
snow(i)=snow(i)-dhsic*sno_den/ice_den |
861 |
|
|
! available energy (freeze sea water + bring to tice) |
862 |
|
|
e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ & |
863 |
|
|
ice_cap/2.*(t_freeze-tice(ki))) |
864 |
|
|
! update ice temperature |
865 |
|
|
tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki)) |
866 |
|
|
END IF |
867 |
|
|
END DO |
868 |
|
|
|
869 |
|
|
! New albedo |
870 |
|
|
DO i=1,knon |
871 |
|
|
ki=knindex(i) |
872 |
|
|
! snow albedo: update snow age |
873 |
|
|
IF (snow(i).GT.0.0001) THEN |
874 |
|
|
agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)& |
875 |
|
|
* EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.) |
876 |
|
|
ELSE |
877 |
|
|
agesno(i)=0.0 |
878 |
|
|
END IF |
879 |
|
|
! snow albedo |
880 |
|
|
alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.) |
881 |
|
|
! ice albedo (varies with ice tkickness and temp) |
882 |
|
|
alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) |
883 |
|
|
IF (tice(ki).GT.t_freeze-0.01) THEN |
884 |
|
|
alb_ice=MIN(alb_ice,alb_ice_wet) |
885 |
|
|
ELSE |
886 |
|
|
alb_ice=MIN(alb_ice,alb_ice_dry) |
887 |
|
|
END IF |
888 |
|
|
! pond albedo |
889 |
|
|
alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) |
890 |
|
|
! pond fraction |
891 |
|
|
frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) |
892 |
|
|
! snow fraction |
893 |
|
|
frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min)) |
894 |
|
|
! ice fraction |
895 |
|
|
frac_ice=MAX(0.0,1.-frac_pond-frac_snow) |
896 |
|
|
! total albedo |
897 |
|
|
alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond |
898 |
|
|
END DO |
899 |
|
|
alb2_new(:) = alb1_new(:) |
900 |
|
|
|
901 |
|
|
!**************************************************************************************** |
902 |
|
|
! 3) Recalculate new ocean temperature (add fluxes below ice) |
903 |
|
|
! Melt / freeze from below |
904 |
|
|
!***********************************************o***************************************** |
905 |
|
|
!cumul fluxes |
906 |
|
|
bilg_cum(:)=bilg_cum(:)+slab_bilg(:) |
907 |
|
|
IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction |
908 |
|
|
! Add cumulated surface fluxes |
909 |
|
|
tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime |
910 |
|
|
DO i=1,knon |
911 |
|
|
ki=knindex(i) |
912 |
|
|
! split lateral/top melt-freeze |
913 |
|
|
frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin))) |
914 |
|
|
IF (tslab(ki,1).LE.t_freeze) THEN |
915 |
|
|
! ****** Form new ice from below ******* |
916 |
|
|
! quantity of new ice |
917 |
|
|
e_melt=(t_freeze-tslab(ki,1))/cyang & |
918 |
|
|
/(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) |
919 |
|
|
! first increase height to h_thin |
920 |
|
|
dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki))) |
921 |
|
|
seaice(ki)=dhsic+seaice(ki) |
922 |
|
|
e_melt=e_melt-fsic(ki)*dhsic |
923 |
|
|
IF (e_melt.GT.0.) THEN |
924 |
|
|
! frac_mf fraction used for lateral increase |
925 |
|
|
dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki)) |
926 |
|
|
fsic(ki)=fsic(ki)+dfsic |
927 |
|
|
e_melt=e_melt-dfsic*seaice(ki) |
928 |
|
|
! rest used to increase height |
929 |
|
|
seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki)) |
930 |
|
|
END IF |
931 |
|
|
tslab(ki,1)=t_freeze |
932 |
|
|
ELSE ! slab temperature above freezing |
933 |
|
|
! ****** melt ice from below ******* |
934 |
|
|
! quantity of melted ice |
935 |
|
|
e_melt=(tslab(ki,1)-t_freeze)/cyang & |
936 |
|
|
/(ice_lat+ice_cap/2.*(tice(ki)-t_freeze)) |
937 |
|
|
! first decrease height to h_thick |
938 |
|
|
dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki))) |
939 |
|
|
seaice(ki)=seaice(ki)-dhsic |
940 |
|
|
e_melt=e_melt-fsic(ki)*dhsic |
941 |
|
|
IF (e_melt.GT.0) THEN |
942 |
|
|
! frac_mf fraction used for height decrease |
943 |
|
|
dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki))) |
944 |
|
|
seaice(ki)=seaice(ki)-dhsic |
945 |
|
|
e_melt=e_melt-fsic(ki)*dhsic |
946 |
|
|
! rest used to decrease fraction (up to 0!) |
947 |
|
|
dfsic=MIN(fsic(ki),e_melt/seaice(ki)) |
948 |
|
|
! keep remaining in ocean |
949 |
|
|
e_melt=e_melt-dfsic*seaice(ki) |
950 |
|
|
END IF |
951 |
|
|
tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang |
952 |
|
|
fsic(ki)=fsic(ki)-dfsic |
953 |
|
|
END IF |
954 |
|
|
END DO |
955 |
|
|
bilg_cum(:)=0. |
956 |
|
|
END IF ! coupling time |
957 |
|
|
|
958 |
|
|
!tests ice fraction |
959 |
|
|
WHERE (fsic.LT.ice_frac_min) |
960 |
|
|
tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang |
961 |
|
|
tice=t_melt |
962 |
|
|
fsic=0. |
963 |
|
|
seaice=0. |
964 |
|
|
END WHERE |
965 |
|
|
|
966 |
|
|
END SUBROUTINE ocean_slab_ice |
967 |
|
|
! |
968 |
|
|
!**************************************************************************************** |
969 |
|
|
! |
970 |
|
|
SUBROUTINE ocean_slab_final |
971 |
|
|
|
972 |
|
|
!**************************************************************************************** |
973 |
|
|
! Deallocate module variables |
974 |
|
|
!**************************************************************************************** |
975 |
|
|
IF (ALLOCATED(tslab)) DEALLOCATE(tslab) |
976 |
|
|
IF (ALLOCATED(fsic)) DEALLOCATE(fsic) |
977 |
|
|
IF (ALLOCATED(tice)) DEALLOCATE(tice) |
978 |
|
|
IF (ALLOCATED(seaice)) DEALLOCATE(seaice) |
979 |
|
|
IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg) |
980 |
|
|
IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum) |
981 |
|
|
IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum) |
982 |
|
|
IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum) |
983 |
|
|
IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum) |
984 |
|
|
IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman) |
985 |
|
|
IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff) |
986 |
|
|
IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm) |
987 |
|
|
IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux) |
988 |
|
|
|
989 |
|
|
END SUBROUTINE ocean_slab_final |
990 |
|
|
|
991 |
|
|
END MODULE ocean_slab_mod |