Context Navigation

testphys1d.F @ 1146

Last change on this file since 1146 was 1130, checked in by emillour, 12 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
dyn3d/control.h is now module control_mod.F90
rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

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

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

Download in other formats:

Original Format