Context Navigation

testphys1d.F @ 3026

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

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

Download in other formats:

Original Format