Context Navigation

testphys1d.F @ 2980

Last change on this file since 2980 was 2980, checked in by romain.vande, 19 months ago

Mars PEM :

Adapt PEM to 1d runs.
Cleaning of names and unused variables.
Correct minor errors.
Adapt and correct reshape_xios_output utilitary for 1d diagfi output.

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)
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