Context Navigation

ocean_slab_mod.F90 @ 1308

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

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

Download in other formats:

Original Format