Context Navigation

testphys1d.F @ 2953

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

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

Download in other formats:

Original Format