Context Navigation

hydrol.F90 @ 815

Last change on this file since 815 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

Property svn:executable set to *

File size: 10.1 KB

Rev	Line
[787]	1	subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf, &
[253]	2	qsurf,dqsurf,dqs_hyd,pcapcal, &
	3	albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
	4
[650]	5	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
[787]	6	USE surfdat_h
	7	use comdiurn_h
	8	USE comgeomfi_h
	9	USE tracer_h
[253]	10
	11	implicit none
	12
	13	!==================================================================
	14	!
	15	! Purpose
	16	! -------
	17	! Calculate the surface hydrology and albedo changes.
	18	!
	19	! Authors
	20	! -------
	21	! Adapted from LMDTERRE by B. Charnay (2010). Further
	22	! modifications by R. Wordsworth (2010).
	23	!
	24	! Called by
	25	! ---------
	26	! physiq.F
	27	!
	28	! Calls
	29	! -----
	30	! none
	31	!
	32	! Notes
	33	! -----
	34	! rnat is terrain type: 0-ocean; 1-continent
	35	!
	36	!==================================================================
	37
	38	#include "dimensions.h"
	39	#include "dimphys.h"
	40	#include "comcstfi.h"
	41	#include "callkeys.h"
	42
[787]	43	integer ngrid,nq
	44
[253]	45	! Inputs
	46	! ------
	47	real albedoocean
	48	parameter (albedoocean=0.07)
	49	real albedoice
	50	save albedoice
	51
	52	real snowlayer
	53	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
	54	real oceantime
[305]	55	parameter (oceantime=10243600)
[253]	56
	57	logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
	58	logical activerunoff ! enable simple runoff scheme?
	59	logical oceanalbvary ! ocean albedo varies with the diurnal cycle?
	60
	61	! Arguments
	62	! ---------
[787]	63	integer rnat(ngrid) ! I changed this to integer (RW)
	64	real,dimension(:),allocatable,save :: runoff
[253]	65	real totalrunoff, tsea, oceanarea
	66	save oceanarea
	67
	68	real ptimestep
[787]	69	real mu0(ngrid)
	70	real qsurf(ngrid,nq), tsurf(ngrid)
	71	real dqsurf(ngrid,nq), pdtsurf(ngrid)
	72	real hice(ngrid)
	73	real albedo0(ngrid), albedo(ngrid)
[253]	74
	75	real oceanarea2
	76
	77	! Output
	78	! ------
[787]	79	real dqs_hyd(ngrid,nq)
	80	real pdtsurf_hyd(ngrid)
[253]	81
	82	! Local
	83	! -----
	84	real a,b,E
	85	integer ig,iq, icap ! wld like to remove icap
	86	real fsnoi, subli, fauxo
[787]	87	real twater(ngrid)
	88	real pcapcal(ngrid)
	89	real hicebis(ngrid)
	90	real zqsurf(ngrid,nq)
	91	real ztsurf(ngrid)
[253]	92
	93	integer ivap, iliq, iice
	94	save ivap, iliq, iice
	95
	96	logical firstcall
	97	save firstcall
	98
	99	data firstcall /.true./
	100
	101	oceanbulkavg=.false.
[368]	102	activerunoff=.false.
[253]	103	oceanalbvary=.false.
	104
	105	if(firstcall)then
	106
[787]	107	ALLOCATE(runoff(ngrid))
	108
[253]	109	ivap=igcm_h2o_vap
	110	iliq=igcm_h2o_vap
	111	iice=igcm_h2o_ice
	112
	113	write(,) "hydrol: ivap=",ivap
	114	write(,) " iliq=",iliq
	115	write(,) " iice=",iice
	116
	117	! Here's the deal: iice is used in place of igcm_h2o_ice both on the
	118	! surface and in the atmosphere. ivap is used in
	119	! place of igcm_h2o_vap ONLY in the atmosphere, while
	120	! iliq is used in place of igcm_h2o_vap ONLY on the
	121	! surface.
	122
	123	! Soon to be extended to the entire water cycle...
	124
	125	! Ice albedo = snow albedo for now
	126	albedoice=albedosnow
	127
	128	! Total ocean surface area
	129	oceanarea=0.
[787]	130	do ig=1,ngrid
[253]	131	if(rnat(ig).eq.0)then
	132	oceanarea=oceanarea+area(ig)
	133	endif
	134	enddo
	135
	136	if(oceanbulkavg.and.(oceanarea.le.0.))then
	137	print*,'How are we supposed to average the ocean'
	138	print*,'temperature, when there are no oceans?'
	139	call abort
	140	endif
	141
	142	if(activerunoff.and.(oceanarea.le.0.))then
	143	print*,'You have enabled runoff, but you have no oceans.'
	144	print*,'Where did you think the water was going to go?'
	145	call abort
	146	endif
	147
	148	firstcall = .false.
	149	endif
	150
	151	! add physical tendencies already calculated
	152	! ------------------------------------------
	153
[787]	154	do ig=1,ngrid
[253]	155	ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
	156	pdtsurf_hyd(ig)=0.0
[787]	157	do iq=1,nq
[253]	158	zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
	159	enddo
	160	enddo
	161
[787]	162	do ig=1,ngrid
	163	do iq=1,nq
[253]	164	dqs_hyd(ig,iq) = 0.0
	165	enddo
	166	enddo
	167
[787]	168	do ig = 1, ngrid
[253]	169
	170	! Ocean
	171	! -----
	172	if(rnat(ig).eq.0)then
	173
	174	! re-calculate oceanic albedo
	175	if(diurnal.and.oceanalbvary)then
	176	fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
	177	albedo(ig) = 1.1( .03 + .630/( 1. + fauxofauxo))
	178	albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
	179	else
	180	albedo(ig) = albedoocean ! modif Benjamin
	181	end if
	182
	183	! calculate oceanic ice height including the latent heat of ice formation
	184	! hice is the height of oceanic ice with a maximum of maxicethick.
	185	hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
	186
	187	! twater(ig) = tsurf(ig) + ptimestep*zdtsurf(ig) &
	188	twater(ig) = ztsurf(ig) - hice(ig)RLFTTrhowater/pcapcal(ig)
	189	! this is temperature water would have if we melted entire ocean ice layer
	190	hicebis(ig) = hice(ig)
	191	hice(ig) = 0.
	192
[650]	193	if(twater(ig) .lt. T_h2O_ice_liq)then
	194	E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))pcapcal(ig),RLFTTrhowater*maxicethick)
[253]	195	hice(ig) = E/(RLFTT*rhowater)
	196	hice(ig) = max(hice(ig),0.0)
	197	hice(ig) = min(hice(ig),maxicethick)
[773]	198	pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))RLFTTrhowater/pcapcal(ig)/ptimestep
[253]	199	albedo(ig) = albedoice
	200
	201	! if (zqsurf(ig,iice).ge.snowlayer) then
	202	! albedo(ig) = albedoice
	203	! else
	204	! albedo(ig) = albedoocean &
	205	! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
	206	! endif
	207
	208	else
	209
	210	pdtsurf_hyd(ig) = -hicebis(ig)RLFTTrhowater/pcapcal(ig)/ptimestep
	211	albedo(ig) = albedoocean
	212
	213	endif
	214
	215	zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
	216	zqsurf(ig,iice) = hice(ig)*rhowater
	217
	218
	219	! Continent
	220	! ---------
	221	elseif (rnat(ig).eq.1) then
	222
	223	! melt the snow
[650]	224	if(ztsurf(ig).gt.T_h2O_ice_liq)then
[253]	225	if(zqsurf(ig,iice).gt.1.0e-8)then
	226
[650]	227	a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
[253]	228	b = zqsurf(ig,iice)
	229	fsnoi = min(a,b)
	230
	231	zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi
	232	zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi
	233
	234	! thermal effects
	235	pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep
	236
	237	endif
	238	else
	239
	240	! freeze the water
	241	if(zqsurf(ig,iliq).gt.1.0e-8)then
	242
[650]	243	a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
[253]	244	b = zqsurf(ig,iliq)
	245
	246	fsnoi = min(a,b)
	247
	248	zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi
	249	zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi
	250
	251	! thermal effects
	252	pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep
	253
	254	endif
	255	endif
	256
	257	! deal with runoff
	258	if(activerunoff)then
	259
	260	runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
[787]	261	if(ngrid.gt.1)then ! runoff only exists in 3D
[253]	262	if(runoff(ig).ne.0.0)then
	263	zqsurf(ig,iliq) = mx_eau_sol
	264	! runoff is added to ocean at end
	265	endif
	266	end if
	267
	268	endif
	269
	270	! re-calculate continental albedo
	271	albedo(ig) = albedo0(ig) ! albedo0 = base values
	272	if (zqsurf(ig,iice).ge.snowlayer) then
	273	albedo(ig) = albedosnow
	274	else
	275	albedo(ig) = albedo0(ig) &
	276	+ (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer
	277	endif
	278
	279	else
	280
	281	print*,'Surface type not recognised in hydrol.F!'
	282	print*,'Exiting...'
	283	call abort
	284
	285	endif
	286
[787]	287	end do ! ig=1,ngrid
[253]	288
	289	! perform crude bulk averaging of temperature in ocean
	290	! ----------------------------------------------------
	291	if(oceanbulkavg)then
	292
	293	oceanarea2=0.
[787]	294	DO ig=1,ngrid
[253]	295	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
	296	oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
	297	end if
	298	END DO
	299
	300	tsea=0.
[787]	301	DO ig=1,ngrid
[253]	302	if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
	303	tsea=tsea+ztsurf(ig)area(ig)pcapcal(ig)/oceanarea2
	304	end if
	305	END DO
	306
[787]	307	DO ig=1,ngrid
[253]	308	if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
	309	pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
	310	end if
	311	END DO
	312
[305]	313	print*,'Mean ocean temperature = ',tsea,' K'
[253]	314
	315	endif
	316
	317	! shove all the runoff water into the ocean
	318	! -----------------------------------------
	319	if(activerunoff)then
	320
	321	totalrunoff=0.
[787]	322	do ig=1,ngrid
[253]	323	if (rnat(ig).eq.1) then
	324	totalrunoff = totalrunoff + area(ig)*runoff(ig)
	325	endif
	326	enddo
	327
[787]	328	do ig=1,ngrid
[253]	329	if (rnat(ig).eq.0) then
	330	zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
	331	totalrunoff/oceanarea
	332	endif
	333	enddo
	334
	335	endif
	336
	337	! Re-add the albedo effects of CO2 ice if necessary
	338	! -------------------------------------------------
	339	if(co2cond)then
	340
	341	icap=1
[787]	342	do ig=1,ngrid
[253]	343	if (qsurf(ig,igcm_co2_ice).gt.0) then
	344	albedo(ig) = albedice(icap)
	345	endif
	346	enddo
	347
	348	endif
	349
	350
[787]	351	do ig=1,ngrid
[253]	352	dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep
	353	dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep
	354	enddo
	355
[787]	356	call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
[253]	357
	358	return
	359	end subroutine hydrol

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

Download in other formats:

Original Format