Context Navigation

phyetat0_mod.F90 @ 3026

Last change on this file since 3026 was 2999, checked in by llange, 18 months ago

Mars PCM
Include perenial_co2ice (equivalent of watercap) to distinguich between CO2 frost and perenial CO2 ice for paleoclimate studies.
When no frost is present and we dig into perenial ice, the surface albedo is changed. The albedo for seasonal ice is set to 0.65, and the perenial ice albedo can be fixed in the callphys.def. I recommand values between 0.8 and 0.9.
To use this, paleoclimate must be set to True and TESalbedo to false in the callphys.def. Else, it runs as usual with TES albedo
LL

File size: 26.7 KB

Rev	Line
[1944]	1	module phyetat0_mod
	2
	3	implicit none
[2545]	4	real,save :: tab_cntrl_mod(100)
[1944]	5
[2578]	6	!$OMP THREADPRIVATE(tab_cntrl_mod)
	7
[1944]	8	contains
	9
[1130]	10	subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
[2826]	11	day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf, &
[2999]	12	tauscaling,totcloudfrac,wstar,watercap,perenial_co2ice, &
	13	def_slope,def_slope_mean,subslope_dist)
[1922]	14
[1621]	15	use tracer_mod, only: noms ! tracer names
[1130]	16	use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
[2884]	17	zmea, zstd, zsig, zgam, zthe, hmons, summit, base,&
	18	watercaptag
[1130]	19	use iostart, only: nid_start, open_startphy, close_startphy, &
	20	get_field, get_var, inquire_field, &
	21	inquire_dimension, inquire_dimension_length
[2220]	22	use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
[2265]	23	use compute_dtau_mod, only: dtau
[2417]	24	use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
	25	use dust_param_mod, only: dustscaling_mode
[2281]	26	USE ioipsl_getin_p_mod, ONLY : getin_p
[2887]	27	use comsoil_h, only: flux_geo
[2896]	28	USE comslope_mod, ONLY: nslope, major_slope
[2994]	29	USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice
[2999]	30	USE comcstfi_h, only: pi
	31	use geometry_mod, only: latitude
[1130]	32	implicit none
[1974]	33
	34	include "callkeys.h"
[1130]	35	!======================================================================
	36	! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
[2562]	37	! Adaptation � Mars : Yann Wanherdrick
[1130]	38	! Objet: Lecture de l etat initial pour la physique
	39	! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
	40	! r4 or r8 restarts independently of having compiled
	41	! the GCM in r4 or r8)
	42	! June 2013 TN : Possibility to read files with a time axis
	43	! November 2013 EM : Enabeling parallel, using iostart module
[2281]	44	! March 2020 AD: Enabling initialization of physics without startfi
	45	! flag: startphy_file
[1130]	46	!======================================================================
	47	INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
	48	PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
	49	!======================================================================
	50	! Arguments:
	51	! ---------
	52	! inputs:
[2281]	53	! logical,intent(in) :: startphy_file ! .true. if reading start file
[1130]	54	character(),intent(in) :: fichnom ! "startfi.nc" file
	55	integer,intent(in) :: tab0
	56	integer,intent(in) :: Lmodif
	57	integer,intent(in) :: nsoil ! # of soil layers
	58	integer,intent(in) :: ngrid ! # of atmospheric columns
	59	integer,intent(in) :: nlay ! # of atmospheric layers
	60	integer,intent(in) :: nq
	61	integer :: day_ini
	62	real :: time0
	63
	64	! outputs:
[2913]	65	real,intent(out) :: tsurf(ngrid,nslope) ! surface temperature
	66	real,intent(out) :: tsoil(ngrid,nsoil,nslope) ! soil temperature
	67	real,intent(out) :: albedo(ngrid,2,nslope) ! surface albedo
	68	real,intent(out) :: emis(ngrid,nslope) ! surface emissivity
[1130]	69	real,intent(out) :: q2(ngrid,nlay+1) !
[2913]	70	real,intent(out) :: qsurf(ngrid,nq,nslope) ! tracers on surface
[1208]	71	real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
[1711]	72	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
[1944]	73	real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
[2913]	74	real,intent(out) :: watercap(ngrid,nslope) ! h2o_ice_cover
[2999]	75	real,intent(out) :: perenial_co2ice(ngrid,nslope) ! perenial co2 ice(kg/m^2)
	76	real,intent(out) :: def_slope(nslope+1) !boundaries for bining of the slopes
	77	real,intent(out) :: def_slope_mean(nslope)
	78	real,intent(out) :: subslope_dist(ngrid,nslope) !undermesh statistics
[1130]	79	!======================================================================
	80	! Local variables:
	81
	82	real surffield(ngrid) ! to temporarily store a surface field
	83	real xmin,xmax ! to display min and max of a field
	84	!
[2896]	85	INTEGER ig,iq,lmax,islope
[1130]	86	INTEGER nid, nvarid
	87	INTEGER ierr, i, nsrf
	88	! integer isoil
	89	! INTEGER length
	90	! PARAMETER (length=100)
	91	CHARACTER*7 str7
	92	CHARACTER*2 str2
	93	CHARACTER*1 yes
	94	!
	95	REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
	96	INTEGER nqold
	97
	98	! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
	99	logical :: oldtracernames=.false.
	100	integer :: count
	101	character(len=30) :: txt ! to store some text
	102
	103	! specific for time
	104	REAL,ALLOCATABLE :: time(:) ! times stored in start
	105	INTEGER timelen ! number of times stored in the file
	106	INTEGER indextime ! index of selected time
	107
	108	INTEGER :: edges(3),corner(3)
	109	LOGICAL :: found
	110
[1525]	111	REAL :: timestart ! to pick which initial state to start from
[2281]	112	REAL :: surfemis ! constant emissivity when no startfi
[2630]	113	REAL :: surfalbedo ! constant albedo when no startfi
[2884]	114
	115	REAL :: watercaptag_tmp(ngrid)
	116
[2896]	117	! Sub-grid scale slopes
	118	LOGICAL :: startphy_slope !to be retrocompatible and add the nslope dimension
	119	REAL, ALLOCATABLE :: default_def_slope(:)
	120	REAL :: sum_dist
	121	REAL :: current_max !var to find max distrib slope
[2999]	122	! Variables for CO2 index
	123	INTEGER :: igcm_co2_tmp
[2896]	124
[2281]	125	CHARACTER(len=5) :: modname="phyetat0"
[1525]	126
[2281]	127	write(,) "phyetat0: startphy_file", startphy_file
[1130]	128
[2281]	129	if (startphy_file) then
	130	! open physics initial state file:
	131	call open_startphy(fichnom)
	132	! possibility to modify tab_cntrl in tabfi
	133	write(,)
	134	write(,) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
	135	call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
	136	p_omeg,p_g,p_mugaz,p_daysec,time0)
	137	else ! "academic" initialization of planetary parameters via tabfi
	138	call tabfi (0,0,0,day_ini,lmax,p_rad, &
	139	p_omeg,p_g,p_mugaz,p_daysec,time0)
	140	endif ! of if (startphy_file)
[1130]	141
[2281]	142	if (startphy_file) then
[2896]	143	call get_var("def_slope",def_slope,found)
	144	if(.not.found) then
	145	startphy_slope=.false.
	146	write(,)'slope_settings: Problem while reading <def_slope>'
[2909]	147	write(,)'Checking that nslope=1'
	148	if(nslope.eq.1) then
	149	write(,)'We take default def_slope and subslope dist'
	150	allocate(default_def_slope(nslope+1))
[2913]	151	default_def_slope(1) = -90.
	152	default_def_slope(2) = 90.
[2909]	153	subslope_dist(:,:)=1.
[2913]	154	def_slope(:)=default_def_slope(:)
[2909]	155	else
	156	write(,)'Variable def_slope is not present in the start and'
	157	write(,)'you are trying to run with nslope!=1'
	158	write(,)'This is not possible, you have to go through newstart'
	159	endif
[2896]	160	else
	161	startphy_slope=.true.
	162	call get_field("subslope_dist",subslope_dist,found,indextime)
	163	if(.not.found) then
	164	write(,)'slope_settings: Problem while reading <subslope_dist>'
[2909]	165	write(,)'We have to abort.'
	166	write(,)'Please check that nslope is coherent, you can modify this parameter with newstart'
	167	call abort_physic(modname, &
	168	"phyetat0: Failed loading <subslope_dist>",1)
[2896]	169	endif
	170	endif
[2909]	171	else ! (no startphy_file)
	172	if(nslope.eq.1) then
	173	allocate(default_def_slope(nslope+1))
	174	default_def_slope(1) = 0.
	175	default_def_slope(2) = 0.
	176	subslope_dist(:,:)=1.
[2913]	177	def_slope(:)=default_def_slope(:)
[2909]	178	else
	179	write(,)'Without starfi, lets first run with nslope=1'
	180	call abort_physic(modname, &
	181	"phyetat0: No startfi and nslope!=1",1)
	182	endif
[2896]	183	endif
	184
	185	do islope=1,nslope
	186	def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
	187	enddo
	188
	189	DO ig = 1,ngrid
	190	sum_dist = 0.
	191	DO islope = 1,nslope
	192	sum_dist = sum_dist + subslope_dist(ig,islope)
	193	ENDDO
	194	DO islope = 1,nslope
	195	subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
	196	ENDDO
	197	ENDDO
	198
	199	!Now determine the major subslope, ie. the maximal distribution
	200
	201	DO ig=1,ngrid
	202	major_slope(ig)=1
	203	current_max=subslope_dist(ig,1)
	204	DO islope=2,nslope
	205	if(subslope_dist(ig,islope).GT.current_max) then
	206	major_slope(ig)=islope
	207	current_max=subslope_dist(ig,islope)
	208	ENDIF
	209	ENDDO
	210	ENDDO
	211
	212	if (startphy_file) then
[2281]	213	! Load surface geopotential:
	214	call get_field("phisfi",phisfi,found)
	215	if (.not.found) then
	216	call abort_physic(modname, &
	217	"phyetat0: Failed loading <phisfi>",1)
	218	endif
[1130]	219	else
[2281]	220	phisfi(:)=0.
	221	endif ! of if (startphy_file)
	222	write(,) "phyetat0: surface geopotential <phisfi> range:", &
	223	minval(phisfi), maxval(phisfi)
[1130]	224
	225
[2281]	226	if (startphy_file) then
	227	! Load bare ground albedo:
	228	call get_field("albedodat",albedodat,found)
	229	if (.not.found) then
	230	call abort_physic(modname, &
	231	"phyetat0: Failed loading <albedodat>",1)
	232	endif
	233	else ! If no startfi file, use parameter surfalbedo in def file
[2294]	234	surfalbedo=0.1
	235	call getin_p("surfalbedo_without_startfi",surfalbedo)
	236	print*,"surfalbedo_without_startfi",surfalbedo
[2281]	237	albedodat(:)=surfalbedo
	238	endif ! of if (startphy_file)
	239	write(,) "phyetat0: Bare ground albedo <albedodat> range:", &
[1130]	240	minval(albedodat), maxval(albedodat)
	241
	242	! ZMEA
[2281]	243	if (startphy_file) then
	244	call get_field("ZMEA",zmea,found)
	245	if (.not.found) then
	246	call abort_physic(modname, &
	247	"phyetat0: Failed loading <ZMEA>",1)
	248	endif
[1130]	249	else
[2281]	250	zmea(:)=0.
	251	endif ! of if (startphy_file)
	252	write(,) "phyetat0: <ZMEA> range:", &
	253	minval(zmea), maxval(zmea)
[1130]	254
	255
	256	! ZSTD
[2281]	257	if (startphy_file) then
	258	call get_field("ZSTD",zstd,found)
	259	if (.not.found) then
	260	call abort_physic(modname, &
	261	"phyetat0: Failed loading <ZSTD>",1)
	262	endif
[1130]	263	else
[2281]	264	zstd(:)=0.
	265	endif ! of if (startphy_file)
	266	write(,) "phyetat0: <ZSTD> range:", &
	267	minval(zstd), maxval(zstd)
[1130]	268
	269
	270	! ZSIG
[2281]	271	if (startphy_file) then
	272	call get_field("ZSIG",zsig,found)
	273	if (.not.found) then
	274	call abort_physic(modname, &
	275	"phyetat0: Failed loading <ZSIG>",1)
	276	endif
[1130]	277	else
[2281]	278	zsig(:)=0.
	279	endif ! of if (startphy_file)
	280	write(,) "phyetat0: <ZSIG> range:", &
	281	minval(zsig), maxval(zsig)
[1130]	282
	283
	284	! ZGAM
[2281]	285	if (startphy_file) then
	286	call get_field("ZGAM",zgam,found)
	287	if (.not.found) then
	288	call abort_physic(modname, &
	289	"phyetat0: Failed loading <ZGAM>",1)
	290	endif
[1130]	291	else
[2281]	292	zgam(:)=0.
	293	endif ! of if (startphy_file)
	294	write(,) "phyetat0: <ZGAM> range:", &
	295	minval(zgam), maxval(zgam)
[1130]	296
	297
	298	! ZTHE
[2281]	299	if (startphy_file) then
	300	call get_field("ZTHE",zthe,found)
	301	if (.not.found) then
	302	call abort_physic(modname, &
	303	"phyetat0: Failed loading <ZTHE>",1)
	304	endif
[1130]	305	else
[2281]	306	zthe(:)=0.
	307	endif ! of if (startphy_file)
	308	write(,) "phyetat0: <ZTHE> range:", &
[1130]	309	minval(zthe), maxval(zthe)
	310
[2281]	311	! hmons
	312	if (startphy_file) then
	313	call get_field("hmons",hmons,found)
	314	if (.not.found) then
	315	write(,) "WARNING: phyetat0: Failed loading <hmons>"
	316	if (rdstorm) then
	317	call abort_physic(modname, &
	318	"phyetat0: Failed loading <hmons>",1)
	319	else
	320	write(,) "will continue anyway..."
	321	write(,) "because you may not need it."
	322	hmons(:)=0.
	323	end if ! if (rdstorm)
	324	else
	325	do ig=1,ngrid
	326	if (hmons(ig) .eq. -999999.) hmons(ig)=0.
	327	enddo
	328	endif ! (.not.found)
[1974]	329	else
[2281]	330	hmons(:)=0.
	331	endif ! if (startphy_file)
	332	write(,) "phyetat0: <hmons> range:", &
	333	minval(hmons), maxval(hmons)
[2079]	334
[2281]	335
	336	! summit
	337	if (startphy_file) then
	338	call get_field("summit",summit,found)
	339	if (.not.found) then
	340	write(,) "WARNING: phyetat0: Failed loading <summit>"
	341	if (rdstorm) then
	342	call abort_physic(modname, &
	343	"phyetat0: Failed loading <summit>",1)
	344	else
	345	write(,) "will continue anyway..."
	346	write(,) "because you may not need it."
	347	summit(:)=0.
	348	end if
	349	else
	350	do ig=1,ngrid
	351	if (summit(ig) .eq. -999999.) summit(ig)=0.
	352	enddo
	353	endif ! if (.not.found)
[2079]	354	else
[2281]	355	summit(:)=0.
	356	endif ! if (startphy_file)
	357	write(,) "phyetat0: <summit> range:", &
	358	minval(summit), maxval(summit)
[2079]	359
[2281]	360
	361	! base
	362	if (startphy_file) then
	363	call get_field("base",base,found)
	364	if (.not.found) then
	365	write(,) "WARNING: phyetat0: Failed loading <base>"
	366	if (rdstorm) then
	367	call abort_physic(modname, &
	368	"phyetat0: Failed loading <base>",1)
	369	else
	370	write(,) "will continue anyway..."
	371	write(,) "because you may not need it."
	372	base(:)=0.
	373	end if
	374	else
	375	do ig=1,ngrid
	376	if (base(ig) .eq. -999999.) base(ig)=0.
	377	enddo
	378	endif ! if(.not.found)
[2079]	379	else
[2281]	380	base(:)=0.
	381	endif ! if (startphy_file)
	382	write(,) "phyetat0: <base> range:", &
	383	minval(base), maxval(base)
	384
[1130]	385	! Time axis
[1525]	386	! obtain timestart from run.def
	387	timestart=-9999 ! default value
[2281]	388	call getin_p("timestart",timestart)
	389	if (startphy_file) then
	390	found=inquire_dimension("Time")
	391	if (.not.found) then
	392	indextime = 1
	393	write(,) "phyetat0: No time axis found in "//trim(fichnom)
	394	else
	395	write(,) "phyetat0: Time axis found in "//trim(fichnom)
	396	timelen=inquire_dimension_length("Time")
	397	allocate(time(timelen))
	398	! load "Time" array:
	399	call get_var("Time",time,found)
	400	if (.not.found) then
	401	call abort_physic(modname, &
	402	"phyetat0: Failed loading <Time>",1)
	403	endif
	404	! seclect the desired time index
	405	IF (timestart .lt. 0) THEN ! default: we use the last time value
	406	indextime = timelen
	407	ELSE ! else we look for the desired value in the time axis
	408	indextime = 0
	409	DO i=1,timelen
	410	IF (abs(time(i) - timestart) .lt. 0.01) THEN
	411	indextime = i
	412	EXIT
	413	ENDIF
	414	ENDDO
	415	IF (indextime .eq. 0) THEN
	416	PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
	417	PRINT*, "Stored times are:"
	418	DO i=1,timelen
	419	PRINT*, time(i)
	420	ENDDO
	421	call abort_physic(modname,"phyetat0: Time error",1)
	422	ENDIF
	423	ENDIF ! of IF (timestart .lt. 0)
	424	! In startfi the absolute date is day_ini + time0 + time
	425	! For now on, in the GCM physics, it is day_ini + time0
	426	time0 = time(indextime) + time0
	427	day_ini = day_ini + INT(time0)
	428	time0 = time0 - INT(time0)
	429	PRINT*, "phyetat0: Selected time ",time(indextime), &
	430	" at index ",indextime
	431	DEALLOCATE(time)
	432	endif ! of if Time not found in file
[2545]	433
	434	call ini_tab_controle_dyn_xios(day_ini)
	435
[2281]	436	else
[1130]	437	indextime = 1
[2281]	438	endif ! if (startphy_file)
[1130]	439
[1208]	440	! Dust conversion factor
[2281]	441	if (startphy_file) then
	442	call get_field("tauscaling",tauscaling,found,indextime)
	443	if (.not.found) then
	444	write(,) "phyetat0: <tauscaling> not in file"
	445	tauscaling(:) = 1
	446	endif
[1208]	447	else
[2281]	448	tauscaling(:) = 1
	449	endif ! if (startphy_file)
	450	write(,) "phyetat0: dust conversion factor <tauscaling> range:", &
	451	minval(tauscaling), maxval(tauscaling)
[1208]	452
[2417]	453	! dust_rad_adjust_* for radiative rescaling of dust
	454	if (dustscaling_mode==2) then
	455	if (startphy_file) then
	456	call get_field("dust_rad_adjust_prev",dust_rad_adjust_prev,found,indextime)
	457	if (.not.found) then
	458	write(,) "phyetat0: <dust_rad_adjust_prev> not in file; set to 1"
	459	dust_rad_adjust_prev(:) = 1
	460	endif
	461	call get_field("dust_rad_adjust_next",dust_rad_adjust_next,found,indextime)
	462	if (.not.found) then
	463	write(,) "phyetat0: <dust_rad_adjust_next> not in file; set to 1"
	464	dust_rad_adjust_next(:) = 1
	465	endif
	466	else
	467	dust_rad_adjust_prev(:)= 0
	468	dust_rad_adjust_next(:)= 0
	469	endif ! if (startphy_file)
	470	write(,) "phyetat0: radiative scaling coeff <dust_rad_adjust_prev> range:", &
	471	minval(dust_rad_adjust_prev), maxval(dust_rad_adjust_prev)
	472	write(,) "phyetat0: radiative scaling coeff <dust_rad_adjust_next> range:", &
	473	minval(dust_rad_adjust_next), maxval(dust_rad_adjust_next)
	474	endif ! of if (dustscaling_mode==2)
[2281]	475
[2265]	476	! dtau: opacity difference between GCM and dust scenario
[2281]	477	if (startphy_file) then
	478	call get_field("dtau",dtau,found,indextime)
	479	if (.not.found) then
	480	write(,) "phyetat0: <dtau> not in file; set to zero"
	481	dtau(:) = 0
	482	endif
[2265]	483	else
[2281]	484	dtau(:)= 0
	485	endif ! if (startphy_file)
	486	write(,) "phyetat0: opacity diff wrt scenario <dtau> range:", &
	487	minval(dtau), maxval(dtau)
[2265]	488
	489
[1711]	490	! Sub-grid cloud fraction
[2281]	491	if (startphy_file) then
	492	call get_field("totcloudfrac",totcloudfrac,found,indextime)
	493	if (.not.found) then
	494	write(,) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
	495	totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
	496	endif
[1711]	497	else
[2281]	498	totcloudfrac(:)=1.0
	499	endif ! if (startphy_file)
	500	write(,) "phyetat0: total cloud fraction <totcloudfrac> range:", &
	501	minval(totcloudfrac), maxval(totcloudfrac)
[1208]	502
[2281]	503
[1944]	504	! Max vertical velocity in thermals
[2281]	505	if (startphy_file) then
	506	call get_field("wstar",wstar,found,indextime)
	507	if (.not.found) then
	508	write(,) "phyetat0: <wstar> not in file! Set to zero"
	509	wstar(:)=0
	510	endif
[1944]	511	else
[2281]	512	wstar(:)=0
	513	endif ! if (startphy_file)
	514	write(,) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
	515	minval(wstar),maxval(wstar)
[1944]	516
[2281]	517
[1130]	518	! Surface temperature :
[2281]	519	if (startphy_file) then !tsurf
	520	call get_field("tsurf",tsurf,found,indextime)
	521	if (.not.found) then
	522	call abort_physic(modname, &
	523	"phyetat0: Failed loading <tsurf>",1)
	524	endif
[1130]	525	else
[2913]	526	tsurf(:,:)=0. ! will be updated afterwards in physiq !
[2281]	527	endif ! of if (startphy_file)
	528	write(,) "phyetat0: Surface temperature <tsurf> range:", &
	529	minval(tsurf), maxval(tsurf)
[1130]	530
[1944]	531	! Surface albedo
[2281]	532	if (startphy_file) then
[2913]	533	call get_field("albedo",albedo(:,1,:),found,indextime)
[2281]	534	if (.not.found) then
	535	write(,) "phyetat0: Failed loading <albedo>"
[2913]	536	do islope=1,nslope
	537	albedo(:,1,islope)=albedodat(:)
	538	enddo
[2281]	539	endif
[1944]	540	else
[2913]	541	do islope=1,nslope
	542	albedo(:,1,islope)=albedodat(:)
	543	enddo
[2281]	544	endif ! of if (startphy_file)
	545	write(,) "phyetat0: Surface albedo <albedo> range:", &
[2913]	546	minval(albedo(:,1,:)), maxval(albedo(:,1,:))
	547	albedo(:,2,:)=albedo(:,1,:)
[1944]	548
[1130]	549	! Surface emissivity
[2281]	550	if (startphy_file) then
	551	call get_field("emis",emis,found,indextime)
	552	if (.not.found) then
	553	call abort_physic(modname, &
	554	"phyetat0: Failed loading <emis>",1)
	555	endif
[1130]	556	else
[2281]	557	! If no startfi file, use parameter surfemis in def file
[2294]	558	surfemis=0.95
	559	call getin_p("surfemis_without_startfi",surfemis)
	560	print*,"surfemis_without_startfi",surfemis
[2913]	561	emis(:,:)=surfemis
[2281]	562	endif ! of if (startphy_file)
	563	write(,) "phyetat0: Surface emissivity <emis> range:", &
	564	minval(emis), maxval(emis)
[1130]	565
	566
	567	! surface roughness length (NB: z0 is a common in surfdat_h)
[2281]	568	if (startphy_file) then
	569	call get_field("z0",z0,found)
	570	if (.not.found) then
	571	write(,) "phyetat0: Failed loading <z0>"
	572	write(,) 'will use constant value of z0_default:',z0_default
	573	z0(:)=z0_default
	574	endif
[1130]	575	else
[2281]	576	z0(:)=z0_default
	577	endif ! of if (startphy_file)
	578	write(,) "phyetat0: Surface roughness <z0> range:", &
	579	minval(z0), maxval(z0)
[1130]	580
	581
	582	! pbl wind variance
[2281]	583	if (startphy_file) then
	584	call get_field("q2",q2,found,indextime)
	585	if (.not.found) then
	586	call abort_physic(modname, &
	587	"phyetat0: Failed loading <q2>",1)
	588	endif
[1130]	589	else
[2281]	590	q2(:,:)=0.
	591	endif ! of if (startphy_file)
	592	write(,) "phyetat0: PBL wind variance <q2> range:", &
	593	minval(q2), maxval(q2)
[1130]	594
[2281]	595	! Non-orographic gravity waves
	596	if (startphy_file) then
	597	call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
	598	if (.not.found) then
	599	write(,) "phyetat0: <du_nonoro_gwd> not in file"
	600	du_nonoro_gwd(:,:)=0.
	601	endif
[2220]	602	else
[2281]	603	du_nonoro_gwd(:,:)=0.
	604	endif ! if (startphy_file)
	605	write(,) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
	606	write(,) " <du_nonoro_gwd> range:", &
	607	minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
[1130]	608
[2281]	609	if (startphy_file) then
	610	call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
	611	if (.not.found) then
	612	write(,) "phyetat0: <dv_nonoro_gwd> not in file"
	613	dv_nonoro_gwd(:,:)=0.
	614	endif
	615	else ! ! if (startphy_file)
	616	dv_nonoro_gwd(:,:)=0.
	617	endif ! if (startphy_file)
	618	write(,) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
	619	write(,) " <dv_nonoro_gwd> range:", &
	620	minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
[2220]	621
[1130]	622	! tracer on surface
	623	if (nq.ge.1) then
	624	do iq=1,nq
[1621]	625	txt=noms(iq)
[1130]	626	if (txt.eq."h2o_vap") then
	627	! There is no surface tracer for h2o_vap;
	628	! "h2o_ice" should be loaded instead
	629	txt="h2o_ice"
	630	write(,) 'phyetat0: loading surface tracer', &
	631	' h2o_ice instead of h2o_vap'
[2260]	632	write(,) 'iq = ', iq
[1130]	633	endif
[2312]	634
	635	if (hdo) then
	636	if (txt.eq."hdo_vap") then
	637	txt="hdo_ice"
	638	write(,) 'phyetat0: loading surface tracer', &
	639	' hdo_ice instead of hdo_vap'
	640	endif
	641	endif !hdo
	642
[2281]	643	if (startphy_file) then
[2826]	644	if (txt.eq."co2") then
	645	! We first check if co2ice exist in the startfi file (old way)
	646	! CO2 ice cover
[2913]	647	call get_field("co2ice",qsurf(:,iq,:),found,indextime)
[2826]	648	! If not, we load the corresponding tracer in qsurf
	649	if (.not.found) then
[2913]	650	call get_field(txt,qsurf(:,iq,:),found,indextime)
[2826]	651	if (.not.found) then
	652	call abort_physic(modname, &
	653	"phyetat0: Failed loading co2ice. there is neither the variable co2ice nor qsurf",1)
	654	endif
	655	endif
	656	else ! (not the co2 tracer)
[2913]	657	call get_field(txt,qsurf(:,iq,:),found,indextime)
[2826]	658	if (.not.found) then
	659	write(,) "phyetat0: Failed loading <",trim(txt),">"
	660	write(,) " ",trim(txt)," is set to zero"
[2913]	661	qsurf(:,iq,:)=0.
[2826]	662	endif
	663	endif !endif co2
	664	else !(not startphy_file)
[2913]	665	qsurf(:,iq,:)=0.
[2281]	666	endif ! of if (startphy_file)
	667	write(,) "phyetat0: Surface tracer <",trim(txt),"> range:", &
[2913]	668	minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
[2281]	669	enddo ! of do iq=1,nq
[2826]	670
	671	if (txt.eq."co2") then
	672	! We first check if co2ice exist in the startfi file (old way)
	673	! CO2 ice cover
	674	if (startphy_file) then
[2913]	675	call get_field("co2ice",qsurf(:,iq,:),found,indextime)
[2826]	676	! If not, we load the corresponding tracer in qsurf
	677	if (.not.found) then
[2913]	678	call get_field(txt,qsurf(:,iq,:),found,indextime)
[2826]	679	if (.not.found) then
	680	call abort_physic(modname, &
	681	"phyetat0: Failed loading co2ice",1)
	682	endif
	683	endif
	684	else
	685	! If we run without startfile, co2ice is set to 0
[2913]	686	qsurf(:,iq,:)=0.
[2826]	687	endif !if (startphy_file)
	688	write(,) "phyetat0: CO2 ice cover <co2ice> range:", &
[2913]	689	minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
[2826]	690	endif
	691
[1130]	692	endif ! of if (nq.ge.1)
	693
[2281]	694	if (startphy_file) then
	695	call get_field("watercap",watercap,found,indextime)
	696	if (.not.found) then
	697	write(,) "phyetat0: Failed loading <watercap> : ", &
	698	"<watercap> is set to zero"
[2913]	699	watercap(:,:)=0
[2281]	700	write(,) 'Now transfer negative surface water ice to', &
	701	' watercap !'
	702	if (nq.ge.1) then
	703	do iq=1,nq
	704	txt=noms(iq)
	705	if (txt.eq."h2o_ice") then
	706	do ig=1,ngrid
[2913]	707	do islope=1,nslope
	708	if (qsurf(ig,iq,islope).lt.0.0) then
	709	watercap(ig,islope) = qsurf(ig,iq,islope)
	710	qsurf(ig,iq,islope) = 0.0
	711	end if
	712	enddo
[2281]	713	end do
	714	endif
[2378]	715
	716	if (txt.eq."hdo_ice") then
	717	do ig=1,ngrid
[2913]	718	do islope=1,nslope
	719	if (qsurf(ig,iq,islope).lt.0.0) then
	720	qsurf(ig,iq,islope) = 0.0
	721	end if
	722	enddo
[2378]	723	end do
	724	endif
	725
[2281]	726	enddo
	727	endif ! of if (nq.ge.1)
[2892]	728	endif ! of if (.not.found)
[2264]	729	else
[2913]	730	watercap(:,:)=0
[2281]	731	endif ! of if (startphy_file)
	732	write(,) "phyetat0: Surface water ice <watercap> range:", &
	733	minval(watercap), maxval(watercap)
[2260]	734
[2281]	735	if (startphy_file) then
	736	! Call to soil_settings, in order to read soil temperatures,
	737	! as well as thermal inertia and volumetric heat capacity
	738	call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
[2887]	739	else
[2919]	740	flux_geo(:,:) = 0.
[2281]	741	endif ! of if (startphy_file)
[2884]	742
	743	if (startphy_file) then
	744	call get_field("watercaptag",watercaptag_tmp,found,indextime)
	745	if (.not.found) then
	746	write(,) "phyetat0: Failed loading <watercaptag> : ", &
	747	"<watercaptag> is set as defined by icelocationmode in surfini.F"
	748	watercaptag(:)=.false.
	749	else
	750	do ig=1,ngrid
	751	if(watercaptag_tmp(ig).lt.0.5) then
	752	watercaptag(ig)=.false.
	753	else
	754	watercaptag(ig)=.true.
	755	endif
	756	enddo
	757	endif
	758	endif !startphy_file
	759
[2994]	760
	761	if (paleoclimate) then
[2999]	762	do iq=1,nq
	763	txt=noms(iq)
	764	if (txt.eq."co2") igcm_co2_tmp = iq
	765	enddo
[2994]	766	if (startphy_file) then
[2999]	767	! Depth of H2O lag
[2994]	768	call get_field("lag_h2o_ice",lag_h2o_ice,found,indextime)
	769	if (.not.found) then
	770	write(,) "phyetat0: Failed loading <lag_h2o_ice> : ", &
	771	"<lag_h2o_ice> is set as -1 (no subsurface ice)"
	772	lag_h2o_ice(:,:) = -1.
	773	endif
	774
[2999]	775	! Depth of CO2 lag
[2994]	776	call get_field("lag_co2_ice",lag_co2_ice,found,indextime)
	777	if (.not.found) then
	778	write(,) "phyetat0: Failed loading <lag_co2_ice> : ", &
	779	"<lag_co2_ice> is set as -1 (no subsurface ice)"
	780	lag_co2_ice(:,:) = -1.
	781	endif
[2999]	782
	783	! Perenial CO2 ice
	784	call get_field("perenial_co2ice",perenial_co2ice,found,indextime)
	785	if (.not.found) then
	786	write(,) "phyetat0: Failed loading <perenial_co2ice> : ", &
	787	"<perenial_co2ice> is set as 10m at the South Pole"
	788	perenial_co2ice(:,:) = 0.
	789	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
	790	DO islope = 1,nslope
	791	perenial_co2ice(ngrid,islope) = 10*1.6e3 ! 10m which is convert to kg/m^2
	792	qsurf(ngrid,igcm_co2_tmp,islope) = qsurf(ngrid-1,igcm_co2_tmp,islope) + & ! perenial ice + frost
	793	perenial_co2ice(ngrid,islope)
	794	ENDDO
	795	endif
	796	endif ! notfound
	797	else ! no startfiphyle
[2994]	798	lag_h2o_ice(:,:) = -1.
	799	lag_co2_ice(:,:) = -1.
[2999]	800	perenial_co2ice(:,:) = 0.
	801	if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
	802	DO islope = 1,nslope
	803	perenial_co2ice(ngrid,islope) = 10*1.6e3 ! 10m which is convert to kg/m^2
	804	qsurf(ngrid,igcm_co2_tmp,islope) = qsurf(ngrid-1,igcm_co2_tmp,islope) + & ! perenial ice + frost
	805	perenial_co2ice(ngrid,islope)
	806	ENDDO
	807	endif
[2994]	808	endif !startphy_file
	809	else
	810	lag_h2o_ice(:,:) = -1.
	811	lag_co2_ice(:,:) = -1.
[2999]	812	perenial_co2ice(:,:) = 0.
[2994]	813	endif !paleoclimate
	814
[1130]	815	! close file:
	816	!
[2281]	817	if (startphy_file) call close_startphy
[1130]	818
[1944]	819	end subroutine phyetat0
	820
[2545]	821
	822	subroutine ini_tab_controle_dyn_xios(idayref)
	823
	824	USE comcstfi_h, only: g, mugaz, omeg, rad, rcp
	825	USE time_phylmdz_mod, ONLY: hour_ini, daysec, dtphys
	826	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
	827	IMPLICIT NONE
	828
	829
	830	INTEGER*4,intent(in) :: idayref ! date (initial date for this run)
	831
	832
	833	INTEGER length,l
	834	parameter (length = 100)
	835	REAL tab_cntrl(length) ! run parameters are stored in this array
	836
	837	DO l=1,length
	838	tab_cntrl(l)=0.
	839	ENDDO
	840	tab_cntrl(1) = real(nbp_lon)
	841	tab_cntrl(2) = real(nbp_lat-1)
	842	tab_cntrl(3) = real(nbp_lev)
	843	tab_cntrl(4) = real(idayref)
	844	tab_cntrl(5) = rad
	845	tab_cntrl(6) = omeg
	846	tab_cntrl(7) = g
	847	tab_cntrl(8) = mugaz
	848	tab_cntrl(9) = rcp
	849	tab_cntrl(10) = daysec
	850	tab_cntrl(11) = dtphys
	851
	852	tab_cntrl(27) = hour_ini
	853
	854	tab_cntrl_mod = tab_cntrl
	855
	856	end subroutine ini_tab_controle_dyn_xios
	857
[1944]	858	end module phyetat0_mod

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

Download in other formats:

Original Format