Context Navigation

ocean_slab_mod.F90 @ 2569

Last change on this file since 2569 was 2408, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2298:2396 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 28.7 KB

Rev	Line
[781]	1	!
	2	MODULE ocean_slab_mod
	3	!
	4	! This module is used for both surface ocean and sea-ice when using the slab ocean,
	5	! "ocean=slab".
	6	!
[2073]	7
	8	USE dimphy
	9	USE indice_sol_mod
[2220]	10	USE surface_data
[2073]	11
[781]	12	IMPLICIT NONE
[996]	13	PRIVATE
[2220]	14	PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
[781]	15
[2220]	16	!****************************************************************************************
	17	! Global saved variables
	18	!****************************************************************************************
	19
[2073]	20	INTEGER, PRIVATE, SAVE :: cpl_pas
	21	!$OMP THREADPRIVATE(cpl_pas)
	22	REAL, PRIVATE, SAVE :: cyang
	23	!$OMP THREADPRIVATE(cyang)
	24	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh
	25	!$OMP THREADPRIVATE(slabh)
	26	REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab
	27	!$OMP THREADPRIVATE(tslab)
[2220]	28	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic
	29	!$OMP THREADPRIVATE(fsic)
	30	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice
	31	!$OMP THREADPRIVATE(tice)
	32	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice
	33	!$OMP THREADPRIVATE(seaice)
[2073]	34	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bils
	35	!$OMP THREADPRIVATE(slab_bils)
	36	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum
	37	!$OMP THREADPRIVATE(bils_cum)
[2220]	38	REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg
	39	!$OMP THREADPRIVATE(slab_bilg)
	40	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum
	41	!$OMP THREADPRIVATE(bilg_cum)
[2073]	42
[2220]	43	!****************************************************************************************
	44	! Parameters (could be read in def file: move to slab_init)
	45	!****************************************************************************************
	46	! snow and ice physical characteristics:
	47	REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
	48	REAL, PARAMETER :: t_melt=273.15 ! melting ice temp
	49	REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
	50	REAL, PARAMETER :: ice_den=917.
	51	REAL, PARAMETER :: sea_den=1025.
	52	REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
	53	REAL, PARAMETER :: sno_cond=0.31*sno_den
	54	REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice
	55	REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
	56
	57	! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
	58	REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
	59	REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
	60	REAL, PARAMETER :: ice_frac_min=0.001
	61	REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction
	62	REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness
	63	REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness
	64	! below ice_thin, priority is melt lateral / grow height
	65	! ice_thin is also height of new ice
	66	REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness
	67	! above ice_thick, priority is melt height / grow lateral
	68	REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice
	69	REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height
	70
	71	! albedo and radiation parameters
	72	REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
	73	REAL, PARAMETER :: alb_sno_del=0.3 !max snow albedo = min + del
	74	REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
	75	REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
	76	REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
	77	REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
	78
	79	!****************************************************************************************
	80
[781]	81	CONTAINS
	82	!
	83	!****************************************************************************************
	84	!
[2073]	85	SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
	86	!, seaice_rst etc
[781]	87
[2073]	88	use IOIPSL
	89
	90	! For ok_xxx vars (Ekman...)
	91	INCLUDE "clesphys.h"
	92
	93	! Input variables
	94	!****************************************************************************************
	95	REAL, INTENT(IN) :: dtime
	96	! Variables read from restart file
	97	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
	98
	99	! Local variables
	100	!****************************************************************************************
	101	INTEGER :: error
	102	CHARACTER (len = 80) :: abort_message
	103	CHARACTER (len = 20) :: modname = 'ocean_slab_intit'
	104
	105	!****************************************************************************************
	106	! Allocate surface fraction read from restart file
	107	!****************************************************************************************
[2220]	108	ALLOCATE(fsic(klon), stat = error)
[2073]	109	IF (error /= 0) THEN
	110	abort_message='Pb allocation tmp_pctsrf_slab'
[2408]	111	CALL abort_physic(modname,abort_message,1)
[2073]	112	ENDIF
[2220]	113	fsic(:)=0.
	114	WHERE (1.-zmasq(:)>EPSFRA)
	115	fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
	116	END WHERE
[2073]	117
	118	!****************************************************************************************
[2220]	119	! Allocate saved variables
[2073]	120	!****************************************************************************************
[2220]	121	ALLOCATE(tslab(klon,nslay), stat=error)
[2408]	122	IF (error /= 0) CALL abort_physic &
[2220]	123	(modname,'pb allocation tslab', 1)
	124
[2073]	125	ALLOCATE(slab_bils(klon), stat = error)
	126	IF (error /= 0) THEN
	127	abort_message='Pb allocation slab_bils'
[2408]	128	CALL abort_physic(modname,abort_message,1)
[2073]	129	ENDIF
	130	slab_bils(:) = 0.0
	131	ALLOCATE(bils_cum(klon), stat = error)
	132	IF (error /= 0) THEN
	133	abort_message='Pb allocation slab_bils_cum'
[2408]	134	CALL abort_physic(modname,abort_message,1)
[2073]	135	ENDIF
	136	bils_cum(:) = 0.0
	137
[2220]	138	IF (version_ocean=='sicINT') THEN
	139	ALLOCATE(slab_bilg(klon), stat = error)
	140	IF (error /= 0) THEN
	141	abort_message='Pb allocation slab_bilg'
[2408]	142	CALL abort_physic(modname,abort_message,1)
[2220]	143	ENDIF
	144	slab_bilg(:) = 0.0
	145	ALLOCATE(bilg_cum(klon), stat = error)
	146	IF (error /= 0) THEN
	147	abort_message='Pb allocation slab_bilg_cum'
[2408]	148	CALL abort_physic(modname,abort_message,1)
[2220]	149	ENDIF
	150	bilg_cum(:) = 0.0
	151	ALLOCATE(tice(klon), stat = error)
	152	IF (error /= 0) THEN
	153	abort_message='Pb allocation slab_tice'
[2408]	154	CALL abort_physic(modname,abort_message,1)
[2220]	155	ENDIF
	156	ALLOCATE(seaice(klon), stat = error)
	157	IF (error /= 0) THEN
	158	abort_message='Pb allocation slab_seaice'
[2408]	159	CALL abort_physic(modname,abort_message,1)
[2220]	160	ENDIF
	161	END IF
	162
	163	!****************************************************************************************
	164	! Define some parameters
	165	!****************************************************************************************
[2073]	166	! Layer thickness
	167	ALLOCATE(slabh(nslay), stat = error)
	168	IF (error /= 0) THEN
	169	abort_message='Pb allocation slabh'
[2408]	170	CALL abort_physic(modname,abort_message,1)
[2073]	171	ENDIF
	172	slabh(1)=50.
	173	! cyang = 1/heat capacity of top layer (rho.c.H)
	174	cyang=1/(slabh(1)*4.228e+06)
	175
	176	! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
	177	! pour un calcul à chaque pas de temps, cpl_pas=1
	178	cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
	179	CALL getin('cpl_pas',cpl_pas)
	180	print *,'cpl_pas',cpl_pas
[2220]	181
[2073]	182	END SUBROUTINE ocean_slab_init
	183	!
	184	!****************************************************************************************
	185	!
	186	SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
	187
[996]	188	USE limit_read_mod
[1795]	189
[996]	190	! INCLUDE "clesphys.h"
[781]	191
[996]	192	! Arguments
[781]	193	!****************************************************************************************
[996]	194	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
	195	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
	196	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
[2073]	197	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction
[996]	198	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
[781]	199
	200	! Local variables
	201	!****************************************************************************************
	202
[2073]	203	IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN
	204	CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
[996]	205	ELSE
[2220]	206	pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
	207	pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
[2073]	208	is_modified=.TRUE.
[996]	209	END IF
[781]	210
[996]	211	END SUBROUTINE ocean_slab_frac
[781]	212	!
	213	!****************************************************************************************
	214	!
	215	SUBROUTINE ocean_slab_noice( &
[996]	216	itime, dtime, jour, knon, knindex, &
[2298]	217	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, &
[1067]	218	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	219	AcoefU, AcoefV, BcoefU, BcoefV, &
[2298]	220	ps, u1, v1, gustiness, tsurf_in, &
[2220]	221	radsol, snow, &
[1067]	222	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
[2073]	223	tsurf_new, dflux_s, dflux_l, qflux)
[1067]	224
	225	USE calcul_fluxs_mod
[1795]	226
[2298]	227	INCLUDE "clesphys.h"
[781]	228
	229	! Input arguments
	230	!****************************************************************************************
[996]	231	INTEGER, INTENT(IN) :: itime
	232	INTEGER, INTENT(IN) :: jour
[781]	233	INTEGER, INTENT(IN) :: knon
	234	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
	235	REAL, INTENT(IN) :: dtime
	236	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
[2298]	237	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm
[781]	238	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
	239	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
[1067]	240	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
	241	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
[781]	242	REAL, DIMENSION(klon), INTENT(IN) :: ps
[2298]	243	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
[996]	244	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
[2220]	245	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
[781]	246
	247	! In/Output arguments
	248	!****************************************************************************************
	249	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
	250
	251	! Output arguments
	252	!****************************************************************************************
	253	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
	254	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
[1067]	255	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
[781]	256	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
	257	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
[2073]	258	REAL, DIMENSION(klon), INTENT(OUT) :: qflux
[781]	259
	260	! Local variables
	261	!****************************************************************************************
[2073]	262	INTEGER :: i,ki
[2220]	263	! surface fluxes
[996]	264	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
[2220]	265	! flux correction
	266	REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
	267	! surface wind stress
[1067]	268	REAL, DIMENSION(klon) :: u0, v0
	269	REAL, DIMENSION(klon) :: u1_lay, v1_lay
[2220]	270	! ice creation
	271	REAL :: e_freeze, h_new, dfsic
[781]	272
	273	!****************************************************************************************
[996]	274	! 1) Flux calculation
	275	!
	276	!****************************************************************************************
	277	cal(:) = 0.
	278	beta(:) = 1.
	279	dif_grnd(:) = 0.
[781]	280
[1067]	281	! Suppose zero surface speed
	282	u0(:)=0.0
	283	v0(:)=0.0
	284	u1_lay(:) = u1(:) - u0(:)
	285	v1_lay(:) = v1(:) - v0(:)
	286
[781]	287	CALL calcul_fluxs(knon, is_oce, dtime, &
[2298]	288	tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
[781]	289	precip_rain, precip_snow, snow, qsurf, &
[2298]	290	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
	291	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
[781]	292	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
	293
[1067]	294	! - Flux calculation at first modele level for U and V
	295	CALL calcul_flux_wind(knon, dtime, &
[2298]	296	u0, v0, u1, v1, gustiness, cdragm, &
[1067]	297	AcoefU, AcoefV, BcoefU, BcoefV, &
	298	p1lay, temp_air, &
	299	flux_u1, flux_v1)
	300
[2073]	301	! Accumulate total fluxes locally
	302	slab_bils(:)=0.
	303	DO i=1,knon
	304	ki=knindex(i)
[2220]	305	slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
	306	-precip_snow(i)ice_lat(1.+snow_wfact*fsic(ki)))
[2073]	307	bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
	308	! Also taux, tauy, saved vars...
	309	END DO
	310
[781]	311	!****************************************************************************************
[2073]	312	! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
[781]	313	!
	314	!****************************************************************************************
[2073]	315	lmt_bils(:)=0.
[2220]	316	CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
	317	! lmt_bils and diff_sst,siv saved by limit_slab
	318	qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)ice_denice_lat/86400.
[2073]	319	! qflux = total QFlux correction (in W/m2)
[781]	320
	321	!****************************************************************************************
[996]	322	! 3) Recalculate new temperature
[781]	323	!
[2073]	324	!*********************************************o***************************************
	325	tsurf_new=tsurf_in
	326	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
	327	! Compute transport
	328	! Add read QFlux and SST tendency
	329	tslab(:,1)=tslab(:,1)+qflux(:)cyangdtime*cpl_pas
	330	! Add cumulated surface fluxes
	331	tslab(:,1)=tslab(:,1)+bils_cum(:)cyangdtime
	332	! Update surface temperature
	333	SELECT CASE(version_ocean)
	334	CASE('sicNO')
	335	DO i=1,knon
	336	ki=knindex(i)
	337	tsurf_new(i)=tslab(ki,1)
	338	END DO
[2220]	339	CASE('sicOBS') ! check for sea ice or tslab below freezing
[2073]	340	DO i=1,knon
	341	ki=knindex(i)
[2220]	342	IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
[2073]	343	tslab(ki,1)=t_freeze
	344	END IF
[2220]	345	tsurf_new(i)=tslab(ki,1)
[2073]	346	END DO
	347	CASE('sicINT')
	348	DO i=1,knon
	349	ki=knindex(i)
[2220]	350	IF (fsic(ki).LT.epsfra) THEN ! Free of ice
	351	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
	352	! quantity of new ice formed
	353	e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
	354	! new ice
	355	tice(ki)=t_freeze
	356	fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
	357	IF (fsic(ki).GT.ice_frac_min) THEN
	358	seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
	359	tslab(ki,1)=t_freeze
	360	ELSE
	361	fsic(ki)=0.
	362	END IF
	363	tsurf_new(i)=t_freeze
	364	ELSE
[2073]	365	tsurf_new(i)=tslab(ki,1)
[2220]	366	END IF
	367	ELSE ! ice present
[2073]	368	tsurf_new(i)=t_freeze
[2220]	369	IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
	370	! quantity of new ice formed over open ocean
	371	e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
	372	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
	373	! new ice height and fraction
	374	h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
	375	dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
	376	h_new=MIN(e_freeze/dfsic,h_ice_max)
	377	! update tslab to freezing over open ocean only
	378	tslab(ki,1)=tslab(ki,1)fsic(ki)+t_freeze(1.-fsic(ki))
	379	! update sea ice
	380	seaice(ki)=(h_newdfsic+seaice(ki)fsic(ki)) &
	381	/(dfsic+fsic(ki))
	382	fsic(ki)=fsic(ki)+dfsic
	383	! update snow?
	384	END IF !freezing
	385	END IF ! sea ice present
[2073]	386	END DO
	387	END SELECT
	388	bils_cum(:)=0.0! clear cumulated fluxes
	389	END IF ! coupling time
	390	END SUBROUTINE ocean_slab_noice
	391	!
[781]	392	!****************************************************************************************
[2220]	393
	394	SUBROUTINE ocean_slab_ice( &
	395	itime, dtime, jour, knon, knindex, &
	396	tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
	397	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	398	AcoefU, AcoefV, BcoefU, BcoefV, &
[2298]	399	ps, u1, v1, gustiness, &
[2220]	400	radsol, snow, qsurf, qsol, agesno, &
	401	alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
	402	tsurf_new, dflux_s, dflux_l, swnet)
	403
	404	USE calcul_fluxs_mod
	405
	406	INCLUDE "YOMCST.h"
[2298]	407	INCLUDE "clesphys.h"
[2220]	408
	409	! Input arguments
[2073]	410	!****************************************************************************************
[2220]	411	INTEGER, INTENT(IN) :: itime, jour, knon
	412	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
	413	REAL, INTENT(IN) :: dtime
	414	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
	415	REAL, DIMENSION(klon), INTENT(IN) :: p1lay
	416	REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
	417	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
	418	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
	419	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
	420	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
	421	REAL, DIMENSION(klon), INTENT(IN) :: ps
[2298]	422	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
[2220]	423	REAL, DIMENSION(klon), INTENT(IN) :: swnet
	424
	425	! In/Output arguments
	426	!****************************************************************************************
	427	REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
	428	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
	429	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
	430
	431	! Output arguments
	432	!****************************************************************************************
	433	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
	434	REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
	435	REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
	436	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
	437	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
	438	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
	439	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
	440
	441	! Local variables
	442	!****************************************************************************************
	443	INTEGER :: i,ki
	444	REAL, DIMENSION(klon) :: cal, beta, dif_grnd
	445	REAL, DIMENSION(klon) :: u0, v0
	446	REAL, DIMENSION(klon) :: u1_lay, v1_lay
	447	! intermediate heat fluxes:
	448	REAL :: f_cond, f_swpen
	449	! for snow/ice albedo:
	450	REAL :: alb_snow, alb_ice, alb_pond
	451	REAL :: frac_snow, frac_ice, frac_pond
	452	! for ice melt / freeze
	453	REAL :: e_melt, snow_evap, h_test
	454	! dhsic, dfsic change in ice mass, fraction.
	455	REAL :: dhsic, dfsic, frac_mf
	456
	457	!****************************************************************************************
[2073]	458	! 1) Flux calculation
	459	!****************************************************************************************
[2220]	460	! Suppose zero surface speed
	461	u0(:)=0.0
	462	v0(:)=0.0
	463	u1_lay(:) = u1(:) - u0(:)
	464	v1_lay(:) = v1(:) - v0(:)
	465
	466	! set beta, cal, compute conduction fluxes inside ice/snow
	467	slab_bilg(:)=0.
	468	dif_grnd(:)=0.
	469	beta(:) = 1.
	470	DO i=1,knon
	471	ki=knindex(i)
	472	IF (snow(i).GT.snow_min) THEN
	473	! snow-layer heat capacity
	474	cal(i)=2.RCPD/(snow(i)ice_cap)
	475	! snow conductive flux
	476	f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
	477	! all shortwave flux absorbed
	478	f_swpen=0.
	479	! bottom flux (ice conduction)
	480	slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
	481	! update ice temperature
	482	tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
	483	(slab_bilg(ki)+f_cond)dtime
	484	ELSE ! bare ice
	485	! ice-layer heat capacity
	486	cal(i)=2.RCPD/(seaice(ki)ice_cap)
	487	! conductive flux
	488	f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
	489	! penetrative shortwave flux...
	490	f_swpen=swnet(i)pen_fracexp(-pen_ext*seaice(ki)/ice_den)
	491	slab_bilg(ki)=f_swpen-f_cond
	492	END IF
	493	radsol(i)=radsol(i)+f_cond-f_swpen
	494	END DO
	495	! weight fluxes to ocean by sea ice fraction
	496	slab_bilg(:)=slab_bilg(:)*fsic(:)
	497
[2073]	498	! calcul_fluxs (sens, lat etc)
[2220]	499	CALL calcul_fluxs(knon, is_sic, dtime, &
[2298]	500	tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
[2220]	501	precip_rain, precip_snow, snow, qsurf, &
[2298]	502	radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
	503	f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
[2220]	504	tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
	505	DO i=1,knon
	506	IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
	507	END DO
	508
[2073]	509	! calcul_flux_wind
[2220]	510	CALL calcul_flux_wind(knon, dtime, &
[2298]	511	u0, v0, u1, v1, gustiness, cdragm, &
[2220]	512	AcoefU, AcoefV, BcoefU, BcoefV, &
	513	p1lay, temp_air, &
	514	flux_u1, flux_v1)
[781]	515
[2073]	516	!****************************************************************************************
[2220]	517	! 2) Update snow and ice surface
[2073]	518	!****************************************************************************************
[2220]	519	! snow precip
	520	DO i=1,knon
	521	ki=knindex(i)
	522	IF (precip_snow(i) > 0.) THEN
	523	snow(i) = snow(i)+precip_snow(i)dtime(1.-snow_wfact*(1.-fsic(ki)))
	524	END IF
	525	! snow and ice sublimation
	526	IF (evap(i) > 0.) THEN
	527	snow_evap = MIN (snow(i) / dtime, evap(i))
	528	snow(i) = snow(i) - snow_evap * dtime
	529	snow(i) = MAX(0.0, snow(i))
	530	seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
	531	ENDIF
	532	! Melt / Freeze from above if Tsurf>0
	533	IF (tsurf_new(i).GT.t_melt) THEN
	534	! energy available for melting snow (in kg/m2 of snow)
	535	e_melt = MIN(MAX(snow(i)(tsurf_new(i)-t_melt)ice_cap/2. &
	536	/(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
	537	! remove snow
	538	IF (snow(i).GT.e_melt) THEN
	539	snow(i)=snow(i)-e_melt
	540	tsurf_new(i)=t_melt
	541	ELSE ! all snow is melted
	542	! add remaining heat flux to ice
	543	e_melt=e_melt-snow(i)
	544	tice(ki)=tice(ki)+e_meltice_lat2./(ice_cap*seaice(ki))
	545	tsurf_new(i)=tice(ki)
	546	END IF
	547	END IF
	548	! melt ice from above if Tice>0
	549	IF (tice(ki).GT.t_melt) THEN
	550	! quantity of ice melted (kg/m2)
	551	e_melt=MAX(seaice(ki)(tice(ki)-t_melt)ice_cap/2. &
	552	/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
	553	! melt from above, height only
	554	dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
	555	e_melt=e_melt-dhsic
	556	IF (e_melt.GT.0) THEN
	557	! lateral melt if ice too thin
	558	dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
	559	! if all melted add remaining heat to ocean
	560	e_melt=MAX(0.,e_meltfsic(ki)-dfsich_ice_min)
	561	slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
	562	! update height and fraction
	563	fsic(ki)=fsic(ki)-dfsic
	564	END IF
	565	seaice(ki)=seaice(ki)-dhsic
	566	! surface temperature at melting point
	567	tice(ki)=t_melt
	568	tsurf_new(i)=t_melt
	569	END IF
	570	! convert snow to ice if below floating line
	571	h_test=(seaice(ki)+snow(i))ice_den-seaice(ki)sea_den
	572	IF (h_test.GT.0.) THEN !snow under water
	573	! extra snow converted to ice (with added frozen sea water)
	574	dhsic=h_test/(sea_den-ice_den+sno_den)
	575	seaice(ki)=seaice(ki)+dhsic
	576	snow(i)=snow(i)-dhsic*sno_den/ice_den
	577	! available energy (freeze sea water + bring to tice)
	578	e_melt=dhsic(1.-sno_den/ice_den)(ice_lat+ &
	579	ice_cap/2.*(t_freeze-tice(ki)))
	580	! update ice temperature
	581	tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
	582	END IF
	583	END DO
	584
[2073]	585	! New albedo
[2220]	586	DO i=1,knon
	587	ki=knindex(i)
	588	! snow albedo: update snow age
	589	IF (snow(i).GT.0.0001) THEN
	590	agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
	591	* EXP(-1.MAX(0.0,precip_snow(i))dtime/5.)
	592	ELSE
	593	agesno(i)=0.0
	594	END IF
	595	! snow albedo
	596	alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
	597	! ice albedo (varies with ice tkickness and temp)
	598	alb_ice=MAX(0.0,0.13LOG(100.seaice(ki)/ice_den)+0.1)
	599	IF (tice(ki).GT.t_freeze-0.01) THEN
	600	alb_ice=MIN(alb_ice,alb_ice_wet)
	601	ELSE
	602	alb_ice=MIN(alb_ice,alb_ice_dry)
	603	END IF
	604	! pond albedo
	605	alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
	606	! pond fraction
	607	frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
	608	! snow fraction
	609	frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
	610	! ice fraction
	611	frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
	612	! total albedo
	613	alb1_new(i)=alb_snowfrac_snow+alb_icefrac_ice+alb_pond*frac_pond
	614	END DO
	615	alb2_new(:) = alb1_new(:)
[2073]	616
	617	!****************************************************************************************
[2220]	618	! 3) Recalculate new ocean temperature (add fluxes below ice)
[2073]	619	! Melt / freeze from below
	620	!*********************************************o***************************************
[2220]	621	!cumul fluxes
	622	bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
	623	IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
	624	! Add cumulated surface fluxes
	625	tslab(:,1)=tslab(:,1)+bilg_cum(:)cyangdtime
	626	DO i=1,knon
	627	ki=knindex(i)
	628	! split lateral/top melt-freeze
	629	frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
	630	IF (tslab(ki,1).LE.t_freeze) THEN
	631	! **** Form new ice from below *****
	632	! quantity of new ice
	633	e_melt=(t_freeze-tslab(ki,1))/cyang &
	634	/(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
	635	! first increase height to h_thin
	636	dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
	637	seaice(ki)=dhsic+seaice(ki)
	638	e_melt=e_melt-fsic(ki)*dhsic
	639	IF (e_melt.GT.0.) THEN
	640	! frac_mf fraction used for lateral increase
	641	dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
	642	fsic(ki)=fsic(ki)+dfsic
	643	e_melt=e_melt-dfsic*seaice(ki)
	644	! rest used to increase height
	645	seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
	646	END IF
	647	tslab(ki,1)=t_freeze
	648	ELSE ! slab temperature above freezing
	649	! **** melt ice from below *****
	650	! quantity of melted ice
	651	e_melt=(tslab(ki,1)-t_freeze)/cyang &
	652	/(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
	653	! first decrease height to h_thick
	654	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
	655	seaice(ki)=seaice(ki)-dhsic
	656	e_melt=e_melt-fsic(ki)*dhsic
	657	IF (e_melt.GT.0) THEN
	658	! frac_mf fraction used for height decrease
	659	dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
	660	seaice(ki)=seaice(ki)-dhsic
	661	e_melt=e_melt-fsic(ki)*dhsic
	662	! rest used to decrease fraction (up to 0!)
	663	dfsic=MIN(fsic(ki),e_melt/seaice(ki))
	664	! keep remaining in ocean
	665	e_melt=e_melt-dfsic*seaice(ki)
	666	END IF
	667	tslab(ki,1)=t_freeze+e_meltice_latcyang
	668	fsic(ki)=fsic(ki)-dfsic
	669	END IF
	670	END DO
	671	bilg_cum(:)=0.
	672	END IF ! coupling time
	673
	674	!tests fraction glace
	675	WHERE (fsic.LT.ice_frac_min)
	676	tslab(:,1)=tslab(:,1)-fsicseaiceice_lat*cyang
	677	tice=t_melt
	678	fsic=0.
	679	seaice=0.
	680	END WHERE
[2073]	681
[2220]	682	END SUBROUTINE ocean_slab_ice
[781]	683	!
	684	!****************************************************************************************
	685	!
[2073]	686	SUBROUTINE ocean_slab_final
	687
	688	!****************************************************************************************
	689	! Deallocate module variables
	690	!****************************************************************************************
	691	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
[2220]	692	IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
	693	IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
	694	IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
	695	IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
	696	IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
	697	IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
[2073]	698
	699	END SUBROUTINE ocean_slab_final
	700
[781]	701	END MODULE ocean_slab_mod

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

Download in other formats:

Original Format