Context Navigation

surf_ocean_mod.F90 @ 5209

Last change on this file since 5209 was 5202, checked in by Laurent Fairhead, 24 hours ago
Updating cirrus branch to trunk revision 5171
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:keywords set to `Author Date Id Revision`
File size: 16.3 KB

Rev	Line
[781]	1	!
[2538]	2	! $Id: surf_ocean_mod.F90 5202 2024-09-20 10:32:04Z fairhead $
	3	!
[781]	4	MODULE surf_ocean_mod
	5
	6	IMPLICIT NONE
	7
	8	CONTAINS
[3395]	9	!
	10	!******************************************************************************
	11	!
[888]	12	SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
[2243]	13	windsp, rmu0, fder, tsurf_in, &
[781]	14	itime, dtime, jour, knon, knindex, &
[4523]	15	p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
[1067]	16	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	17	AcoefU, AcoefV, BcoefU, BcoefV, &
[2240]	18	ps, u1, v1, gustiness, rugoro, pctsrf, &
[888]	19	snow, qsurf, agesno, &
[2243]	20	z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &
[1067]	21	tsurf_new, dflux_s, dflux_l, lmt_bils, &
[4370]	22	flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, &
[5202]	23	dt_ds, tkt, tks, taur, sss &
	24	#ifdef ISO
	25	& ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
	26	& xtsnow,xtevap,h1 &
	27	#endif
	28	& )
[1067]	29
[3395]	30	use albedo, only: alboc, alboc_cd
[3815]	31	use bulk_flux_m, only: bulk_flux
[3395]	32	USE dimphy, ONLY: klon, zmasq
	33	USE surface_data, ONLY : type_ocean
	34	USE ocean_forced_mod, ONLY : ocean_forced_noice
	35	USE ocean_slab_mod, ONLY : ocean_slab_noice
	36	USE ocean_cpl_mod, ONLY : ocean_cpl_noice
	37	USE indice_sol_mod, ONLY : nbsrf, is_oce
[5202]	38	#ifdef ISO
	39	USE infotrac_phy, ONLY : ntraciso=>ntiso,niso
	40	#ifdef ISOVERIF
	41	USE isotopes_mod, ONLY: iso_eau,ridicule
	42	USE isotopes_verif_mod
	43	#endif
	44	#endif
[3395]	45	USE limit_read_mod
[5202]	46	USE config_ocean_skin_m, ONLY: activate_ocean_skin
[3395]	47	!
	48	! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force,
	49	! slab or couple). The calculations of albedo and rugosity for the ocean surface are
	50	! done in here because they are identical for the different modes of ocean.
[1785]	51
	52
[793]	53	INCLUDE "YOMCST.h"
[781]	54
[2178]	55	include "clesphys.h"
[2455]	56	! for cycle_diurne and for iflag_z0_oce==-1 (prescribed z0)
[2178]	57
[3395]	58	! Input variables
	59	!******************************************************************************
[781]	60	INTEGER, INTENT(IN) :: itime, jour, knon
	61	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
	62	REAL, INTENT(IN) :: dtime
	63	REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
[888]	64	REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiation at surface
	65	REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiation at surface
	66	REAL, DIMENSION(klon), INTENT(IN) :: alb1 ! albedo in visible SW interval
[3815]	67	REAL, DIMENSION(klon), INTENT(IN) :: windsp ! wind at 10 m, in m s-1
[781]	68	REAL, DIMENSION(klon), INTENT(IN) :: rmu0
	69	REAL, DIMENSION(klon), INTENT(IN) :: fder
[3906]	70	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! defined only for subscripts 1:knon
[2254]	71	REAL, DIMENSION(klon), INTENT(IN) :: p1lay,z1lay ! pression (Pa) et altitude (m) du premier niveau
[1067]	72	REAL, DIMENSION(klon), INTENT(IN) :: cdragh
	73	REAL, DIMENSION(klon), INTENT(IN) :: cdragm
[4523]	74	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow, precip_bs
[781]	75	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
[1067]	76	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
	77	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
[781]	78	REAL, DIMENSION(klon), INTENT(IN) :: ps
[2240]	79	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
[781]	80	REAL, DIMENSION(klon), INTENT(IN) :: rugoro
	81	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
[5202]	82	#ifdef ISO
	83	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
	84	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtspechum
	85	#endif
[781]	86
[3395]	87	! In/Output variables
	88	!******************************************************************************
[888]	89	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
	90	REAL, DIMENSION(klon), INTENT(INOUT) :: qsurf
[781]	91	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
[3906]	92	REAL, DIMENSION(klon), INTENT(inOUT) :: z0h
[5202]	93	#ifdef ISO
	94	REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsnow
	95	REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce
	96	#endif
[781]	97
[3815]	98	REAL, intent(inout):: delta_sst(:) ! (knon)
	99	! Ocean-air interface temperature minus bulk SST, in K. Defined
	100	! only if activate_ocean_skin >= 1.
	101
	102	real, intent(inout):: delta_sal(:) ! (knon)
	103	! Ocean-air interface salinity minus bulk salinity, in ppt. Defined
	104	! only if activate_ocean_skin >= 1.
	105
	106	REAL, intent(inout):: ds_ns(:) ! (knon)
	107	! "delta salinity near surface". Salinity variation in the
	108	! near-surface turbulent layer. That is subskin salinity minus
	109	! foundation salinity. In ppt.
	110
	111	REAL, intent(inout):: dt_ns(:) ! (knon)
	112	! "delta temperature near surface". Temperature variation in the
	113	! near-surface turbulent layer. That is subskin temperature
	114	! minus foundation temperature. (Can be negative.) In K.
	115
[4370]	116	REAL, intent(inout):: dter(:) ! (knon)
	117	! Temperature variation in the diffusive microlayer, that is
	118	! ocean-air interface temperature minus subskin temperature. In
	119	! K.
	120
	121	REAL, intent(inout):: dser(:) ! (knon)
	122	! Salinity variation in the diffusive microlayer, that is
	123	! ocean-air interface salinity minus subskin salinity. In ppt.
	124
	125	real, intent(inout):: dt_ds(:) ! (knon)
	126	! (tks / tkt) * dTer, in K
	127
[3395]	128	! Output variables
[3815]	129	!**************************************************************************
[3389]	130	REAL, DIMENSION(klon), INTENT(OUT) :: z0m
[3395]	131	!albedo SB >>>
[3906]	132	! REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
	133	! REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
	134	REAL, DIMENSION(6), INTENT(IN) :: SFRWL
	135	REAL, DIMENSION(klon,nsw), INTENT(OUT) :: alb_dir_new,alb_dif_new
[3395]	136	!albedo SB <<<
[781]	137	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
[3906]	138	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! sea surface temperature, in K
[781]	139	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
[996]	140	REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils
[1067]	141	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
[781]	142
[3815]	143	REAL, intent(out):: tkt(:) ! (knon)
	144	! épaisseur (m) de la couche de diffusion thermique (microlayer)
	145	! cool skin thickness
	146
	147	REAL, intent(out):: tks(:) ! (knon)
	148	! épaisseur (m) de la couche de diffusion de masse (microlayer)
	149
	150	REAL, intent(out):: taur(:) ! (knon)
	151	! momentum flux due to rain, in Pa
	152
	153	real, intent(out):: sss(:) ! (klon)
	154	! Bulk salinity of the surface layer of the ocean, in ppt. (Only
	155	! defined for subscripts 1:knon, but we have to declare it with
	156	! size klon because of the coupling machinery.)
	157
[5202]	158	#ifdef ISO
	159	REAL, DIMENSION(ntraciso,klon), INTENT(out) :: xtevap ! isotopes in surface evaporation flux
	160	REAL, DIMENSION(klon), INTENT(out) :: h1 ! just a diagnostic, not useful for the simulation
	161	#endif
	162
[3395]	163	! Local variables
[3815]	164	!*************************************************************************
[2227]	165	INTEGER :: i, k
[1146]	166	REAL :: tmp
	167	REAL, PARAMETER :: cepdu2=(0.1)**2
[3002]	168	REAL, DIMENSION(klon) :: alb_eau, z0_lim
[888]	169	REAL, DIMENSION(klon) :: radsol
[2254]	170	REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
[4523]	171	REAL, DIMENSION(klon) :: precip_totsnow
[2391]	172	CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
[5202]	173	REAL rhoa(knon) ! density of moist air (kg / m3)
[3815]	174	REAL sens_prec_liq(knon)
[781]	175
[3815]	176	REAL t_int(knon) ! ocean-air interface temperature, in K
[5202]	177	REAL s_int(knon) ! ocean-air interface salinity, in ppt
[888]	178
[3815]	179	!**************************************************************************
[888]	180
[5202]	181	#ifdef ISO
	182	#ifdef ISOVERIF
	183	DO i = 1, knon
	184	IF (iso_eau > 0) THEN
	185	CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), &
	186	& spechum(i),'surf_ocean_mod 117', &
	187	& errmax,errmaxrel)
	188	CALL iso_verif_egalite_choix(xtsnow(iso_eau,i), &
	189	& snow(i),'surf_ocean_mod 127', &
	190	& errmax,errmaxrel)
	191	ENDIF !IF (iso_eau > 0) then
	192	ENDDO !DO i=1,klon
	193	#endif
	194	#endif
[3815]	195
[3395]	196	!******************************************************************************
	197	! Calculate total net radiance at surface
	198	!
	199	!******************************************************************************
[2719]	200	radsol(1:klon) = 0.0 ! initialisation a priori inutile
[888]	201	radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
	202
[4523]	203
	204	!****************************************************************************************
	205	!Total solid precip
	206
	207	IF (ok_bs) THEN
[4526]	208	precip_totsnow(:)=precip_snow(:)+precip_bs(:)
[4523]	209	ELSE
[4526]	210	precip_totsnow(:)=precip_snow(:)
[4523]	211	ENDIF
	212
	213
[3395]	214	!******************************************************************************
	215	! Cdragq computed from cdrag
	216	! The difference comes only from a factor (f_z0qh_oce) on z0, so that
	217	! it can be computed inside surf_ocean
	218	! More complicated appraches may require the propagation through
	219	! pbl_surface of an independant cdragq variable.
	220	!******************************************************************************
[2254]	221
	222	IF ( f_z0qh_oce .ne. 1.) THEN
[3395]	223	! Si on suit les formulations par exemple de Tessel, on
	224	! a z0h=0.4nu/u, z0q=0.62nu/u, d'ou f_z0qh_oce=0.62/0.4=1.55
[2718]	225	cdragq(1:knon)=cdragh(1:knon)* &
[3395]	226	log(z1lay(1:knon)/z0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*z0h(1:knon)))
[2254]	227	ELSE
[2718]	228	cdragq(1:knon)=cdragh(1:knon)
[2254]	229	ENDIF
	230
[3815]	231	rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon)))
[3395]	232	!******************************************************************************
	233	! Switch according to type of ocean (couple, slab or forced)
	234	!******************************************************************************
[996]	235	SELECT CASE(type_ocean)
[781]	236	CASE('couple')
[888]	237	CALL ocean_cpl_noice( &
	238	swnet, lwnet, alb1, &
[1067]	239	windsp, fder, &
[781]	240	itime, dtime, knon, knindex, &
[4523]	241	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow,temp_air,spechum,&
[1067]	242	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	243	AcoefU, AcoefV, BcoefU, BcoefV, &
[3815]	244	ps, u1, v1, gustiness, tsurf_in, &
[888]	245	radsol, snow, agesno, &
[1067]	246	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
[3815]	247	tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, delta_sal, rhoa, &
[4370]	248	delta_sst, dTer, dSer, dt_ds)
[781]	249
	250	CASE('slab')
[888]	251	CALL ocean_slab_noice( &
[996]	252	itime, dtime, jour, knon, knindex, &
[4523]	253	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, temp_air, spechum,&
[1067]	254	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	255	AcoefU, AcoefV, BcoefU, BcoefV, &
[2240]	256	ps, u1, v1, gustiness, tsurf_in, &
[2209]	257	radsol, snow, &
[1067]	258	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
[996]	259	tsurf_new, dflux_s, dflux_l, lmt_bils)
[3395]	260
[781]	261	CASE('force')
[888]	262	CALL ocean_forced_noice( &
	263	itime, dtime, jour, knon, knindex, &
[4523]	264	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, &
[781]	265	temp_air, spechum, &
[1067]	266	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	267	AcoefU, AcoefV, BcoefU, BcoefV, &
[3815]	268	ps, u1, v1, gustiness, tsurf_in, &
[888]	269	radsol, snow, agesno, &
[1067]	270	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
[5202]	271	tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
	272	#ifdef ISO
	273	,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
	274	xtsnow,xtevap,h1 &
	275	#endif
	276	)
[781]	277	END SELECT
	278
[3395]	279	!******************************************************************************
	280	! fcodron: compute lmt_bils forced case (same as wfbils_oce / 1.-contfracatm)
	281	!******************************************************************************
[2057]	282	IF (type_ocean.NE.'slab') THEN
[3395]	283	lmt_bils(1:klon)=0.
	284	DO i=1,knon
	285	lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
	286	*pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i)))
	287	END DO
[2057]	288	END IF
	289
[3395]	290	!******************************************************************************
	291	! Calculate ocean surface albedo
	292	!******************************************************************************
	293	!albedo SB >>>
	294	IF (iflag_albedo==0) THEN
	295	!--old parametrizations of ocean surface albedo
	296	!
	297	IF (iflag_cycle_diurne.GE.1) THEN
	298	!
	299	CALL alboc_cd(rmu0,alb_eau)
	300	!
	301	!--ad-hoc correction for model radiative balance tuning
	302	!--now outside alboc_cd routine
	303	alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
	304	alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.0),1.0)
	305	!
	306	ELSE
	307	!
	308	CALL alboc(REAL(jour),rlat,alb_eau)
	309	!--ad-hoc correction for model radiative balance tuning
	310	!--now outside alboc routine
	311	alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
	312	alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.04),0.60)
	313	!
	314	ENDIF
	315	!
	316	DO i =1, knon
	317	DO k=1,nsw
	318	alb_dir_new(i,k) = alb_eau(knindex(i))
	319	ENDDO
	320	ENDDO
	321	!IM 09122015 next line corresponds to the old way of doing in LMDZ5A/IPSLCM5A versions
	322	!albedo for diffuse radiation is taken the same as for direct radiation
	323	alb_dif_new(1:knon,:)=alb_dir_new(1:knon,:)
	324	!IM 09122015 end
	325	!
	326	ELSE IF (iflag_albedo==1) THEN
	327	!--new parametrization of ocean surface albedo by Sunghye Baek
	328	!--albedo for direct and diffuse radiation are different
	329	!
	330	CALL ocean_albedo(knon,rmu0,knindex,windsp,SFRWL,alb_dir_new,alb_dif_new)
	331	!
	332	!--ad-hoc correction for model radiative balance tuning
	333	alb_dir_new(1:knon,:) = fmagic*alb_dir_new(1:knon,:) + pmagic
	334	alb_dif_new(1:knon,:) = fmagic*alb_dif_new(1:knon,:) + pmagic
	335	alb_dir_new(1:knon,:)=MIN(MAX(alb_dir_new(1:knon,:),0.0),1.0)
	336	alb_dif_new(1:knon,:)=MIN(MAX(alb_dif_new(1:knon,:),0.0),1.0)
	337	!
	338	ELSE IF (iflag_albedo==2) THEN
	339	! F. Codron albedo read from limit.nc
	340	CALL limit_read_rug_alb(itime, dtime, jour,&
	341	knon, knindex, z0_lim, alb_eau)
	342	DO i =1, knon
	343	DO k=1,nsw
	344	alb_dir_new(i,k) = alb_eau(i)
	345	ENDDO
	346	ENDDO
	347	alb_dif_new=alb_dir_new
[781]	348	ENDIF
[3395]	349	!albedo SB <<<
[2227]	350
[3395]	351	!******************************************************************************
	352	! Calculate the rugosity
	353	!******************************************************************************
	354	IF (iflag_z0_oce==0) THEN
	355	DO i = 1, knon
	356	tmp = MAX(cepdu2,gustiness(i)+u1(i)2+v1(i)2)
	357	z0m(i) = 0.018cdragm(i) (gustiness(i)+u1(i)2+v1(i)2)/RG &
	358	+ 0.1114e-6 / SQRT(cdragm(i) tmp)
	359	z0m(i) = MAX(1.5e-05,z0m(i))
	360	ENDDO
	361	z0h(1:knon)=z0m(1:knon) ! En attendant mieux
[2243]	362
[3395]	363	ELSE IF (iflag_z0_oce==1) THEN
	364	DO i = 1, knon
	365	tmp = MAX(cepdu2,gustiness(i)+u1(i)2+v1(i)2)
	366	z0m(i) = 0.018cdragm(i) (gustiness(i)+u1(i)2+v1(i)2)/RG &
	367	+ 0.1114e-6 / SQRT(cdragm(i) tmp)
	368	z0m(i) = MAX(1.5e-05,z0m(i))
	369	z0h(i)=0.414e-6 / SQRT(cdragm(i) tmp)
	370	ENDDO
	371	ELSE IF (iflag_z0_oce==-1) THEN
	372	DO i = 1, knon
	373	z0m(i) = z0min
	374	z0h(i) = z0min
	375	ENDDO
	376	ELSE
[2391]	377	CALL abort_physic(modname,'version non prevue',1)
[3395]	378	ENDIF
[3815]	379
	380	if (activate_ocean_skin >= 1) then
	381	if (type_ocean /= 'couple') sss(:knon) = 35.
	382	call bulk_flux(tkt, tks, taur, dter, dser, t_int, s_int, ds_ns, dt_ns, &
	383	u = windsp(:knon), t_ocean_1 = tsurf_new(:knon), s1 = sss(:knon), &
[4523]	384	rain = precip_rain(:knon) + precip_totsnow(:knon), &
[3815]	385	hf = - fluxsens(:knon), hlb = - fluxlat(:knon), &
	386	rnl = - lwnet(:knon), &
	387	tau = sqrt(flux_u1(:knon)2 + flux_v1(:knon)2), rhoa = rhoa, &
	388	xlv = [(rlvtt, i = 1, knon)], rf = - sens_prec_liq, dtime = dtime, &
	389	rns = swnet(:knon))
	390	delta_sst = t_int - tsurf_new(:knon)
	391	delta_sal = s_int - sss(:knon)
[4370]	392
	393	if (activate_ocean_skin == 2) then
	394	tsurf_new(:knon) = t_int
	395	if (type_ocean == 'couple') dt_ds = (tks / tkt) * dter
	396	end if
[3815]	397	end if
	398
[781]	399	END SUBROUTINE surf_ocean
[3815]	400	!****************************************************************************
[3395]	401	!
[781]	402	END MODULE surf_ocean_mod

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

Download in other formats:

Original Format