Context Navigation

ocean_slab_mod.F90 @ 2613

Last change on this file since 2613 was 1682, checked in by emillour, 8 years ago

All GCMs: set things up to enable pluging physics with dynamico

dyn3d
gcm.F90 : move I/O initialization (dates) to be done before physics

initialization

dyn3dpar
gcm.F : move I/O initialization (dates) to be done before physics

initialization

dynphy_lonlat:
inigeomphy_mod.F90 : add ind_cell_glo computation and transfer

to init_geometry

phy_common:
geometry_mod.F90 : add ind_cell_glo module variable to store global

column index

print_control_mod.F90 : make initialization occur via init_print_control_mod

to avoid circular module dependencies

init_print_control_mod.F90 : added to initialize print_control_mod module

variables

mod_phys_lmdz_mpi_data.F90 : use print_control_mod (rather than iniprint.h)
mod_phys_lmdz_para.F90 : use print_control_mod (rather than iniprint.h)
mod_phys_lmdz_omp_data.F90 : add is_omp_master (alias of is_omp_root) module

variable and use print_control_mod (rather than
iniprint.h)

physics_distribution_mod.F90 : add call to init_dimphy in

init_physics_distribution

xios_writefield.F90 : generic routine to output field with XIOS (for debug)

misc:
handle_err_m.F90 : call abort_physic, rather than abort_gcm
wxios.F90 : updates to enable unstructured grids

set module variable g_ctx_name to "LMDZ"
wxios_init(): remove call to wxios_context_init
wxios_context_init(): call xios_context_initialize with COMM_LMDZ_PHY
add routine wxios_set_context() to get handle and set context to XIOS
wxios_domain_param(): change arguments and generate the domain in-place
add wxios_domain_param_unstructured(): generate domain for unstructured case

NB: access is via "domain group" (whereas it is via "domain" in

wxios_domain_param)

dynphy_lonlat/phy[std|mars|venus|titan]:
iniphysiq_mod.F90 : Remove call to init_dimphy (which is now done in

phy_common/physics_distribution_mod.F90)

Property svn:executable set to *

File size: 26.5 KB

Rev	Line
[1298]	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
[1315]	22	!!$OMP THREADPRIVATE(ok_slab_sic,ok_slab_heaT_h2O_ice_liqBS)
[1298]	23	!INTEGER, PRIVATE, SAVE :: slab_ekman, slab_cadj
[1315]	24	!!$OMP THREADPRIVATE(slab_ekman,slab_cadj)
[1298]	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'
[1682]	96	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	102	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	109	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	116	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	124	CALL abort_physic(modname,abort_message,1)
[1298]	125	ENDIF
	126
	127	ALLOCATE(tmp_flux_o(ngrid), stat = error)
	128	IF (error /= 0) THEN
	129	abort_message='Pb allocation tmp_flux_o'
[1682]	130	CALL abort_physic(modname,abort_message,1)
[1298]	131	ENDIF
	132
	133	ALLOCATE(tmp_flux_g(ngrid), stat = error)
	134	IF (error /= 0) THEN
	135	abort_message='Pb allocation tmp_flux_g'
[1682]	136	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	143	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	150	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	157	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	165	CALL abort_physic(modname,abort_message,1)
[1298]	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'
[1682]	172	CALL abort_physic(modname,abort_message,1)
[1298]	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
[1397]	199	use callkeys_mod, only: ok_slab_sic
[1298]	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
[1397]	435	use callkeys_mod, only: ok_slab_sic,ok_slab_heat_transp
[1298]	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
[1308]	515	CALL divgrad_phy(ngrid,noceanmx,tmp_tslab_loc,dtdiff_loc)
[1298]	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
[1308]	525	CALL slab_ekman2(ngrid,taux_slab,tauy_slab,tslab,dtekman_loc)
[1298]	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