Context Navigation

testphys1d.F @ 2997

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

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

Download in other formats:

Original Format