Context Navigation

testphys1d.F @ 2981

Last change on this file since 2981 was 2981, checked in by evos, 19 months ago

EV 1D Mars

Created an option to output the profiles in 1D using write_prof (default is false)
for the several years run mostly the PEM

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

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

Download in other formats:

Original Format