Context Navigation

testphys1d.F @ 1229

Last change on this file since 1229 was 1229, checked in by emillour, 11 years ago

Mars GCM:
Some minor fixes and cleanup.

allocation of qsurf() was missing in start2archive
'tauscaling' was missing from argument list when calling physdem1 in testphys1d

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