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

Line
1
2	PROGRAM testphys1d
3	! to use 'getin'
4	USE ioipsl_getincom, only: getin
5	use dimphy, only : init_dimphy
6	use mod_grid_phy_lmdz, only : regular_lonlat
7	use infotrac, only: nqtot, tname, nqperes,nqfils
8	use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx
9	use comgeomfi_h, only: sinlat, ini_fillgeom
10	use surfdat_h, only: albedodat, z0_default, emissiv, emisice,
11	& albedice, iceradius, dtemisice, z0,
12	& zmea, zstd, zsig, zgam, zthe, phisfi,
13	& watercaptag, watercap, hmons, summit, base
14	use slope_mod, only: theta_sl, psi_sl
15	use phyredem, only: physdem0,physdem1
16	use geometry_mod, only: init_geometry
17	use planete_h, only: year_day, periheli, aphelie, peri_day,
18	& obliquit, emin_turb, lmixmin
19	use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp
20	use time_phylmdz_mod, only: daysec, dtphys, day_step,
21	& ecritphy, iphysiq
22	use dimradmars_mod, only: tauscaling,tauvis,totcloudfrac
23	use co2cloud_mod, only: mem_Mccn_co2, mem_Mh2o_co2,
24	& mem_Nccn_co2
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
28	USE logic_mod, ONLY: hybrid
29	use physics_distribution_mod, only: init_physics_distribution
30	use regular_lonlat_mod, only: init_regular_lonlat
31	use mod_interface_dyn_phys, only: init_interface_dyn_phys
32	USE phys_state_var_init_mod, ONLY: phys_state_var_init
33	USE physiq_mod, ONLY: physiq
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
55	include "dimensions.h"
56	integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm)
57	integer, parameter :: nlayer = llm
58	!#include "dimradmars.h"
59	!#include "comgeomfi.h"
60	!#include "surfdat.h"
61	!#include "slope.h"
62	!#include "comsoil.h"
63	!#include "comdiurn.h"
64	include "callkeys.h"
65	!#include "comsaison.h"
66	!#include "control.h"
67	include "netcdf.inc"
68	include "comg1d.h"
69	!#include "advtrac.h"
70
71	c --------------------------------------------------------------
72	c Declarations
73	c --------------------------------------------------------------
74	c
75	INTEGER unitstart ! unite d'ecriture de "startfi"
76	INTEGER nlevel,nsoil,ndt
77	INTEGER ilayer,ilevel,isoil,idt,iq
78	LOGICAl firstcall,lastcall
79	c
80	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
81	c
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)
85	REAL play(nlayer) ! Pressure at the middle of the layers (Pa)
86	REAL plev(nlayer+1) ! intermediate pressure levels (pa)
87	REAL psurf,tsurf(1)
88	REAL u(nlayer),v(nlayer) ! zonal, meridional wind
89	REAL gru,grv ! prescribed "geostrophic" background wind
90	REAL temp(nlayer) ! temperature at the middle of the layers
91	REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
92	REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
93	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
94	REAL co2ice(1) ! co2ice layer (kg.m-2)
95	REAL emis(1) ! surface layer
96	REAL albedo(1,1) ! surface albedo
97	REAL :: wstar(1)=0. ! Thermals vertical velocity
98	REAL q2(nlayer+1) ! Turbulent Kinetic Energy
99	REAL zlay(nlayer) ! altitude estimee dans les couches (km)
100
101	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
102	REAL du(nlayer),dv(nlayer),dtemp(nlayer)
103	REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer)
104	REAL dpsurf(1)
105	REAL,ALLOCATABLE :: dq(:,:)
106	REAL,ALLOCATABLE :: dqdyn(:,:)
107
108	c Various intermediate variables
109	INTEGER thermo
110	REAL zls
111	REAL phi(nlayer),h(nlayer),s(nlayer)
112	REAL pks, ptif, w(nlayer)
113	REAL qtotinit,qtot
114	real,allocatable :: mqtot(:)
115	INTEGER ierr, aslun
116	REAL tmp1(0:nlayer),tmp2(0:nlayer)
117	integer :: nq=1 ! number of tracers
118	real :: latitude(1), longitude(1), cell_area(1)
119	integer iqh2ovap
120	integer iqh2oice
121
122	character*2 str2
123	character (len=7) :: str7
124	character(len=44) :: txt
125
126	c New flag to compute paleo orbital configurations + few variables JN
127	REAL halfaxe, excentric, Lsperi
128	Logical paleomars
129
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
137	c=======================================================================
138
139	c=======================================================================
140	c INITIALISATION
141	c=======================================================================
142	! initialize "serial/parallel" related stuff
143	! CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
144	! CALL init_phys_lmdz(1,1,llm,1,(/1/))
145	! call initcomgeomphy
146
147	c ------------------------------------------------------
148	c Prescribed constants to be set here
149	c ------------------------------------------------------
150
151	pi=2.E+0*asin(1.E+0)
152
153	c Mars planetary constants
154	c ----------------------------
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
160	rcp=.256793 ! = r/cp ~0.256793
161	r= 8.314511E+0 *1000.E+0/mugaz
162	cpp= r/rcp
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
166	halfaxe = 227.94 ! demi-grand axe de l'ellipse
167	peri_day = 485. ! perihelion date (sols since N. Spring)
168	obliquit = 25.2 ! Obliquity (deg) ~25.2
169	excentric = 0.0934 ! Eccentricity (0.0934)
170
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
176
177	c cap properties and surface emissivities
178	c ----------------------------------------------------
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
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
189	c mesh surface (not a very usefull quantity in 1D)
190	c ----------------------------------------------------
191	cell_area(1)=1.E+0
192
193	c ------------------------------------------------------
194	c Loading run parameters from "run.def" file
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
221	! and process it.
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
235	nqtot=nq ! set value of nqtot (in infotrac module) as nq
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
242	! allocate arrays:
243	allocate(tname(nq))
244	allocate(q(nlayer,nq))
245	allocate(qsurf(nq))
246	allocate(dq(nlayer,nq))
247	allocate(dqdyn(nlayer,nq))
248	allocate(mqtot(nq))
249	allocate(tnom_transp(nq))
250
251	! read tracer names from file traceur.def
252	do iq=1,nq
253	read(90,'(80a)',iostat=ierr) line ! store the line from traceur.def
254	if (ierr.ne.0) then
255	write(,) 'testphys1d: error reading tracer names...'
256	stop
257	endif
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
265	enddo
266	close(90)
267
268	! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
269	ALLOCATE(nqfils(nqtot))
270	nqperes=0
271	nqfils(:)=0
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))
290	nqfils(ipere)=nqfils(ipere)+1
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
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)
313	do iq=1,nq
314	txt=""
315	write(txt,"(a)") tname(iq)
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)
326	do ilayer=1,nlayer
327	read(91,*) q(ilayer,iq)
328	enddo
329	endif
330	close(91)
331	endif ! of if (txt.eq."co2")
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)
341	do ilayer=1,nlayer
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
350	! WATER VAPOUR
351	if (txt.eq."h2o_vap") then
352	iqh2ovap=iq !remember index for water vap
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)
358	do ilayer=1,nlayer
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
368	iqh2oice = iq !remember index for water ice
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)
374	do ilayer=1,nlayer
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")
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
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)
434	do ilayer=1,nlayer
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)
450	do ilayer=1,nlayer
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")
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)
466	do ilayer=1,nlayer
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)
481	do ilayer=1,nlayer
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")
489	enddo ! of do iq=1,nq
490
491	else
492	! we still need to set (dummy) tracer number and names for physdem1
493	nq=1
494	nqtot=nq ! set value of nqtot (in infotrac module) as nq
495	! allocate arrays:
496	allocate(tname(nq))
497	allocate(q(nlayer,nq))
498	allocate(qsurf(nq))
499	allocate(dq(nlayer,nq))
500	allocate(dqdyn(nlayer,nq))
501	allocate(mqtot(nq))
502	do iq=1,nq
503	write(str7,'(a1,i2.2)')'t',iq
504	tname(iq)=str7
505	enddo
506	! and just to be clean, also initialize tracers to zero for physdem1
507	q(:,:)=0
508	qsurf(:)=0
509	endif ! of if (tracer)
510
511	!write(,) "testphys1d q", q(1,:)
512	!write(,) "testphys1d qsurf", qsurf
513
514	c Date and local time at beginning of run
515	c ---------------------------------------
516	c Date (in sols since spring solstice) at beginning of run
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
522	c Local time at beginning of run
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
529	c Discretization (Definition of grid and time steps)
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
540	ecritphy=day_step ! default value for ecritphy, output every time step
541
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
549
550	c Imposed surface pressure
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
557	c Reference pressures
558	pa=20. ! transition pressure (for hybrid coord.)
559	preff=610. ! reference surface pressure
560
561	c Aerosol properties
562	c --------------------------------
563	tauvis=0.2 ! default value for tauvis (dust opacity)
564	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
565	call getin("tauvis",tauvis)
566	write(,) " tauvis = ",tauvis
567
568	c Orbital parameters
569	c ------------------
570	paleomars=.false. ! Default: no water ice reservoir
571	call getin("paleomars",paleomars)
572	if (paleomars.eqv..true.) then
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
597
598	print *,'Max. distance Sun-Mars (Mkm)?'
599	call getin("aphelie",aphelie)
600	write(,) " aphelie = ",aphelie
601
602	print *,'Day of perihelion?'
603	call getin("periday",peri_day)
604	write(,) " periday = ",peri_day
605
606	print *,'Obliquity?'
607	call getin("obliquit",obliquit)
608	write(,) " obliquit = ",obliquit
609	end if
610
611	c latitude/longitude
612	c ------------------
613	latitude(1)=0 ! default value for latitude
614	PRINT *,'latitude (in degrees) ?'
615	call getin("latitude",latitude(1))
616	write(,) " latitude = ",latitude
617	latitude=latitude*pi/180.E+0
618	longitude=0.E+0
619	longitude=longitude*pi/180.E+0
620
621	! some initializations (some of which have already been
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!!!
625	! Initializations below should mimick what is done in iniphysiq for 3D GCM
626	call init_physics_distribution(regular_lonlat,4,
627	& 1,1,1,nlayer,1)
628	call init_interface_dyn_phys
629	call init_regular_lonlat(1,1,longitude,latitude,
630	& (/0.,0./),(/0.,0./))
631	call init_geometry(1,longitude,latitude,
632	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
633	& cell_area)
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)
637	call init_dimphy(1,nlayer) ! Initialize dimphy module
638	call phys_state_var_init(1,llm,nq,tname,
639	. day0,time,daysec,dtphys,rad,g,r,cpp,
640	. nqperes,nqfils)! MVals: variables isotopes
641	call ini_fillgeom(1,latitude,longitude,(/1.0/))
642	call conf_phys(1,llm,nq)
643
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
650
651	c Initialize albedo / soil thermal inertia
652	c ----------------------------------------
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)
658	albedo(1,1)=albedodat(1)
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)
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)
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
680	c
681	c for the gravity wave scheme
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
688	zthe(1)=0.E+0
689	c
690	c for the slope wind scheme
691	c ---------------------------------
692	c
693	hmons(1)=0.E+0
694	PRINT *,'hmons is initialized to ',hmons(1)
695	summit(1)=0.E+0
696	PRINT *,'summit is initialized to ',summit(1)
697	base(1)=0.E+0
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
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
709	phisfi(1)=0.E+0
710
711	c Initialization to take into account prescribed winds
712	c ------------------------------------------------------
713	ptif=2.E+0omegsinlat(1)
714
715	c geostrophic wind
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
726	c Initialize winds for first time step
727	DO ilayer=1,nlayer
728	u(ilayer)=gru
729	v(ilayer)=grv
730	w(ilayer)=0 ! default: no vertical wind
731	ENDDO
732
733	c Initialize turbulente kinetic energy
734	DO ilevel=1,nlevel
735	q2(ilevel)=0.E+0
736	ENDDO
737
738	c CO2 ice on the surface
739	c -------------------
740	co2ice(1)=0.E+0 ! default value for co2ice
741	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
742	call getin("co2ice",co2ice)
743	write(,) " co2ice = ",co2ice
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
748	c
749	c emissivity
750	c ----------
751	emis=emissiv
752	IF (co2ice(1).eq.1.E+0) THEN
753	emis=emisice(1) ! northern hemisphere
754	IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
755	ENDIF
756
757
758
759	c Compute pressures and altitudes of atmospheric levels
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
769	CALL disvert_noterre
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)
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
788	c Initialize temperature profile
789	c --------------------------------------
790	pks=psurf**rcp
791
792	c altitude in km in profile: divide zlay by 1000
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
812	tsoil(isoil)=tsurf(1) ! soil temperature
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
824	c Initialize traceurs
825	c ---------------------------
826
827	if (photochem.or.callthermos) then
828	write(,) 'Initializing chemical species'
829	! thermo=0: initialize over all atmospheric layers
830	thermo=0
831	call inichim_newstart(q,psurf,sig,nq,latitude,longitude,
832	$ cell_area,thermo,qsurf)
833	endif
834
835	c Check if the surface is a water ice reservoir
836	c --------------------------------------------------
837	watercap(1)=0 ! Initialize watercap
838	watercaptag(1)=.false. ! Default: no water ice reservoir
839	print *,'Water ice cap on ground ?'
840	call getin("watercaptag",watercaptag)
841	write(,) " watercaptag = ",watercaptag
842
843
844	c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
845	c ----------------------------------------------------------------
846	c (output done in "writeg1d", typically called by "physiq.F")
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
856	c Write a "startfi" file
857	c --------------------
858	c This file will be read during the first call to "physiq".
859	c It is needed to transfert physics variables to "physiq"...
860
861	call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
862	& nq,dtphys,float(day0),0.,cell_area,
863	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
864	& hmons,summit,base)
865	call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
866	& dtphys,time,
867	& tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,tauscaling,
868	& totcloudfrac,wstar,
869	& mem_Mccn_co2,mem_Nccn_co2,
870	& mem_Mh2o_co2,watercap)
871
872	c=======================================================================
873	c 1D MODEL TIME STEPPING LOOP
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
885	write(103,) 'Lat=', latitude(1)180./pi
886	write(103,) 'Tau=', tauvis/odprefpsurf
887	write(103,*) 'RunEnd - Atmos. Temp. File'
888	write(103,*) 'RunEnd - Atmos. Temp. File'
889	write(104,) 'Ls=',zls180./pi
890	write(104,*) 'Lat=', latitude(1)
891	write(104,) 'Tau=', tauvis/odprefpsurf
892	write(104,*) 'RunEnd - Atmos. Temp. File'
893	ENDIF
894
895	c compute geopotential
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
909	c call physics
910	c --------------------
911	! write(,) "testphys1d avant q", q(1,:)
912	CALL physiq (1,llm,nq,
913	, firstcall,lastcall,
914	, day,time,dtphys,
915	, plev,play,phi,
916	, u, v,temp, q,
917	, w,
918	C - outputs
919	s du, dv, dtemp, dq,dpsurf)
920	! write(,) "testphys1d apres q", q(1,:)
921
922
923	c wind increment : specific for 1D
924	c --------------------------------
925
926	c The physics compute the tendencies on u and v,
927	c here we just add Coriolos effect
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
934	c For some tests : No coriolis force at equator
935	c if(latitude(1).eq.0.) then
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
943	c Compute time for next time step
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
952	c compute winds and temperature for next time step
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
961	c compute pressure for next time step
962	c ----------------------------------------------------------
963
964	psurf=psurf+dtphys*dpsurf(1) ! surface pressure change
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
972	! increment tracers
973	DO iq = 1, nq
974	DO ilayer=1,nlayer
975	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
976	ENDDO
977	ENDDO
978
979	ENDDO ! of idt=1,ndt ! end of time stepping loop
980
981	c ========================================================
982	c OUTPUTS
983	c ========================================================
984
985	c finalize and close grads files "g1d.dat" and "g1d.ctl"
986
987	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
988	CALL endg1d(1,nlayer,zlay/1000.,ndt)
989
990	write(,) "testphys1d: Everything is cool."
991
992	END
993
994	c***********************************************************************
995	c***********************************************************************
996	c Dummy subroutines used only in 3D, but required to
997	c compile testphys1d (to cleanly use writediagfi)
998
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
1015
1016	c***********************************************************************
1017	c***********************************************************************
1018

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

Download in other formats:

Original Format