Context Navigation

testphys1d.F @ 1224

Last change on this file since 1224 was 1224, checked in by aslmd, 11 years ago

LMDZ.MARS. If not major, a quite important commit.

No more SAVE,ALLOCATABLE arrays outside modules.

This is important to solve the nesting conundrum in MESOSCALE.
And overall this is good for the harmony of the universe.
(Joke apart, this is good for any interfacing task. And compliant with a F90 spirit).
Note that bit-to-bit compatibility of results in debug mode was checked.

inifis is split in two : phys_state_var_init + conf_phys

This makes interfacing with MESOSCALE more transparent.
This is also clearer for LMDZ.MARS.
Before, inifis has two very different tasks to do.

a bit of cleaning as far as modules and saves are concerned

Point 1

Removed SAVE,ALLOCATABLE arrays from

physiq, aeropacity, updatereffrad, soil

and put those in

dimradmars_mod, surfdat_h, tracer_mod, comsoil_h

and changed accordingly the initialization subroutines associated to each module.
Allocating these arrays is thus done at initialization.

Point 2

Created a subroutine phys_state_var_init which does all the allocation / initialization work for modules. This was previously done in inifis.
Replaced inifis which was then (after the previous modification) just about reading callphys.def and setting a few constants by conf_phys. This mimics the new LMDZ terminology (cf. LMDZ.VENUS for instance)
Bye bye inifis.

Point 3

Removed comdiurn and put everything in comgeomfi
Created a turb_mod module for turbulence variables (e.g. l0 in yamada4)
dryness had nothing to do in tracer_h, put it in surfdat_h (like watercaptag)
topdust0 does not need to be SAVE in aeropacity. better use sinlat.
emisref does not need to be SAVE in newcondens. made it automatic array.
Removed useless co2ice argument in initracer.

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

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

Download in other formats:

Original Format