Context Navigation

testphys1d.F @ 706

Last change on this file since 706 was 627, checked in by emillour, 13 years ago

Mars GCM:

Some cleanup on messages and comments in code about the reference pressure

for dust opacity which is now 610 Pa.

File size: 24.9 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"
	31	#include "comsoil.h"
	32	#include "comdiurn.h"
	33	#include "callkeys.h"
	34	#include "comcstfi.h"
	35	#include "planete.h"
	36	#include "comsaison.h"
	37	#include "yomaer.h"
	38	#include "control.h"
	39	#include "comvert.h"
	40	#include "netcdf.inc"
	41	#include "comg1d.h"
	42	#include "logic.h"
	43	#include "advtrac.h"
	44
	45	c --------------------------------------------------------------
	46	c Declarations
	47	c --------------------------------------------------------------
	48	c
	49	INTEGER unitstart ! unite d'ecriture de "startfi"
	50	INTEGER nlayer,nlevel,nsoil,ndt
	51	INTEGER ilayer,ilevel,isoil,idt,iq
	52	LOGICAl firstcall,lastcall
	53	c
[627]	54	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
	55	c
[38]	56	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
	57	REAL day ! date durant le run
	58	REAL time ! time (0<time<1 ; time=0.5 a midi)
	59	REAL play(nlayermx) ! Pressure at the middle of the layers (Pa)
	60	REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
	61	REAL psurf,tsurf
	62	REAL u(nlayermx),v(nlayermx) ! zonal, meridional wind
	63	REAL gru,grv ! prescribed "geostrophic" background wind
	64	REAL temp(nlayermx) ! temperature at the middle of the layers
	65	REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
	66	REAL qsurf(nqmx) ! tracer surface budget (e.g. kg.m-2)
	67	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
	68	REAL co2ice ! co2ice layer (kg.m-2)
	69	REAL emis ! surface layer
	70	REAL q2(nlayermx+1) ! Turbulent Kinetic Energy
	71	REAL zlay(nlayermx) ! altitude estimee dans les couches (km)
	72
	73	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
	74	REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
	75	REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
	76	REAL dpsurf
	77	REAL dq(nlayermx,nqmx)
	78	REAL dqdyn(nlayermx,nqmx)
	79
	80	c Various intermediate variables
	81	INTEGER thermo
	82	REAL zls
	83	REAL phi(nlayermx),h(nlayermx),s(nlayermx)
	84	REAL pks, ptif, w(nlayermx)
	85	REAL qtotinit, mqtot(nqmx),qtot
	86	INTEGER ierr, aslun
	87	REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
	88	Logical tracerdyn
	89	integer :: nq=1 ! number of tracers
	90
	91	character*2 str2
	92	character (len=7) :: str7
	93	character(len=44) :: txt
	94
	95	c=======================================================================
	96
	97	c=======================================================================
	98	c INITIALISATION
	99	c=======================================================================
	100
	101	c ------------------------------------------------------
	102	c Constantes prescrites ICI
	103	c ------------------------------------------------------
	104
	105	pi=2.E+0*asin(1.E+0)
	106
	107	c Constante de la Planete Mars
	108	c ----------------------------
	109	rad=3397200. ! rayon de mars (m) ~3397200 m
	110	daysec=88775. ! duree du sol (s) ~88775 s
	111	omeg=4.*asin(1.)/(daysec) ! vitesse de rotation (rad.s-1)
	112	g=3.72 ! gravite (m.s-2) ~3.72
	113	mugaz=43.49 ! Masse molaire de l'atm (g.mol-1) ~43.49
	114	rcp=.256793 ! = r/cp ~0.256793
	115	r= 8.314511E+0 *1000.E+0/mugaz
	116	cpp= r/rcp
	117	year_day = 669 ! duree de l'annee (sols) ~668.6
	118	periheli = 206.66 ! dist.min. soleil-mars (Mkm) ~206.66
	119	aphelie = 249.22 ! dist.max. soleil-mars (Mkm) ~249.22
	120	peri_day = 485. ! date du perihelie (sols depuis printemps)
	121	obliquit = 25.2 ! Obliquite de la planete (deg) ~25.2
	122
	123	c Parametres Couche limite et Turbulence
	124	c --------------------------------------
[224]	125	z0_default = 1.e-2 ! surface roughness (m) ~0.01
[38]	126	emin_turb = 1.e-6 ! energie minimale ~1.e-8
	127	lmixmin = 30 ! longueur de melange ~100
	128
	129	c propriete optiques des calottes et emissivite du sol
	130	c ----------------------------------------------------
	131	emissiv= 0.95 ! Emissivite du sol martien ~.95
	132	emisice(1)=0.95 ! Emissivite calotte nord
	133	emisice(2)=0.95 ! Emissivite calotte sud
	134	albedice(1)=0.5 ! Albedo calotte nord
	135	albedice(2)=0.5 ! Albedo calotte sud
	136	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
	137	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
	138	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
	139	dtemisice(2) = 2. ! time scale for snow metamorphism (south
	140	hybrid=.false.
	141
	142	c ------------------------------------------------------
	143	c Lecture des parametres dans "run.def"
	144	c ------------------------------------------------------
	145
	146
	147	! check if 'run.def' file is around (otherwise reading parameters
	148	! from callphys.def via getin() routine won't work.
	149	open(99,file='run.def',status='old',form='formatted',
	150	& iostat=ierr)
	151	if (ierr.ne.0) then
	152	write(,) 'Cannot find required file "run.def"'
	153	write(,) ' (which should contain some input parameters'
	154	write(,) ' along with the following line:'
	155	write(,) ' INCLUDEDEF=callphys.def'
	156	write(,) ' )'
	157	write(,) ' ... might as well stop here ...'
	158	stop
	159	else
	160	close(99)
	161	endif
	162
	163	! check if we are going to run with or without tracers
	164	write(,) "Run with or without tracer transport ?"
	165	tracer=.false. ! default value
	166	call getin("tracer",tracer)
	167	write(,) " tracer = ",tracer
	168
	169	! while we're at it, check if there is a 'traceur.def' file
	170	! and preocess it, if necessary. Otherwise initialize tracer names
	171	if (tracer) then
	172	! load tracer names from file 'traceur.def'
	173	open(90,file='traceur.def',status='old',form='formatted',
	174	& iostat=ierr)
	175	if (ierr.ne.0) then
	176	write(,) 'Cannot find required file "traceur.def"'
	177	write(,) ' If you want to run with tracers, I need it'
	178	write(,) ' ... might as well stop here ...'
	179	stop
	180	else
	181	write(,) "testphys1d: Reading file traceur.def"
	182	! read number of tracers:
	183	read(90,*,iostat=ierr) nq
	184	if (ierr.ne.0) then
	185	write(,) "testphys1d: error reading number of tracers"
	186	write(,) " (first line of traceur.def) "
	187	stop
	188	else
	189	! check that the number of tracers is indeed nqmx
	190	if (nq.ne.nqmx) then
	191	write(,) "testphys1d: error, wrong number of tracers:"
	192	write(,) "nq=",nq," whereas nqmx=",nqmx
	193	stop
	194	endif
	195	endif
	196	endif
	197	! read tracer names from file traceur.def
	198	do iq=1,nqmx
	199	read(90,*,iostat=ierr) tnom(iq)
	200	if (ierr.ne.0) then
	201	write(,) 'testphys1d: error reading tracer names...'
	202	stop
	203	endif
	204	enddo
	205	close(90)
	206
	207	! initialize tracers here:
	208	write(,) "testphys1d: initializing tracers"
	209	q(:,:)=0 ! default, set everything to zero
	210	qsurf(:)=0
	211	! "smarter" initialization of some tracers
	212	! (get values from "profile_*" files, if these are available)
	213	do iq=1,nqmx
	214	txt=""
	215	write(txt,"(a)") tnom(iq)
	216	write(,)" tracer:",trim(txt)
	217	! CO2
	218	if (txt.eq."co2") then
	219	q(:,iq)=0.95 ! kg /kg of atmosphere
	220	qsurf(iq)=0. ! kg/m2 (not used for CO2)
	221	! even better, look for a "profile_co2" input file
	222	open(91,file='profile_co2',status='old',
	223	& form='formatted',iostat=ierr)
	224	if (ierr.eq.0) then
	225	read(91,*) qsurf(iq)
	226	do ilayer=1,nlayermx
	227	read(91,*) q(ilayer,iq)
	228	enddo
	229	endif
	230	close(91)
	231	endif ! of if (txt.eq."co2")
[161]	232	! Allow for an initial profile of argon
	233	! Can also be used to introduce a decaying tracer
	234	! in the 1D (TBD) to study thermals
	235	if (txt.eq."ar") then
	236	!look for a "profile_ar" input file
	237	open(91,file='profile_ar',status='old',
	238	& form='formatted',iostat=ierr)
	239	if (ierr.eq.0) then
	240	read(91,*) qsurf(iq)
	241	do ilayer=1,nlayermx
	242	read(91,*) q(ilayer,iq)
	243	enddo
	244	else
	245	write(,) "No profile_ar file!"
	246	endif
	247	close(91)
	248	endif ! of if (txt.eq."ar")
	249
[38]	250	! WATER VAPOUR
	251	if (txt.eq."h2o_vap") then
	252	!look for a "profile_h2o_vap" input file
	253	open(91,file='profile_h2o_vap',status='old',
	254	& form='formatted',iostat=ierr)
	255	if (ierr.eq.0) then
	256	read(91,*) qsurf(iq)
	257	do ilayer=1,nlayermx
	258	read(91,*) q(ilayer,iq)
	259	enddo
	260	else
	261	write(,) "No profile_h2o_vap file!"
	262	endif
	263	close(91)
	264	endif ! of if (txt.eq."h2o_ice")
	265	! WATER ICE
	266	if (txt.eq."h2o_ice") then
	267	!look for a "profile_h2o_vap" input file
	268	open(91,file='profile_h2o_ice',status='old',
	269	& form='formatted',iostat=ierr)
	270	if (ierr.eq.0) then
	271	read(91,*) qsurf(iq)
	272	do ilayer=1,nlayermx
	273	read(91,*) q(ilayer,iq)
	274	enddo
	275	else
	276	write(,) "No profile_h2o_ice file!"
	277	endif
	278	close(91)
	279	endif ! of if (txt.eq."h2o_ice")
	280	! DUST
	281	!if (txt(1:4).eq."dust") then
	282	! q(:,iq)=0.4 ! kg/kg of atmosphere
	283	! qsurf(iq)=100 ! kg/m2
	284	!endif
	285	! DUST MMR
	286	if (txt.eq."dust_mass") then
	287	!look for a "profile_dust_mass" input file
	288	open(91,file='profile_dust_mass',status='old',
	289	& form='formatted',iostat=ierr)
	290	if (ierr.eq.0) then
	291	read(91,*) qsurf(iq)
	292	do ilayer=1,nlayermx
	293	read(91,*) q(ilayer,iq)
	294	! write(,) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
	295	enddo
	296	else
	297	write(,) "No profile_dust_mass file!"
	298	endif
	299	close(91)
	300	endif ! of if (txt.eq."dust_mass")
	301	! DUST NUMBER
	302	if (txt.eq."dust_number") then
	303	!look for a "profile_dust_number" input file
	304	open(91,file='profile_dust_number',status='old',
	305	& form='formatted',iostat=ierr)
	306	if (ierr.eq.0) then
	307	read(91,*) qsurf(iq)
	308	do ilayer=1,nlayermx
	309	read(91,*) q(ilayer,iq)
	310	enddo
	311	else
	312	write(,) "No profile_dust_number file!"
	313	endif
	314	close(91)
	315	endif ! of if (txt.eq."dust_number")
[358]	316	! NB: some more initializations (chemistry) is done later
	317	! CCN MASS
	318	if (txt.eq."ccn_mass") then
	319	!look for a "profile_ccn_mass" input file
	320	open(91,file='profile_ccn_mass',status='old',
	321	& form='formatted',iostat=ierr)
	322	if (ierr.eq.0) then
	323	read(91,*) qsurf(iq)
	324	do ilayer=1,nlayermx
	325	read(91,*) q(ilayer,iq)
	326	enddo
	327	else
	328	write(,) "No profile_ccn_mass file!"
	329	endif
	330	close(91)
	331	endif ! of if (txt.eq."ccn_mass")
	332	! CCN NUMBER
	333	if (txt.eq."ccn_number") then
	334	!look for a "profile_ccn_number" input file
	335	open(91,file='profile_ccn_number',status='old',
	336	& form='formatted',iostat=ierr)
	337	if (ierr.eq.0) then
	338	read(91,*) qsurf(iq)
	339	do ilayer=1,nlayermx
	340	read(91,*) q(ilayer,iq)
	341	enddo
	342	else
	343	write(,) "No profile_ccn_number file!"
	344	endif
	345	close(91)
	346	endif ! of if (txt.eq."ccn_number")
[38]	347	enddo ! of do iq=1,nqmx
	348
	349	else
	350	! we still need to set (dummy) tracer names for physdem1
	351	nq=nqmx
	352	do iq=1,nq
	353	write(str7,'(a1,i2.2)')'q',iq
	354	tnom(iq)=str7
	355	enddo
	356	endif ! of if (tracer)
[358]	357
	358	!write(,) "testphys1d q", q(1,:)
	359	!write(,) "testphys1d qsurf", qsurf
[38]	360
	361	c Date et heure locale du debut du run
	362	c ------------------------------------
	363	c Date (en sols depuis le solstice de printemps) du debut du run
	364	day0 = 0 ! default value for day0
	365	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
	366	call getin("day0",day0)
	367	day=float(day0)
	368	write(,) " day0 = ",day0
	369	c Heure de demarrage
	370	time=0 ! default value for time
	371	write(,)'Initial local time (in hours, between 0 and 24)?'
	372	call getin("time",time)
	373	write(,)" time = ",time
	374	time=time/24.E+0 ! convert time (hours) to fraction of sol
	375
	376	c Discretisation (Definition de la grille et des pas de temps)
	377	c --------------
	378	c
	379	nlayer=nlayermx
	380	nlevel=nlayer+1
	381	nsoil=nsoilmx
	382
	383	day_step=48 ! default value for day_step
	384	PRINT *,'Number of time steps per sol ?'
	385	call getin("day_step",day_step)
	386	write(,) " day_step = ",day_step
	387
	388	ndt=10 ! default value for ndt
	389	PRINT *,'Number of sols to run ?'
	390	call getin("ndt",ndt)
	391	write(,) " ndt = ",ndt
	392
	393	ndt=ndt*day_step
	394	dtphys=daysec/day_step
	395	c Pression de surface sur la planete
	396	c ------------------------------------
	397	c
	398	psurf=610. ! default value for psurf
	399	PRINT *,'Surface pressure (Pa) ?'
	400	call getin("psurf",psurf)
	401	write(,) " psurf = ",psurf
	402	c Pression de reference
	403	pa=20.
	404	preff=610.
	405
	406	c Proprietes des poussiere aerosol
	407	c --------------------------------
	408	tauvis=0.2 ! default value for tauvis
[627]	409	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
[38]	410	call getin("tauvis",tauvis)
	411	write(,) " tauvis = ",tauvis
	412
	413	c latitude/longitude
	414	c ------------------
	415	lati(1)=0 ! default value for lati(1)
	416	PRINT *,'latitude (in degrees) ?'
	417	call getin("latitude",lati(1))
	418	write(,) " latitude = ",lati(1)
	419	lati(1)=lati(1)*pi/180.E+0
	420	long(1)=0.E+0
	421	long(1)=long(1)*pi/180.E+0
	422
	423	c Initialisation albedo / inertie du sol
	424	c --------------------------------------
	425	c
	426	albedodat(1)=0.2 ! default value for albedodat
	427	PRINT *,'Albedo of bare ground ?'
	428	call getin("albedo",albedodat(1))
	429	write(,) " albedo = ",albedodat(1)
	430
	431	inertiedat(1,1)=400 ! default value for inertiedat
	432	PRINT *,'Soil thermal inertia (SI) ?'
	433	call getin("inertia",inertiedat(1,1))
	434	write(,) " inertia = ",inertiedat(1,1)
[224]	435
	436	z0(1)=z0_default ! default value for roughness
	437	write(,) 'Surface roughness length z0 (m)?'
	438	call getin("z0",z0(1))
	439	write(,) " z0 = ",z0(1)
[38]	440	c
	441	c pour le schema d'ondes de gravite
	442	c ---------------------------------
	443	c
	444	zmea(1)=0.E+0
	445	zstd(1)=0.E+0
	446	zsig(1)=0.E+0
	447	zgam(1)=0.E+0
	448	zthe(1)=0.E+0
	449
	450
	451
	452	c Initialisation speciales "physiq"
	453	c ---------------------------------
	454	c la surface de chaque maille est inutile en 1D --->
	455	area(1)=1.E+0
	456
	457	c le geopotentiel au sol est inutile en 1D car tout est controle
	458	c par la pression de surface --->
	459	phisfi(1)=0.E+0
	460
	461	c "inifis" reproduit un certain nombre d'initialisations deja faites
	462	c + lecture des clefs de callphys.def
	463	c + calcul de la frequence d'appel au rayonnement
	464	c + calcul des sinus et cosinus des longitude latitude
	465
	466	!Mars possible matter with dtphys in input and include!!!
	467	CALL inifis(1,llm,day0,daysec,dtphys,
	468	. lati,long,area,rad,g,r,cpp)
	469	c Initialisation pour prendre en compte les vents en 1-D
	470	c ------------------------------------------------------
	471	ptif=2.E+0omegsinlat(1)
	472
	473	c vent geostrophique
	474	gru=10. ! default value for gru
	475	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
	476	call getin("u",gru)
	477	write(,) " u = ",gru
	478	grv=0. !default value for grv
	479	PRINT *,'meridional northward component of the geostrophic',
	480	&' wind (m/s) ?'
	481	call getin("v",grv)
	482	write(,) " v = ",grv
	483
	484	c Initialisation des vents au premier pas de temps
	485	DO ilayer=1,nlayer
	486	u(ilayer)=gru
	487	v(ilayer)=grv
	488	ENDDO
	489
	490	c energie cinetique turbulente
	491	DO ilevel=1,nlevel
	492	q2(ilevel)=0.E+0
	493	ENDDO
	494
	495	c glace de CO2 au sol
	496	c -------------------
	497	co2ice=0.E+0 ! default value for co2ice
	498	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
	499	call getin("co2ice",co2ice)
	500	write(,) " co2ice = ",co2ice
	501
	502	c
	503	c emissivite
	504	c ----------
	505	emis=emissiv
	506	IF (co2ice.eq.1.E+0) THEN
	507	emis=emisice(1) ! northern hemisphere
	508	IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
	509	ENDIF
	510
	511
	512
	513	c calcul des pressions et altitudes en utilisant les niveaux sigma
	514	c ----------------------------------------------------------------
	515
	516	c Vertical Coordinates
	517	c """"""""""""""""""""
	518	hybrid=.true.
	519	PRINT *,'Hybrid coordinates ?'
	520	call getin("hybrid",hybrid)
	521	write(,) " hybrid = ", hybrid
	522
	523	CALL disvert
	524
	525	DO ilevel=1,nlevel
	526	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	527	ENDDO
	528
	529	DO ilayer=1,nlayer
	530	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	531	ENDDO
	532
	533	DO ilayer=1,nlayer
	534	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
	535	& /g
	536	ENDDO
	537
	538
	539	c profil de temperature au premier appel
	540	c --------------------------------------
	541	pks=psurf**rcp
	542
	543	c altitude en km dans profile: on divise zlay par 1000
	544	tmp1(0)=0.E+0
	545	DO ilayer=1,nlayer
	546	tmp1(ilayer)=zlay(ilayer)/1000.E+0
	547	ENDDO
	548
	549	call profile(nlayer+1,tmp1,tmp2)
	550
	551	tsurf=tmp2(0)
	552	DO ilayer=1,nlayer
	553	temp(ilayer)=tmp2(ilayer)
	554	ENDDO
	555
	556
	557
	558	! Initialize soil properties and temperature
	559	! ------------------------------------------
	560	volcapa=1.e6 ! volumetric heat capacity
	561	DO isoil=1,nsoil
	562	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
	563	tsoil(isoil)=tsurf ! soil temperature
	564	ENDDO
	565
	566	! Initialize depths
	567	! -----------------
	568	do isoil=0,nsoil-1
	569	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
	570	enddo
	571	do isoil=1,nsoil
	572	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
	573	enddo
	574
	575	c Initialisation des traceurs
	576	c ---------------------------
	577
	578	if (photochem.or.callthermos) then
	579	write(,) 'Especes chimiques initialisees'
	580	! thermo=0: initialisation pour toutes les couches
	581	thermo=0
	582	call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
	583	$ thermo,qsurf)
	584	endif
[520]	585
	586	c Regarder si le sol est un reservoir de glace d'eau
	587	c --------------------------------------------------
	588	watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol
	589	print *,'Water ice cap on ground ?'
	590	call getin("watercaptag",watercaptag)
	591	write(,) " watercaptag = ",watercaptag
[38]	592
	593
	594	c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
	595	c ----------------------------------------------------------------
	596	c (effectuee avec "writeg1d" a partir de "physiq.F"
	597	c ou tout autre subroutine physique, y compris celle ci).
	598
	599	g1d_nlayer=nlayer
	600	g1d_nomfich='g1d.dat'
	601	g1d_unitfich=40
	602	g1d_nomctl='g1d.ctl'
	603	g1d_unitctl=41
	604	g1d_premier=.true.
	605	g2d_premier=.true.
	606
	607	c Ecriture de "startfi"
	608	c --------------------
	609	c (Ce fichier sera aussitot relu au premier
	610	c appel de "physiq", mais il est necessaire pour passer
	611	c les variables purement physiques a "physiq"...
	612
	613	call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
[283]	614	. dtphys,float(day0),time,tsurf,
	615	. tsoil,co2ice,emis,q2,qsurf,area,albedodat,
	616	. inertiedat,zmea,zstd,zsig,zgam,zthe)
[38]	617	c=======================================================================
	618	c BOUCLE TEMPORELLE DU MODELE 1D
	619	c=======================================================================
	620	c
	621	firstcall=.true.
	622	lastcall=.false.
	623
	624	DO idt=1,ndt
	625	c IF (idt.eq.ndt) lastcall=.true.
	626	IF (idt.eq.ndt-day_step-1) then !test
	627	lastcall=.true.
	628	call solarlong(day*1.0,zls)
	629	write(103,) 'Ls=',zls180./pi
	630	write(103,) 'Lat=', lati(1)180./pi
[627]	631	write(103,) 'Tau=', tauvis/odprefpsurf
[38]	632	write(103,*) 'RunEnd - Atmos. Temp. File'
	633	write(103,*) 'RunEnd - Atmos. Temp. File'
	634	write(104,) 'Ls=',zls180./pi
	635	write(104,*) 'Lat=', lati(1)
[627]	636	write(104,) 'Tau=', tauvis/odprefpsurf
[38]	637	write(104,*) 'RunEnd - Atmos. Temp. File'
	638	ENDIF
	639
	640	c calcul du geopotentiel
	641	c ~~~~~~~~~~~~~~~~~~~~~
	642	DO ilayer=1,nlayer
	643	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	644	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	645	ENDDO
	646	phi(1)=pksh(1)(1.E+0-s(1))
	647	DO ilayer=2,nlayer
	648	phi(ilayer)=phi(ilayer-1)+
	649	& pks(h(ilayer-1)+h(ilayer)).5E+0
	650	& *(s(ilayer-1)-s(ilayer))
	651
	652	ENDDO
	653
	654	c appel de la physique
	655	c --------------------
[358]	656	! write(,) "testphys1d avant q", q(1,:)
[38]	657	CALL physiq (1,llm,nqmx,
	658	, firstcall,lastcall,
	659	, day,time,dtphys,
	660	, plev,play,phi,
	661	, u, v,temp, q,
	662	, w,
	663	C - sorties
	664	s du, dv, dtemp, dq,dpsurf,tracerdyn)
[358]	665	! write(,) "testphys1d apres q", q(1,:)
[38]	666
[358]	667
[38]	668	c evolution du vent : modele 1D
	669	c -----------------------------
	670
	671	c la physique calcule les derivees temporelles de u et v.
	672	c on y rajoute betement un effet Coriolis.
	673	c
	674	c DO ilayer=1,nlayer
	675	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
	676	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
	677	c ENDDO
	678
	679	c Pour certain test : pas de coriolis a l'equateur
	680	c if(lati(1).eq.0.) then
	681	DO ilayer=1,nlayer
	682	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
	683	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
	684	ENDDO
	685	c end if
	686	c
	687	c
	688	c Calcul du temps au pas de temps suivant
	689	c ---------------------------------------
	690	firstcall=.false.
	691	time=time+dtphys/daysec
	692	IF (time.gt.1.E+0) then
	693	time=time-1.E+0
	694	day=day+1
	695	ENDIF
	696
	697	c calcul des vitesses et temperature au pas de temps suivant
	698	c ----------------------------------------------------------
	699
	700	DO ilayer=1,nlayer
	701	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
	702	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
	703	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
	704	ENDDO
	705
	706	c calcul des pressions au pas de temps suivant
	707	c ----------------------------------------------------------
	708
	709	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
	710	DO ilevel=1,nlevel
	711	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	712	ENDDO
	713	DO ilayer=1,nlayer
	714	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	715	ENDDO
	716
	717	! increment tracer values
	718	DO iq = 1, nqmx
	719	DO ilayer=1,nlayer
	720	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
	721	ENDDO
	722	ENDDO
	723
	724	ENDDO ! fin de la boucle temporelle
	725
	726	c ========================================================
	727	c GESTION DES SORTIE
	728	c ========================================================
	729
	730	c fermeture pour conclure les sorties format grads dans "g1d.dat"
	731	c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
	732	c ou tout autre subroutine physique, y compris celle ci).
	733
	734	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
	735	CALL endg1d(1,nlayer,zlay/1000.,ndt)
	736
	737	c ========================================================
	738	END
	739
	740	c***********************************************************************
	741	c***********************************************************************
	742	c Subroutines Bidons utilise seulement en 3D, mais
	743	c necessaire a la compilation de testphys1d en 1-D
	744
	745	subroutine gr_fi_dyn
	746	RETURN
	747	END
	748
	749	c***********************************************************************
	750	c***********************************************************************
	751
	752	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format