Context Navigation

surf_ocean_mod.F90 @ 4441

Last change on this file since 4441 was 4374, checked in by lguez, 19 months ago

Report modifications from phylmd into phylmdiso

Report modifications of revision 4370 from phylmd into phylmdiso.

File size: 16.1 KB

Line
1	!
2	! $Id: surf_ocean_mod.F90 3395 2018-09-25 15:22:13Z lguez $
3	!
4	MODULE surf_ocean_mod
5
6	IMPLICIT NONE
7
8	CONTAINS
9	!
10	!******************************************************************************
11	!
12	SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
13	windsp, rmu0, fder, tsurf_in, &
14	itime, dtime, jour, knon, knindex, &
15	p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
16	AcoefH, AcoefQ, BcoefH, BcoefQ, &
17	AcoefU, AcoefV, BcoefU, BcoefV, &
18	ps, u1, v1, gustiness, rugoro, pctsrf, &
19	snow, qsurf, agesno, &
20	z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &
21	tsurf_new, dflux_s, dflux_l, lmt_bils, &
22	flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, &
23	dt_ds, tkt, tks, taur, sss &
24	#ifdef ISO
25	& ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
26	& xtsnow,xtevap,h1 &
27	#endif
28	& )
29
30	use albedo, only: alboc, alboc_cd
31	use bulk_flux_m, only: bulk_flux
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
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
45	USE limit_read_mod
46	use config_ocean_skin_m, only: activate_ocean_skin
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.
51
52
53	INCLUDE "YOMCST.h"
54
55	include "clesphys.h"
56	! for cycle_diurne and for iflag_z0_oce==-1 (prescribed z0)
57
58	! Input variables
59	!******************************************************************************
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
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
67	REAL, DIMENSION(klon), INTENT(IN) :: windsp ! wind at 10 m, in m s-1
68	REAL, DIMENSION(klon), INTENT(IN) :: rmu0
69	REAL, DIMENSION(klon), INTENT(IN) :: fder
70	REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! defined only for subscripts 1:knon
71	REAL, DIMENSION(klon), INTENT(IN) :: p1lay,z1lay ! pression (Pa) et altitude (m) du premier niveau
72	REAL, DIMENSION(klon), INTENT(IN) :: cdragh
73	REAL, DIMENSION(klon), INTENT(IN) :: cdragm
74	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
75	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
76	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
77	REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
78	REAL, DIMENSION(klon), INTENT(IN) :: ps
79	REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
80	REAL, DIMENSION(klon), INTENT(IN) :: rugoro
81	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
82	#ifdef ISO
83	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
84	REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtspechum
85	#endif
86
87	! In/Output variables
88	!******************************************************************************
89	REAL, DIMENSION(klon), INTENT(INOUT) :: snow
90	REAL, DIMENSION(klon), INTENT(INOUT) :: qsurf
91	REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
92	REAL, DIMENSION(klon), INTENT(inOUT) :: z0h
93	#ifdef ISO
94	REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsnow
95	REAL, DIMENSION(niso,klon), INTENT(INOUT) :: Roce
96	#endif
97
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
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
128	! Output variables
129	!**************************************************************************
130	REAL, DIMENSION(klon), INTENT(OUT) :: z0m
131	!albedo SB >>>
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
136	!albedo SB <<<
137	REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
138	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! sea surface temperature, in K
139	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
140	REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils
141	REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
142
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
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
163	! Local variables
164	!*************************************************************************
165	INTEGER :: i, k
166	REAL :: tmp
167	REAL, PARAMETER :: cepdu2=(0.1)**2
168	REAL, DIMENSION(klon) :: alb_eau, z0_lim
169	REAL, DIMENSION(klon) :: radsol
170	REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
171	CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
172	real rhoa(knon) ! density of moist air (kg / m3)
173	REAL sens_prec_liq(knon)
174
175	REAL t_int(knon) ! ocean-air interface temperature, in K
176	real s_int(knon) ! ocean-air interface salinity, in ppt
177
178	!**************************************************************************
179
180	#ifdef ISO
181	#ifdef ISOVERIF
182	do i=1,knon
183	if (iso_eau.gt.0) then
184	call iso_verif_egalite_choix(xtspechum(iso_eau,i), &
185	& spechum(i),'surf_ocean_mod 117', &
186	& errmax,errmaxrel)
187	call iso_verif_egalite_choix(xtsnow(iso_eau,i), &
188	& snow(i),'surf_ocean_mod 127', &
189	& errmax,errmaxrel)
190	endif !if (iso_eau.gt.0) then
191	enddo !do i=1,klon
192	#endif
193	#endif
194
195	!******************************************************************************
196	! Calculate total net radiance at surface
197	!
198	!******************************************************************************
199	radsol(1:klon) = 0.0 ! initialisation a priori inutile
200	radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
201
202	!******************************************************************************
203	! Cdragq computed from cdrag
204	! The difference comes only from a factor (f_z0qh_oce) on z0, so that
205	! it can be computed inside surf_ocean
206	! More complicated appraches may require the propagation through
207	! pbl_surface of an independant cdragq variable.
208	!******************************************************************************
209
210	IF ( f_z0qh_oce .ne. 1.) THEN
211	! Si on suit les formulations par exemple de Tessel, on
212	! a z0h=0.4nu/u, z0q=0.62nu/u, d'ou f_z0qh_oce=0.62/0.4=1.55
213	cdragq(1:knon)=cdragh(1:knon)* &
214	log(z1lay(1:knon)/z0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*z0h(1:knon)))
215	ELSE
216	cdragq(1:knon)=cdragh(1:knon)
217	ENDIF
218
219	rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon)))
220	!******************************************************************************
221	! Switch according to type of ocean (couple, slab or forced)
222	!******************************************************************************
223	SELECT CASE(type_ocean)
224	CASE('couple')
225	CALL ocean_cpl_noice( &
226	swnet, lwnet, alb1, &
227	windsp, fder, &
228	itime, dtime, knon, knindex, &
229	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow,temp_air,spechum,&
230	AcoefH, AcoefQ, BcoefH, BcoefQ, &
231	AcoefU, AcoefV, BcoefU, BcoefV, &
232	ps, u1, v1, gustiness, tsurf_in, &
233	radsol, snow, agesno, &
234	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
235	tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, delta_sal, rhoa, &
236	delta_sst, dTer, dSer, dt_ds)
237
238	CASE('slab')
239	CALL ocean_slab_noice( &
240	itime, dtime, jour, knon, knindex, &
241	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum,&
242	AcoefH, AcoefQ, BcoefH, BcoefQ, &
243	AcoefU, AcoefV, BcoefU, BcoefV, &
244	ps, u1, v1, gustiness, tsurf_in, &
245	radsol, snow, &
246	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
247	tsurf_new, dflux_s, dflux_l, lmt_bils)
248
249	CASE('force')
250	CALL ocean_forced_noice( &
251	itime, dtime, jour, knon, knindex, &
252	p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
253	temp_air, spechum, &
254	AcoefH, AcoefQ, BcoefH, BcoefQ, &
255	AcoefU, AcoefV, BcoefU, BcoefV, &
256	ps, u1, v1, gustiness, tsurf_in, &
257	radsol, snow, agesno, &
258	qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
259	tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
260	#ifdef ISO
261	,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
262	xtsnow,xtevap,h1 &
263	#endif
264	)
265	END SELECT
266
267	!******************************************************************************
268	! fcodron: compute lmt_bils forced case (same as wfbils_oce / 1.-contfracatm)
269	!******************************************************************************
270	IF (type_ocean.NE.'slab') THEN
271	lmt_bils(1:klon)=0.
272	DO i=1,knon
273	lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
274	*pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i)))
275	END DO
276	END IF
277
278	!******************************************************************************
279	! Calculate ocean surface albedo
280	!******************************************************************************
281	!albedo SB >>>
282	IF (iflag_albedo==0) THEN
283	!--old parametrizations of ocean surface albedo
284	!
285	IF (iflag_cycle_diurne.GE.1) THEN
286	!
287	CALL alboc_cd(rmu0,alb_eau)
288	!
289	!--ad-hoc correction for model radiative balance tuning
290	!--now outside alboc_cd routine
291	alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
292	alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.0),1.0)
293	!
294	ELSE
295	!
296	CALL alboc(REAL(jour),rlat,alb_eau)
297	!--ad-hoc correction for model radiative balance tuning
298	!--now outside alboc routine
299	alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
300	alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.04),0.60)
301	!
302	ENDIF
303	!
304	DO i =1, knon
305	DO k=1,nsw
306	alb_dir_new(i,k) = alb_eau(knindex(i))
307	ENDDO
308	ENDDO
309	!IM 09122015 next line corresponds to the old way of doing in LMDZ5A/IPSLCM5A versions
310	!albedo for diffuse radiation is taken the same as for direct radiation
311	alb_dif_new(1:knon,:)=alb_dir_new(1:knon,:)
312	!IM 09122015 end
313	!
314	ELSE IF (iflag_albedo==1) THEN
315	!--new parametrization of ocean surface albedo by Sunghye Baek
316	!--albedo for direct and diffuse radiation are different
317	!
318	CALL ocean_albedo(knon,rmu0,knindex,windsp,SFRWL,alb_dir_new,alb_dif_new)
319	!
320	!--ad-hoc correction for model radiative balance tuning
321	alb_dir_new(1:knon,:) = fmagic*alb_dir_new(1:knon,:) + pmagic
322	alb_dif_new(1:knon,:) = fmagic*alb_dif_new(1:knon,:) + pmagic
323	alb_dir_new(1:knon,:)=MIN(MAX(alb_dir_new(1:knon,:),0.0),1.0)
324	alb_dif_new(1:knon,:)=MIN(MAX(alb_dif_new(1:knon,:),0.0),1.0)
325	!
326	ELSE IF (iflag_albedo==2) THEN
327	! F. Codron albedo read from limit.nc
328	CALL limit_read_rug_alb(itime, dtime, jour,&
329	knon, knindex, z0_lim, alb_eau)
330	DO i =1, knon
331	DO k=1,nsw
332	alb_dir_new(i,k) = alb_eau(i)
333	ENDDO
334	ENDDO
335	alb_dif_new=alb_dir_new
336	ENDIF
337	!albedo SB <<<
338
339	!******************************************************************************
340	! Calculate the rugosity
341	!******************************************************************************
342	IF (iflag_z0_oce==0) THEN
343	DO i = 1, knon
344	tmp = MAX(cepdu2,gustiness(i)+u1(i)2+v1(i)2)
345	z0m(i) = 0.018cdragm(i) (gustiness(i)+u1(i)2+v1(i)2)/RG &
346	+ 0.1114e-6 / SQRT(cdragm(i) tmp)
347	z0m(i) = MAX(1.5e-05,z0m(i))
348	ENDDO
349	z0h(1:knon)=z0m(1:knon) ! En attendant mieux
350
351	ELSE IF (iflag_z0_oce==1) THEN
352	DO i = 1, knon
353	tmp = MAX(cepdu2,gustiness(i)+u1(i)2+v1(i)2)
354	z0m(i) = 0.018cdragm(i) (gustiness(i)+u1(i)2+v1(i)2)/RG &
355	+ 0.1114e-6 / SQRT(cdragm(i) tmp)
356	z0m(i) = MAX(1.5e-05,z0m(i))
357	z0h(i)=0.414e-6 / SQRT(cdragm(i) tmp)
358	ENDDO
359	ELSE IF (iflag_z0_oce==-1) THEN
360	DO i = 1, knon
361	z0m(i) = z0min
362	z0h(i) = z0min
363	ENDDO
364	ELSE
365	CALL abort_physic(modname,'version non prevue',1)
366	ENDIF
367
368	if (activate_ocean_skin >= 1) then
369	if (type_ocean /= 'couple') sss(:knon) = 35.
370	call bulk_flux(tkt, tks, taur, dter, dser, t_int, s_int, ds_ns, dt_ns, &
371	u = windsp(:knon), t_ocean_1 = tsurf_new(:knon), s1 = sss(:knon), &
372	rain = precip_rain(:knon) + precip_snow(:knon), &
373	hf = - fluxsens(:knon), hlb = - fluxlat(:knon), &
374	rnl = - lwnet(:knon), &
375	tau = sqrt(flux_u1(:knon)2 + flux_v1(:knon)2), rhoa = rhoa, &
376	xlv = [(rlvtt, i = 1, knon)], rf = - sens_prec_liq, dtime = dtime, &
377	rns = swnet(:knon))
378	delta_sst = t_int - tsurf_new(:knon)
379	delta_sal = s_int - sss(:knon)
380
381	if (activate_ocean_skin == 2) then
382	tsurf_new(:knon) = t_int
383	if (type_ocean == 'couple') dt_ds = (tks / tkt) * dter
384	end if
385	end if
386
387	END SUBROUTINE surf_ocean
388	!****************************************************************************
389	!
390	END MODULE surf_ocean_mod

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

Download in other formats:

Original Format