Context Navigation

testphys1d.F @ 2829

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

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

Download in other formats:

Original Format