Context Navigation

ocean_slab_mod.F90 @ 1351

Last change on this file since 1351 was 1315, checked in by milmd, 10 years ago
LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.
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