Context Navigation

testphys1d.F @ 1524

Last change on this file since 1524 was 1524, checked in by emillour, 9 years ago

All GCMS:
More updates to enforce dynamics/physics separation:

get rid of references to "temps_mod" from physics packages;
make a "time_phylmdz_mod.F90" module to store that
information and fill it via "iniphysiq".

File size: 27.5 KB

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

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

Download in other formats:

Original Format