Context Navigation

testphys1d.F @ 1240

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

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

Download in other formats:

Original Format