Context Navigation

testphys1d.F @ 2790

Last change on this file since 2790 was 2790, checked in by jnaar, 2 years ago
Minor change of 1d flag "atm_wat_profile" which can now be used to prescribe directly the desired MMR profile for water vapor. JN
File size: 32.1 KB

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

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

Download in other formats:

Original Format