Context Navigation

testphys1d.F @ 2316

Last change on this file since 2316 was 2312, checked in by lrossi, 5 years ago

MARS GCM:
Implementation of HDO in the physics
New tracers hdo_vap and hdo_ice are added. Cf. traceur.def.isotopes in deftank for exemple of traceur.def.
All HDO related computations are activated by flag hdo=.true. in callphys.def. (see callphys.def.hdo in deftank for an example)
Additional option hdofrac (true by default) allows for turning on/off fractionation (for tests).
For now, only functional with simplified cloud scheme (so microphys=.false. and activice=.false. also recommended).
Initialisation of start files can be done with option 'inihdo' in newstart.
LR

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

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

Download in other formats:

Original Format