Context Navigation

testphys1d.F @ 2947

Last change on this file since 2947 was 2947, checked in by llange, 20 months ago
Mars PCM Following r-2942: Fix a bug in the 1D: inertiesoil was not written in the startfile I introduce the possibility to prescribe subsurface ice (when no startfi file): In the callphys.def, set subsurface_ice_depth to the depth where you want ice (set value <0 for no ice). It will assign for these depth a thermal inertia of 2100. LL
File size: 35.7 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
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	c
90	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
91	c
92	INTEGER day0,dayn ! date initial (sol ; =0 a Ls=0) and final
93	REAL day ! date durant le run
94	REAL time ! time (0<time<1 ; time=0.5 a midi)
95	REAL play(nlayer) ! Pressure at the middle of the layers (Pa)
96	REAL plev(nlayer+1) ! intermediate pressure levels (pa)
97	REAL psurf,tsurf(1)
98	REAL u(nlayer),v(nlayer) ! zonal, meridional wind
99	REAL gru,grv ! prescribed "geostrophic" background wind
100	REAL temp(nlayer) ! temperature at the middle of the layers
101	REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
102	REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
103	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
104	REAL emis(1) ! surface layer
105	REAL albedo(1,1) ! surface albedo
106	REAL :: wstar(1)=0. ! Thermals vertical velocity
107	REAL q2(nlayer+1) ! Turbulent Kinetic Energy
108	REAL zlay(nlayer) ! altitude estimee dans les couches (km)
109
110	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
111	REAL du(nlayer),dv(nlayer),dtemp(nlayer)
112	REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer)
113	REAL dpsurf(1)
114	REAL,ALLOCATABLE :: dq(:,:)
115	REAL,ALLOCATABLE :: dqdyn(:,:)
116
117	c Various intermediate variables
118	INTEGER flagthermo,flagh2o
119	REAL zls
120	REAL phi(nlayer),h(nlayer),s(nlayer)
121	REAL pks, ptif, w(nlayer)
122	REAL qtotinit,qtot
123	real,allocatable :: mqtot(:)
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(mqtot(nq))
289	allocate(tnom_transp(nq))
290
291	! read tracer names from file traceur.def
292	do iq=1,nq
293	read(90,'(80a)',iostat=ierr) line ! store the line from traceur.def
294	if (ierr.ne.0) then
295	write(,) 'testphys1d: error reading tracer names...'
296	stop
297	endif
298	! if format is tnom_0, tnom_transp (isotopes)
299	read(line,*,iostat=ierr0) tname(iq),tnom_transp(iq)
300	if (ierr0.ne.0) then
301	read(line,*) tname(iq)
302	tnom_transp(iq)='air'
303	endif
304
305	enddo
306	close(90)
307
308	! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
309	ALLOCATE(nqfils(nqtot))
310	nqperes=0
311	nqfils(:)=0
312	DO iq=1,nqtot
313	if (tnom_transp(iq) == 'air') then
314	! ceci est un traceur père
315	WRITE(,) 'Le traceur',iq,', appele ',
316	& trim(tname(iq)),', est un pere'
317	nqperes=nqperes+1
318	else !if (tnom_transp(iq) == 'air') then
319	! ceci est un fils. Qui est son père?
320	WRITE(,) 'Le traceur',iq,', appele ',
321	& trim(tname(iq)),', est un fils'
322	continu=.true.
323	ipere=1
324	do while (continu)
325	if (tnom_transp(iq) .eq. tname(ipere)) then
326	! Son père est ipere
327	WRITE(,) 'Le traceur',iq,'appele ',
328	& trim(tname(iq)),' est le fils de ',
329	& ipere,'appele ',trim(tname(ipere))
330	nqfils(ipere)=nqfils(ipere)+1
331	continu=.false.
332	else !if (tnom_transp(iq) == tnom_0(ipere)) then
333	ipere=ipere+1
334	if (ipere.gt.nqtot) then
335	WRITE(,) 'Le traceur',iq,'appele ',
336	& trim(tname(iq)),', est orpelin.'
337	CALL abort_gcm('infotrac_init',
338	& 'Un traceur est orphelin',1)
339	endif !if (ipere.gt.nqtot) then
340	endif !if (tnom_transp(iq) == tnom_0(ipere)) then
341	enddo !do while (continu)
342	endif !if (tnom_transp(iq) == 'air') then
343	enddo !DO iq=1,nqtot
344	WRITE(,) 'nqperes=',nqperes
345	WRITE(,) 'nqfils=',nqfils
346
347	! initialize tracers here:
348	write(,) "testphys1d: initializing tracers"
349	call read_profile(nq, nlayer, qsurf, q)
350
351	call init_physics_distribution(regular_lonlat,4,
352	& 1,1,1,nlayer,1)
353
354	if(.not. startfile_1D ) then
355
356	c Date and local time at beginning of run
357	c ---------------------------------------
358	c Date (in sols since spring solstice) at beginning of run
359	day0 = 0 ! default value for day0
360	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
361	call getin("day0",day0)
362	day=float(day0)
363	write(,) " day0 = ",day0
364	c Local time at beginning of run
365	time=0 ! default value for time
366	write(,)'Initial local time (in hours, between 0 and 24)?'
367	call getin("time",time)
368	write(,)" time = ",time
369	time=time/24.E+0 ! convert time (hours) to fraction of sol
370
371	else
372	call open_startphy("startfi.nc")
373	call get_var("controle",tab_cntrl,found)
374	if (.not.found) then
375	call abort_physic("open_startphy",
376	& "tabfi: Failed reading <controle> array",1)
377	else
378	write(,)'tabfi: tab_cntrl',tab_cntrl
379	endif
380	day0 = tab_cntrl(3)
381	day=float(day0)
382
383	call get_var("Time",time,found)
384
385	call close_startphy
386
387	endif !startfile_1D
388
389	c Discretization (Definition of grid and time steps)
390	c --------------
391	c
392	nlevel=nlayer+1
393	nsoil=nsoilmx
394
395	day_step=48 ! default value for day_step
396	PRINT *,'Number of time steps per sol ?'
397	call getin("day_step",day_step)
398	write(,) " day_step = ",day_step
399
400	ecritphy=day_step ! default value for ecritphy, output every time step
401
402	ndt=10 ! default value for ndt
403	PRINT *,'Number of sols to run ?'
404	call getin("ndt",ndt)
405	write(,) " ndt = ",ndt
406
407	dayn=day0+ndt
408	ndt=ndt*day_step
409	dtphys=daysec/day_step
410
411	c Imposed surface pressure
412	c ------------------------------------
413	c
414
415	psurf=610. ! default value for psurf
416	PRINT *,'Surface pressure (Pa) ?'
417	call getin("psurf",psurf)
418	write(,) " psurf = ",psurf
419
420	c Reference pressures
421	pa=20. ! transition pressure (for hybrid coord.)
422	preff=610. ! reference surface pressure
423
424	c Aerosol properties
425	c --------------------------------
426	tauvis=0.2 ! default value for tauvis (dust opacity)
427	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
428	call getin("tauvis",tauvis)
429	write(,) " tauvis = ",tauvis
430
431	c Orbital parameters
432	c ------------------
433
434	if(.not. startfile_1D ) then
435
436	paleomars=.false. ! Default: no water ice reservoir
437	call getin("paleomars",paleomars)
438	if (paleomars.eqv..true.) then
439	write(,) "paleomars=", paleomars
440	write(,) "Orbital parameters from callphys.def"
441	write(,) "Enter eccentricity & Lsperi"
442	print *, 'Martian eccentricity (0<e<1) ?'
443	call getin('excentric ',excentric)
444	write(,)"excentric =",excentric
445	print *, 'Solar longitude of perihelion (0<Ls<360) ?'
446	call getin('Lsperi',Lsperi )
447	write(,)"Lsperi=",Lsperi
448	Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
449	periheli = halfaxe*(1-excentric)
450	aphelie = halfaxe*(1+excentric)
451	call call_dayperi(Lsperi,excentric,peri_day,year_day)
452	write(,) "Corresponding orbital params for GCM"
453	write(,) " periheli = ",periheli
454	write(,) " aphelie = ",aphelie
455	write(,) "date of perihelion (sol)",peri_day
456	else
457	write(,) "paleomars=", paleomars
458	write(,) "Default present-day orbital parameters"
459	write(,) "Unless specified otherwise"
460	print *,'Min. distance Sun-Mars (Mkm)?'
461	call getin("periheli",periheli)
462	write(,) " periheli = ",periheli
463
464	print *,'Max. distance Sun-Mars (Mkm)?'
465	call getin("aphelie",aphelie)
466	write(,) " aphelie = ",aphelie
467
468	print *,'Day of perihelion?'
469	call getin("periday",peri_day)
470	write(,) " periday = ",peri_day
471
472	print *,'Obliquity?'
473	call getin("obliquit",obliquit)
474	write(,) " obliquit = ",obliquit
475	end if
476
477	endif !(.not. startfile_1D )
478
479	c latitude/longitude
480	c ------------------
481	latitude(1)=0 ! default value for latitude
482	PRINT *,'latitude (in degrees) ?'
483	call getin("latitude",latitude(1))
484	write(,) " latitude = ",latitude
485	latitude=latitude*pi/180.E+0
486	longitude=0.E+0
487	longitude=longitude*pi/180.E+0
488
489	! some initializations (some of which have already been
490	! done above!) and loads parameters set in callphys.def
491	! and allocates some arrays
492	!Mars possible matter with dtphys in input and include!!!
493	! Initializations below should mimick what is done in iniphysiq for 3D GCM
494	call init_interface_dyn_phys
495	call init_regular_lonlat(1,1,longitude,latitude,
496	& (/0.,0./),(/0.,0./))
497	call init_geometry(1,longitude,latitude,
498	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
499	& cell_area)
500	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
501	! call init_vertical_layers(nlayer,preff,scaleheight,
502	! & ap,bp,aps,bps,presnivs,pseudoalt)
503	call init_dimphy(1,nlayer) ! Initialize dimphy module
504	call phys_state_var_init(1,llm,nq,tname,
505	. day0,dayn,time,
506	. daysec,dtphys,
507	. rad,g,r,cpp,
508	. nqperes,nqfils)! MVals: variables isotopes
509	call ini_fillgeom(1,latitude,longitude,(/1.0/))
510	call conf_phys(1,llm,nq)
511
512	! in 1D model physics are called every time step
513	! ovverride iphysiq value that has been set by conf_phys
514	if (iphysiq/=1) then
515	write(,) "testphys1d: setting iphysiq=1"
516	iphysiq=1
517	endif
518
519	c Initialize albedo / soil thermal inertia
520	c ----------------------------------------
521	c
522
523	if(.not. startfile_1D ) then
524
525	albedodat(1)=0.2 ! default value for albedodat
526	PRINT *,'Albedo of bare ground ?'
527	call getin("albedo",albedodat(1))
528	write(,) " albedo = ",albedodat(1)
529	albedo(1,1)=albedodat(1)
530
531	inertiedat(1,1)=400 ! default value for inertiedat
532	PRINT *,'Soil thermal inertia (SI) ?'
533	call getin("inertia",inertiedat(1,1))
534	write(,) " inertia = ",inertiedat(1,1)
535
536	ice_depth = -1 ! default value: no ice
537	CALL getin("subsurface_ice_depth",ice_depth)
538
539	z0(1)=z0_default ! default value for roughness
540	write(,) 'Surface roughness length z0 (m)?'
541	call getin("z0",z0(1))
542	write(,) " z0 = ",z0(1)
543
544	endif !(.not. startfile_1D )
545
546	! Initialize local slope parameters (only matters if "callslope"
547	! is .true. in callphys.def)
548	! slope inclination angle (deg) 0: horizontal, 90: vertical
549	theta_sl(1)=0.0 ! default: no inclination
550	call getin("slope_inclination",theta_sl(1))
551	! slope orientation (deg)
552	! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
553	psi_sl(1)=0.0 ! default value
554	call getin("slope_orientation",psi_sl(1))
555
556	! sub-slopes parameters (assuming no sub-slopes distribution for now).
557	def_slope(1)=-90 ! minimum slope angle
558	def_slope(2)=90 ! maximum slope angle
559	subslope_dist(1,1)=1 ! fraction of subslopes in mesh
560	c
561	c for the gravity wave scheme
562	c ---------------------------------
563	c
564	zmea(1)=0.E+0
565	zstd(1)=0.E+0
566	zsig(1)=0.E+0
567	zgam(1)=0.E+0
568	zthe(1)=0.E+0
569	c
570	c for the slope wind scheme
571	c ---------------------------------
572	c
573	hmons(1)=0.E+0
574	PRINT *,'hmons is initialized to ',hmons(1)
575	summit(1)=0.E+0
576	PRINT *,'summit is initialized to ',summit(1)
577	base(1)=0.E+0
578	c
579	c Default values initializing the coefficients calculated later
580	c ---------------------------------
581	c
582	tauscaling(1)=1. ! calculated in aeropacity_mod.F
583	totcloudfrac(1)=1. ! calculated in watercloud_mod.F
584
585	c Specific initializations for "physiq"
586	c -------------------------------------
587	c surface geopotential is not used (or useful) since in 1D
588	c everything is controled by surface pressure
589	phisfi(1)=0.E+0
590
591	c Initialization to take into account prescribed winds
592	c ------------------------------------------------------
593	ptif=2.E+0omegsinlat(1)
594
595	c geostrophic wind
596	gru=10. ! default value for gru
597	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
598	call getin("u",gru)
599	write(,) " u = ",gru
600	grv=0. !default value for grv
601	PRINT *,'meridional northward component of the geostrophic',
602	&' wind (m/s) ?'
603	call getin("v",grv)
604	write(,) " v = ",grv
605
606	c Initialize winds for first time step
607	DO ilayer=1,nlayer
608	u(ilayer)=gru
609	v(ilayer)=grv
610	w(ilayer)=0 ! default: no vertical wind
611	ENDDO
612
613	c Initialize turbulente kinetic energy
614	DO ilevel=1,nlevel
615	q2(ilevel)=0.E+0
616	ENDDO
617
618	c CO2 ice on the surface
619	c -------------------
620	! get the index of co2 tracer (not known at this stage)
621	igcm_co2=0
622	do iq=1,nq
623	if (trim(tname(iq))=="co2") then
624	igcm_co2=iq
625	endif
626	enddo
627	if (igcm_co2==0) then
628	write(,) "testphys1d error, missing co2 tracer!"
629	stop
630	endif
631
632	if(.not. startfile_1D ) then
633	qsurf(igcm_co2)=0.E+0 ! default value for co2ice
634	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
635	call getin("co2ice",qsurf(igcm_co2))
636	write(,) " co2ice = ",qsurf(igcm_co2)
637	endif !(.not. startfile_1D )
638
639	c
640	c emissivity
641	c ----------
642	if(.not. startfile_1D ) then
643	emis=emissiv
644	IF (qsurf(igcm_co2).eq.1.E+0) THEN
645	emis=emisice(1) ! northern hemisphere
646	IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
647	ENDIF
648	endif !(.not. startfile_1D )
649
650
651	c Compute pressures and altitudes of atmospheric levels
652	c ----------------------------------------------------------------
653
654	c Vertical Coordinates
655	c """"""""""""""""""""
656	hybrid=.true.
657	PRINT *,'Hybrid coordinates ?'
658	call getin("hybrid",hybrid)
659	write(,) " hybrid = ", hybrid
660
661	CALL disvert_noterre
662	! now that disvert has been called, initialize module vertical_layers_mod
663	call init_vertical_layers(nlayer,preff,scaleheight,
664	& ap,bp,aps,bps,presnivs,pseudoalt)
665
666	DO ilevel=1,nlevel
667	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
668	ENDDO
669
670	DO ilayer=1,nlayer
671	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
672	ENDDO
673
674	DO ilayer=1,nlayer
675	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
676	& /g
677	ENDDO
678
679
680	c Initialize temperature profile
681	c --------------------------------------
682	pks=psurf**rcp
683
684	c altitude in km in profile: divide zlay by 1000
685	tmp1(0)=0.E+0
686	DO ilayer=1,nlayer
687	tmp1(ilayer)=zlay(ilayer)/1000.E+0
688	ENDDO
689
690	call profile(nlayer+1,tmp1,tmp2)
691
692	tsurf=tmp2(0)
693	DO ilayer=1,nlayer
694	temp(ilayer)=tmp2(ilayer)
695	ENDDO
696
697
698	! Initialize soil properties and temperature
699	! ------------------------------------------
700	volcapa=1.e6 ! volumetric heat capacity
701
702	if(.not. startfile_1D ) then
703
704	! Initialize depths
705	! -----------------
706	do isoil=1,nsoil
707	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
708	enddo
709
710	! Creating the new soil inertia table if there is subsurface ice :
711	IF (ice_depth.gt.0) THEN
712	iref = 1 ! ice/regolith boundary index
713	IF (ice_depth.lt.layer(1)) THEN
714	inertiedat(1,1) = sqrt( layer(1) /
715	& ( (ice_depth/inertiedat(1,1)**2) +
716	& ((layer(1)-ice_depth)/inertieice**2) ) )
717	DO isoil=1,nsoil
718	inertiedat(1,isoil) = inertieice
719	ENDDO
720	ELSE ! searching for the ice/regolith boundary :
721	DO isoil=1,nsoil
722	IF ((ice_depth.ge.layer(isoil)).and.
723	& (ice_depth.lt.layer(isoil+1))) THEN
724	iref=isoil+1
725	EXIT
726	ENDIF
727	ENDDO
728	! We then change the soil inertia table :
729	DO isoil=1,iref-1
730	inertiedat(1,isoil) = inertiedat(1,1)
731	ENDDO
732	! We compute the transition in layer(iref)
733	inertiedat(1,iref) = sqrt( (layer(iref)-layer(iref-1)) /
734	& ( ((ice_depth-layer(iref-1))/inertiedat(1,1)**2) +
735	& ((layer(iref)-ice_depth)/inertieice**2) ) )
736	! Finally, we compute the underlying ice :
737	DO isoil=iref+1,nsoil
738	inertiedat(1,isoil) = inertieice
739	ENDDO
740	ENDIF ! (ice_depth.lt.layer(1))
741	ELSE ! ice_depth <0 all is set to surface thermal inertia
742	DO isoil=1,nsoil
743	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
744	ENDDO
745	ENDIF ! ice_depth.gt.0
746
747	inertiesoil(1,:,1) = inertiedat(1,:)
748
749	DO isoil = 1,nsoil
750	tsoil(isoil)=tsurf(1) ! soil temperature
751	ENDDO
752
753	endif !(.not. startfile_1D )
754
755	flux_geo_tmp=0.
756	call getin("flux_geo",flux_geo_tmp)
757	flux_geo(:,:) = flux_geo_tmp
758
759	! Initialize depths
760	! -----------------
761	do isoil=0,nsoil-1
762	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
763	enddo
764	do isoil=1,nsoil
765	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
766	enddo
767
768	c Initialize traceurs
769	c ---------------------------
770
771	if (photochem.or.callthermos) then
772	write(,) 'Initializing chemical species'
773	! flagthermo=0: initialize over all atmospheric layers
774	flagthermo=0
775	! check if "h2o_vap" has already been initialized
776	! (it has been if there is a "profile_h2o_vap" file around)
777	inquire(file="profile_h2o_vap",exist=present)
778	if (present) then
779	flagh2o=0 ! 0: do not initialize h2o_vap
780	else
781	flagh2o=1 ! 1: initialize h2o_vap in inichim_newstart
782	endif
783
784	! hack to accomodate that inichim_newstart assumes that
785	! q & psurf arrays are on the dynamics scalar grid
786	allocate(qdyn(2,1,llm,nq),psdyn(2,1))
787	qdyn(1,1,1:llm,1:nq)=q(1:llm,1:nq)
788	psdyn(1:2,1)=psurf
789	call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn,
790	$ flagh2o,flagthermo)
791	q(1:llm,1:nq)=qdyn(1,1,1:llm,1:nq)
792	endif
793
794	c Check if the surface is a water ice reservoir
795	c --------------------------------------------------
796	if(.not. startfile_1D ) then
797	watercap(1,:)=0 ! Initialize watercap
798	endif !(.not. startfile_1D )
799	watercaptag(1)=.false. ! Default: no water ice reservoir
800	print *,'Water ice cap on ground ?'
801	call getin("watercaptag",watercaptag)
802	write(,) " watercaptag = ",watercaptag
803
804	c Check if the atmospheric water profile is specified
805	c ---------------------------------------------------
806	! Adding an option to force atmospheric water values JN
807	atm_wat_profile=-1 ! Default: free atm wat profile
808	print *,'Force atmospheric water vapor profile ?'
809	call getin("atm_wat_profile",atm_wat_profile)
810	write(,) "atm_wat_profile = ", atm_wat_profile
811	if (atm_wat_profile.EQ.-1) then
812	write(,) "Free atmospheric water vapor profile"
813	write(,) "Total water is conserved on the column"
814	else if (atm_wat_profile.EQ.0) then
815	write(,) "Dry atmospheric water vapor profile"
816	else if ((atm_wat_profile.GT.0).and.(atm_wat_profile.LE.1)) then
817	write(,) "Prescribed atmospheric water vapor MMR profile"
818	write(,) "Unless it reaches saturation (maximal value)"
819	write(,) "MMR chosen : ", atm_wat_profile
820	endif
821
822
823	c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
824	c ----------------------------------------------------------------
825	c (output done in "writeg1d", typically called by "physiq.F")
826
827	g1d_nlayer=nlayer
828	g1d_nomfich='g1d.dat'
829	g1d_unitfich=40
830	g1d_nomctl='g1d.ctl'
831	g1d_unitctl=41
832	g1d_premier=.true.
833	g2d_premier=.true.
834
835	c Write a "startfi" file
836	c --------------------
837	c This file will be read during the first call to "physiq".
838	c It is needed to transfert physics variables to "physiq"...
839
840	if(.not. startfile_1D ) then
841
842	call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,
843	& llm,nq,dtphys,float(day0),0.,cell_area,
844	& albedodat,inertiedat,def_slope,subslope_dist)
845	call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
846	& dtphys,time,
847	& tsurf,tsoil,inertiesoil,albedo,emis,
848	& q2,qsurf,tauscaling,
849	& totcloudfrac,wstar,watercap)
850	endif !(.not. startfile_1D )
851
852	c=======================================================================
853	c 1D MODEL TIME STEPPING LOOP
854	c=======================================================================
855	c
856	firstcall=.true.
857	lastcall=.false.
858	DO idt=1,ndt
859	IF (idt.eq.ndt) lastcall=.true.
860	c IF (idt.eq.ndt-day_step-1) then !test
861	c lastcall=.true.
862	c call solarlong(day*1.0,zls)
863	c write(103,) 'Ls=',zls180./pi
864	c write(103,) 'Lat=', latitude(1)180./pi
865	c write(103,) 'Tau=', tauvis/odprefpsurf
866	c write(103,*) 'RunEnd - Atmos. Temp. File'
867	c write(103,*) 'RunEnd - Atmos. Temp. File'
868	c write(104,) 'Ls=',zls180./pi
869	c write(104,*) 'Lat=', latitude(1)
870	c write(104,) 'Tau=', tauvis/odprefpsurf
871	c write(104,*) 'RunEnd - Atmos. Temp. File'
872	c ENDIF
873
874	c compute geopotential
875	c ~~~~~~~~~~~~~~~~~~~~~
876	DO ilayer=1,nlayer
877	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
878	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
879	ENDDO
880	phi(1)=pksh(1)(1.E+0-s(1))
881	DO ilayer=2,nlayer
882	phi(ilayer)=phi(ilayer-1)+
883	& pks(h(ilayer-1)+h(ilayer)).5E+0
884	& *(s(ilayer-1)-s(ilayer))
885
886	ENDDO
887
888	c Force atmospheric water if asked
889	c Added "atm_wat_profile" flag (JN)
890	c -------------------------------------
891	call watersat(nlayer,temp,play,zqsat)
892	DO iq = 1, nq
893	IF ((iq.eq.igcm_h2o_vap).and.(atm_wat_profile.eq.0)) THEN
894	DO ilayer=1,nlayer
895	q(ilayer,igcm_h2o_vap)=0.
896	c write(,) "atm_wat_profile dry"
897	ENDDO! ilayer=1,nlayer
898	ELSE IF((iq.eq.igcm_h2o_vap).and.(atm_wat_profile.gt.0)
899	& .and.(atm_wat_profile.le.1)) THEN
900	DO ilayer=1,nlayer
901	q(ilayer,igcm_h2o_vap)=min(zqsat(ilayer),atm_wat_profile)
902	c write(,) "atm_wat_profile wet"
903	ENDDO! ilayer=1,nlayer
904	ELSE IF ((iq.eq.igcm_h2o_ice).and.
905	& (atm_wat_profile.ne.-1)) THEN
906	DO ilayer=1,nlayer
907	q(ilayer,igcm_h2o_ice)=0.
908	c write(,) "atm_wat_profile : reset ice"
909	ENDDO! ilayer=1,nlayer
910	ENDIF !((iq.eq.igcm_h2o_vap).and.(atm_wat_profile.eq.2)) THEN
911	ENDDO
912	CALL write_output('qsat',
913	& 'qsat',
914	& 'kg/kg',zqsat(:))
915
916
917
918	! write(,) "testphys1d avant q", q(1,:)
919	c call physics
920	c --------------------
921	CALL physiq (1,llm,nq,
922	, firstcall,lastcall,
923	, day,time,dtphys,
924	, plev,play,phi,
925	, u, v,temp, q,
926	, w,
927	C - outputs
928	s du, dv, dtemp, dq,dpsurf)
929	! write(,) "testphys1d apres q", q(1,:)
930
931
932	c wind increment : specific for 1D
933	c --------------------------------
934
935	c The physics compute the tendencies on u and v,
936	c here we just add Coriolos effect
937	c
938	c DO ilayer=1,nlayer
939	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
940	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
941	c ENDDO
942
943	c For some tests : No coriolis force at equator
944	c if(latitude(1).eq.0.) then
945	DO ilayer=1,nlayer
946	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
947	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
948	ENDDO
949	c end if
950	c
951	c
952	c Compute time for next time step
953	c ---------------------------------------
954	firstcall=.false.
955	time=time+dtphys/daysec
956	IF (time.gt.1.E+0) then
957	time=time-1.E+0
958	day=day+1
959	ENDIF
960
961	c compute winds and temperature for next time step
962	c ----------------------------------------------------------
963
964	DO ilayer=1,nlayer
965	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
966	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
967	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
968	ENDDO
969
970	c compute pressure for next time step
971	c ----------------------------------------------------------
972
973	psurf=psurf+dtphys*dpsurf(1) ! surface pressure change
974	DO ilevel=1,nlevel
975	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
976	ENDDO
977	DO ilayer=1,nlayer
978	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
979	ENDDO
980
981	! increment tracers
982	DO iq = 1, nq
983	DO ilayer=1,nlayer
984	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
985	ENDDO
986	ENDDO
987
988	ENDDO ! of idt=1,ndt ! end of time stepping loop
989
990	c ========================================================
991	c OUTPUTS
992	c ========================================================
993
994	c finalize and close grads files "g1d.dat" and "g1d.ctl"
995
996	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
997	CALL endg1d(1,nlayer,zlay/1000.,ndt)
998
999	write(,) "testphys1d: Everything is cool."
1000
1001	END
1002
1003	c***********************************************************************
1004	c***********************************************************************
1005	c Dummy subroutines used only in 3D, but required to
1006	c compile testphys1d (to cleanly use writediagfi)
1007
1008	subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
1009
1010	IMPLICIT NONE
1011
1012	INTEGER im,jm,ngrid,nfield
1013	REAL pdyn(im,jm,nfield)
1014	REAL pfi(ngrid,nfield)
1015
1016	if (ngrid.ne.1) then
1017	write(,) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
1018	stop
1019	endif
1020
1021	pdyn(1,1,1:nfield)=pfi(1,1:nfield)
1022
1023	end
1024
1025	c***********************************************************************
1026	c***********************************************************************
1027

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

Download in other formats:

Original Format