Context Navigation

testphys1d.F @ 1009

Last change on this file since 1009 was 999, checked in by tnavarro, 11 years ago
Possibility to store multiple initial states in one start/startfi. This is RETROCOMPATIBLE. New option ecrithist in run.def to write data in start/startfi every ecrithist dynamical timestep. New option timestart in run.def to initialize the GCM with the time timestart stored in start
File size: 26.0 KB

Rev	Line
[38]	1
	2	PROGRAM testphys1d
	3	! to use 'getin'
	4	USE ioipsl_getincom
	5	IMPLICIT NONE
	6
	7	c=======================================================================
	8	c subject:
	9	c --------
	10	c PROGRAM useful to run physical part of the martian GCM in a 1D column
	11	c
	12	c Can be compiled with a command like (e.g. for 25 layers)
	13	c "makegcm -p mars -d 25 testphys1d"
	14	c It requires the files "testphys1d.def" "callphys.def"
	15	c and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
	16	c and a file describing the sigma layers (e.g. "z2sig.def")
	17	c
	18	c author: Frederic Hourdin, R.Fournier,F.Forget
	19	c -------
	20	c
	21	c update: 12/06/2003 including chemistry (S. Lebonnois)
	22	c and water ice (F. Montmessin)
	23	c
	24	c=======================================================================
	25
	26	#include "dimensions.h"
	27	#include "dimphys.h"
	28	#include "dimradmars.h"
	29	#include "comgeomfi.h"
	30	#include "surfdat.h"
[899]	31	#include "slope.h"
[38]	32	#include "comsoil.h"
	33	#include "comdiurn.h"
	34	#include "callkeys.h"
	35	#include "comcstfi.h"
	36	#include "planete.h"
	37	#include "comsaison.h"
	38	#include "yomaer.h"
	39	#include "control.h"
	40	#include "comvert.h"
	41	#include "netcdf.inc"
	42	#include "comg1d.h"
	43	#include "logic.h"
	44	#include "advtrac.h"
	45
	46	c --------------------------------------------------------------
	47	c Declarations
	48	c --------------------------------------------------------------
	49	c
	50	INTEGER unitstart ! unite d'ecriture de "startfi"
	51	INTEGER nlayer,nlevel,nsoil,ndt
	52	INTEGER ilayer,ilevel,isoil,idt,iq
	53	LOGICAl firstcall,lastcall
	54	c
[627]	55	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
	56	c
[38]	57	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
	58	REAL day ! date durant le run
	59	REAL time ! time (0<time<1 ; time=0.5 a midi)
	60	REAL play(nlayermx) ! Pressure at the middle of the layers (Pa)
	61	REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
	62	REAL psurf,tsurf
	63	REAL u(nlayermx),v(nlayermx) ! zonal, meridional wind
	64	REAL gru,grv ! prescribed "geostrophic" background wind
	65	REAL temp(nlayermx) ! temperature at the middle of the layers
	66	REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
	67	REAL qsurf(nqmx) ! tracer surface budget (e.g. kg.m-2)
	68	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
	69	REAL co2ice ! co2ice layer (kg.m-2)
	70	REAL emis ! surface layer
	71	REAL q2(nlayermx+1) ! Turbulent Kinetic Energy
	72	REAL zlay(nlayermx) ! altitude estimee dans les couches (km)
	73
	74	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
	75	REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
	76	REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
	77	REAL dpsurf
	78	REAL dq(nlayermx,nqmx)
	79	REAL dqdyn(nlayermx,nqmx)
	80
	81	c Various intermediate variables
	82	INTEGER thermo
	83	REAL zls
	84	REAL phi(nlayermx),h(nlayermx),s(nlayermx)
	85	REAL pks, ptif, w(nlayermx)
	86	REAL qtotinit, mqtot(nqmx),qtot
	87	INTEGER ierr, aslun
	88	REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
	89	Logical tracerdyn
	90	integer :: nq=1 ! number of tracers
	91
	92	character*2 str2
	93	character (len=7) :: str7
	94	character(len=44) :: txt
	95
	96	c=======================================================================
	97
	98	c=======================================================================
	99	c INITIALISATION
	100	c=======================================================================
	101
	102	c ------------------------------------------------------
[899]	103	c Prescribed constants to be set here
[38]	104	c ------------------------------------------------------
	105
	106	pi=2.E+0*asin(1.E+0)
	107
[899]	108	c Mars planetary constants
[38]	109	c ----------------------------
[899]	110	rad=3397200. ! mars radius (m) ~3397200 m
	111	daysec=88775. ! length of a sol (s) ~88775 s
	112	omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1)
	113	g=3.72 ! gravity (m.s-2) ~3.72
	114	mugaz=43.49 ! atmosphere mola mass (g.mol-1) ~43.49
[38]	115	rcp=.256793 ! = r/cp ~0.256793
	116	r= 8.314511E+0 *1000.E+0/mugaz
	117	cpp= r/rcp
[899]	118	year_day = 669 ! lenght of year (sols) ~668.6
	119	periheli = 206.66 ! minimum sun-mars distance (Mkm) ~206.66
	120	aphelie = 249.22 ! maximum sun-mars distance (Mkm) ~249.22
	121	peri_day = 485. ! perihelion date (sols since N. Spring)
	122	obliquit = 25.2 ! Obliquity (deg) ~25.2
[38]	123
[899]	124	c Planetary Boundary Layer and Turbulence parameters
	125	c --------------------------------------------------
	126	z0_default = 1.e-2 ! surface roughness (m) ~0.01
	127	emin_turb = 1.e-6 ! minimal turbulent energy ~1.e-8
	128	lmixmin = 30 ! mixing length ~100
[38]	129
[899]	130	c cap properties and surface emissivities
[38]	131	c ----------------------------------------------------
[899]	132	emissiv= 0.95 ! Bare ground emissivity ~.95
	133	emisice(1)=0.95 ! Northern cap emissivity
	134	emisice(2)=0.95 ! Southern cap emisssivity
	135	albedice(1)=0.5 ! Northern cap albedo
	136	albedice(2)=0.5 ! Southern cap albedo
[38]	137	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
	138	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
	139	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
	140	dtemisice(2) = 2. ! time scale for snow metamorphism (south
	141
[899]	142
[38]	143	c ------------------------------------------------------
[899]	144	c Loading run parameters from "run.def" file
[38]	145	c ------------------------------------------------------
	146
	147
	148	! check if 'run.def' file is around (otherwise reading parameters
	149	! from callphys.def via getin() routine won't work.
	150	open(99,file='run.def',status='old',form='formatted',
	151	& iostat=ierr)
	152	if (ierr.ne.0) then
	153	write(,) 'Cannot find required file "run.def"'
	154	write(,) ' (which should contain some input parameters'
	155	write(,) ' along with the following line:'
	156	write(,) ' INCLUDEDEF=callphys.def'
	157	write(,) ' )'
	158	write(,) ' ... might as well stop here ...'
	159	stop
	160	else
	161	close(99)
	162	endif
	163
	164	! check if we are going to run with or without tracers
	165	write(,) "Run with or without tracer transport ?"
	166	tracer=.false. ! default value
	167	call getin("tracer",tracer)
	168	write(,) " tracer = ",tracer
	169
	170	! while we're at it, check if there is a 'traceur.def' file
	171	! and preocess it, if necessary. Otherwise initialize tracer names
	172	if (tracer) then
	173	! load tracer names from file 'traceur.def'
	174	open(90,file='traceur.def',status='old',form='formatted',
	175	& iostat=ierr)
	176	if (ierr.ne.0) then
	177	write(,) 'Cannot find required file "traceur.def"'
	178	write(,) ' If you want to run with tracers, I need it'
	179	write(,) ' ... might as well stop here ...'
	180	stop
	181	else
	182	write(,) "testphys1d: Reading file traceur.def"
	183	! read number of tracers:
	184	read(90,*,iostat=ierr) nq
	185	if (ierr.ne.0) then
	186	write(,) "testphys1d: error reading number of tracers"
	187	write(,) " (first line of traceur.def) "
	188	stop
	189	else
	190	! check that the number of tracers is indeed nqmx
	191	if (nq.ne.nqmx) then
	192	write(,) "testphys1d: error, wrong number of tracers:"
	193	write(,) "nq=",nq," whereas nqmx=",nqmx
	194	stop
	195	endif
	196	endif
	197	endif
	198	! read tracer names from file traceur.def
	199	do iq=1,nqmx
	200	read(90,*,iostat=ierr) tnom(iq)
	201	if (ierr.ne.0) then
	202	write(,) 'testphys1d: error reading tracer names...'
	203	stop
	204	endif
	205	enddo
	206	close(90)
	207
	208	! initialize tracers here:
	209	write(,) "testphys1d: initializing tracers"
	210	q(:,:)=0 ! default, set everything to zero
	211	qsurf(:)=0
	212	! "smarter" initialization of some tracers
	213	! (get values from "profile_*" files, if these are available)
	214	do iq=1,nqmx
	215	txt=""
	216	write(txt,"(a)") tnom(iq)
	217	write(,)" tracer:",trim(txt)
	218	! CO2
	219	if (txt.eq."co2") then
	220	q(:,iq)=0.95 ! kg /kg of atmosphere
	221	qsurf(iq)=0. ! kg/m2 (not used for CO2)
	222	! even better, look for a "profile_co2" input file
	223	open(91,file='profile_co2',status='old',
	224	& form='formatted',iostat=ierr)
	225	if (ierr.eq.0) then
	226	read(91,*) qsurf(iq)
	227	do ilayer=1,nlayermx
	228	read(91,*) q(ilayer,iq)
	229	enddo
	230	endif
	231	close(91)
	232	endif ! of if (txt.eq."co2")
[161]	233	! Allow for an initial profile of argon
	234	! Can also be used to introduce a decaying tracer
	235	! in the 1D (TBD) to study thermals
	236	if (txt.eq."ar") then
	237	!look for a "profile_ar" input file
	238	open(91,file='profile_ar',status='old',
	239	& form='formatted',iostat=ierr)
	240	if (ierr.eq.0) then
	241	read(91,*) qsurf(iq)
	242	do ilayer=1,nlayermx
	243	read(91,*) q(ilayer,iq)
	244	enddo
	245	else
	246	write(,) "No profile_ar file!"
	247	endif
	248	close(91)
	249	endif ! of if (txt.eq."ar")
	250
[38]	251	! WATER VAPOUR
	252	if (txt.eq."h2o_vap") then
	253	!look for a "profile_h2o_vap" input file
	254	open(91,file='profile_h2o_vap',status='old',
	255	& form='formatted',iostat=ierr)
	256	if (ierr.eq.0) then
	257	read(91,*) qsurf(iq)
	258	do ilayer=1,nlayermx
	259	read(91,*) q(ilayer,iq)
	260	enddo
	261	else
	262	write(,) "No profile_h2o_vap file!"
	263	endif
	264	close(91)
	265	endif ! of if (txt.eq."h2o_ice")
	266	! WATER ICE
	267	if (txt.eq."h2o_ice") then
	268	!look for a "profile_h2o_vap" input file
	269	open(91,file='profile_h2o_ice',status='old',
	270	& form='formatted',iostat=ierr)
	271	if (ierr.eq.0) then
	272	read(91,*) qsurf(iq)
	273	do ilayer=1,nlayermx
	274	read(91,*) q(ilayer,iq)
	275	enddo
	276	else
	277	write(,) "No profile_h2o_ice file!"
	278	endif
	279	close(91)
	280	endif ! of if (txt.eq."h2o_ice")
	281	! DUST
	282	!if (txt(1:4).eq."dust") then
	283	! q(:,iq)=0.4 ! kg/kg of atmosphere
	284	! qsurf(iq)=100 ! kg/m2
	285	!endif
	286	! DUST MMR
	287	if (txt.eq."dust_mass") then
	288	!look for a "profile_dust_mass" input file
	289	open(91,file='profile_dust_mass',status='old',
	290	& form='formatted',iostat=ierr)
	291	if (ierr.eq.0) then
	292	read(91,*) qsurf(iq)
	293	do ilayer=1,nlayermx
	294	read(91,*) q(ilayer,iq)
	295	! write(,) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
	296	enddo
	297	else
	298	write(,) "No profile_dust_mass file!"
	299	endif
	300	close(91)
	301	endif ! of if (txt.eq."dust_mass")
	302	! DUST NUMBER
	303	if (txt.eq."dust_number") then
	304	!look for a "profile_dust_number" input file
	305	open(91,file='profile_dust_number',status='old',
	306	& form='formatted',iostat=ierr)
	307	if (ierr.eq.0) then
	308	read(91,*) qsurf(iq)
	309	do ilayer=1,nlayermx
	310	read(91,*) q(ilayer,iq)
	311	enddo
	312	else
	313	write(,) "No profile_dust_number file!"
	314	endif
	315	close(91)
	316	endif ! of if (txt.eq."dust_number")
[358]	317	! NB: some more initializations (chemistry) is done later
	318	! CCN MASS
	319	if (txt.eq."ccn_mass") then
	320	!look for a "profile_ccn_mass" input file
	321	open(91,file='profile_ccn_mass',status='old',
	322	& form='formatted',iostat=ierr)
	323	if (ierr.eq.0) then
	324	read(91,*) qsurf(iq)
	325	do ilayer=1,nlayermx
	326	read(91,*) q(ilayer,iq)
	327	enddo
	328	else
	329	write(,) "No profile_ccn_mass file!"
	330	endif
	331	close(91)
	332	endif ! of if (txt.eq."ccn_mass")
	333	! CCN NUMBER
	334	if (txt.eq."ccn_number") then
	335	!look for a "profile_ccn_number" input file
	336	open(91,file='profile_ccn_number',status='old',
	337	& form='formatted',iostat=ierr)
	338	if (ierr.eq.0) then
	339	read(91,*) qsurf(iq)
	340	do ilayer=1,nlayermx
	341	read(91,*) q(ilayer,iq)
	342	enddo
	343	else
	344	write(,) "No profile_ccn_number file!"
	345	endif
	346	close(91)
	347	endif ! of if (txt.eq."ccn_number")
[38]	348	enddo ! of do iq=1,nqmx
	349
	350	else
	351	! we still need to set (dummy) tracer names for physdem1
	352	nq=nqmx
	353	do iq=1,nq
	354	write(str7,'(a1,i2.2)')'q',iq
	355	tnom(iq)=str7
	356	enddo
[899]	357	! and just to be clean, also initialize tracers to zero for physdem1
	358	q(:,:)=0
	359	qsurf(:)=0
[38]	360	endif ! of if (tracer)
[358]	361
	362	!write(,) "testphys1d q", q(1,:)
	363	!write(,) "testphys1d qsurf", qsurf
[38]	364
[899]	365	c Date and local time at beginning of run
	366	c ---------------------------------------
	367	c Date (in sols since spring solstice) at beginning of run
[38]	368	day0 = 0 ! default value for day0
	369	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
	370	call getin("day0",day0)
	371	day=float(day0)
	372	write(,) " day0 = ",day0
[899]	373	c Local time at beginning of run
[38]	374	time=0 ! default value for time
	375	write(,)'Initial local time (in hours, between 0 and 24)?'
	376	call getin("time",time)
	377	write(,)" time = ",time
	378	time=time/24.E+0 ! convert time (hours) to fraction of sol
	379
[899]	380	c Discretization (Definition of grid and time steps)
[38]	381	c --------------
	382	c
	383	nlayer=nlayermx
	384	nlevel=nlayer+1
	385	nsoil=nsoilmx
	386
	387	day_step=48 ! default value for day_step
	388	PRINT *,'Number of time steps per sol ?'
	389	call getin("day_step",day_step)
	390	write(,) " day_step = ",day_step
	391
	392	ndt=10 ! default value for ndt
	393	PRINT *,'Number of sols to run ?'
	394	call getin("ndt",ndt)
	395	write(,) " ndt = ",ndt
	396
	397	ndt=ndt*day_step
	398	dtphys=daysec/day_step
[899]	399
	400	c Imposed surface pressure
[38]	401	c ------------------------------------
	402	c
	403	psurf=610. ! default value for psurf
	404	PRINT *,'Surface pressure (Pa) ?'
	405	call getin("psurf",psurf)
	406	write(,) " psurf = ",psurf
[899]	407	c Reference pressures
	408	pa=20. ! transition pressure (for hybrid coord.)
	409	preff=610. ! reference surface pressure
[38]	410
[899]	411	c Aerosol properties
[38]	412	c --------------------------------
[899]	413	tauvis=0.2 ! default value for tauvis (dust opacity)
[627]	414	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
[38]	415	call getin("tauvis",tauvis)
	416	write(,) " tauvis = ",tauvis
[720]	417
	418	c Orbital parameters
	419	c ------------------
	420	print *,'Min. distance Sun-Mars (Mkm)?'
	421	call getin("periheli",periheli)
	422	write(,) " periheli = ",periheli
	423
	424	print *,'Max. distance Sun-Mars (Mkm)?'
	425	call getin("aphelie",aphelie)
	426	write(,) " aphelie = ",aphelie
	427
	428	print *,'Day of perihelion?'
	429	call getin("periday",peri_day)
	430	write(,) " periday = ",peri_day
	431
	432	print *,'Obliquity?'
	433	call getin("obliquit",obliquit)
	434	write(,) " obliquit = ",obliquit
[38]	435
	436	c latitude/longitude
	437	c ------------------
	438	lati(1)=0 ! default value for lati(1)
	439	PRINT *,'latitude (in degrees) ?'
	440	call getin("latitude",lati(1))
	441	write(,) " latitude = ",lati(1)
	442	lati(1)=lati(1)*pi/180.E+0
	443	long(1)=0.E+0
	444	long(1)=long(1)*pi/180.E+0
	445
[899]	446	c Initialize albedo / soil thermal inertia
	447	c ----------------------------------------
[38]	448	c
	449	albedodat(1)=0.2 ! default value for albedodat
	450	PRINT *,'Albedo of bare ground ?'
	451	call getin("albedo",albedodat(1))
	452	write(,) " albedo = ",albedodat(1)
	453
	454	inertiedat(1,1)=400 ! default value for inertiedat
	455	PRINT *,'Soil thermal inertia (SI) ?'
	456	call getin("inertia",inertiedat(1,1))
	457	write(,) " inertia = ",inertiedat(1,1)
[224]	458
	459	z0(1)=z0_default ! default value for roughness
	460	write(,) 'Surface roughness length z0 (m)?'
	461	call getin("z0",z0(1))
	462	write(,) " z0 = ",z0(1)
[899]	463
	464	! Initialize local slope parameters (only matters if "callslope"
	465	! is .true. in callphys.def)
	466	! slope inclination angle (deg) 0: horizontal, 90: vertical
	467	theta_sl(1)=0.0 ! default: no inclination
	468	call getin("slope_inclination",theta_sl(1))
	469	! slope orientation (deg)
	470	! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
	471	psi_sl(1)=0.0 ! default value
	472	call getin("slope_orientation",psi_sl(1))
	473
[38]	474	c
[899]	475	c for the gravity wave scheme
[38]	476	c ---------------------------------
	477	c
	478	zmea(1)=0.E+0
	479	zstd(1)=0.E+0
	480	zsig(1)=0.E+0
	481	zgam(1)=0.E+0
	482	zthe(1)=0.E+0
	483
	484
[899]	485	c Specific initializations for "physiq"
	486	c -------------------------------------
	487	c mesh surface (not a very usefull quantity in 1D)
[38]	488	area(1)=1.E+0
	489
[899]	490	c surface geopotential is not used (or useful) since in 1D
	491	c everything is controled by surface pressure
[38]	492	phisfi(1)=0.E+0
	493
[899]	494	c "inifis" does some initializations (some of which have already been
	495	c done above!) and loads parameters set in callphys.def
[38]	496
	497	!Mars possible matter with dtphys in input and include!!!
	498	CALL inifis(1,llm,day0,daysec,dtphys,
	499	. lati,long,area,rad,g,r,cpp)
[899]	500
	501	c Initialization to take into account prescribed winds
[38]	502	c ------------------------------------------------------
	503	ptif=2.E+0omegsinlat(1)
	504
[899]	505	c geostrophic wind
[38]	506	gru=10. ! default value for gru
	507	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
	508	call getin("u",gru)
	509	write(,) " u = ",gru
	510	grv=0. !default value for grv
	511	PRINT *,'meridional northward component of the geostrophic',
	512	&' wind (m/s) ?'
	513	call getin("v",grv)
	514	write(,) " v = ",grv
	515
[899]	516	c Initialize winds for first time step
[38]	517	DO ilayer=1,nlayer
	518	u(ilayer)=gru
	519	v(ilayer)=grv
	520	ENDDO
	521
[899]	522	c Initialize turbulente kinetic energy
[38]	523	DO ilevel=1,nlevel
	524	q2(ilevel)=0.E+0
	525	ENDDO
	526
[899]	527	c CO2 ice on the surface
[38]	528	c -------------------
	529	co2ice=0.E+0 ! default value for co2ice
	530	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
	531	call getin("co2ice",co2ice)
	532	write(,) " co2ice = ",co2ice
	533
	534	c
[899]	535	c emissivity
[38]	536	c ----------
	537	emis=emissiv
	538	IF (co2ice.eq.1.E+0) THEN
	539	emis=emisice(1) ! northern hemisphere
	540	IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
	541	ENDIF
	542
	543
	544
[899]	545	c Compute pressures and altitudes of atmospheric levels
[38]	546	c ----------------------------------------------------------------
	547
	548	c Vertical Coordinates
	549	c """"""""""""""""""""
	550	hybrid=.true.
	551	PRINT *,'Hybrid coordinates ?'
	552	call getin("hybrid",hybrid)
	553	write(,) " hybrid = ", hybrid
	554
	555	CALL disvert
	556
	557	DO ilevel=1,nlevel
	558	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	559	ENDDO
	560
	561	DO ilayer=1,nlayer
	562	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	563	ENDDO
	564
	565	DO ilayer=1,nlayer
	566	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
	567	& /g
	568	ENDDO
	569
	570
[899]	571	c Initialize temperature profile
[38]	572	c --------------------------------------
	573	pks=psurf**rcp
	574
[899]	575	c altitude in km in profile: divide zlay by 1000
[38]	576	tmp1(0)=0.E+0
	577	DO ilayer=1,nlayer
	578	tmp1(ilayer)=zlay(ilayer)/1000.E+0
	579	ENDDO
	580
	581	call profile(nlayer+1,tmp1,tmp2)
	582
	583	tsurf=tmp2(0)
	584	DO ilayer=1,nlayer
	585	temp(ilayer)=tmp2(ilayer)
	586	ENDDO
	587
	588
	589
	590	! Initialize soil properties and temperature
	591	! ------------------------------------------
	592	volcapa=1.e6 ! volumetric heat capacity
	593	DO isoil=1,nsoil
	594	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
	595	tsoil(isoil)=tsurf ! soil temperature
	596	ENDDO
	597
	598	! Initialize depths
	599	! -----------------
	600	do isoil=0,nsoil-1
	601	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
	602	enddo
	603	do isoil=1,nsoil
	604	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
	605	enddo
	606
[899]	607	c Initialize traceurs
[38]	608	c ---------------------------
	609
	610	if (photochem.or.callthermos) then
[899]	611	write(,) 'Initializing chemical species'
	612	! thermo=0: initialize over all atmospheric layers
[38]	613	thermo=0
	614	call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
	615	$ thermo,qsurf)
	616	endif
[520]	617
[899]	618	c Check if the surface is a water ice reservoir
[520]	619	c --------------------------------------------------
[899]	620	watercaptag(ngridmx)=.false. ! Default: no water ice reservoir
[520]	621	print *,'Water ice cap on ground ?'
	622	call getin("watercaptag",watercaptag)
	623	write(,) " watercaptag = ",watercaptag
[38]	624
	625
[899]	626	c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
[38]	627	c ----------------------------------------------------------------
[899]	628	c (output done in "writeg1d", typically called by "physiq.F")
[38]	629
	630	g1d_nlayer=nlayer
	631	g1d_nomfich='g1d.dat'
	632	g1d_unitfich=40
	633	g1d_nomctl='g1d.ctl'
	634	g1d_unitctl=41
	635	g1d_premier=.true.
	636	g2d_premier=.true.
	637
[899]	638	c Write a "startfi" file
[38]	639	c --------------------
[899]	640	c This file will be read during the first call to "physiq".
	641	c It is needed to transfert physics variables to "physiq"...
[38]	642
[999]	643	call physdem0("startfi.nc",long,lati,nsoilmx,nqmx,
	644	. dtphys,float(day0),time,area,
	645	. albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
	646	call physdem1("startfi.nc",nsoilmx,nqmx,
	647	. dtphys,time,
	648	. tsurf,tsoil,co2ice,emis,q2,qsurf)
[899]	649
[38]	650	c=======================================================================
[899]	651	c 1D MODEL TIME STEPPING LOOP
[38]	652	c=======================================================================
	653	c
	654	firstcall=.true.
	655	lastcall=.false.
	656
	657	DO idt=1,ndt
	658	c IF (idt.eq.ndt) lastcall=.true.
	659	IF (idt.eq.ndt-day_step-1) then !test
	660	lastcall=.true.
	661	call solarlong(day*1.0,zls)
	662	write(103,) 'Ls=',zls180./pi
	663	write(103,) 'Lat=', lati(1)180./pi
[627]	664	write(103,) 'Tau=', tauvis/odprefpsurf
[38]	665	write(103,*) 'RunEnd - Atmos. Temp. File'
	666	write(103,*) 'RunEnd - Atmos. Temp. File'
	667	write(104,) 'Ls=',zls180./pi
	668	write(104,*) 'Lat=', lati(1)
[627]	669	write(104,) 'Tau=', tauvis/odprefpsurf
[38]	670	write(104,*) 'RunEnd - Atmos. Temp. File'
	671	ENDIF
	672
[899]	673	c compute geopotential
[38]	674	c ~~~~~~~~~~~~~~~~~~~~~
	675	DO ilayer=1,nlayer
	676	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	677	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	678	ENDDO
	679	phi(1)=pksh(1)(1.E+0-s(1))
	680	DO ilayer=2,nlayer
	681	phi(ilayer)=phi(ilayer-1)+
	682	& pks(h(ilayer-1)+h(ilayer)).5E+0
	683	& *(s(ilayer-1)-s(ilayer))
	684
	685	ENDDO
	686
[899]	687	c call physics
[38]	688	c --------------------
[358]	689	! write(,) "testphys1d avant q", q(1,:)
[38]	690	CALL physiq (1,llm,nqmx,
	691	, firstcall,lastcall,
	692	, day,time,dtphys,
	693	, plev,play,phi,
	694	, u, v,temp, q,
	695	, w,
[899]	696	C - outputs
[38]	697	s du, dv, dtemp, dq,dpsurf,tracerdyn)
[358]	698	! write(,) "testphys1d apres q", q(1,:)
[38]	699
[358]	700
[899]	701	c wind increment : specific for 1D
	702	c --------------------------------
[38]	703
[899]	704	c The physics compute the tendencies on u and v,
	705	c here we just add Coriolos effect
[38]	706	c
	707	c DO ilayer=1,nlayer
	708	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
	709	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
	710	c ENDDO
	711
[899]	712	c For some tests : No coriolis force at equator
[38]	713	c if(lati(1).eq.0.) then
	714	DO ilayer=1,nlayer
	715	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
	716	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
	717	ENDDO
	718	c end if
	719	c
	720	c
[899]	721	c Compute time for next time step
[38]	722	c ---------------------------------------
	723	firstcall=.false.
	724	time=time+dtphys/daysec
	725	IF (time.gt.1.E+0) then
	726	time=time-1.E+0
	727	day=day+1
	728	ENDIF
	729
[899]	730	c compute winds and temperature for next time step
[38]	731	c ----------------------------------------------------------
	732
	733	DO ilayer=1,nlayer
	734	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
	735	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
	736	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
	737	ENDDO
	738
[899]	739	c compute pressure for next time step
[38]	740	c ----------------------------------------------------------
	741
[899]	742	psurf=psurf+dtphys*dpsurf ! surface pressure change
[38]	743	DO ilevel=1,nlevel
	744	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	745	ENDDO
	746	DO ilayer=1,nlayer
	747	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	748	ENDDO
	749
[899]	750	! increment tracers
[38]	751	DO iq = 1, nqmx
	752	DO ilayer=1,nlayer
	753	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
	754	ENDDO
	755	ENDDO
	756
[899]	757	ENDDO ! of idt=1,ndt ! end of time stepping loop
[38]	758
	759	c ========================================================
[899]	760	c OUTPUTS
[38]	761	c ========================================================
	762
[899]	763	c finalize and close grads files "g1d.dat" and "g1d.ctl"
[38]	764
	765	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
	766	CALL endg1d(1,nlayer,zlay/1000.,ndt)
	767
	768	c ========================================================
	769	END
	770
	771	c***********************************************************************
	772	c***********************************************************************
[899]	773	c Dummy subroutines used only in 3D, but required to
	774	c compile testphys1d (to cleanly use writediagfi)
[38]	775
[899]	776	subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
	777
	778	IMPLICIT NONE
	779
	780	INTEGER im,jm,ngrid,nfield
	781	REAL pdyn(im,jm,nfield)
	782	REAL pfi(ngrid,nfield)
	783
	784	if (ngrid.ne.1) then
	785	write(,) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
	786	stop
	787	endif
	788
	789	pdyn(1,1,1:nfield)=pfi(1,1:nfield)
	790
	791	end
[38]	792
	793	c***********************************************************************
	794	c***********************************************************************
	795
	796	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format