Context Navigation

testphys1d.F @ 2347

Last change on this file since 2347 was 2332, checked in by mvals, 5 years ago

Mars GCM:
Follow-up of the last commit for the transport of the isotopic ratio: simplification of the transmission of variables from the dynamics to the
physics:

libf/dynphy_lonlat/phymars/: iniphysiq_mod.F90: transmission of the content of 2 variables describing the isotopes instead of 4 (nqperes: number of tracers "peres", nqfils:

number of tracers "fils")

libf/phymars/: phys_state_var_init_mod.F90, tracer_mod.F: idem callsedim_mod.F: idem co2condens_mod.F: idem
libf/phymars/dyn1d: testphys1d.F: idem (the reading interface for traceur.def has been completed to fill the variables nqperes and nqfils).

File size: 35.5 KB

Rev	Line
[38]	1
	2	PROGRAM testphys1d
	3	! to use 'getin'
[1036]	4	USE ioipsl_getincom, only: getin
[1543]	5	use dimphy, only : init_dimphy
	6	use mod_grid_phy_lmdz, only : regular_lonlat
[2332]	7	use infotrac, only: nqtot, tname, nqperes,nqfils
[1130]	8	use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx
[1541]	9	use comgeomfi_h, only: sinlat, ini_fillgeom
[1047]	10	use surfdat_h, only: albedodat, z0_default, emissiv, emisice,
	11	& albedice, iceradius, dtemisice, z0,
	12	& zmea, zstd, zsig, zgam, zthe, phisfi,
[2293]	13	& watercaptag, watercap, hmons, summit, base
[1047]	14	use slope_mod, only: theta_sl, psi_sl
[1130]	15	use phyredem, only: physdem0,physdem1
[1543]	16	use geometry_mod, only: init_geometry
[1229]	17	use planete_h, only: year_day, periheli, aphelie, peri_day,
	18	& obliquit, emin_turb, lmixmin
[1524]	19	use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp
[1535]	20	use time_phylmdz_mod, only: daysec, dtphys, day_step,
	21	& ecritphy, iphysiq
[1919]	22	use dimradmars_mod, only: tauscaling,tauvis,totcloudfrac
[1944]	23	use co2cloud_mod, only: mem_Mccn_co2, mem_Mh2o_co2,
	24	& mem_Nccn_co2
[1621]	25	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig,
	26	& presnivs,pseudoalt,scaleheight
	27	USE vertical_layers_mod, ONLY: init_vertical_layers
[1462]	28	USE logic_mod, ONLY: hybrid
[1543]	29	use physics_distribution_mod, only: init_physics_distribution
[1535]	30	use regular_lonlat_mod, only: init_regular_lonlat
[1543]	31	use mod_interface_dyn_phys, only: init_interface_dyn_phys
[1524]	32	USE phys_state_var_init_mod, ONLY: phys_state_var_init
[1549]	33	USE physiq_mod, ONLY: physiq
[38]	34	IMPLICIT NONE
	35
	36	c=======================================================================
	37	c subject:
	38	c --------
	39	c PROGRAM useful to run physical part of the martian GCM in a 1D column
	40	c
	41	c Can be compiled with a command like (e.g. for 25 layers)
	42	c "makegcm -p mars -d 25 testphys1d"
	43	c It requires the files "testphys1d.def" "callphys.def"
	44	c and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
	45	c and a file describing the sigma layers (e.g. "z2sig.def")
	46	c
	47	c author: Frederic Hourdin, R.Fournier,F.Forget
	48	c -------
	49	c
	50	c update: 12/06/2003 including chemistry (S. Lebonnois)
	51	c and water ice (F. Montmessin)
	52	c
	53	c=======================================================================
	54
[1621]	55	include "dimensions.h"
[1130]	56	integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm)
[1266]	57	integer, parameter :: nlayer = llm
[1047]	58	!#include "dimradmars.h"
	59	!#include "comgeomfi.h"
	60	!#include "surfdat.h"
	61	!#include "slope.h"
	62	!#include "comsoil.h"
	63	!#include "comdiurn.h"
[1621]	64	include "callkeys.h"
[1047]	65	!#include "comsaison.h"
[1130]	66	!#include "control.h"
[1621]	67	include "netcdf.inc"
	68	include "comg1d.h"
[1036]	69	!#include "advtrac.h"
[38]	70
	71	c --------------------------------------------------------------
	72	c Declarations
	73	c --------------------------------------------------------------
	74	c
	75	INTEGER unitstart ! unite d'ecriture de "startfi"
[1273]	76	INTEGER nlevel,nsoil,ndt
[38]	77	INTEGER ilayer,ilevel,isoil,idt,iq
	78	LOGICAl firstcall,lastcall
	79	c
[627]	80	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
	81	c
[38]	82	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
	83	REAL day ! date durant le run
	84	REAL time ! time (0<time<1 ; time=0.5 a midi)
[1266]	85	REAL play(nlayer) ! Pressure at the middle of the layers (Pa)
	86	REAL plev(nlayer+1) ! intermediate pressure levels (pa)
[1130]	87	REAL psurf,tsurf(1)
[1266]	88	REAL u(nlayer),v(nlayer) ! zonal, meridional wind
[38]	89	REAL gru,grv ! prescribed "geostrophic" background wind
[1266]	90	REAL temp(nlayer) ! temperature at the middle of the layers
[1036]	91	REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
	92	REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
[38]	93	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
[1130]	94	REAL co2ice(1) ! co2ice layer (kg.m-2)
	95	REAL emis(1) ! surface layer
[1944]	96	REAL albedo(1,1) ! surface albedo
	97	REAL :: wstar(1)=0. ! Thermals vertical velocity
[1266]	98	REAL q2(nlayer+1) ! Turbulent Kinetic Energy
	99	REAL zlay(nlayer) ! altitude estimee dans les couches (km)
[38]	100
	101	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
[1266]	102	REAL du(nlayer),dv(nlayer),dtemp(nlayer)
	103	REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer)
[1549]	104	REAL dpsurf(1)
[1036]	105	REAL,ALLOCATABLE :: dq(:,:)
	106	REAL,ALLOCATABLE :: dqdyn(:,:)
[38]	107
	108	c Various intermediate variables
	109	INTEGER thermo
	110	REAL zls
[1266]	111	REAL phi(nlayer),h(nlayer),s(nlayer)
	112	REAL pks, ptif, w(nlayer)
[1036]	113	REAL qtotinit,qtot
	114	real,allocatable :: mqtot(:)
[38]	115	INTEGER ierr, aslun
[1266]	116	REAL tmp1(0:nlayer),tmp2(0:nlayer)
[38]	117	integer :: nq=1 ! number of tracers
[1543]	118	real :: latitude(1), longitude(1), cell_area(1)
[2312]	119	integer iqh2ovap
	120	integer iqh2oice
[38]	121
	122	character*2 str2
	123	character (len=7) :: str7
	124	character(len=44) :: txt
	125
[2228]	126	c New flag to compute paleo orbital configurations + few variables JN
	127	REAL halfaxe, excentric, Lsperi
	128	Logical paleomars
	129
[2322]	130	c MVals: isotopes as in the dynamics (CRisi)
	131	INTEGER :: ifils,ipere,generation
	132	CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name
	133	CHARACTER(len=80) :: line ! to store a line of text
	134	INTEGER ierr0
	135	LOGICAL :: continu
	136
[38]	137	c=======================================================================
	138
	139	c=======================================================================
	140	c INITIALISATION
	141	c=======================================================================
[1130]	142	! initialize "serial/parallel" related stuff
	143	! CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
[1543]	144	! CALL init_phys_lmdz(1,1,llm,1,(/1/))
	145	! call initcomgeomphy
[38]	146
	147	c ------------------------------------------------------
[899]	148	c Prescribed constants to be set here
[38]	149	c ------------------------------------------------------
	150
	151	pi=2.E+0*asin(1.E+0)
	152
[899]	153	c Mars planetary constants
[38]	154	c ----------------------------
[899]	155	rad=3397200. ! mars radius (m) ~3397200 m
	156	daysec=88775. ! length of a sol (s) ~88775 s
	157	omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1)
	158	g=3.72 ! gravity (m.s-2) ~3.72
	159	mugaz=43.49 ! atmosphere mola mass (g.mol-1) ~43.49
[38]	160	rcp=.256793 ! = r/cp ~0.256793
	161	r= 8.314511E+0 *1000.E+0/mugaz
	162	cpp= r/rcp
[899]	163	year_day = 669 ! lenght of year (sols) ~668.6
	164	periheli = 206.66 ! minimum sun-mars distance (Mkm) ~206.66
	165	aphelie = 249.22 ! maximum sun-mars distance (Mkm) ~249.22
[2228]	166	halfaxe = 227.94 ! demi-grand axe de l'ellipse
[899]	167	peri_day = 485. ! perihelion date (sols since N. Spring)
	168	obliquit = 25.2 ! Obliquity (deg) ~25.2
[2228]	169	excentric = 0.0934 ! Eccentricity (0.0934)
[38]	170
[899]	171	c Planetary Boundary Layer and Turbulence parameters
	172	c --------------------------------------------------
	173	z0_default = 1.e-2 ! surface roughness (m) ~0.01
	174	emin_turb = 1.e-6 ! minimal turbulent energy ~1.e-8
	175	lmixmin = 30 ! mixing length ~100
[38]	176
[899]	177	c cap properties and surface emissivities
[38]	178	c ----------------------------------------------------
[899]	179	emissiv= 0.95 ! Bare ground emissivity ~.95
	180	emisice(1)=0.95 ! Northern cap emissivity
	181	emisice(2)=0.95 ! Southern cap emisssivity
	182	albedice(1)=0.5 ! Northern cap albedo
	183	albedice(2)=0.5 ! Southern cap albedo
[38]	184	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
	185	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
	186	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
	187	dtemisice(2) = 2. ! time scale for snow metamorphism (south
	188
[1920]	189	c mesh surface (not a very usefull quantity in 1D)
	190	c ----------------------------------------------------
	191	cell_area(1)=1.E+0
	192
[38]	193	c ------------------------------------------------------
[899]	194	c Loading run parameters from "run.def" file
[38]	195	c ------------------------------------------------------
	196
	197
	198	! check if 'run.def' file is around (otherwise reading parameters
	199	! from callphys.def via getin() routine won't work.
	200	open(99,file='run.def',status='old',form='formatted',
	201	& iostat=ierr)
	202	if (ierr.ne.0) then
	203	write(,) 'Cannot find required file "run.def"'
	204	write(,) ' (which should contain some input parameters'
	205	write(,) ' along with the following line:'
	206	write(,) ' INCLUDEDEF=callphys.def'
	207	write(,) ' )'
	208	write(,) ' ... might as well stop here ...'
	209	stop
	210	else
	211	close(99)
	212	endif
	213
	214	! check if we are going to run with or without tracers
	215	write(,) "Run with or without tracer transport ?"
	216	tracer=.false. ! default value
	217	call getin("tracer",tracer)
	218	write(,) " tracer = ",tracer
	219
	220	! while we're at it, check if there is a 'traceur.def' file
[1036]	221	! and process it.
[38]	222	if (tracer) then
	223	! load tracer names from file 'traceur.def'
	224	open(90,file='traceur.def',status='old',form='formatted',
	225	& iostat=ierr)
	226	if (ierr.ne.0) then
	227	write(,) 'Cannot find required file "traceur.def"'
	228	write(,) ' If you want to run with tracers, I need it'
	229	write(,) ' ... might as well stop here ...'
	230	stop
	231	else
	232	write(,) "testphys1d: Reading file traceur.def"
	233	! read number of tracers:
	234	read(90,*,iostat=ierr) nq
[1036]	235	nqtot=nq ! set value of nqtot (in infotrac module) as nq
[38]	236	if (ierr.ne.0) then
	237	write(,) "testphys1d: error reading number of tracers"
	238	write(,) " (first line of traceur.def) "
	239	stop
	240	endif
	241	endif
[1036]	242	! allocate arrays:
[1130]	243	allocate(tname(nq))
[1266]	244	allocate(q(nlayer,nq))
[1036]	245	allocate(qsurf(nq))
[1266]	246	allocate(dq(nlayer,nq))
	247	allocate(dqdyn(nlayer,nq))
[1036]	248	allocate(mqtot(nq))
[2322]	249	allocate(tnom_transp(nq))
[1036]	250
[38]	251	! read tracer names from file traceur.def
[1036]	252	do iq=1,nq
[2322]	253	read(90,'(80a)',iostat=ierr) line ! store the line from traceur.def
[38]	254	if (ierr.ne.0) then
	255	write(,) 'testphys1d: error reading tracer names...'
	256	stop
	257	endif
[2322]	258	! if format is tnom_0, tnom_transp (isotopes)
	259	read(line,*,iostat=ierr0) tname(iq),tnom_transp(iq)
	260	if (ierr0.ne.0) then
	261	read(line,*) tname(iq)
	262	tnom_transp(iq)='air'
	263	endif
	264
[38]	265	enddo
	266	close(90)
	267
[2322]	268	! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
[2332]	269	ALLOCATE(nqfils(nqtot))
[2322]	270	nqperes=0
[2332]	271	nqfils(:)=0
[2322]	272	DO iq=1,nqtot
	273	if (tnom_transp(iq) == 'air') then
	274	! ceci est un traceur père
	275	WRITE(,) 'Le traceur',iq,', appele ',
	276	& trim(tname(iq)),', est un pere'
	277	nqperes=nqperes+1
	278	else !if (tnom_transp(iq) == 'air') then
	279	! ceci est un fils. Qui est son père?
	280	WRITE(,) 'Le traceur',iq,', appele ',
	281	& trim(tname(iq)),', est un fils'
	282	continu=.true.
	283	ipere=1
	284	do while (continu)
	285	if (tnom_transp(iq) .eq. tname(ipere)) then
	286	! Son père est ipere
	287	WRITE(,) 'Le traceur',iq,'appele ',
	288	& trim(tname(iq)),' est le fils de ',
	289	& ipere,'appele ',trim(tname(ipere))
[2332]	290	nqfils(ipere)=nqfils(ipere)+1
[2322]	291	continu=.false.
	292	else !if (tnom_transp(iq) == tnom_0(ipere)) then
	293	ipere=ipere+1
	294	if (ipere.gt.nqtot) then
	295	WRITE(,) 'Le traceur',iq,'appele ',
	296	& trim(tname(iq)),', est orpelin.'
	297	CALL abort_gcm('infotrac_init',
	298	& 'Un traceur est orphelin',1)
	299	endif !if (ipere.gt.nqtot) then
	300	endif !if (tnom_transp(iq) == tnom_0(ipere)) then
	301	enddo !do while (continu)
	302	endif !if (tnom_transp(iq) == 'air') then
	303	enddo !DO iq=1,nqtot
	304	WRITE(,) 'nqperes=',nqperes
	305	WRITE(,) 'nqfils=',nqfils
	306
[38]	307	! initialize tracers here:
	308	write(,) "testphys1d: initializing tracers"
	309	q(:,:)=0 ! default, set everything to zero
	310	qsurf(:)=0
	311	! "smarter" initialization of some tracers
	312	! (get values from "profile_*" files, if these are available)
[1036]	313	do iq=1,nq
[38]	314	txt=""
[1130]	315	write(txt,"(a)") tname(iq)
[38]	316	write(,)" tracer:",trim(txt)
	317	! CO2
	318	if (txt.eq."co2") then
	319	q(:,iq)=0.95 ! kg /kg of atmosphere
	320	qsurf(iq)=0. ! kg/m2 (not used for CO2)
	321	! even better, look for a "profile_co2" input file
	322	open(91,file='profile_co2',status='old',
	323	& form='formatted',iostat=ierr)
	324	if (ierr.eq.0) then
	325	read(91,*) qsurf(iq)
[1266]	326	do ilayer=1,nlayer
[38]	327	read(91,*) q(ilayer,iq)
	328	enddo
	329	endif
	330	close(91)
	331	endif ! of if (txt.eq."co2")
[161]	332	! Allow for an initial profile of argon
	333	! Can also be used to introduce a decaying tracer
	334	! in the 1D (TBD) to study thermals
	335	if (txt.eq."ar") then
	336	!look for a "profile_ar" input file
	337	open(91,file='profile_ar',status='old',
	338	& form='formatted',iostat=ierr)
	339	if (ierr.eq.0) then
	340	read(91,*) qsurf(iq)
[1266]	341	do ilayer=1,nlayer
[161]	342	read(91,*) q(ilayer,iq)
	343	enddo
	344	else
	345	write(,) "No profile_ar file!"
	346	endif
	347	close(91)
	348	endif ! of if (txt.eq."ar")
	349
[38]	350	! WATER VAPOUR
	351	if (txt.eq."h2o_vap") then
[2312]	352	iqh2ovap=iq !remember index for water vap
[38]	353	!look for a "profile_h2o_vap" input file
	354	open(91,file='profile_h2o_vap',status='old',
	355	& form='formatted',iostat=ierr)
	356	if (ierr.eq.0) then
	357	read(91,*) qsurf(iq)
[1266]	358	do ilayer=1,nlayer
[38]	359	read(91,*) q(ilayer,iq)
	360	enddo
	361	else
	362	write(,) "No profile_h2o_vap file!"
	363	endif
	364	close(91)
	365	endif ! of if (txt.eq."h2o_ice")
	366	! WATER ICE
	367	if (txt.eq."h2o_ice") then
[2312]	368	iqh2oice = iq !remember index for water ice
[38]	369	!look for a "profile_h2o_vap" input file
	370	open(91,file='profile_h2o_ice',status='old',
	371	& form='formatted',iostat=ierr)
	372	if (ierr.eq.0) then
	373	read(91,*) qsurf(iq)
[1266]	374	do ilayer=1,nlayer
[38]	375	read(91,*) q(ilayer,iq)
	376	enddo
	377	else
	378	write(,) "No profile_h2o_ice file!"
	379	endif
	380	close(91)
	381	endif ! of if (txt.eq."h2o_ice")
[2312]	382
	383	! HDO VAPOUR
	384	if (txt.eq."hdo_vap") then
	385	!look for a "profile_hdo_vap" input file
	386	open(91,file='profile_hdo_vap',status='old',
	387	& form='formatted',iostat=ierr)
	388	if (ierr.eq.0) then
	389	read(91,*) qsurf(iq)
	390	do ilayer=1,nlayer
	391	read(91,*) q(ilayer,iq)
	392	enddo
	393	else
	394	write(,) "No profile_hdo_vap file!"
	395	do ilayer=1,nlayer
	396	q(ilayer,iq) = q(ilayer,iqh2ovap)2155.76e-6*5
	397	enddo
	398	endif
	399	close(91)
	400	endif ! of if (txt.eq."hdo_ice")
	401	! HDO ICE
	402	if (txt.eq."hdo_ice") then
	403	!look for a "profile_hdo_vap" input file
	404	open(91,file='profile_hdo_ice',status='old',
	405	& form='formatted',iostat=ierr)
	406	if (ierr.eq.0) then
	407	read(91,*) qsurf(iq)
	408	do ilayer=1,nlayer
	409	read(91,*) q(ilayer,iq)
	410	enddo
	411	else
	412	write(,) "No profile_hdo_ice file!"
	413	qsurf(iq) = qsurf(iqh2oice) * 2155.76e-65
	414	do ilayer=1,nlayer
	415	q(ilayer,iq) = q(ilayer,iqh2oice) * 2155.76e-65
	416	enddo
	417	endif
	418	close(91)
	419	endif ! of if (txt.eq."hdo_ice")
	420
	421
[38]	422	! DUST
	423	!if (txt(1:4).eq."dust") then
	424	! q(:,iq)=0.4 ! kg/kg of atmosphere
	425	! qsurf(iq)=100 ! kg/m2
	426	!endif
	427	! DUST MMR
	428	if (txt.eq."dust_mass") then
	429	!look for a "profile_dust_mass" input file
	430	open(91,file='profile_dust_mass',status='old',
	431	& form='formatted',iostat=ierr)
	432	if (ierr.eq.0) then
	433	read(91,*) qsurf(iq)
[1266]	434	do ilayer=1,nlayer
[38]	435	read(91,*) q(ilayer,iq)
	436	! write(,) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
	437	enddo
	438	else
	439	write(,) "No profile_dust_mass file!"
	440	endif
	441	close(91)
	442	endif ! of if (txt.eq."dust_mass")
	443	! DUST NUMBER
	444	if (txt.eq."dust_number") then
	445	!look for a "profile_dust_number" input file
	446	open(91,file='profile_dust_number',status='old',
	447	& form='formatted',iostat=ierr)
	448	if (ierr.eq.0) then
	449	read(91,*) qsurf(iq)
[1266]	450	do ilayer=1,nlayer
[38]	451	read(91,*) q(ilayer,iq)
	452	enddo
	453	else
	454	write(,) "No profile_dust_number file!"
	455	endif
	456	close(91)
	457	endif ! of if (txt.eq."dust_number")
[358]	458	! NB: some more initializations (chemistry) is done later
	459	! CCN MASS
	460	if (txt.eq."ccn_mass") then
	461	!look for a "profile_ccn_mass" input file
	462	open(91,file='profile_ccn_mass',status='old',
	463	& form='formatted',iostat=ierr)
	464	if (ierr.eq.0) then
	465	read(91,*) qsurf(iq)
[1266]	466	do ilayer=1,nlayer
[358]	467	read(91,*) q(ilayer,iq)
	468	enddo
	469	else
	470	write(,) "No profile_ccn_mass file!"
	471	endif
	472	close(91)
	473	endif ! of if (txt.eq."ccn_mass")
	474	! CCN NUMBER
	475	if (txt.eq."ccn_number") then
	476	!look for a "profile_ccn_number" input file
	477	open(91,file='profile_ccn_number',status='old',
	478	& form='formatted',iostat=ierr)
	479	if (ierr.eq.0) then
	480	read(91,*) qsurf(iq)
[1266]	481	do ilayer=1,nlayer
[358]	482	read(91,*) q(ilayer,iq)
	483	enddo
	484	else
	485	write(,) "No profile_ccn_number file!"
	486	endif
	487	close(91)
	488	endif ! of if (txt.eq."ccn_number")
[1036]	489	enddo ! of do iq=1,nq
[38]	490
	491	else
[1036]	492	! we still need to set (dummy) tracer number and names for physdem1
	493	nq=1
[1380]	494	nqtot=nq ! set value of nqtot (in infotrac module) as nq
[1036]	495	! allocate arrays:
[1130]	496	allocate(tname(nq))
[1266]	497	allocate(q(nlayer,nq))
[1036]	498	allocate(qsurf(nq))
[1266]	499	allocate(dq(nlayer,nq))
	500	allocate(dqdyn(nlayer,nq))
[1036]	501	allocate(mqtot(nq))
[38]	502	do iq=1,nq
[1380]	503	write(str7,'(a1,i2.2)')'t',iq
[1130]	504	tname(iq)=str7
[38]	505	enddo
[899]	506	! and just to be clean, also initialize tracers to zero for physdem1
	507	q(:,:)=0
	508	qsurf(:)=0
[38]	509	endif ! of if (tracer)
[358]	510
	511	!write(,) "testphys1d q", q(1,:)
	512	!write(,) "testphys1d qsurf", qsurf
[38]	513
[899]	514	c Date and local time at beginning of run
	515	c ---------------------------------------
	516	c Date (in sols since spring solstice) at beginning of run
[38]	517	day0 = 0 ! default value for day0
	518	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
	519	call getin("day0",day0)
	520	day=float(day0)
	521	write(,) " day0 = ",day0
[899]	522	c Local time at beginning of run
[38]	523	time=0 ! default value for time
	524	write(,)'Initial local time (in hours, between 0 and 24)?'
	525	call getin("time",time)
	526	write(,)" time = ",time
	527	time=time/24.E+0 ! convert time (hours) to fraction of sol
	528
[899]	529	c Discretization (Definition of grid and time steps)
[38]	530	c --------------
	531	c
	532	nlevel=nlayer+1
	533	nsoil=nsoilmx
	534
	535	day_step=48 ! default value for day_step
	536	PRINT *,'Number of time steps per sol ?'
	537	call getin("day_step",day_step)
	538	write(,) " day_step = ",day_step
	539
[1535]	540	ecritphy=day_step ! default value for ecritphy, output every time step
	541
[38]	542	ndt=10 ! default value for ndt
	543	PRINT *,'Number of sols to run ?'
	544	call getin("ndt",ndt)
	545	write(,) " ndt = ",ndt
	546
	547	ndt=ndt*day_step
	548	dtphys=daysec/day_step
[899]	549
	550	c Imposed surface pressure
[38]	551	c ------------------------------------
	552	c
	553	psurf=610. ! default value for psurf
	554	PRINT *,'Surface pressure (Pa) ?'
	555	call getin("psurf",psurf)
	556	write(,) " psurf = ",psurf
[899]	557	c Reference pressures
	558	pa=20. ! transition pressure (for hybrid coord.)
	559	preff=610. ! reference surface pressure
[38]	560
[899]	561	c Aerosol properties
[38]	562	c --------------------------------
[899]	563	tauvis=0.2 ! default value for tauvis (dust opacity)
[627]	564	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
[38]	565	call getin("tauvis",tauvis)
	566	write(,) " tauvis = ",tauvis
[720]	567
	568	c Orbital parameters
	569	c ------------------
[2228]	570	paleomars=.false. ! Default: no water ice reservoir
	571	call getin("paleomars",paleomars)
[2293]	572	if (paleomars.eqv..true.) then
[2228]	573	write(,) "paleomars=", paleomars
	574	write(,) "Orbital parameters from callphys.def"
	575	write(,) "Enter eccentricity & Lsperi"
	576	print *, 'Martian eccentricity (0<e<1) ?'
	577	call getin('excentric ',excentric)
	578	write(,)"excentric =",excentric
	579	print *, 'Solar longitude of perihelion (0<Ls<360) ?'
	580	call getin('Lsperi',Lsperi )
	581	write(,)"Lsperi=",Lsperi
	582	Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
	583	periheli = halfaxe*(1-excentric)
	584	aphelie = halfaxe*(1+excentric)
	585	call call_dayperi(Lsperi,excentric,peri_day,year_day)
	586	write(,) "Corresponding orbital params for GCM"
	587	write(,) " periheli = ",periheli
	588	write(,) " aphelie = ",aphelie
	589	write(,) "date of perihelion (sol)",peri_day
	590	else
	591	write(,) "paleomars=", paleomars
	592	write(,) "Default present-day orbital parameters"
	593	write(,) "Unless specified otherwise"
	594	print *,'Min. distance Sun-Mars (Mkm)?'
	595	call getin("periheli",periheli)
	596	write(,) " periheli = ",periheli
[720]	597
[2228]	598	print *,'Max. distance Sun-Mars (Mkm)?'
	599	call getin("aphelie",aphelie)
	600	write(,) " aphelie = ",aphelie
[720]	601
[2228]	602	print *,'Day of perihelion?'
	603	call getin("periday",peri_day)
	604	write(,) " periday = ",peri_day
[720]	605
[2228]	606	print *,'Obliquity?'
	607	call getin("obliquit",obliquit)
	608	write(,) " obliquit = ",obliquit
	609	end if
[38]	610
	611	c latitude/longitude
	612	c ------------------
[1311]	613	latitude(1)=0 ! default value for latitude
[38]	614	PRINT *,'latitude (in degrees) ?'
[1311]	615	call getin("latitude",latitude(1))
[1047]	616	write(,) " latitude = ",latitude
	617	latitude=latitude*pi/180.E+0
	618	longitude=0.E+0
	619	longitude=longitude*pi/180.E+0
[38]	620
[1224]	621	! some initializations (some of which have already been
[1047]	622	! done above!) and loads parameters set in callphys.def
	623	! and allocates some arrays
	624	!Mars possible matter with dtphys in input and include!!!
[1535]	625	! Initializations below should mimick what is done in iniphysiq for 3D GCM
[1543]	626	call init_physics_distribution(regular_lonlat,4,
	627	& 1,1,1,nlayer,1)
	628	call init_interface_dyn_phys
[1535]	629	call init_regular_lonlat(1,1,longitude,latitude,
	630	& (/0.,0./),(/0.,0./))
[1543]	631	call init_geometry(1,longitude,latitude,
	632	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
	633	& cell_area)
[1621]	634	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
	635	! call init_vertical_layers(nlayer,preff,scaleheight,
	636	! & ap,bp,aps,bps,presnivs,pseudoalt)
[1543]	637	call init_dimphy(1,nlayer) ! Initialize dimphy module
[1621]	638	call phys_state_var_init(1,llm,nq,tname,
[2322]	639	. day0,time,daysec,dtphys,rad,g,r,cpp,
[2332]	640	. nqperes,nqfils)! MVals: variables isotopes
[1311]	641	call ini_fillgeom(1,latitude,longitude,(/1.0/))
[1246]	642	call conf_phys(1,llm,nq)
[1047]	643
[1550]	644	! in 1D model physics are called every time step
	645	! ovverride iphysiq value that has been set by conf_phys
	646	if (iphysiq/=1) then
	647	write(,) "testphys1d: setting iphysiq=1"
	648	iphysiq=1
	649	endif
[1047]	650
[899]	651	c Initialize albedo / soil thermal inertia
	652	c ----------------------------------------
[38]	653	c
	654	albedodat(1)=0.2 ! default value for albedodat
	655	PRINT *,'Albedo of bare ground ?'
	656	call getin("albedo",albedodat(1))
	657	write(,) " albedo = ",albedodat(1)
[1944]	658	albedo(1,1)=albedodat(1)
[38]	659
	660	inertiedat(1,1)=400 ! default value for inertiedat
	661	PRINT *,'Soil thermal inertia (SI) ?'
	662	call getin("inertia",inertiedat(1,1))
	663	write(,) " inertia = ",inertiedat(1,1)
[224]	664
	665	z0(1)=z0_default ! default value for roughness
	666	write(,) 'Surface roughness length z0 (m)?'
	667	call getin("z0",z0(1))
	668	write(,) " z0 = ",z0(1)
[899]	669
	670	! Initialize local slope parameters (only matters if "callslope"
	671	! is .true. in callphys.def)
	672	! slope inclination angle (deg) 0: horizontal, 90: vertical
	673	theta_sl(1)=0.0 ! default: no inclination
	674	call getin("slope_inclination",theta_sl(1))
	675	! slope orientation (deg)
	676	! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
	677	psi_sl(1)=0.0 ! default value
	678	call getin("slope_orientation",psi_sl(1))
	679
[38]	680	c
[899]	681	c for the gravity wave scheme
[38]	682	c ---------------------------------
	683	c
	684	zmea(1)=0.E+0
	685	zstd(1)=0.E+0
	686	zsig(1)=0.E+0
	687	zgam(1)=0.E+0
[2199]	688	zthe(1)=0.E+0
	689	c
	690	c for the slope wind scheme
	691	c ---------------------------------
	692	c
[1974]	693	hmons(1)=0.E+0
[2199]	694	PRINT *,'hmons is initialized to ',hmons(1)
[2079]	695	summit(1)=0.E+0
[2199]	696	PRINT *,'summit is initialized to ',summit(1)
[2079]	697	base(1)=0.E+0
[2199]	698	c
	699	c Default values initializing the coefficients calculated later
	700	c ---------------------------------
	701	c
	702	tauscaling(1)=1. ! calculated in aeropacity_mod.F
	703	totcloudfrac(1)=1. ! calculated in watercloud_mod.F
	704
[899]	705	c Specific initializations for "physiq"
	706	c -------------------------------------
	707	c surface geopotential is not used (or useful) since in 1D
	708	c everything is controled by surface pressure
[38]	709	phisfi(1)=0.E+0
	710
[899]	711	c Initialization to take into account prescribed winds
[38]	712	c ------------------------------------------------------
	713	ptif=2.E+0omegsinlat(1)
	714
[899]	715	c geostrophic wind
[38]	716	gru=10. ! default value for gru
	717	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
	718	call getin("u",gru)
	719	write(,) " u = ",gru
	720	grv=0. !default value for grv
	721	PRINT *,'meridional northward component of the geostrophic',
	722	&' wind (m/s) ?'
	723	call getin("v",grv)
	724	write(,) " v = ",grv
	725
[899]	726	c Initialize winds for first time step
[38]	727	DO ilayer=1,nlayer
	728	u(ilayer)=gru
	729	v(ilayer)=grv
[1541]	730	w(ilayer)=0 ! default: no vertical wind
[38]	731	ENDDO
	732
[899]	733	c Initialize turbulente kinetic energy
[38]	734	DO ilevel=1,nlevel
	735	q2(ilevel)=0.E+0
	736	ENDDO
	737
[899]	738	c CO2 ice on the surface
[38]	739	c -------------------
[1130]	740	co2ice(1)=0.E+0 ! default value for co2ice
[38]	741	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
	742	call getin("co2ice",co2ice)
	743	write(,) " co2ice = ",co2ice
[1944]	744	! Initialization for CO2 clouds (could be improved to read initial profiles)
	745	mem_Mccn_co2(:,:)=0
	746	mem_Mh2o_co2(:,:)=0
	747	mem_Nccn_co2(:,:)=0
[38]	748	c
[899]	749	c emissivity
[38]	750	c ----------
	751	emis=emissiv
[1130]	752	IF (co2ice(1).eq.1.E+0) THEN
[38]	753	emis=emisice(1) ! northern hemisphere
[1541]	754	IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
[38]	755	ENDIF
	756
	757
	758
[899]	759	c Compute pressures and altitudes of atmospheric levels
[38]	760	c ----------------------------------------------------------------
	761
	762	c Vertical Coordinates
	763	c """"""""""""""""""""
	764	hybrid=.true.
	765	PRINT *,'Hybrid coordinates ?'
	766	call getin("hybrid",hybrid)
	767	write(,) " hybrid = ", hybrid
	768
[2167]	769	CALL disvert_noterre
[1621]	770	! now that disvert has been called, initialize module vertical_layers_mod
	771	call init_vertical_layers(nlayer,preff,scaleheight,
	772	& ap,bp,aps,bps,presnivs,pseudoalt)
[38]	773
	774	DO ilevel=1,nlevel
	775	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	776	ENDDO
	777
	778	DO ilayer=1,nlayer
	779	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	780	ENDDO
	781
	782	DO ilayer=1,nlayer
	783	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
	784	& /g
	785	ENDDO
	786
	787
[899]	788	c Initialize temperature profile
[38]	789	c --------------------------------------
	790	pks=psurf**rcp
	791
[899]	792	c altitude in km in profile: divide zlay by 1000
[38]	793	tmp1(0)=0.E+0
	794	DO ilayer=1,nlayer
	795	tmp1(ilayer)=zlay(ilayer)/1000.E+0
	796	ENDDO
	797
	798	call profile(nlayer+1,tmp1,tmp2)
	799
	800	tsurf=tmp2(0)
	801	DO ilayer=1,nlayer
	802	temp(ilayer)=tmp2(ilayer)
	803	ENDDO
	804
	805
	806
	807	! Initialize soil properties and temperature
	808	! ------------------------------------------
	809	volcapa=1.e6 ! volumetric heat capacity
	810	DO isoil=1,nsoil
	811	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
[1130]	812	tsoil(isoil)=tsurf(1) ! soil temperature
[38]	813	ENDDO
	814
	815	! Initialize depths
	816	! -----------------
	817	do isoil=0,nsoil-1
	818	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
	819	enddo
	820	do isoil=1,nsoil
	821	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
	822	enddo
	823
[899]	824	c Initialize traceurs
[38]	825	c ---------------------------
	826
	827	if (photochem.or.callthermos) then
[899]	828	write(,) 'Initializing chemical species'
	829	! thermo=0: initialize over all atmospheric layers
[38]	830	thermo=0
[1541]	831	call inichim_newstart(q,psurf,sig,nq,latitude,longitude,
	832	$ cell_area,thermo,qsurf)
[38]	833	endif
[520]	834
[899]	835	c Check if the surface is a water ice reservoir
[520]	836	c --------------------------------------------------
[2293]	837	watercap(1)=0 ! Initialize watercap
[1047]	838	watercaptag(1)=.false. ! Default: no water ice reservoir
[520]	839	print *,'Water ice cap on ground ?'
	840	call getin("watercaptag",watercaptag)
	841	write(,) " watercaptag = ",watercaptag
[38]	842
	843
[899]	844	c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
[38]	845	c ----------------------------------------------------------------
[899]	846	c (output done in "writeg1d", typically called by "physiq.F")
[38]	847
	848	g1d_nlayer=nlayer
	849	g1d_nomfich='g1d.dat'
	850	g1d_unitfich=40
	851	g1d_nomctl='g1d.ctl'
	852	g1d_unitctl=41
	853	g1d_premier=.true.
	854	g2d_premier=.true.
	855
[899]	856	c Write a "startfi" file
[38]	857	c --------------------
[899]	858	c This file will be read during the first call to "physiq".
	859	c It is needed to transfert physics variables to "physiq"...
[38]	860
[1541]	861	call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
[2214]	862	& nq,dtphys,float(day0),0.,cell_area,
[2079]	863	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
	864	& hmons,summit,base)
[1112]	865	call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
[1944]	866	& dtphys,time,
	867	& tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,tauscaling,
	868	& totcloudfrac,wstar,
	869	& mem_Mccn_co2,mem_Nccn_co2,
[2293]	870	& mem_Mh2o_co2,watercap)
[899]	871
[38]	872	c=======================================================================
[899]	873	c 1D MODEL TIME STEPPING LOOP
[38]	874	c=======================================================================
	875	c
	876	firstcall=.true.
	877	lastcall=.false.
	878
	879	DO idt=1,ndt
	880	c IF (idt.eq.ndt) lastcall=.true.
	881	IF (idt.eq.ndt-day_step-1) then !test
	882	lastcall=.true.
	883	call solarlong(day*1.0,zls)
	884	write(103,) 'Ls=',zls180./pi
[1541]	885	write(103,) 'Lat=', latitude(1)180./pi
[627]	886	write(103,) 'Tau=', tauvis/odprefpsurf
[38]	887	write(103,*) 'RunEnd - Atmos. Temp. File'
	888	write(103,*) 'RunEnd - Atmos. Temp. File'
	889	write(104,) 'Ls=',zls180./pi
[1541]	890	write(104,*) 'Lat=', latitude(1)
[627]	891	write(104,) 'Tau=', tauvis/odprefpsurf
[38]	892	write(104,*) 'RunEnd - Atmos. Temp. File'
	893	ENDIF
	894
[899]	895	c compute geopotential
[38]	896	c ~~~~~~~~~~~~~~~~~~~~~
	897	DO ilayer=1,nlayer
	898	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	899	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	900	ENDDO
	901	phi(1)=pksh(1)(1.E+0-s(1))
	902	DO ilayer=2,nlayer
	903	phi(ilayer)=phi(ilayer-1)+
	904	& pks(h(ilayer-1)+h(ilayer)).5E+0
	905	& *(s(ilayer-1)-s(ilayer))
	906
	907	ENDDO
	908
[899]	909	c call physics
[38]	910	c --------------------
[358]	911	! write(,) "testphys1d avant q", q(1,:)
[1036]	912	CALL physiq (1,llm,nq,
[38]	913	, firstcall,lastcall,
	914	, day,time,dtphys,
	915	, plev,play,phi,
	916	, u, v,temp, q,
	917	, w,
[899]	918	C - outputs
[1576]	919	s du, dv, dtemp, dq,dpsurf)
[358]	920	! write(,) "testphys1d apres q", q(1,:)
[38]	921
[358]	922
[899]	923	c wind increment : specific for 1D
	924	c --------------------------------
[38]	925
[899]	926	c The physics compute the tendencies on u and v,
	927	c here we just add Coriolos effect
[38]	928	c
	929	c DO ilayer=1,nlayer
	930	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
	931	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
	932	c ENDDO
	933
[899]	934	c For some tests : No coriolis force at equator
[1541]	935	c if(latitude(1).eq.0.) then
[38]	936	DO ilayer=1,nlayer
	937	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
	938	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
	939	ENDDO
	940	c end if
	941	c
	942	c
[899]	943	c Compute time for next time step
[38]	944	c ---------------------------------------
	945	firstcall=.false.
	946	time=time+dtphys/daysec
	947	IF (time.gt.1.E+0) then
	948	time=time-1.E+0
	949	day=day+1
	950	ENDIF
	951
[899]	952	c compute winds and temperature for next time step
[38]	953	c ----------------------------------------------------------
	954
	955	DO ilayer=1,nlayer
	956	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
	957	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
	958	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
	959	ENDDO
	960
[899]	961	c compute pressure for next time step
[38]	962	c ----------------------------------------------------------
	963
[1549]	964	psurf=psurf+dtphys*dpsurf(1) ! surface pressure change
[38]	965	DO ilevel=1,nlevel
	966	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	967	ENDDO
	968	DO ilayer=1,nlayer
	969	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	970	ENDDO
	971
[899]	972	! increment tracers
[1036]	973	DO iq = 1, nq
[38]	974	DO ilayer=1,nlayer
	975	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
	976	ENDDO
	977	ENDDO
	978
[899]	979	ENDDO ! of idt=1,ndt ! end of time stepping loop
[38]	980
	981	c ========================================================
[899]	982	c OUTPUTS
[38]	983	c ========================================================
	984
[899]	985	c finalize and close grads files "g1d.dat" and "g1d.ctl"
[38]	986
	987	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
	988	CALL endg1d(1,nlayer,zlay/1000.,ndt)
	989
[1380]	990	write(,) "testphys1d: Everything is cool."
	991
[38]	992	END
	993
	994	c***********************************************************************
	995	c***********************************************************************
[899]	996	c Dummy subroutines used only in 3D, but required to
	997	c compile testphys1d (to cleanly use writediagfi)
[38]	998
[899]	999	subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
	1000
	1001	IMPLICIT NONE
	1002
	1003	INTEGER im,jm,ngrid,nfield
	1004	REAL pdyn(im,jm,nfield)
	1005	REAL pfi(ngrid,nfield)
	1006
	1007	if (ngrid.ne.1) then
	1008	write(,) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
	1009	stop
	1010	endif
	1011
	1012	pdyn(1,1,1:nfield)=pfi(1,1:nfield)
	1013
	1014	end
[38]	1015
	1016	c***********************************************************************
	1017	c***********************************************************************
	1018

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

Download in other formats:

Original Format