Context Navigation

testphys1d.F @ 2606

Last change on this file since 2606 was 2565, checked in by emillour, 4 years ago
Mars GCM: Some follow-up of commit r2562; update newstart, start2archive and testphys1d accordingly. EM
File size: 29.0 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: tauvis,totcloudfrac
23	use dust_param_mod, only: tauscaling
24	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig,
25	& presnivs,pseudoalt,scaleheight
26	USE vertical_layers_mod, ONLY: init_vertical_layers
27	USE logic_mod, ONLY: hybrid
28	use physics_distribution_mod, only: init_physics_distribution
29	use regular_lonlat_mod, only: init_regular_lonlat
30	use mod_interface_dyn_phys, only: init_interface_dyn_phys
31	USE phys_state_var_init_mod, ONLY: phys_state_var_init
32	USE physiq_mod, ONLY: physiq
33	USE read_profile_mod, ONLY: read_profile
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,dayn ! date initial (sol ; =0 a Ls=0) and final
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
120	character*2 str2
121	character (len=7) :: str7
122	character(len=44) :: txt
123
124	c New flag to compute paleo orbital configurations + few variables JN
125	REAL halfaxe, excentric, Lsperi
126	Logical paleomars
127
128	c MVals: isotopes as in the dynamics (CRisi)
129	INTEGER :: ifils,ipere,generation
130	CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name
131	CHARACTER(len=80) :: line ! to store a line of text
132	INTEGER ierr0
133	LOGICAL :: continu
134
135	c=======================================================================
136
137	c=======================================================================
138	c INITIALISATION
139	c=======================================================================
140	! initialize "serial/parallel" related stuff
141	! CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
142	! CALL init_phys_lmdz(1,1,llm,1,(/1/))
143	! call initcomgeomphy
144
145	c ------------------------------------------------------
146	c Prescribed constants to be set here
147	c ------------------------------------------------------
148
149	pi=2.E+0*asin(1.E+0)
150
151	c Mars planetary constants
152	c ----------------------------
153	rad=3397200. ! mars radius (m) ~3397200 m
154	daysec=88775. ! length of a sol (s) ~88775 s
155	omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1)
156	g=3.72 ! gravity (m.s-2) ~3.72
157	mugaz=43.49 ! atmosphere mola mass (g.mol-1) ~43.49
158	rcp=.256793 ! = r/cp ~0.256793
159	r= 8.314511E+0 *1000.E+0/mugaz
160	cpp= r/rcp
161	year_day = 669 ! lenght of year (sols) ~668.6
162	periheli = 206.66 ! minimum sun-mars distance (Mkm) ~206.66
163	aphelie = 249.22 ! maximum sun-mars distance (Mkm) ~249.22
164	halfaxe = 227.94 ! demi-grand axe de l'ellipse
165	peri_day = 485. ! perihelion date (sols since N. Spring)
166	obliquit = 25.2 ! Obliquity (deg) ~25.2
167	excentric = 0.0934 ! Eccentricity (0.0934)
168
169	c Planetary Boundary Layer and Turbulence parameters
170	c --------------------------------------------------
171	z0_default = 1.e-2 ! surface roughness (m) ~0.01
172	emin_turb = 1.e-6 ! minimal turbulent energy ~1.e-8
173	lmixmin = 30 ! mixing length ~100
174
175	c cap properties and surface emissivities
176	c ----------------------------------------------------
177	emissiv= 0.95 ! Bare ground emissivity ~.95
178	emisice(1)=0.95 ! Northern cap emissivity
179	emisice(2)=0.95 ! Southern cap emisssivity
180	albedice(1)=0.5 ! Northern cap albedo
181	albedice(2)=0.5 ! Southern cap albedo
182	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
183	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
184	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
185	dtemisice(2) = 2. ! time scale for snow metamorphism (south
186
187	c mesh surface (not a very usefull quantity in 1D)
188	c ----------------------------------------------------
189	cell_area(1)=1.E+0
190
191	c ------------------------------------------------------
192	c Loading run parameters from "run.def" file
193	c ------------------------------------------------------
194
195
196	! check if 'run.def' file is around (otherwise reading parameters
197	! from callphys.def via getin() routine won't work.
198	open(99,file='run.def',status='old',form='formatted',
199	& iostat=ierr)
200	if (ierr.ne.0) then
201	write(,) 'Cannot find required file "run.def"'
202	write(,) ' (which should contain some input parameters'
203	write(,) ' along with the following line:'
204	write(,) ' INCLUDEDEF=callphys.def'
205	write(,) ' )'
206	write(,) ' ... might as well stop here ...'
207	stop
208	else
209	close(99)
210	endif
211
212	! check if we are going to run with or without tracers
213	write(,) "Run with or without tracer transport ?"
214	tracer=.false. ! default value
215	call getin("tracer",tracer)
216	write(,) " tracer = ",tracer
217
218	! while we're at it, check if there is a 'traceur.def' file
219	! and process it.
220	if (tracer) then
221	! load tracer names from file 'traceur.def'
222	open(90,file='traceur.def',status='old',form='formatted',
223	& iostat=ierr)
224	if (ierr.ne.0) then
225	write(,) 'Cannot find required file "traceur.def"'
226	write(,) ' If you want to run with tracers, I need it'
227	write(,) ' ... might as well stop here ...'
228	stop
229	else
230	write(,) "testphys1d: Reading file traceur.def"
231	! read number of tracers:
232	read(90,*,iostat=ierr) nq
233	nqtot=nq ! set value of nqtot (in infotrac module) as nq
234	if (ierr.ne.0) then
235	write(,) "testphys1d: error reading number of tracers"
236	write(,) " (first line of traceur.def) "
237	stop
238	endif
239	endif
240	! allocate arrays:
241	allocate(tname(nq))
242	allocate(q(nlayer,nq))
243	allocate(qsurf(nq))
244	allocate(dq(nlayer,nq))
245	allocate(dqdyn(nlayer,nq))
246	allocate(mqtot(nq))
247	allocate(tnom_transp(nq))
248
249	! read tracer names from file traceur.def
250	do iq=1,nq
251	read(90,'(80a)',iostat=ierr) line ! store the line from traceur.def
252	if (ierr.ne.0) then
253	write(,) 'testphys1d: error reading tracer names...'
254	stop
255	endif
256	! if format is tnom_0, tnom_transp (isotopes)
257	read(line,*,iostat=ierr0) tname(iq),tnom_transp(iq)
258	if (ierr0.ne.0) then
259	read(line,*) tname(iq)
260	tnom_transp(iq)='air'
261	endif
262
263	enddo
264	close(90)
265
266	! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
267	ALLOCATE(nqfils(nqtot))
268	nqperes=0
269	nqfils(:)=0
270	DO iq=1,nqtot
271	if (tnom_transp(iq) == 'air') then
272	! ceci est un traceur père
273	WRITE(,) 'Le traceur',iq,', appele ',
274	& trim(tname(iq)),', est un pere'
275	nqperes=nqperes+1
276	else !if (tnom_transp(iq) == 'air') then
277	! ceci est un fils. Qui est son père?
278	WRITE(,) 'Le traceur',iq,', appele ',
279	& trim(tname(iq)),', est un fils'
280	continu=.true.
281	ipere=1
282	do while (continu)
283	if (tnom_transp(iq) .eq. tname(ipere)) then
284	! Son père est ipere
285	WRITE(,) 'Le traceur',iq,'appele ',
286	& trim(tname(iq)),' est le fils de ',
287	& ipere,'appele ',trim(tname(ipere))
288	nqfils(ipere)=nqfils(ipere)+1
289	continu=.false.
290	else !if (tnom_transp(iq) == tnom_0(ipere)) then
291	ipere=ipere+1
292	if (ipere.gt.nqtot) then
293	WRITE(,) 'Le traceur',iq,'appele ',
294	& trim(tname(iq)),', est orpelin.'
295	CALL abort_gcm('infotrac_init',
296	& 'Un traceur est orphelin',1)
297	endif !if (ipere.gt.nqtot) then
298	endif !if (tnom_transp(iq) == tnom_0(ipere)) then
299	enddo !do while (continu)
300	endif !if (tnom_transp(iq) == 'air') then
301	enddo !DO iq=1,nqtot
302	WRITE(,) 'nqperes=',nqperes
303	WRITE(,) 'nqfils=',nqfils
304
305	! initialize tracers here:
306	write(,) "testphys1d: initializing tracers"
307	call read_profile(nq, nlayer, qsurf, q)
308	else
309	! we still need to set (dummy) tracer number and names for physdem1
310	nq=1
311	nqtot=nq ! set value of nqtot (in infotrac module) as nq
312	! allocate arrays:
313	allocate(tname(nq))
314	allocate(q(nlayer,nq))
315	allocate(qsurf(nq))
316	allocate(dq(nlayer,nq))
317	allocate(dqdyn(nlayer,nq))
318	allocate(mqtot(nq))
319	do iq=1,nq
320	write(str7,'(a1,i2.2)')'t',iq
321	tname(iq)=str7
322	enddo
323	! and just to be clean, also initialize tracers to zero for physdem1
324	q(:,:)=0
325	qsurf(:)=0
326	endif ! of if (tracer)
327
328	!write(,) "testphys1d q", q(1,:)
329	!write(,) "testphys1d qsurf", qsurf
330
331	c Date and local time at beginning of run
332	c ---------------------------------------
333	c Date (in sols since spring solstice) at beginning of run
334	day0 = 0 ! default value for day0
335	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
336	call getin("day0",day0)
337	day=float(day0)
338	write(,) " day0 = ",day0
339	c Local time at beginning of run
340	time=0 ! default value for time
341	write(,)'Initial local time (in hours, between 0 and 24)?'
342	call getin("time",time)
343	write(,)" time = ",time
344	time=time/24.E+0 ! convert time (hours) to fraction of sol
345
346	c Discretization (Definition of grid and time steps)
347	c --------------
348	c
349	nlevel=nlayer+1
350	nsoil=nsoilmx
351
352	day_step=48 ! default value for day_step
353	PRINT *,'Number of time steps per sol ?'
354	call getin("day_step",day_step)
355	write(,) " day_step = ",day_step
356
357	ecritphy=day_step ! default value for ecritphy, output every time step
358
359	ndt=10 ! default value for ndt
360	PRINT *,'Number of sols to run ?'
361	call getin("ndt",ndt)
362	write(,) " ndt = ",ndt
363
364	dayn=day0+ndt
365	ndt=ndt*day_step
366	dtphys=daysec/day_step
367
368	c Imposed surface pressure
369	c ------------------------------------
370	c
371	psurf=610. ! default value for psurf
372	PRINT *,'Surface pressure (Pa) ?'
373	call getin("psurf",psurf)
374	write(,) " psurf = ",psurf
375	c Reference pressures
376	pa=20. ! transition pressure (for hybrid coord.)
377	preff=610. ! reference surface pressure
378
379	c Aerosol properties
380	c --------------------------------
381	tauvis=0.2 ! default value for tauvis (dust opacity)
382	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
383	call getin("tauvis",tauvis)
384	write(,) " tauvis = ",tauvis
385
386	c Orbital parameters
387	c ------------------
388	paleomars=.false. ! Default: no water ice reservoir
389	call getin("paleomars",paleomars)
390	if (paleomars.eqv..true.) then
391	write(,) "paleomars=", paleomars
392	write(,) "Orbital parameters from callphys.def"
393	write(,) "Enter eccentricity & Lsperi"
394	print *, 'Martian eccentricity (0<e<1) ?'
395	call getin('excentric ',excentric)
396	write(,)"excentric =",excentric
397	print *, 'Solar longitude of perihelion (0<Ls<360) ?'
398	call getin('Lsperi',Lsperi )
399	write(,)"Lsperi=",Lsperi
400	Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
401	periheli = halfaxe*(1-excentric)
402	aphelie = halfaxe*(1+excentric)
403	call call_dayperi(Lsperi,excentric,peri_day,year_day)
404	write(,) "Corresponding orbital params for GCM"
405	write(,) " periheli = ",periheli
406	write(,) " aphelie = ",aphelie
407	write(,) "date of perihelion (sol)",peri_day
408	else
409	write(,) "paleomars=", paleomars
410	write(,) "Default present-day orbital parameters"
411	write(,) "Unless specified otherwise"
412	print *,'Min. distance Sun-Mars (Mkm)?'
413	call getin("periheli",periheli)
414	write(,) " periheli = ",periheli
415
416	print *,'Max. distance Sun-Mars (Mkm)?'
417	call getin("aphelie",aphelie)
418	write(,) " aphelie = ",aphelie
419
420	print *,'Day of perihelion?'
421	call getin("periday",peri_day)
422	write(,) " periday = ",peri_day
423
424	print *,'Obliquity?'
425	call getin("obliquit",obliquit)
426	write(,) " obliquit = ",obliquit
427	end if
428
429	c latitude/longitude
430	c ------------------
431	latitude(1)=0 ! default value for latitude
432	PRINT *,'latitude (in degrees) ?'
433	call getin("latitude",latitude(1))
434	write(,) " latitude = ",latitude
435	latitude=latitude*pi/180.E+0
436	longitude=0.E+0
437	longitude=longitude*pi/180.E+0
438
439	! some initializations (some of which have already been
440	! done above!) and loads parameters set in callphys.def
441	! and allocates some arrays
442	!Mars possible matter with dtphys in input and include!!!
443	! Initializations below should mimick what is done in iniphysiq for 3D GCM
444	call init_physics_distribution(regular_lonlat,4,
445	& 1,1,1,nlayer,1)
446	call init_interface_dyn_phys
447	call init_regular_lonlat(1,1,longitude,latitude,
448	& (/0.,0./),(/0.,0./))
449	call init_geometry(1,longitude,latitude,
450	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
451	& cell_area)
452	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
453	! call init_vertical_layers(nlayer,preff,scaleheight,
454	! & ap,bp,aps,bps,presnivs,pseudoalt)
455	call init_dimphy(1,nlayer) ! Initialize dimphy module
456	call phys_state_var_init(1,llm,nq,tname,
457	. day0,dayn,time,
458	. daysec,dtphys,
459	. rad,g,r,cpp,
460	. nqperes,nqfils)! MVals: variables isotopes
461	call ini_fillgeom(1,latitude,longitude,(/1.0/))
462	call conf_phys(1,llm,nq)
463
464	! in 1D model physics are called every time step
465	! ovverride iphysiq value that has been set by conf_phys
466	if (iphysiq/=1) then
467	write(,) "testphys1d: setting iphysiq=1"
468	iphysiq=1
469	endif
470
471	c Initialize albedo / soil thermal inertia
472	c ----------------------------------------
473	c
474	albedodat(1)=0.2 ! default value for albedodat
475	PRINT *,'Albedo of bare ground ?'
476	call getin("albedo",albedodat(1))
477	write(,) " albedo = ",albedodat(1)
478	albedo(1,1)=albedodat(1)
479
480	inertiedat(1,1)=400 ! default value for inertiedat
481	PRINT *,'Soil thermal inertia (SI) ?'
482	call getin("inertia",inertiedat(1,1))
483	write(,) " inertia = ",inertiedat(1,1)
484
485	z0(1)=z0_default ! default value for roughness
486	write(,) 'Surface roughness length z0 (m)?'
487	call getin("z0",z0(1))
488	write(,) " z0 = ",z0(1)
489
490	! Initialize local slope parameters (only matters if "callslope"
491	! is .true. in callphys.def)
492	! slope inclination angle (deg) 0: horizontal, 90: vertical
493	theta_sl(1)=0.0 ! default: no inclination
494	call getin("slope_inclination",theta_sl(1))
495	! slope orientation (deg)
496	! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
497	psi_sl(1)=0.0 ! default value
498	call getin("slope_orientation",psi_sl(1))
499
500	c
501	c for the gravity wave scheme
502	c ---------------------------------
503	c
504	zmea(1)=0.E+0
505	zstd(1)=0.E+0
506	zsig(1)=0.E+0
507	zgam(1)=0.E+0
508	zthe(1)=0.E+0
509	c
510	c for the slope wind scheme
511	c ---------------------------------
512	c
513	hmons(1)=0.E+0
514	PRINT *,'hmons is initialized to ',hmons(1)
515	summit(1)=0.E+0
516	PRINT *,'summit is initialized to ',summit(1)
517	base(1)=0.E+0
518	c
519	c Default values initializing the coefficients calculated later
520	c ---------------------------------
521	c
522	tauscaling(1)=1. ! calculated in aeropacity_mod.F
523	totcloudfrac(1)=1. ! calculated in watercloud_mod.F
524
525	c Specific initializations for "physiq"
526	c -------------------------------------
527	c surface geopotential is not used (or useful) since in 1D
528	c everything is controled by surface pressure
529	phisfi(1)=0.E+0
530
531	c Initialization to take into account prescribed winds
532	c ------------------------------------------------------
533	ptif=2.E+0omegsinlat(1)
534
535	c geostrophic wind
536	gru=10. ! default value for gru
537	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
538	call getin("u",gru)
539	write(,) " u = ",gru
540	grv=0. !default value for grv
541	PRINT *,'meridional northward component of the geostrophic',
542	&' wind (m/s) ?'
543	call getin("v",grv)
544	write(,) " v = ",grv
545
546	c Initialize winds for first time step
547	DO ilayer=1,nlayer
548	u(ilayer)=gru
549	v(ilayer)=grv
550	w(ilayer)=0 ! default: no vertical wind
551	ENDDO
552
553	c Initialize turbulente kinetic energy
554	DO ilevel=1,nlevel
555	q2(ilevel)=0.E+0
556	ENDDO
557
558	c CO2 ice on the surface
559	c -------------------
560	co2ice(1)=0.E+0 ! default value for co2ice
561	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
562	call getin("co2ice",co2ice)
563	write(,) " co2ice = ",co2ice
564
565	c
566	c emissivity
567	c ----------
568	emis=emissiv
569	IF (co2ice(1).eq.1.E+0) THEN
570	emis=emisice(1) ! northern hemisphere
571	IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
572	ENDIF
573
574
575
576	c Compute pressures and altitudes of atmospheric levels
577	c ----------------------------------------------------------------
578
579	c Vertical Coordinates
580	c """"""""""""""""""""
581	hybrid=.true.
582	PRINT *,'Hybrid coordinates ?'
583	call getin("hybrid",hybrid)
584	write(,) " hybrid = ", hybrid
585
586	CALL disvert_noterre
587	! now that disvert has been called, initialize module vertical_layers_mod
588	call init_vertical_layers(nlayer,preff,scaleheight,
589	& ap,bp,aps,bps,presnivs,pseudoalt)
590
591	DO ilevel=1,nlevel
592	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
593	ENDDO
594
595	DO ilayer=1,nlayer
596	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
597	ENDDO
598
599	DO ilayer=1,nlayer
600	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
601	& /g
602	ENDDO
603
604
605	c Initialize temperature profile
606	c --------------------------------------
607	pks=psurf**rcp
608
609	c altitude in km in profile: divide zlay by 1000
610	tmp1(0)=0.E+0
611	DO ilayer=1,nlayer
612	tmp1(ilayer)=zlay(ilayer)/1000.E+0
613	ENDDO
614
615	call profile(nlayer+1,tmp1,tmp2)
616
617	tsurf=tmp2(0)
618	DO ilayer=1,nlayer
619	temp(ilayer)=tmp2(ilayer)
620	ENDDO
621
622
623
624	! Initialize soil properties and temperature
625	! ------------------------------------------
626	volcapa=1.e6 ! volumetric heat capacity
627	DO isoil=1,nsoil
628	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
629	tsoil(isoil)=tsurf(1) ! soil temperature
630	ENDDO
631
632	! Initialize depths
633	! -----------------
634	do isoil=0,nsoil-1
635	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
636	enddo
637	do isoil=1,nsoil
638	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
639	enddo
640
641	c Initialize traceurs
642	c ---------------------------
643
644	if (photochem.or.callthermos) then
645	write(,) 'Initializing chemical species'
646	! thermo=0: initialize over all atmospheric layers
647	thermo=0
648	call inichim_newstart(q,psurf,sig,nq,latitude,longitude,
649	$ cell_area,thermo,qsurf)
650	endif
651
652	c Check if the surface is a water ice reservoir
653	c --------------------------------------------------
654	watercap(1)=0 ! Initialize watercap
655	watercaptag(1)=.false. ! Default: no water ice reservoir
656	print *,'Water ice cap on ground ?'
657	call getin("watercaptag",watercaptag)
658	write(,) " watercaptag = ",watercaptag
659
660
661	c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
662	c ----------------------------------------------------------------
663	c (output done in "writeg1d", typically called by "physiq.F")
664
665	g1d_nlayer=nlayer
666	g1d_nomfich='g1d.dat'
667	g1d_unitfich=40
668	g1d_nomctl='g1d.ctl'
669	g1d_unitctl=41
670	g1d_premier=.true.
671	g2d_premier=.true.
672
673	c Write a "startfi" file
674	c --------------------
675	c This file will be read during the first call to "physiq".
676	c It is needed to transfert physics variables to "physiq"...
677
678	call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
679	& nq,dtphys,float(day0),0.,cell_area,
680	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
681	& hmons,summit,base)
682	call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
683	& dtphys,time,
684	& tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,tauscaling,
685	& totcloudfrac,wstar,watercap)
686
687	c=======================================================================
688	c 1D MODEL TIME STEPPING LOOP
689	c=======================================================================
690	c
691	firstcall=.true.
692	lastcall=.false.
693
694	DO idt=1,ndt
695	c IF (idt.eq.ndt) lastcall=.true.
696	IF (idt.eq.ndt-day_step-1) then !test
697	lastcall=.true.
698	call solarlong(day*1.0,zls)
699	write(103,) 'Ls=',zls180./pi
700	write(103,) 'Lat=', latitude(1)180./pi
701	write(103,) 'Tau=', tauvis/odprefpsurf
702	write(103,*) 'RunEnd - Atmos. Temp. File'
703	write(103,*) 'RunEnd - Atmos. Temp. File'
704	write(104,) 'Ls=',zls180./pi
705	write(104,*) 'Lat=', latitude(1)
706	write(104,) 'Tau=', tauvis/odprefpsurf
707	write(104,*) 'RunEnd - Atmos. Temp. File'
708	ENDIF
709
710	c compute geopotential
711	c ~~~~~~~~~~~~~~~~~~~~~
712	DO ilayer=1,nlayer
713	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
714	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
715	ENDDO
716	phi(1)=pksh(1)(1.E+0-s(1))
717	DO ilayer=2,nlayer
718	phi(ilayer)=phi(ilayer-1)+
719	& pks(h(ilayer-1)+h(ilayer)).5E+0
720	& *(s(ilayer-1)-s(ilayer))
721
722	ENDDO
723
724	c call physics
725	c --------------------
726	! write(,) "testphys1d avant q", q(1,:)
727	CALL physiq (1,llm,nq,
728	, firstcall,lastcall,
729	, day,time,dtphys,
730	, plev,play,phi,
731	, u, v,temp, q,
732	, w,
733	C - outputs
734	s du, dv, dtemp, dq,dpsurf)
735	! write(,) "testphys1d apres q", q(1,:)
736
737
738	c wind increment : specific for 1D
739	c --------------------------------
740
741	c The physics compute the tendencies on u and v,
742	c here we just add Coriolos effect
743	c
744	c DO ilayer=1,nlayer
745	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
746	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
747	c ENDDO
748
749	c For some tests : No coriolis force at equator
750	c if(latitude(1).eq.0.) then
751	DO ilayer=1,nlayer
752	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
753	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
754	ENDDO
755	c end if
756	c
757	c
758	c Compute time for next time step
759	c ---------------------------------------
760	firstcall=.false.
761	time=time+dtphys/daysec
762	IF (time.gt.1.E+0) then
763	time=time-1.E+0
764	day=day+1
765	ENDIF
766
767	c compute winds and temperature for next time step
768	c ----------------------------------------------------------
769
770	DO ilayer=1,nlayer
771	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
772	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
773	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
774	ENDDO
775
776	c compute pressure for next time step
777	c ----------------------------------------------------------
778
779	psurf=psurf+dtphys*dpsurf(1) ! surface pressure change
780	DO ilevel=1,nlevel
781	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
782	ENDDO
783	DO ilayer=1,nlayer
784	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
785	ENDDO
786
787	! increment tracers
788	DO iq = 1, nq
789	DO ilayer=1,nlayer
790	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
791	ENDDO
792	ENDDO
793
794	ENDDO ! of idt=1,ndt ! end of time stepping loop
795
796	c ========================================================
797	c OUTPUTS
798	c ========================================================
799
800	c finalize and close grads files "g1d.dat" and "g1d.ctl"
801
802	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
803	CALL endg1d(1,nlayer,zlay/1000.,ndt)
804
805	write(,) "testphys1d: Everything is cool."
806
807	END
808
809	c***********************************************************************
810	c***********************************************************************
811	c Dummy subroutines used only in 3D, but required to
812	c compile testphys1d (to cleanly use writediagfi)
813
814	subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
815
816	IMPLICIT NONE
817
818	INTEGER im,jm,ngrid,nfield
819	REAL pdyn(im,jm,nfield)
820	REAL pfi(ngrid,nfield)
821
822	if (ngrid.ne.1) then
823	write(,) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
824	stop
825	endif
826
827	pdyn(1,1,1:nfield)=pfi(1,1:nfield)
828
829	end
830
831	c***********************************************************************
832	c***********************************************************************
833

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

Download in other formats:

Original Format