Context Navigation

testphys1d.F @ 1884

Last change on this file since 1884 was 1711, checked in by mvals, 8 years ago

Mars GCM

implementation of a sub-grid water cloud fraction scheme (by A. Pottier): If CLFvarying is set to true in callphys.def, then the sub-grid cloud fraction routine is applied in watercloud.F and aeropacity.F.
accordingly modified files are: watercloud.F, aeropacity.F, callradite.F, conf_phys.F, phyetat0.F, phyredem.F90, physiq_mod.F, testphys1d.F, callkeys.h, newtsart.F, start2archive.F, lect_start_archive.F
added file: tcondwater.F90, used by watercloud.F to calculate the condensation temperature of water
watercloud.F, aeropacity.F, callradite.F are converted to module files as watercloud_mod.F, aeropacity_mod.F, callradite_mod.F

File size: 29.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
8	use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx
9	use comgeomfi_h, only: sinlat, ini_fillgeom
10	use surfdat_h, only: albedodat, z0_default, emissiv, emisice,
11	& albedice, iceradius, dtemisice, z0,
12	& zmea, zstd, zsig, zgam, zthe, phisfi,
13	& watercaptag
14	use slope_mod, only: theta_sl, psi_sl
15	use phyredem, only: physdem0,physdem1
16	use geometry_mod, only: init_geometry
17	use planete_h, only: year_day, periheli, aphelie, peri_day,
18	& obliquit, emin_turb, lmixmin
19	use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp
20	use time_phylmdz_mod, only: daysec, dtphys, day_step,
21	& ecritphy, iphysiq
22	use dimradmars_mod, only: tauscaling,tauvis
23	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig,
24	& presnivs,pseudoalt,scaleheight
25	USE vertical_layers_mod, ONLY: init_vertical_layers
26	USE logic_mod, ONLY: hybrid
27	use physics_distribution_mod, only: init_physics_distribution
28	use regular_lonlat_mod, only: init_regular_lonlat
29	use mod_interface_dyn_phys, only: init_interface_dyn_phys
30	USE phys_state_var_init_mod, ONLY: phys_state_var_init
31	USE physiq_mod, ONLY: physiq
32	IMPLICIT NONE
33
34	c=======================================================================
35	c subject:
36	c --------
37	c PROGRAM useful to run physical part of the martian GCM in a 1D column
38	c
39	c Can be compiled with a command like (e.g. for 25 layers)
40	c "makegcm -p mars -d 25 testphys1d"
41	c It requires the files "testphys1d.def" "callphys.def"
42	c and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
43	c and a file describing the sigma layers (e.g. "z2sig.def")
44	c
45	c author: Frederic Hourdin, R.Fournier,F.Forget
46	c -------
47	c
48	c update: 12/06/2003 including chemistry (S. Lebonnois)
49	c and water ice (F. Montmessin)
50	c
51	c=======================================================================
52
53	include "dimensions.h"
54	integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm)
55	integer, parameter :: nlayer = llm
56	!#include "dimradmars.h"
57	!#include "comgeomfi.h"
58	!#include "surfdat.h"
59	!#include "slope.h"
60	!#include "comsoil.h"
61	!#include "comdiurn.h"
62	include "callkeys.h"
63	!#include "comsaison.h"
64	!#include "control.h"
65	include "netcdf.inc"
66	include "comg1d.h"
67	!#include "advtrac.h"
68
69	c --------------------------------------------------------------
70	c Declarations
71	c --------------------------------------------------------------
72	c
73	INTEGER unitstart ! unite d'ecriture de "startfi"
74	INTEGER nlevel,nsoil,ndt
75	INTEGER ilayer,ilevel,isoil,idt,iq
76	LOGICAl firstcall,lastcall
77	c
78	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
79	c
80	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
81	REAL day ! date durant le run
82	REAL time ! time (0<time<1 ; time=0.5 a midi)
83	REAL play(nlayer) ! Pressure at the middle of the layers (Pa)
84	REAL plev(nlayer+1) ! intermediate pressure levels (pa)
85	REAL psurf,tsurf(1)
86	REAL u(nlayer),v(nlayer) ! zonal, meridional wind
87	REAL gru,grv ! prescribed "geostrophic" background wind
88	REAL temp(nlayer) ! temperature at the middle of the layers
89	REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
90	REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
91	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
92	REAL co2ice(1) ! co2ice layer (kg.m-2)
93	REAL emis(1) ! surface layer
94	REAL q2(nlayer+1) ! Turbulent Kinetic Energy
95	REAL zlay(nlayer) ! altitude estimee dans les couches (km)
96
97	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
98	REAL du(nlayer),dv(nlayer),dtemp(nlayer)
99	REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer)
100	REAL dpsurf(1)
101	REAL,ALLOCATABLE :: dq(:,:)
102	REAL,ALLOCATABLE :: dqdyn(:,:)
103
104	c Various intermediate variables
105	INTEGER thermo
106	REAL zls
107	REAL phi(nlayer),h(nlayer),s(nlayer)
108	REAL pks, ptif, w(nlayer)
109	REAL qtotinit,qtot
110	real,allocatable :: mqtot(:)
111	INTEGER ierr, aslun
112	REAL tmp1(0:nlayer),tmp2(0:nlayer)
113	integer :: nq=1 ! number of tracers
114	real :: latitude(1), longitude(1), cell_area(1)
115
116	character*2 str2
117	character (len=7) :: str7
118	character(len=44) :: txt
119
120	!MV17 nuages alizee
121	! included by AP14 for cloud fraction computations- see BC generic
122	!real,allocatable,dimension(:,:),save :: cloudfrac
123	real,allocatable,dimension(:),save :: totcloudfrac
124	c=======================================================================
125
126	c=======================================================================
127	c INITIALISATION
128	c=======================================================================
129	! initialize "serial/parallel" related stuff
130	! CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
131	! CALL init_phys_lmdz(1,1,llm,1,(/1/))
132	! call initcomgeomphy
133
134	c ------------------------------------------------------
135	c Prescribed constants to be set here
136	c ------------------------------------------------------
137
138	pi=2.E+0*asin(1.E+0)
139
140	c Mars planetary constants
141	c ----------------------------
142	rad=3397200. ! mars radius (m) ~3397200 m
143	daysec=88775. ! length of a sol (s) ~88775 s
144	omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1)
145	g=3.72 ! gravity (m.s-2) ~3.72
146	mugaz=43.49 ! atmosphere mola mass (g.mol-1) ~43.49
147	rcp=.256793 ! = r/cp ~0.256793
148	r= 8.314511E+0 *1000.E+0/mugaz
149	cpp= r/rcp
150	year_day = 669 ! lenght of year (sols) ~668.6
151	periheli = 206.66 ! minimum sun-mars distance (Mkm) ~206.66
152	aphelie = 249.22 ! maximum sun-mars distance (Mkm) ~249.22
153	peri_day = 485. ! perihelion date (sols since N. Spring)
154	obliquit = 25.2 ! Obliquity (deg) ~25.2
155
156	c Planetary Boundary Layer and Turbulence parameters
157	c --------------------------------------------------
158	z0_default = 1.e-2 ! surface roughness (m) ~0.01
159	emin_turb = 1.e-6 ! minimal turbulent energy ~1.e-8
160	lmixmin = 30 ! mixing length ~100
161
162	c cap properties and surface emissivities
163	c ----------------------------------------------------
164	emissiv= 0.95 ! Bare ground emissivity ~.95
165	emisice(1)=0.95 ! Northern cap emissivity
166	emisice(2)=0.95 ! Southern cap emisssivity
167	albedice(1)=0.5 ! Northern cap albedo
168	albedice(2)=0.5 ! Southern cap albedo
169	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
170	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
171	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
172	dtemisice(2) = 2. ! time scale for snow metamorphism (south
173
174
175	c ------------------------------------------------------
176	c Loading run parameters from "run.def" file
177	c ------------------------------------------------------
178
179
180	! check if 'run.def' file is around (otherwise reading parameters
181	! from callphys.def via getin() routine won't work.
182	open(99,file='run.def',status='old',form='formatted',
183	& iostat=ierr)
184	if (ierr.ne.0) then
185	write(,) 'Cannot find required file "run.def"'
186	write(,) ' (which should contain some input parameters'
187	write(,) ' along with the following line:'
188	write(,) ' INCLUDEDEF=callphys.def'
189	write(,) ' )'
190	write(,) ' ... might as well stop here ...'
191	stop
192	else
193	close(99)
194	endif
195
196	! check if we are going to run with or without tracers
197	write(,) "Run with or without tracer transport ?"
198	tracer=.false. ! default value
199	call getin("tracer",tracer)
200	write(,) " tracer = ",tracer
201
202	! while we're at it, check if there is a 'traceur.def' file
203	! and process it.
204	if (tracer) then
205	! load tracer names from file 'traceur.def'
206	open(90,file='traceur.def',status='old',form='formatted',
207	& iostat=ierr)
208	if (ierr.ne.0) then
209	write(,) 'Cannot find required file "traceur.def"'
210	write(,) ' If you want to run with tracers, I need it'
211	write(,) ' ... might as well stop here ...'
212	stop
213	else
214	write(,) "testphys1d: Reading file traceur.def"
215	! read number of tracers:
216	read(90,*,iostat=ierr) nq
217	nqtot=nq ! set value of nqtot (in infotrac module) as nq
218	if (ierr.ne.0) then
219	write(,) "testphys1d: error reading number of tracers"
220	write(,) " (first line of traceur.def) "
221	stop
222	endif
223	endif
224	! allocate arrays:
225	allocate(tname(nq))
226	allocate(q(nlayer,nq))
227	allocate(qsurf(nq))
228	allocate(dq(nlayer,nq))
229	allocate(dqdyn(nlayer,nq))
230	allocate(mqtot(nq))
231	!MV17 nuages alizee
232	!allocate(cloudfrac(ngrid,nlayermx))!essai allocate AP14
233	allocate(totcloudfrac(ngrid))
234
235	! read tracer names from file traceur.def
236	do iq=1,nq
237	read(90,*,iostat=ierr) tname(iq)
238	if (ierr.ne.0) then
239	write(,) 'testphys1d: error reading tracer names...'
240	stop
241	endif
242	enddo
243	close(90)
244
245	! initialize tracers here:
246	write(,) "testphys1d: initializing tracers"
247	q(:,:)=0 ! default, set everything to zero
248	qsurf(:)=0
249	! "smarter" initialization of some tracers
250	! (get values from "profile_*" files, if these are available)
251	do iq=1,nq
252	txt=""
253	write(txt,"(a)") tname(iq)
254	write(,)" tracer:",trim(txt)
255	! CO2
256	if (txt.eq."co2") then
257	q(:,iq)=0.95 ! kg /kg of atmosphere
258	qsurf(iq)=0. ! kg/m2 (not used for CO2)
259	! even better, look for a "profile_co2" input file
260	open(91,file='profile_co2',status='old',
261	& form='formatted',iostat=ierr)
262	if (ierr.eq.0) then
263	read(91,*) qsurf(iq)
264	do ilayer=1,nlayer
265	read(91,*) q(ilayer,iq)
266	enddo
267	endif
268	close(91)
269	endif ! of if (txt.eq."co2")
270	! Allow for an initial profile of argon
271	! Can also be used to introduce a decaying tracer
272	! in the 1D (TBD) to study thermals
273	if (txt.eq."ar") then
274	!look for a "profile_ar" input file
275	open(91,file='profile_ar',status='old',
276	& form='formatted',iostat=ierr)
277	if (ierr.eq.0) then
278	read(91,*) qsurf(iq)
279	do ilayer=1,nlayer
280	read(91,*) q(ilayer,iq)
281	enddo
282	else
283	write(,) "No profile_ar file!"
284	endif
285	close(91)
286	endif ! of if (txt.eq."ar")
287
288	! WATER VAPOUR
289	if (txt.eq."h2o_vap") then
290	!look for a "profile_h2o_vap" input file
291	open(91,file='profile_h2o_vap',status='old',
292	& form='formatted',iostat=ierr)
293	if (ierr.eq.0) then
294	read(91,*) qsurf(iq)
295	do ilayer=1,nlayer
296	read(91,*) q(ilayer,iq)
297	enddo
298	else
299	write(,) "No profile_h2o_vap file!"
300	endif
301	close(91)
302	endif ! of if (txt.eq."h2o_ice")
303	! WATER ICE
304	if (txt.eq."h2o_ice") then
305	!look for a "profile_h2o_vap" input file
306	open(91,file='profile_h2o_ice',status='old',
307	& form='formatted',iostat=ierr)
308	if (ierr.eq.0) then
309	read(91,*) qsurf(iq)
310	do ilayer=1,nlayer
311	read(91,*) q(ilayer,iq)
312	enddo
313	else
314	write(,) "No profile_h2o_ice file!"
315	endif
316	close(91)
317	endif ! of if (txt.eq."h2o_ice")
318	! DUST
319	!if (txt(1:4).eq."dust") then
320	! q(:,iq)=0.4 ! kg/kg of atmosphere
321	! qsurf(iq)=100 ! kg/m2
322	!endif
323	! DUST MMR
324	if (txt.eq."dust_mass") then
325	!look for a "profile_dust_mass" input file
326	open(91,file='profile_dust_mass',status='old',
327	& form='formatted',iostat=ierr)
328	if (ierr.eq.0) then
329	read(91,*) qsurf(iq)
330	do ilayer=1,nlayer
331	read(91,*) q(ilayer,iq)
332	! write(,) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
333	enddo
334	else
335	write(,) "No profile_dust_mass file!"
336	endif
337	close(91)
338	endif ! of if (txt.eq."dust_mass")
339	! DUST NUMBER
340	if (txt.eq."dust_number") then
341	!look for a "profile_dust_number" input file
342	open(91,file='profile_dust_number',status='old',
343	& form='formatted',iostat=ierr)
344	if (ierr.eq.0) then
345	read(91,*) qsurf(iq)
346	do ilayer=1,nlayer
347	read(91,*) q(ilayer,iq)
348	enddo
349	else
350	write(,) "No profile_dust_number file!"
351	endif
352	close(91)
353	endif ! of if (txt.eq."dust_number")
354	! NB: some more initializations (chemistry) is done later
355	! CCN MASS
356	if (txt.eq."ccn_mass") then
357	!look for a "profile_ccn_mass" input file
358	open(91,file='profile_ccn_mass',status='old',
359	& form='formatted',iostat=ierr)
360	if (ierr.eq.0) then
361	read(91,*) qsurf(iq)
362	do ilayer=1,nlayer
363	read(91,*) q(ilayer,iq)
364	enddo
365	else
366	write(,) "No profile_ccn_mass file!"
367	endif
368	close(91)
369	endif ! of if (txt.eq."ccn_mass")
370	! CCN NUMBER
371	if (txt.eq."ccn_number") then
372	!look for a "profile_ccn_number" input file
373	open(91,file='profile_ccn_number',status='old',
374	& form='formatted',iostat=ierr)
375	if (ierr.eq.0) then
376	read(91,*) qsurf(iq)
377	do ilayer=1,nlayer
378	read(91,*) q(ilayer,iq)
379	enddo
380	else
381	write(,) "No profile_ccn_number file!"
382	endif
383	close(91)
384	endif ! of if (txt.eq."ccn_number")
385	enddo ! of do iq=1,nq
386
387	else
388	! we still need to set (dummy) tracer number and names for physdem1
389	nq=1
390	nqtot=nq ! set value of nqtot (in infotrac module) as nq
391	! allocate arrays:
392	allocate(tname(nq))
393	allocate(q(nlayer,nq))
394	allocate(qsurf(nq))
395	allocate(dq(nlayer,nq))
396	allocate(dqdyn(nlayer,nq))
397	allocate(mqtot(nq))
398	do iq=1,nq
399	write(str7,'(a1,i2.2)')'t',iq
400	tname(iq)=str7
401	enddo
402	! and just to be clean, also initialize tracers to zero for physdem1
403	q(:,:)=0
404	qsurf(:)=0
405	endif ! of if (tracer)
406
407	!write(,) "testphys1d q", q(1,:)
408	!write(,) "testphys1d qsurf", qsurf
409
410	c Date and local time at beginning of run
411	c ---------------------------------------
412	c Date (in sols since spring solstice) at beginning of run
413	day0 = 0 ! default value for day0
414	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
415	call getin("day0",day0)
416	day=float(day0)
417	write(,) " day0 = ",day0
418	c Local time at beginning of run
419	time=0 ! default value for time
420	write(,)'Initial local time (in hours, between 0 and 24)?'
421	call getin("time",time)
422	write(,)" time = ",time
423	time=time/24.E+0 ! convert time (hours) to fraction of sol
424
425	c Discretization (Definition of grid and time steps)
426	c --------------
427	c
428	nlevel=nlayer+1
429	nsoil=nsoilmx
430
431	day_step=48 ! default value for day_step
432	PRINT *,'Number of time steps per sol ?'
433	call getin("day_step",day_step)
434	write(,) " day_step = ",day_step
435
436	ecritphy=day_step ! default value for ecritphy, output every time step
437
438	ndt=10 ! default value for ndt
439	PRINT *,'Number of sols to run ?'
440	call getin("ndt",ndt)
441	write(,) " ndt = ",ndt
442
443	ndt=ndt*day_step
444	dtphys=daysec/day_step
445
446	c Imposed surface pressure
447	c ------------------------------------
448	c
449	psurf=610. ! default value for psurf
450	PRINT *,'Surface pressure (Pa) ?'
451	call getin("psurf",psurf)
452	write(,) " psurf = ",psurf
453	c Reference pressures
454	pa=20. ! transition pressure (for hybrid coord.)
455	preff=610. ! reference surface pressure
456
457	c Aerosol properties
458	c --------------------------------
459	tauvis=0.2 ! default value for tauvis (dust opacity)
460	write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
461	call getin("tauvis",tauvis)
462	write(,) " tauvis = ",tauvis
463
464	c Orbital parameters
465	c ------------------
466	print *,'Min. distance Sun-Mars (Mkm)?'
467	call getin("periheli",periheli)
468	write(,) " periheli = ",periheli
469
470	print *,'Max. distance Sun-Mars (Mkm)?'
471	call getin("aphelie",aphelie)
472	write(,) " aphelie = ",aphelie
473
474	print *,'Day of perihelion?'
475	call getin("periday",peri_day)
476	write(,) " periday = ",peri_day
477
478	print *,'Obliquity?'
479	call getin("obliquit",obliquit)
480	write(,) " obliquit = ",obliquit
481
482	c latitude/longitude
483	c ------------------
484	latitude(1)=0 ! default value for latitude
485	PRINT *,'latitude (in degrees) ?'
486	call getin("latitude",latitude(1))
487	write(,) " latitude = ",latitude
488	latitude=latitude*pi/180.E+0
489	longitude=0.E+0
490	longitude=longitude*pi/180.E+0
491
492	! some initializations (some of which have already been
493	! done above!) and loads parameters set in callphys.def
494	! and allocates some arrays
495	!Mars possible matter with dtphys in input and include!!!
496	! Initializations below should mimick what is done in iniphysiq for 3D GCM
497	call init_physics_distribution(regular_lonlat,4,
498	& 1,1,1,nlayer,1)
499	call init_interface_dyn_phys
500	call init_regular_lonlat(1,1,longitude,latitude,
501	& (/0.,0./),(/0.,0./))
502	call init_geometry(1,longitude,latitude,
503	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
504	& cell_area)
505	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
506	! call init_vertical_layers(nlayer,preff,scaleheight,
507	! & ap,bp,aps,bps,presnivs,pseudoalt)
508	call init_dimphy(1,nlayer) ! Initialize dimphy module
509	call phys_state_var_init(1,llm,nq,tname,
510	. day0,time,daysec,dtphys,rad,g,r,cpp)
511	call ini_fillgeom(1,latitude,longitude,(/1.0/))
512	call conf_phys(1,llm,nq)
513
514	! in 1D model physics are called every time step
515	! ovverride iphysiq value that has been set by conf_phys
516	if (iphysiq/=1) then
517	write(,) "testphys1d: setting iphysiq=1"
518	iphysiq=1
519	endif
520
521	c Initialize albedo / soil thermal inertia
522	c ----------------------------------------
523	c
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
529	inertiedat(1,1)=400 ! default value for inertiedat
530	PRINT *,'Soil thermal inertia (SI) ?'
531	call getin("inertia",inertiedat(1,1))
532	write(,) " inertia = ",inertiedat(1,1)
533
534	z0(1)=z0_default ! default value for roughness
535	write(,) 'Surface roughness length z0 (m)?'
536	call getin("z0",z0(1))
537	write(,) " z0 = ",z0(1)
538
539	! Initialize local slope parameters (only matters if "callslope"
540	! is .true. in callphys.def)
541	! slope inclination angle (deg) 0: horizontal, 90: vertical
542	theta_sl(1)=0.0 ! default: no inclination
543	call getin("slope_inclination",theta_sl(1))
544	! slope orientation (deg)
545	! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
546	psi_sl(1)=0.0 ! default value
547	call getin("slope_orientation",psi_sl(1))
548
549	c
550	c for the gravity wave scheme
551	c ---------------------------------
552	c
553	zmea(1)=0.E+0
554	zstd(1)=0.E+0
555	zsig(1)=0.E+0
556	zgam(1)=0.E+0
557	zthe(1)=0.E+0
558
559
560	c Specific initializations for "physiq"
561	c -------------------------------------
562	c mesh surface (not a very usefull quantity in 1D)
563	cell_area(1)=1.E+0
564
565	c surface geopotential is not used (or useful) since in 1D
566	c everything is controled by surface pressure
567	phisfi(1)=0.E+0
568
569	c Initialization to take into account prescribed winds
570	c ------------------------------------------------------
571	ptif=2.E+0omegsinlat(1)
572
573	c geostrophic wind
574	gru=10. ! default value for gru
575	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
576	call getin("u",gru)
577	write(,) " u = ",gru
578	grv=0. !default value for grv
579	PRINT *,'meridional northward component of the geostrophic',
580	&' wind (m/s) ?'
581	call getin("v",grv)
582	write(,) " v = ",grv
583
584	c Initialize winds for first time step
585	DO ilayer=1,nlayer
586	u(ilayer)=gru
587	v(ilayer)=grv
588	w(ilayer)=0 ! default: no vertical wind
589	ENDDO
590
591	c Initialize turbulente kinetic energy
592	DO ilevel=1,nlevel
593	q2(ilevel)=0.E+0
594	ENDDO
595
596	c CO2 ice on the surface
597	c -------------------
598	co2ice(1)=0.E+0 ! default value for co2ice
599	PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
600	call getin("co2ice",co2ice)
601	write(,) " co2ice = ",co2ice
602
603	c
604	c emissivity
605	c ----------
606	emis=emissiv
607	IF (co2ice(1).eq.1.E+0) THEN
608	emis=emisice(1) ! northern hemisphere
609	IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
610	ENDIF
611
612
613
614	c Compute pressures and altitudes of atmospheric levels
615	c ----------------------------------------------------------------
616
617	c Vertical Coordinates
618	c """"""""""""""""""""
619	hybrid=.true.
620	PRINT *,'Hybrid coordinates ?'
621	call getin("hybrid",hybrid)
622	write(,) " hybrid = ", hybrid
623
624	CALL disvert
625	! now that disvert has been called, initialize module vertical_layers_mod
626	call init_vertical_layers(nlayer,preff,scaleheight,
627	& ap,bp,aps,bps,presnivs,pseudoalt)
628
629	DO ilevel=1,nlevel
630	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
631	ENDDO
632
633	DO ilayer=1,nlayer
634	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
635	ENDDO
636
637	DO ilayer=1,nlayer
638	zlay(ilayer)=-200.E+0 rlog(play(ilayer)/plev(1))
639	& /g
640	ENDDO
641
642
643	c Initialize temperature profile
644	c --------------------------------------
645	pks=psurf**rcp
646
647	c altitude in km in profile: divide zlay by 1000
648	tmp1(0)=0.E+0
649	DO ilayer=1,nlayer
650	tmp1(ilayer)=zlay(ilayer)/1000.E+0
651	ENDDO
652
653	call profile(nlayer+1,tmp1,tmp2)
654
655	tsurf=tmp2(0)
656	DO ilayer=1,nlayer
657	temp(ilayer)=tmp2(ilayer)
658	ENDDO
659
660
661
662	! Initialize soil properties and temperature
663	! ------------------------------------------
664	volcapa=1.e6 ! volumetric heat capacity
665	DO isoil=1,nsoil
666	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
667	tsoil(isoil)=tsurf(1) ! soil temperature
668	ENDDO
669
670	! Initialize depths
671	! -----------------
672	do isoil=0,nsoil-1
673	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
674	enddo
675	do isoil=1,nsoil
676	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
677	enddo
678
679	c Initialize traceurs
680	c ---------------------------
681
682	if (photochem.or.callthermos) then
683	write(,) 'Initializing chemical species'
684	! thermo=0: initialize over all atmospheric layers
685	thermo=0
686	call inichim_newstart(q,psurf,sig,nq,latitude,longitude,
687	$ cell_area,thermo,qsurf)
688	endif
689
690	c Check if the surface is a water ice reservoir
691	c --------------------------------------------------
692	watercaptag(1)=.false. ! Default: no water ice reservoir
693	print *,'Water ice cap on ground ?'
694	call getin("watercaptag",watercaptag)
695	write(,) " watercaptag = ",watercaptag
696
697
698	c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
699	c ----------------------------------------------------------------
700	c (output done in "writeg1d", typically called by "physiq.F")
701
702	g1d_nlayer=nlayer
703	g1d_nomfich='g1d.dat'
704	g1d_unitfich=40
705	g1d_nomctl='g1d.ctl'
706	g1d_unitctl=41
707	g1d_premier=.true.
708	g2d_premier=.true.
709
710	c Write a "startfi" file
711	c --------------------
712	c This file will be read during the first call to "physiq".
713	c It is needed to transfert physics variables to "physiq"...
714
715	call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
716	& nq,dtphys,float(day0),time,cell_area,
717	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
718	call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
719	. dtphys,time,
720	. tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
721	. totcloudfrac) !MV17 nuages alizee
722
723	c=======================================================================
724	c 1D MODEL TIME STEPPING LOOP
725	c=======================================================================
726	c
727	firstcall=.true.
728	lastcall=.false.
729
730	DO idt=1,ndt
731	c IF (idt.eq.ndt) lastcall=.true.
732	IF (idt.eq.ndt-day_step-1) then !test
733	lastcall=.true.
734	call solarlong(day*1.0,zls)
735	write(103,) 'Ls=',zls180./pi
736	write(103,) 'Lat=', latitude(1)180./pi
737	write(103,) 'Tau=', tauvis/odprefpsurf
738	write(103,*) 'RunEnd - Atmos. Temp. File'
739	write(103,*) 'RunEnd - Atmos. Temp. File'
740	write(104,) 'Ls=',zls180./pi
741	write(104,*) 'Lat=', latitude(1)
742	write(104,) 'Tau=', tauvis/odprefpsurf
743	write(104,*) 'RunEnd - Atmos. Temp. File'
744	ENDIF
745
746	c compute geopotential
747	c ~~~~~~~~~~~~~~~~~~~~~
748	DO ilayer=1,nlayer
749	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
750	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
751	ENDDO
752	phi(1)=pksh(1)(1.E+0-s(1))
753	DO ilayer=2,nlayer
754	phi(ilayer)=phi(ilayer-1)+
755	& pks(h(ilayer-1)+h(ilayer)).5E+0
756	& *(s(ilayer-1)-s(ilayer))
757
758	ENDDO
759
760	c call physics
761	c --------------------
762	! write(,) "testphys1d avant q", q(1,:)
763	CALL physiq (1,llm,nq,
764	, firstcall,lastcall,
765	, day,time,dtphys,
766	, plev,play,phi,
767	, u, v,temp, q,
768	, w,
769	C - outputs
770	s du, dv, dtemp, dq,dpsurf)
771	! write(,) "testphys1d apres q", q(1,:)
772
773
774	c wind increment : specific for 1D
775	c --------------------------------
776
777	c The physics compute the tendencies on u and v,
778	c here we just add Coriolos effect
779	c
780	c DO ilayer=1,nlayer
781	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
782	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
783	c ENDDO
784
785	c For some tests : No coriolis force at equator
786	c if(latitude(1).eq.0.) then
787	DO ilayer=1,nlayer
788	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
789	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
790	ENDDO
791	c end if
792	c
793	c
794	c Compute time for next time step
795	c ---------------------------------------
796	firstcall=.false.
797	time=time+dtphys/daysec
798	IF (time.gt.1.E+0) then
799	time=time-1.E+0
800	day=day+1
801	ENDIF
802
803	c compute winds and temperature for next time step
804	c ----------------------------------------------------------
805
806	DO ilayer=1,nlayer
807	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
808	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
809	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
810	ENDDO
811
812	c compute pressure for next time step
813	c ----------------------------------------------------------
814
815	psurf=psurf+dtphys*dpsurf(1) ! surface pressure change
816	DO ilevel=1,nlevel
817	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
818	ENDDO
819	DO ilayer=1,nlayer
820	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
821	ENDDO
822
823	! increment tracers
824	DO iq = 1, nq
825	DO ilayer=1,nlayer
826	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
827	ENDDO
828	ENDDO
829
830	ENDDO ! of idt=1,ndt ! end of time stepping loop
831
832	c ========================================================
833	c OUTPUTS
834	c ========================================================
835
836	c finalize and close grads files "g1d.dat" and "g1d.ctl"
837
838	c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
839	CALL endg1d(1,nlayer,zlay/1000.,ndt)
840
841	write(,) "testphys1d: Everything is cool."
842
843	END
844
845	c***********************************************************************
846	c***********************************************************************
847	c Dummy subroutines used only in 3D, but required to
848	c compile testphys1d (to cleanly use writediagfi)
849
850	subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
851
852	IMPLICIT NONE
853
854	INTEGER im,jm,ngrid,nfield
855	REAL pdyn(im,jm,nfield)
856	REAL pfi(ngrid,nfield)
857
858	if (ngrid.ne.1) then
859	write(,) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
860	stop
861	endif
862
863	pdyn(1,1,1:nfield)=pfi(1,1:nfield)
864
865	end
866
867	c***********************************************************************
868	c***********************************************************************
869

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

Download in other formats:

Original Format