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>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<2) THEN |
---|
189 | slab_gm=.FALSE. |
---|
190 | ENDIF |
---|
191 | ! Convective adjustment |
---|
192 | IF (nslay==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>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>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>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>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)==0) THEN ! time to update tslab & fraction |
---|
490 | ! *********************************** |
---|
491 | ! Horizontal transport |
---|
492 | ! *********************************** |
---|
493 | IF (slab_ekman>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>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>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>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>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>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>1).AND.(slab_cadj>0)) THEN |
---|
600 | DO i=1,klon |
---|
601 | IF (tslab(i,2)>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)<t_freeze).OR.(fsic(ki)>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)<epsfra) THEN ! Free of ice |
---|
631 | IF (tslab(ki,1)<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)>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)<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)>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)<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)>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)>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)>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>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>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)>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)>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)==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)<=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>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>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<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 |
---|