Context Navigation

hydrol.F90 @ 3103

Last change on this file since 3103 was 3100, checked in by bhatnags, 2 years ago

Generic-PCM

This commit updates the slab ocean module to a parallelisable dynamic slab ocean module. This is particularly relevant if you want to be able to use oceanic heat transport in parallel mode.

It has the following features:

(a) Computes sea ice creation and evolution.
(b) Snow has thermodynamic properties.
(c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
(d) Can be used in parallel mode.

Required callphys.def flags:
The slab ocean and its dependencies can be activated with the following flags (already added to deftank):
## Ocean options
## ~
# Model slab-ocean (Main flag for slab ocean)
ok_slab_ocean = .true.
# The following flags can only be set to true if ok_slab_ocean is true
# Ekman transport
slab_ekman = .true.
# Ekman zonal advection
slab_ekman_zonadv = .true.
# Horizontal diffusion (default coef_hdiff=25000., can be changed)
slab_hdiff = .true.
# Slab-ocean timestep (in physics timesteps)
cpl_pas = 1
# Gent-McWilliams? Scheme (can only be true if slab_ekman is true)
slab_gm = .true.

Notes:
In the current state, the model crashes if moistadjustment = .true. Unsure whether this is due to the slab or is an inherent issue with moistadj (under investigation).

SB and EM

Property svn:executable set to *

File size: 14.8 KB

Line
1	subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf, &
2	qsurf,dqsurf,dqs_hyd,pcapcal, &
3	albedo,albedo_bareground, &
4	albedo_snow_SPECTV,albedo_co2_ice_SPECTV, &
5	mu0,pdtsurf,pdtsurf_hyd,hice, &
6	pctsrf_sic,sea_ice)
7
8	use ioipsl_getin_p_mod, only: getin_p
9	#ifndef MESOSCALE
10	use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
11	use mod_phys_lmdz_para, only : is_master, gather, scatter
12	#endif
13	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
14	USE surfdat_h
15	use comdiurn_h
16	USE geometry_mod, only: cell_area
17	USE tracer_h
18	! use slab_ice_h
19	USE ocean_slab_mod, ONLY: alb_ice_max
20	USE ocean_slab_mod, ONLY: alb_ice_min
21	USE ocean_slab_mod, ONLY: h_alb_ice
22	use callkeys_mod, only: albedosnow,alb_ocean,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
23	use radinc_h, only : L_NSPECTV
24
25	implicit none
26
27	!==================================================================
28	!
29	! Purpose
30	! -------
31	! Calculate the surface hydrology and albedo changes.
32	!
33	! Authors
34	! -------
35	! Adapted from LMDTERRE by B. Charnay (2010). Further
36	! Modifications by R. Wordsworth (2010).
37	! Spectral albedo by M. Turbet (2015).
38	!
39	! Called by
40	! ---------
41	! physiq.F
42	!
43	! Calls
44	! -----
45	! none
46	!
47	! Notes
48	! -----
49	! rnat is terrain type: 0-ocean; 1-continent
50	!
51	!==================================================================
52
53	integer ngrid,nq
54
55	! Inputs
56	! ------
57	real snowlayer
58	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
59	real oceantime
60	parameter (oceantime=10243600)
61
62	logical,save :: oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
63	logical,save :: activerunoff ! enable simple runoff scheme?
64	logical,save :: oceanalbvary ! ocean albedo varies with the diurnal cycle?
65	!$OMP THREADPRIVATE(oceanbulkavg,activerunoff,oceanalbvary)
66
67	! Arguments
68	! ---------
69	real rnat(ngrid) ! I changed this to integer (RW)
70	real,dimension(:),allocatable,save :: runoff
71	real totalrunoff, tsea, oceanarea
72	save oceanarea
73	!$OMP THREADPRIVATE(runoff,oceanarea)
74
75	real ptimestep
76	real mu0(ngrid)
77	real qsurf(ngrid,nq), tsurf(ngrid)
78	real dqsurf(ngrid,nq), pdtsurf(ngrid)
79	real hice(ngrid)
80	real albedo(ngrid,L_NSPECTV)
81	real albedo_bareground(ngrid)
82	real albedo_snow_SPECTV(L_NSPECTV)
83	real albedo_co2_ice_SPECTV(L_NSPECTV)
84	real pctsrf_sic(ngrid), sea_ice(ngrid)
85
86	real oceanarea2
87
88	! Output
89	! ------
90	real dqs_hyd(ngrid,nq)
91	real pdtsurf_hyd(ngrid)
92
93	! Local
94	! -----
95	real a,b,E
96	integer ig,iq, nw
97	real fsnoi, subli, fauxo
98	real twater(ngrid)
99	real pcapcal(ngrid)
100	real hicebis(ngrid)
101	real zqsurf(ngrid,nq)
102	real ztsurf(ngrid)
103	real albedo_sic, alb_ice
104	real zfra
105
106	integer, save :: ivap, iliq, iice
107	!$OMP THREADPRIVATE(ivap,iliq,iice)
108
109	logical, save :: firstcall=.true.
110	!$OMP THREADPRIVATE(firstcall)
111
112	real :: runoffamount(ngrid)
113	!#ifdef CPP_PARA
114	real :: runoffamount_glo(klon_glo)
115	real :: zqsurf_iliq_glo(klon_glo)
116	real :: rnat_glo(klon_glo)
117	real :: oceanarea_glo
118	real :: cell_area_glo(klon_glo)
119	!#else
120	! real :: runoffamount_glo(ngrid)
121	! real :: zqsurf_iliq_glo(ngrid)
122	!#endif
123
124
125	if(firstcall)then
126
127	oceanbulkavg=.false.
128	oceanalbvary=.false.
129	write(,)"Activate runnoff into oceans?"
130	activerunoff=.false.
131	call getin_p("activerunoff",activerunoff)
132	write(,)" activerunoff = ",activerunoff
133
134
135
136	if (activerunoff) then
137	ALLOCATE(runoff(ngrid))
138	runoff(1:ngrid)=0
139	endif
140
141	ivap=igcm_h2o_vap
142	iliq=igcm_h2o_vap
143	iice=igcm_h2o_ice
144
145	write(,) "hydrol: ivap=",ivap
146	write(,) " iliq=",iliq
147	write(,) " iice=",iice
148
149	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
150	! surface and in the atmosphere. ivap is used in
151	! place of igcm_h2o_vap ONLY in the atmosphere, while
152	! iliq is used in place of igcm_h2o_vap ONLY on the
153	! surface.
154	! Soon to be extended to the entire water cycle...
155
156	! LOCAL ocean surface area
157	oceanarea=0.
158	do ig=1,ngrid
159	if(nint(rnat(ig)).eq.0)then
160	oceanarea=oceanarea+cell_area(ig)
161	endif
162	enddo
163
164	if(oceanbulkavg.and.(oceanarea.le.0.))then
165	print*,'How are we supposed to average the ocean'
166	print*,'temperature, when there are no oceans?'
167	call abort
168	endif
169
170	if(activerunoff.and.(oceanarea.le.0.))then
171	print*,'You have enabled runoff, but you have no oceans.'
172	print*,'Where did you think the water was going to go?'
173	call abort
174	endif
175
176	firstcall = .false.
177	endif
178
179	! call writediagfi(ngrid,"hydrol_rnat"," "," ",2,rnat)
180	!! write (,) "oceanarea", oceanarea
181
182	! add physical tendencies already calculated
183	! ------------------------------------------
184
185	do ig=1,ngrid
186	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
187	pdtsurf_hyd(ig)=0.0
188	do iq=1,nq
189	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
190	enddo
191	enddo
192
193	! call writediagfi(ngrid,"hydrol_zqsurf_iliq_1"," "," ",2,zqsurf(1:ngrid,iliq))
194
195	do ig=1,ngrid
196	do iq=1,nq
197	dqs_hyd(ig,iq) = 0.0
198	enddo
199	enddo
200
201	do ig = 1, ngrid
202
203	! Ocean
204	! -----
205	if(nint(rnat(ig)).eq.0)then
206
207	! re-calculate oceanic albedo
208	! if(diurnal.and.oceanalbvary)then
209	! fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
210	! albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
211	! albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
212	! else
213	do nw=1,L_NSPECTV
214	albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in inifis_mod.F90. Later we could introduce spectral dependency for alb_ocean.
215	enddo
216	! end if
217
218
219	if(ok_slab_ocean) then
220
221	zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0)) ! Snow Fraction (Critical height 45kg/m2~15cm)
222	alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) & ! Ice Albedo
223	*exp(-sea_ice(ig)/h_alb_ice)
224	! Albedo final calculation :
225	do nw=1,L_NSPECTV
226	albedo(ig,nw) = pctsrf_sic(ig)* &
227	(albedo_snow_SPECTV(nw)zfra + alb_ice(1.0-zfra)) &
228	+ (1.-pctsrf_sic(ig))*alb_ocean
229	enddo
230
231	! Oceanic ice height, just for diagnostics
232	hice(ig) = MIN(10.,sea_ice(ig)/rhowater)
233	else !ok_slab_ocean
234
235
236	! calculate oceanic ice height including the latent heat of ice formation
237	! hice is the height of oceanic ice with a maximum of maxicethick.
238	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
239
240	! twater(ig) = tsurf(ig) + ptimestep*zdtsurf(ig) &
241	twater(ig) = ztsurf(ig) - hice(ig)RLFTTrhowater/pcapcal(ig)
242	! this is temperature water would have if we melted entire ocean ice layer
243	hicebis(ig) = hice(ig)
244	hice(ig) = 0.
245
246	if(twater(ig) .lt. T_h2O_ice_liq)then
247	E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))pcapcal(ig),RLFTTrhowater*maxicethick)
248	hice(ig) = E/(RLFTT*rhowater)
249	hice(ig) = max(hice(ig),0.0)
250	hice(ig) = min(hice(ig),maxicethick)
251	pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))RLFTTrhowater/pcapcal(ig)/ptimestep
252	do nw=1,L_NSPECTV
253	albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
254	enddo
255
256	! if (zqsurf(ig,iice).ge.snowlayer) then
257	! albedo(ig) = albedoice
258	! else
259	! albedo(ig) = albedoocean &
260	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
261	! endif
262
263	else
264
265	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
266	DO nw=1,L_NSPECTV
267	albedo(ig,nw) = alb_ocean
268	ENDDO
269
270	endif
271
272	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
273	zqsurf(ig,iice) = hice(ig)*rhowater
274
275	endif!(ok_slab_ocean)
276
277
278	! Continent
279	! ---------
280	elseif (nint(rnat(ig)).eq.1) then
281
282	! melt the snow
283	if(ztsurf(ig).gt.T_h2O_ice_liq)then
284	if(zqsurf(ig,iice).gt.1.0e-8)then
285
286	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
287	b = zqsurf(ig,iice)
288	fsnoi = min(a,b)
289
290	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
291	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
292
293	! thermal effects
294	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
295
296	endif
297	else
298
299	! freeze the water
300	if(zqsurf(ig,iliq).gt.1.0e-8)then
301
302	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
303	b = zqsurf(ig,iliq)
304
305	fsnoi = min(a,b)
306
307	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
308	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
309
310	! thermal effects
311	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
312
313	endif
314	endif
315
316	! deal with runoff
317	if(activerunoff)then
318	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
319	if(ngrid.gt.1)then ! runoff only exists in 3D
320	if(runoff(ig).ne.0.0)then
321	zqsurf(ig,iliq) = mx_eau_sol
322	! runoff is added to ocean at end
323	endif
324	end if
325
326	endif
327
328	! re-calculate continental albedo
329	DO nw=1,L_NSPECTV
330	albedo(ig,nw) = albedo_bareground(ig)
331	ENDDO
332	if (zqsurf(ig,iice).ge.snowlayer) then
333	DO nw=1,L_NSPECTV
334	albedo(ig,nw) = albedo_snow_SPECTV(nw)
335	ENDDO
336	else
337	DO nw=1,L_NSPECTV
338	albedo(ig,nw) = albedo_bareground(ig) &
339	+ (albedo_snow_SPECTV(nw) - albedo_bareground(ig)) &
340	*zqsurf(ig,iice)/snowlayer
341	ENDDO
342	endif
343
344	else
345
346	print*,'Surface type not recognised in hydrol.F!'
347	print*,'Exiting...'
348	call abort
349
350	endif
351
352	end do ! ig=1,ngrid
353
354	! call writediagfi(ngrid,"hydrol_zqsurf_iliq_2"," "," ",2,zqsurf(1:ngrid,iliq))
355
356	! perform crude bulk averaging of temperature in ocean
357	! ----------------------------------------------------
358	if(oceanbulkavg)then
359
360	oceanarea2=0.
361	DO ig=1,ngrid
362	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
363	oceanarea2=oceanarea2+cell_area(ig)*pcapcal(ig)
364	end if
365	END DO
366
367	tsea=0.
368	DO ig=1,ngrid
369	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
370	tsea=tsea+ztsurf(ig)cell_area(ig)pcapcal(ig)/oceanarea2
371	end if
372	END DO
373
374	DO ig=1,ngrid
375	if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then
376	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
377	end if
378	END DO
379
380	print*,'Mean ocean temperature = ',tsea,' K'
381
382	endif
383
384	! shove all the runoff water into the ocean
385	! -----------------------------------------
386	if(activerunoff)then
387
388	! totalrunoff=0.
389	do ig=1,ngrid
390	runoffamount(ig) = cell_area(ig)*runoff(ig)
391	! if (nint(rnat(ig)).eq.1) then
392	! totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
393	! endif
394	enddo
395
396	! collect on the full grid
397	call gather(runoffamount,runoffamount_glo)
398	call gather(zqsurf(1:ngrid,iliq),zqsurf_iliq_glo)
399	call gather(rnat,rnat_glo)
400	call gather(cell_area,cell_area_glo)
401
402	if (is_master) then
403	totalrunoff=0.
404	oceanarea_glo=0.
405	do ig=1,klon_glo
406	if (nint(rnat_glo(ig)).eq.1) then
407	totalrunoff = totalrunoff + runoffamount_glo(ig)
408	endif
409	if (nint(rnat_glo(ig)).eq.0) then
410	oceanarea_glo = oceanarea_glo + cell_area_glo(ig)
411	endif
412	enddo
413
414	do ig=1,klon_glo
415	if (nint(rnat_glo(ig)).eq.0) then
416	zqsurf_iliq_glo(ig) = zqsurf_iliq_glo(ig) + &
417	totalrunoff/oceanarea_glo
418	endif
419	enddo
420
421	endif! is_master
422
423	! scatter the field back on all processes
424	call scatter(zqsurf_iliq_glo,zqsurf(1:ngrid,iliq))
425
426
427
428	! do ig=1,ngrid
429	! if (nint(rnat(ig)).eq.0) then
430	! zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
431	! totalrunoff/oceanarea
432	! endif
433	! enddo
434
435	endif !activerunoff
436
437	! call writediagfi(ngrid,"hydrol_zqsurf_iliq_3"," "," ",2,zqsurf(1:ngrid,iliq))
438
439	! Re-add the albedo effects of CO2 ice if necessary
440	! -------------------------------------------------
441	if(co2cond)then
442
443	do ig=1,ngrid
444	if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
445	DO nw=1,L_NSPECTV
446	albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
447	ENDDO
448	endif
449	enddo ! ngrid
450
451	endif ! co2cond
452
453
454	do ig=1,ngrid ! We calculate here the tracer tendencies. Don't forget that we have to retrieve the dqsurf tendencies we added at the beginning of the routine !
455	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep - dqsurf(ig,iliq)
456	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep - dqsurf(ig,iice)
457	enddo
458
459	! call writediagfi(ngrid,"hydrol_dqs_hyd_iliq"," "," ",2,dqs_hyd(1:ngrid,iliq))
460
461	if (activerunoff) then
462	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
463	endif
464
465	return
466	end subroutine hydrol

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

Download in other formats:

Original Format