Context Navigation

ocean_slab_mod.F90 @ 4257

Last change on this file since 4257 was 2656, checked in by Ehouarn Millour, 8 years ago

Making the slab work:

added a slab_heat_transp_mod module for horizontal diffusion and Ekman transport
added storage and output of relevent variables in phyredem, phyetat0, phy_output_ctrlout_mod, phys_output_write_mod
moved nslay (number of slab layers) out of dimphy into ocean_slab_mod.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 37.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format