Context Navigation

ocean_slab_mod.F90 @ 1477

Last change on this file since 1477 was 1397, checked in by milmd, 10 years ago
In LMDZ.GENERIC replacement of all phystd .h files by module files.
Property svn:executable set to ``*
File size: 26.5 KB

Line
1	!
2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ocean_slab_mod.F90,v 1.3 2008-02-04 16:24:28 fairhead Exp $
3	!
4	MODULE ocean_slab_mod
5	!
6	! This module is used for both surface ocean and sea-ice when using the slab ocean,
7	! "ocean=slab".
8	!
9
10
11	use slab_ice_h
12	use watercommon_h, only: T_h2O_ice_liq
13	use surf_heat_transp_mod
14	implicit none
15
16
17
18
19
20
21	!LOGICAL, PRIVATE, SAVE :: ok_slab_sic,ok_slab_heaT_h2O_ice_liqBS
22	!!$OMP THREADPRIVATE(ok_slab_sic,ok_slab_heaT_h2O_ice_liqBS)
23	!INTEGER, PRIVATE, SAVE :: slab_ekman, slab_cadj
24	!!$OMP THREADPRIVATE(slab_ekman,slab_cadj)
25	INTEGER, PRIVATE, SAVE :: lmt_pas, julien, idayvrai
26	!$OMP THREADPRIVATE(lmt_pas,julien,idayvrai)
27	INTEGER, PRIVATE, SAVE :: cpl_pas
28	!$OMP THREADPRIVATE(cpl_pas)
29	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tmp_seaice
30	!$OMP THREADPRIVATE(tmp_seaice)
31	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tmp_tslab_loc
32	!$OMP THREADPRIVATE(tmp_tslab_loc)
33	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slab_bils
34	!$OMP THREADPRIVATE(slab_bils)
35	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE , SAVE :: lmt_bils
36	!$OMP THREADPRIVATE(lmt_bils)
37	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tmp_pctsrf_slab
38	!$OMP THREADPRIVATE(tmp_pctsrf_slab)
39	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tmp_tslab
40	!$OMP THREADPRIVATE(tmp_tslab)
41	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tmp_radsol
42	!$OMP THREADPRIVATE(tmp_radsol)
43	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: dt_hdiff
44	!$OMP THREADPRIVATE(dt_hdiff)
45	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: dt_ekman
46	!$OMP THREADPRIVATE(dt_ekman)
47	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tmp_flux_o, tmp_flux_g
48	!$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
49	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh
50	!$OMP THREADPRIVATE(slabh)
51	LOGICAL, PRIVATE, SAVE :: check = .FALSE.
52	!$OMP THREADPRIVATE(check)
53
54	CONTAINS
55	!
56	!****************************************************************************************
57	!
58	SUBROUTINE ocean_slab_init(ngrid,dtime, tslab_rst, seaice_rst, pctsrf_rst)
59
60	use slab_ice_h
61
62
63
64
65	! Input variables
66	!****************************************************************************************
67
68	integer,intent(in) :: ngrid ! number of atmospherci columns
69	REAL, INTENT(IN) :: dtime
70	! Variables read from restart file
71	REAL, DIMENSION(ngrid,noceanmx), INTENT(IN) :: tslab_rst
72	REAL, DIMENSION(ngrid), INTENT(IN) :: seaice_rst
73	REAL, DIMENSION(ngrid), INTENT(IN) :: pctsrf_rst
74
75
76	! Local variables
77	!****************************************************************************************
78	INTEGER :: error
79	CHARACTER (len = 80) :: abort_message
80	CHARACTER (len = 20) :: modname = 'ocean_slab_intit'
81
82
83	print,'***********************'
84	print*,'SLAB OCEAN est actif, prenez precautions !'
85	print,'***********************'
86
87	! Lecture des parametres:
88	! CALL getpar('ok_slab_sic',.true.,ok_slab_sic,'glace de mer dans slab')
89	! CALL getpar('slab_ekman',0,slab_ekman,'type transport Ekman pour slab (0=rien)')
90	! CALL getpar('slab_cadj',1,slab_cadj,'type ajustement conv slab 2 couches')
91
92	! Allocate variables initialize from restart fields
93	ALLOCATE(tmp_tslab(ngrid,noceanmx), stat = error)
94	IF (error /= 0) THEN
95	abort_message='Pb allocation tmp_tslab'
96	CALL abort_gcm(modname,abort_message,1)
97	ENDIF
98	tmp_tslab(:,:) = tslab_rst(:,:)
99	ALLOCATE(tmp_tslab_loc(ngrid,noceanmx), stat = error)
100	IF (error /= 0) THEN
101	abort_message='Pb allocation tmp_tslab_loc'
102	CALL abort_gcm(modname,abort_message,1)
103	ENDIF
104	tmp_tslab_loc(:,:) = tslab_rst(:,:)
105
106	ALLOCATE(tmp_seaice(ngrid), stat = error)
107	IF (error /= 0) THEN
108	abort_message='Pb allocation tmp_seaice'
109	CALL abort_gcm(modname,abort_message,1)
110	ENDIF
111	tmp_seaice(:) = seaice_rst(:)
112
113	ALLOCATE(tmp_pctsrf_slab(ngrid), stat = error)
114	IF (error /= 0) THEN
115	abort_message='Pb allocation tmp_pctsrf_slab'
116	CALL abort_gcm(modname,abort_message,1)
117	ENDIF
118	tmp_pctsrf_slab(:) = pctsrf_rst(:)
119
120	! Allocate some other variables internal in module mod_oceanslab
121	ALLOCATE(tmp_radsol(ngrid), stat = error)
122	IF (error /= 0) THEN
123	abort_message='Pb allocation tmp_radsol'
124	CALL abort_gcm(modname,abort_message,1)
125	ENDIF
126
127	ALLOCATE(tmp_flux_o(ngrid), stat = error)
128	IF (error /= 0) THEN
129	abort_message='Pb allocation tmp_flux_o'
130	CALL abort_gcm(modname,abort_message,1)
131	ENDIF
132
133	ALLOCATE(tmp_flux_g(ngrid), stat = error)
134	IF (error /= 0) THEN
135	abort_message='Pb allocation tmp_flux_g'
136	CALL abort_gcm(modname,abort_message,1)
137	ENDIF
138
139	! a mettre un slab_bils aussi en force !!!
140	ALLOCATE(slab_bils(ngrid), stat = error)
141	IF (error /= 0) THEN
142	abort_message='Pb allocation slab_bils'
143	CALL abort_gcm(modname,abort_message,1)
144	ENDIF
145	slab_bils(:) = 0.0
146
147	ALLOCATE(dt_hdiff(ngrid,noceanmx), stat = error)
148	IF (error /= 0) THEN
149	abort_message='Pb allocation dt_hdiff'
150	CALL abort_gcm(modname,abort_message,1)
151	ENDIF
152	dt_hdiff = 0.0
153
154	ALLOCATE(dt_ekman(ngrid,noceanmx), stat = error)
155	IF (error /= 0) THEN
156	abort_message='Pb allocation dt_hdiff'
157	CALL abort_gcm(modname,abort_message,1)
158	ENDIF
159	dt_ekman = 0.0
160
161
162	ALLOCATE(lmt_bils(ngrid), stat = error)
163	IF (error /= 0) THEN
164	abort_message='Pb allocation lmt_bils'
165	CALL abort_gcm(modname,abort_message,1)
166	ENDIF
167	lmt_bils(:) = 0.0
168
169	ALLOCATE(slabh(noceanmx), stat = error)
170	IF (error /= 0) THEN
171	abort_message='Pb allocation slabh'
172	CALL abort_gcm(modname,abort_message,1)
173	ENDIF
174	slabh(1)=50.
175	IF (noceanmx.GE.2) slabh(2)=150.
176	IF (noceanmx.GE.3) slabh(3)=2800.
177
178
179	! CALL init_masquv
180
181
182
183
184	END SUBROUTINE ocean_slab_init
185	!
186	!****************************************************************************************
187	!
188
189	!
190	!****************************************************************************************
191	!
192	SUBROUTINE ocean_slab_ice(dtime, &
193	ngrid, knindex, tsurf, radsol, &
194	precip_snow, snow, evap, &
195	fluxsens,tsurf_new, pctsrf_sic, &
196	taux_slab,tauy_slab,icount)
197
198	use slab_ice_h
199	use callkeys_mod, only: ok_slab_sic
200
201
202	! Input arguments
203	!****************************************************************************************
204	INTEGER, INTENT(IN) :: ngrid
205	INTEGER, DIMENSION(ngrid), INTENT(IN) :: knindex
206	REAL, INTENT(IN) :: dtime
207	REAL, DIMENSION(ngrid), INTENT(IN) :: tsurf
208	REAL, DIMENSION(ngrid), INTENT(IN) :: taux_slab
209	REAL, DIMENSION(ngrid), INTENT(IN) :: tauy_slab
210	REAL, DIMENSION(ngrid), INTENT(IN) :: evap, fluxsens
211	REAL, DIMENSION(ngrid), INTENT(IN) :: precip_snow
212	REAL, DIMENSION(ngrid), INTENT(IN) :: radsol
213	INTEGER, INTENT(IN) :: icount
214
215
216	!In/Output arguments
217	!****************************************************************************************l
218	REAL, DIMENSION(ngrid), INTENT(INOUT) :: snow
219
220	! REAL, DIMENSION(ngridmx), INTENT(INOUT) :: agesno
221
222
223	! Output arguments
224	!****************************************************************************************
225	! REAL, DIMENSION(ngridmx), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
226	! REAL, DIMENSION(ngridmx), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
227	REAL, DIMENSION(ngrid), INTENT(OUT) :: tsurf_new
228	REAL, DIMENSION(ngrid), INTENT(OUT) :: pctsrf_sic
229
230	! Local variables
231	!****************************************************************************************
232	INTEGER :: i
233	REAL, DIMENSION(ngrid) :: cal, beta, dif_grnd, capsol
234	REAL, DIMENSION(ngrid) :: alb_neig, alb_ice!, tsurf_temp
235	REAL, DIMENSION(ngrid) :: soilcap, soilflux
236	REAL, DIMENSION(ngrid) :: zfra
237	REAL :: snow_evap, fonte
238	REAL, DIMENSION(ngrid,noceanmx) :: tslab
239	REAL, DIMENSION(ngrid) :: seaice,tsea_ice ! glace de mer (kg/m2)
240	REAL, DIMENSION(ngrid) :: pctsrf_new
241
242
243	!*****************************************************************************
244
245
246	! Initialization of output variables
247	! alb1_new(:) = 0.0
248
249	!******************************************************************************
250	! F. Codron: 3 cas
251	! -Glace interactive, quantité seaice : T glace suit modèle simple
252	! -Pas de glace: T oce peut descendre en dessous de 0°C sans geler
253	!*****************************************************************************
254
255	pctsrf_new(:)=tmp_pctsrf_slab(:)
256	tmp_radsol(:)=radsol(:)
257	tmp_flux_o(:)=fluxsens(:)
258
259	DO i = 1, ngrid
260	tsurf_new(i) = tsurf(knindex(i))
261	seaice(i)=tmp_seaice(knindex(i))
262
263
264
265	IF (pctsrf_new(knindex(i)) < EPSFRA) THEN
266	snow(i) = 0.0
267	tsurf_new(i) = T_h2O_ice_liq - 1.8
268	!IF (soil_model) tsoil(i,:) = T_h2O_ice_liq -1.8
269	ENDIF
270	ENDDO
271	tmp_flux_g(:) = 0.0
272	tsea_ice(:)=T_h2O_ice_liq-1.8
273	! Calculs flux glace/mer et glace/air
274	IF (ok_slab_sic) THEN
275	!*********************************
276	! Calcul de beta, heat capacity
277	!*********************************
278	! CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
279
280	! IF ((soil_model)) THEN
281
282	! ELSE
283	! dif_grnd = 1.0 / tau_gl
284	! cal = RCPD * calice
285	! WHERE (snow > 0.0) cal = RCPD * calsno
286	! ENDIF
287	! tsurf_temp = tsurf_new
288	! beta = 1.0
289
290
291	! **********************************************
292	! Evolution avec glace interactive:
293	! *************************************
294	! tsurf_new=tsurf_temp
295	DO i = 1, ngrid
296	IF (pctsrf_new(knindex(i)) .GT. epsfra) THEN
297	! *************************************
298	! + Calcul Flux entre glace et océan
299	! *************************************
300	tmp_flux_g(knindex(i))=(tsurf_new(i)-(T_h2O_ice_liq-1.8)) &
301	* ice_cond*ice_den/seaice(i)
302
303	! ****************************************
304	! + Calcul nouvelle température de la glace
305	! ****************************************
306	tsurf_new(i)=tsurf_new(i)+2*(radsol(i)+fluxsens(i) &
307	-tmp_flux_g(knindex(i))) &
308	/(ice_capseaice(i))dtime
309
310	! ***************************************
311	! + Precip and evaporation of snow and ice
312	! ***************************************
313	! Add precip
314	IF (precip_snow(i).GT.0.) THEN
315	snow(i)=snow(i)+precip_snow(i)*dtime
316	ENDIF
317	! Evaporation
318	snow_evap=0.
319	IF (evap(i).GT.0.) THEN
320	snow_evap = MIN (snow(i) / dtime, evap(i))
321	snow(i) = snow(i) - snow_evap*dtime
322	snow(i) = MAX(0.0, snow(i))
323	seaice(i)=MAX(0.0,seaice(i)-(evap(i)-snow_evap)*dtime)
324	ENDIF
325
326	! *****************************************
327	! + Fonte neige & glace par le haut
328	! *****************************************
329	! Snow melt
330	IF ((snow(i) > epsfra) .AND. (tsurf_new(i) >= T_h2O_ice_liq)) THEN
331	fonte=MIN(MAX((tsurf_new(i)-T_h2O_ice_liq)*seaice(i) &
332	/2.*ice_cap/ice_lat,0.0),snow(i))
333	snow(i) = MAX(0.,(snow(i)-fonte))
334	tsurf_new(i)=tsurf_new(i)-fonte2ice_lat/(seaice(i)*ice_cap)
335	ENDIF
336	snow(i)=MIN(snow(i),3000.)
337	! Ice melt
338	IF ((seaice(i) > epsfra) .AND. (tsurf_new(i) >= T_h2O_ice_liq)) THEN
339	fonte=seaice(i)ice_cap(tsurf_new(i)-T_h2O_ice_liq) &
340	/(2ice_lat+ice_cap1.8)
341	CALL icemelt(fonte,pctsrf_new(knindex(i)),seaice(i))
342	tmp_flux_g(knindex(i))=tmp_flux_g(knindex(i)) &
343	+fonte*ice_lat/dtime
344	tsurf_new(i)=T_h2O_ice_liq
345	ENDIF
346	tmp_seaice(knindex(i))=seaice(i)
347	ENDIF!(pctsrf_new(knindex(i)) .GT. epsfra)
348	tsea_ice(knindex(i))=tsurf_new(i)
349
350	ENDDO
351	ENDIF
352	! ******************************************
353	! CALL interfoce:
354	! cumul/calcul nouvelle T_h2O_ice_liqce
355	! fonte/formation de glace en dessous
356	! ******************************************
357	tslab = tmp_tslab
358
359	CALL interfoce_slab(ngrid, dtime, &
360	tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
361	tslab, tsea_ice, pctsrf_new,taux_slab,tauy_slab,icount)
362
363	tmp_pctsrf_slab(:)=pctsrf_new(:)
364	! tmp_pctsrf_slab(:,is_oce)=1.0-tmp_pctsrf_slab(:) &
365	! -tmp_pctsrf_slab(:,is_lic)-tmp_pctsrf_slab(:,is_ter)
366
367	DO i=1, ngrid
368	tmp_tslab(knindex(i),:)=tslab(knindex(i),:)
369	seaice(i)=tmp_seaice(knindex(i))
370	tsurf_new(i)=tsea_ice(knindex(i))
371
372	ENDDO
373
374	! ****************************
375	! calcul new albedo
376	! ****************************
377	! Albedo neige
378	! CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
379	! WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
380	! Fraction neige (hauteur critique 45kg/m2~15cm)
381	! zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/45.0))
382	! Albedo glace
383	! IF (ok_slab_sicOBS) THEN
384	! alb_ice=0.6
385	! ELSE
386	! alb_ice(1:knon)=alb_ice_max-(alb_ice_max-alb_ice_min) &
387	! *exp(-seaice(1:knon)/h_alb_ice)
388	! ENDIF
389	!Albedo final
390	! alb1_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
391	! alb_ice * (1.0-zfra(1:knon))
392	! alb2_new(:) = alb1_new(:)
393
394
395	! ********************************
396	! Return the fraction of sea-ice
397	! ********************************
398	pctsrf_sic(:) = tmp_pctsrf_slab(:)
399
400
401	END SUBROUTINE ocean_slab_ice
402
403
404	!
405	!***************************************************************************
406	!
407	SUBROUTINE interfoce_slab(ngrid, dtime, &
408	radsol, fluxo, fluxg, pctsrf, &
409	tslab, tsea_ice, pctsrf_slab, &
410	taux_slab, tauy_slab,icount)
411	!
412	! Cette routine calcule la temperature d'un slab ocean, la glace de mer
413	! et les pourcentages de la maille couverte par l'ocean libre et/ou
414	! la glace de mer pour un "slab" ocean de 50m
415	!
416	! Conception: Laurent Li
417	! Re-ecriture + adaptation LMDZ4: I. Musat
418	! Transport, nouveau modèle glace, 2 couches: F.Codron
419	!
420	! input:
421	! klon nombre T_h2O_ice_liqtal de points de grille
422	! itap numero du pas de temps
423	! dtime pas de temps de la physique (en s)
424	! ijour jour dans l'annee en cours
425	! radsol rayonnement net au sol (LW + SW)
426	! fluxo flux turbulent (sensible + latent) sur les mailles oceaniques
427	! fluxg flux de conduction entre la surface de la glace de mer et l'ocean
428	! pctsrf tableau des pourcentages de surface de chaque maille
429	! output:
430	! tslab temperature de l'ocean libre
431	! tsea_ice temperature de la glace (surface)
432	! pctsrf_slab "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
433
434	use slab_ice_h
435	use callkeys_mod, only: ok_slab_sic,ok_slab_heat_transp
436
437	! Input arguments
438	!****************************************************************************************
439	INTEGER, INTENT(IN) :: ngrid,icount
440	! INTEGER, INTENT(IN) :: itap
441	REAL, INTENT(IN) :: dtime ! not used
442	! INTEGER, INTENT(IN) :: ijour
443	REAL, DIMENSION(ngrid), INTENT(IN) :: radsol
444	REAL, DIMENSION(ngrid), INTENT(IN) :: fluxo
445	REAL, DIMENSION(ngrid), INTENT(IN) :: fluxg
446	REAL, DIMENSION(ngrid), INTENT(IN) :: pctsrf
447	REAL, DIMENSION(ngrid), INTENT(IN) :: taux_slab
448	REAL, DIMENSION(ngrid), INTENT(IN) :: tauy_slab
449
450	! Output arguments
451	!****************************************************************************************
452	REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: tslab
453	REAL, DIMENSION(ngrid), INTENT(INOUT) :: pctsrf_slab
454	REAL, DIMENSION(ngrid), INTENT(INOUT) :: tsea_ice
455
456	! Local variables
457	!****************************************************************************************
458	REAL :: fonte,t_cadj
459	INTEGER :: i,k,kt,kb
460	REAL :: zz, za, zb
461	REAL :: cyang ! capacite calorifique slab J/(m2 K)
462	REAL, PARAMETER :: tau_conv=432000. ! temps ajust conv (5 jours)
463	REAL, DIMENSION(ngrid,noceanmx) :: dtdiff_loc,dtekman_loc
464
465
466
467	!***************************************************
468	! Capacite calorifique de la couche de surface
469	cyang=capcalocean!slabh(1)*4.228e+06
470	cpl_pas=1!4*iradia
471
472	!
473	! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
474	!! julien = MOD(ijour,360)
475	!! IF (ok_slab_heaT_h2O_ice_liqBS .and. (MOD(itap,lmt_pas) .EQ. 1)) THEN
476	!! ! 1er pas de temps de la journee
477	!! idayvrai = ijour
478	!! CALL condsurf(julien,idayvrai, lmt_bils)
479	!! ENDIF
480
481	! ----------------------------------------------
482	! A chaque pas de temps: cumul flux de chaleur
483	! ----------------------------------------------
484	! bilan du flux de chaleur dans l'ocean :
485	!
486	DO i = 1, ngrid
487	za = radsol(i) + fluxo(i)
488	zb = fluxg(i)
489	! IF(((1-pctsrf(i)).GT.epsfra).OR. &
490	! (pctsrf(i).GT.epsfra)) THEN
491	slab_bils(i)=slab_bils(i)+(za*(1-pctsrf(i)) &
492	+zb*pctsrf(i))/ cpl_pas
493	! ENDIF
494
495	ENDDO !klon
496
497	! ---------------------------------------------
498	! T_h2O_ice_liqus les cpl_pas: update de tslab et seaice
499	! ---------------------------------------------
500
501	IF (mod(icount-1,cpl_pas).eq.0) THEN !fin de journee
502	!
503	! ---------------------------------------------
504	! Transport de chaleur par circulation
505	! Decoupage par pas de temps pour stabilite numérique
506	! (diffusion, schema centre pour advection)
507	! ---------------------------------------------
508	dt_hdiff(:,:)=0.
509	dt_ekman(:,:)=0.
510
511	IF (ok_slab_heat_transp) THEN
512	! DO i=1,cpl_pas
513	! Transport diffusif
514	! IF (ok_soil_hdiff) THEN
515	CALL divgrad_phy(ngrid,noceanmx,tmp_tslab_loc,dtdiff_loc)
516	dtdiff_loc=dtdiff_locsoil_hdiff50./SUM(slabh)!*100.
517	! dtdiff_loc(:,1)=dtdiff_loc(:,1)soil_hdiff50./SUM(slabh)*0.8
518	! dtdiff_loc(:,2)=dtdiff_loc(:,2)soil_hdiff50./SUM(slabh)*0.2
519	dt_hdiff=dt_hdiff+dtdiff_loc
520	tmp_tslab_loc=tmp_tslab_loc+dtdiff_loc*dtime
521	! END IF
522
523	! Calcul de transport par Ekman
524
525	CALL slab_ekman2(ngrid,taux_slab,tauy_slab,tslab,dtekman_loc)
526
527
528	! END SELECT
529	! IF (slab_ekman.GT.0) THEN
530	do k=1,noceanmx
531	dtekman_loc(:,k)=dtekman_loc(:,k)/(slabh(k)1000.)!0.
532	enddo
533	dt_ekman(:,:)=dt_ekman(:,:)+dtekman_loc(:,:)
534	tmp_tslab_loc=tmp_tslab_loc+dtekman_loc*dtime
535	! ENDIF
536	! ENDDO ! time splitting 1,cpl_pas
537	! IF (slab_ekman.GT.0) THEN
538	! taux_slab(:)=0.
539	! tauy_slab(:)=0.
540	! ENDIF
541
542	! print*, 'slab_bils=',slab_bils(1)
543
544
545	ENDIF!(ok_slab_heat_transp)
546
547	DO i = 1, ngrid
548	! IF (((1-pctsrf(i)).GT.epsfra).OR. &
549	! (pctsrf(i).GT.epsfra)) THEN
550	! ---------------------------------------------
551	! Ajout des flux de chaleur dans l'océan
552	! ---------------------------------------
553
554	!print*, 'T_h2O_ice_liqcean1=',tmp_tslab_loc(i,1)
555	tmp_tslab_loc(i,1) = tmp_tslab_loc(i,1) + &
556	slab_bils(i)/cyangdtimecpl_pas
557
558	!print*, 'capcalocean=',capcalocean
559	!print*, 'cyang=',cyang
560	!print, 'dT_h2O_ice_liqcean=',slab_bils(i)/cyangdtime
561	!print*, 'T_h2O_ice_liqcean2=',tmp_tslab_loc(i,1)
562
563
564	! on remet l'accumulation a 0
565	slab_bils(i) = 0.
566	! ---------------------------------------------
567	! Glace interactive
568	! ---------------------------------------------
569	IF (ok_slab_sic) THEN
570	! Fondre la glace si température > 0.
571	! -----------------------------------
572	IF ((tmp_tslab_loc(i,1).GT.T_h2O_ice_liq-1.8) .AND. (tmp_seaice(i).GT.0.0)) THEN
573	fonte=(tmp_tslab_loc(i,1)-T_h2O_ice_liq+1.8)*cyang &
574	/(ice_lat+ice_cap/2*(T_h2O_ice_liq-1.8-tsea_ice(i)))
575	CALL icemelt(fonte,pctsrf_slab(i),tmp_seaice(i))
576	tmp_tslab_loc(i,1)=T_h2O_ice_liq-1.8+fonte*ice_lat/cyang
577	ENDIF
578	! fabriquer de la glace si congelation atteinte:
579	! ----------------------------------------------
580	IF (tmp_tslab_loc(i,1).LT.(T_h2O_ice_liq-1.8)) THEN
581
582	IF (tmp_seaice(i).LT.h_ice_min) THEN
583	! Creation glace nouvelle
584	! IF (pctsrf_slab(i).LT.(epsfra)) THEN
585	fonte=(T_h2O_ice_liq-1.8-tmp_tslab_loc(i,1))*cyang/ice_lat
586	IF (fonte.GT.h_ice_min*ice_frac_min) THEN
587	tmp_seaice(i)=MIN(h_ice_thin,fonte/ice_frac_min)
588	tmp_seaice(i)=MAX(tmp_seaice(i),fonte/ice_frac_max)
589	! IF (fonte.GT.0) THEN
590	! tmp_seaice(i)=fonte
591	tsea_ice(i)=T_h2O_ice_liq-1.8
592	pctsrf_slab(i)=(1-pctsrf_slab(i))*fonte/tmp_seaice(i)
593	! pctsrf_slab(i)=1.
594	tmp_tslab_loc(i,1)=T_h2O_ice_liq-1.8
595	ENDIF
596	ELSE
597	! glace déjà présente
598	! Augmenter epaisseur
599	fonte=(T_h2O_ice_liq-1.8-tmp_tslab_loc(i,1))*cyang &
600	/(ice_lat+ice_cap/2*(T_h2O_ice_liq-1.8-tsea_ice(i)))
601	zz=tmp_seaice(i)
602	tmp_seaice(i)=MAX(zz,MIN(h_ice_thick,fonte+zz))
603	! Augmenter couverture (oce libre et h>h_thick)
604	za=1-pctsrf_slab(i)
605	zb=pctsrf_slab(i)
606	fonte=fonteza+MAX(0.0,fonte+zz-tmp_seaice(i))zb
607	pctsrf_slab(i)=MIN(zb+fonte/tmp_seaice(i), &
608	(za+zb)*ice_frac_max)
609	fonte=MAX(0.0,fonte-(pctsrf_slab(i)-zb)*tmp_seaice(i))
610	! Augmenter epaisseur si couverture complete
611	tmp_seaice(i)=tmp_seaice(i)+fonte/pctsrf_slab(i)
612	tmp_tslab_loc(i,1) = T_h2O_ice_liq-1.8
613	ENDIF ! presence glace
614	ENDIF !congelation
615	! vérifier limites de hauteur de glace
616	IF(tmp_seaice(i).GT.h_ice_min) THEN
617	tmp_seaice(i) = MIN(tmp_seaice(i),h_ice_max)
618	ELSE
619	tmp_seaice(i) = 0.
620	pctsrf_slab(i)=0.
621	ENDIF! (tmp_seaice(i).GT.h_ice_min)
622
623	ENDIF !(ok_slab_sic) Glace interactive
624	! ----------------------------------
625	! Ajustement convectif si > 1 layers
626	! ----------------------------------
627
628	IF ((noceanmx.GT.1)) THEN
629	DO kt=1,noceanmx-1
630	kb=kt
631	DO k=kt+1,noceanmx !look for instability
632	IF (tmp_tslab_loc(i,k).GT.tmp_tslab_loc(i,kt)) THEN
633	kb=k
634	ENDIF
635	ENDDO
636	IF (kb.GT.kt) THEN !ajust conv
637	! IF (slab_cadj.EQ.1) THEN
638	! : t_cadj=SUM(tmp_tslab_loc(i,kt:kb)*slabh(kt:kb))/SUM(slabh(kt:kb))
639	! DO k=kt,kb
640	! tmp_tslab_loc(i,k)=t_cadj
641	! ENDDO
642	! ELSEIF (slab_cadj.EQ.2) THEN
643	t_cadj=(tmp_tslab_loc(i,kb)-tmp_tslab_loc(i,kt))dtime/tau_convcpl_pas
644	tmp_tslab_loc(i,kt)=tmp_tslab_loc(i,kt)+t_cadj
645	tmp_tslab_loc(i,kb)=tmp_tslab_loc(i,kb)-t_cadj*slabh(kt)/slabh(kb)
646	!ENDIF
647	ENDIF
648	ENDDO
649	ENDIF !ajust 2 layers
650	! ENDIF !pctsrf
651	ENDDO !klon
652
653	! On met a jour la temperature et la glace
654	tslab = tmp_tslab_loc
655
656
657	ENDIF !(mod(icount-1,cpl_pas).eq.0)
658
659	END SUBROUTINE interfoce_slab
660	!
661	!******************************************************************************
662	SUBROUTINE icemelt(fonte,pctsrf,seaice)
663
664	use slab_ice_h
665
666
667
668	REAL, INTENT(INOUT) :: pctsrf
669	REAL , INTENT(INOUT) ::fonte,seaice !kg/m2
670	REAL :: hh !auxilliary
671	REAL :: ff !auxilliary
672
673
674	! ice gt h_ice_thick: decrease thickness up T_h2O_ice_liq h_ice_thick
675	IF (seaice.GT.h_ice_thick) THEN
676	hh=seaice
677	! ff=fonte
678	seaice=max(h_ice_thick,hh-fonte)
679	fonte=max(0.0,fonte+h_ice_thick-hh)
680
681	! seaice=max(0.,hh-fonte)
682	! fonte=max(0.0,fonte-(seaice-hh))
683	! IF (seaice.LT.epsfra) THEN
684	! pctsrf=0.
685	! seaice=0.
686	! fonte=ff-hh
687	! ENDIF
688
689	ENDIF
690	! ice gt h_ice_thin: partially decrease thickness
691	IF ((seaice.GE.h_ice_thin).AND.(fonte.GT.0.0)) THEN
692	hh=seaice
693	seaice=MAX(hh-0.6*fonte,h_ice_thin)
694	fonte=MAX(0.0,fonte-hh+seaice)
695	ENDIF
696	! use rest T_h2O_ice_liq decrease area
697	hh=pctsrf
698	pctsrf=MIN(hh,MAX(0.0,hh-fonte/seaice))
699	fonte=MAX(0.0,fonte-(hh-pctsrf)*seaice)
700
701	END SUBROUTINE icemelt
702
703	!****************************************************************************************
704	!
705	SUBROUTINE ocean_slab_final!(tslab_rst, seaice_rst)
706
707	! This subroutine will send T_h2O_ice_liq phyredem the variables concerning the slab
708	! ocean that should be written T_h2O_ice_liq restart file.
709
710	!****************************************************************************************
711
712	! REAL, DIMENSION(ngridmx,noceanmx), INTENT(OUT) :: tslab_rst
713	! REAL, DIMENSION(ngridmx), INTENT(OUT) :: seaice_rst
714
715	!****************************************************************************************
716	! Set the output variables
717	! tslab_rst(:,:) = tmp_tslab(:,:)
718	! tslab_rst(:) = tmp_tslab_loc(:)
719	! seaice_rst(:) = tmp_seaice(:)
720
721	! Deallocation of all variables in module
722	DEALLOCATE(tmp_tslab, tmp_tslab_loc, tmp_pctsrf_slab)
723	DEALLOCATE(tmp_seaice, tmp_radsol, tmp_flux_o, tmp_flux_g)
724	DEALLOCATE(slab_bils, lmt_bils)
725	DEALLOCATE(dt_hdiff)
726
727	END SUBROUTINE ocean_slab_final
728	!
729	!****************************************************************************************
730	!
731	SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_o_loc, flux_g_loc, &
732	dt_hdiff_loc,dt_ekman_loc)
733
734	! "Get some variables from module ocean_slab_mod"
735	! This subroutine prints variables T_h2O_ice_liq a external routine
736
737	INTEGER, INTENT(IN) :: ngrid
738	REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: tslab_loc
739	REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc
740	REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_o_loc
741	REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc
742	REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: dt_hdiff_loc
743	REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: dt_ekman_loc
744	INTEGER :: i
745
746
747	! Set the output variables
748	tslab_loc=0.
749	dt_hdiff_loc=0.
750	dt_ekman_loc=0.
751	tslab_loc = tmp_tslab
752	seaice_loc(:) = tmp_seaice(:)
753	flux_o_loc(:) = tmp_flux_o(:)
754	flux_g_loc(:) = tmp_flux_g(:)
755	DO i=1,noceanmx
756	dt_hdiff_loc(:,i) = dt_hdiff(:,i)slabh(i)1000.*4228. !Convert en W/m2
757	dt_ekman_loc(:,i) = dt_ekman(:,i)slabh(i)1000.*4228.
758	ENDDO
759
760
761
762	END SUBROUTINE ocean_slab_get_vars
763	!
764	!****************************************************************************************
765	!
766	END MODULE ocean_slab_mod

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

Download in other formats:

Original Format