Context Navigation

testphys1d.F @ 1103

Last change on this file since 1103 was 1047, checked in by emillour, 12 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

File size: 26.8 KB

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

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

Download in other formats:

Original Format