Context Navigation

surf_ocean_mod.F90 @ 5210

Last change on this file since 5210 was 5202, checked in by Laurent Fairhead, 7 weeks 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 evignon $
	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