Context Navigation

testphys1d.F @ 1242

Last change on this file since 1242 was 1233, checked in by aslmd, 11 years ago
LMDZ.MARS. Filling geom arrays is now out of phys_var_state_init. Done through a merged function ini_fillgeom within the comgeomfi_h module. Cosmetic changes. New interface with the mesoscale model: lesser amount of dirty MESOSCALE includes.
File size: 27.5 KB

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

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

Download in other formats:

Original Format