Context Navigation

rcm1d.F @ 3580

Last change on this file since 3580 was 3576, checked in by jbclement, 6 months ago

COMMON:
Follow-up of r3574:

Small corrections to make it work with Git;
Addition of the functionality with the programs 'Generic newstart' and 'Generic start2archive';
Improvements of the visualization format;
Diplaying version control information of every sub-folder in the "trunk" instead of only the "trunk" (who can do more can do less).

JBC

Property svn:executable set to *

File size: 44.0 KB

Rev	Line
[253]	1	program rcm1d
	2
	3	! to use 'getin'
[1216]	4	use ioipsl_getincom, only: getin
[1543]	5	use dimphy, only : init_dimphy
	6	use mod_grid_phy_lmdz, only : regular_lonlat
[1216]	7	use infotrac, only: nqtot, tname
[2706]	8	use tracer_h, only: noms, is_condensable
[3335]	9	use radinc_h, only : L_NSPECTV
[2784]	10	use surfdat_h, only: albedodat, phisfi, dryness,
[1216]	11	& zmea, zstd, zsig, zgam, zthe,
[1494]	12	& emissiv, emisice, iceradius,
[1216]	13	& dtemisice
	14	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
	15	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
	16	use phyredem, only: physdem0,physdem1
[1543]	17	use geometry_mod, only: init_geometry
[1318]	18	use planete_mod, only: apoastr,periastr,year_day,peri_day,
[3515]	19	& obliquit,nres,z0,coefvis,coefir,
[1318]	20	& timeperi,e_elips,p_elips
[1524]	21	use comcstfi_mod, only: pi, cpp, rad, g, r,
[1384]	22	& mugaz, rcp, omeg
[1525]	23	use time_phylmdz_mod, only: daysec, dtphys, day_step, ecritphy,
[1531]	24	& nday, iphysiq
[2635]	25	use callkeys_mod, only: tracer,rings_shadow,
[1802]	26	& specOLR,water,pceil,ok_slab_ocean,photochem
[1621]	27	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff, sig,
	28	& presnivs,pseudoalt,scaleheight
	29	USE vertical_layers_mod, ONLY: init_vertical_layers
[2354]	30	USE logic_mod, ONLY: hybrid
[2972]	31	use radinc_h, only: naerkind
[1531]	32	use regular_lonlat_mod, only: init_regular_lonlat
	33	use planete_mod, only: ini_planete_mod
[1543]	34	use physics_distribution_mod, only: init_physics_distribution
	35	use regular_lonlat_mod, only: init_regular_lonlat
	36	use mod_interface_dyn_phys, only: init_interface_dyn_phys
[1524]	37	use inifis_mod, only: inifis
[1883]	38	use phys_state_var_mod, only: phys_state_var_init
[1549]	39	use physiq_mod, only: physiq
[3562]	40	use restart1D_mod, only: writerestart1D
	41	use iostart, only: length
	42	use netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE,
	43	& nf90_strerror,NF90_INQ_VARID, NF90_GET_VAR,
	44	& NF90_CLOSE
[3574]	45	use version_info_mod, only: print_version_info
	46
[253]	47	implicit none
	48
	49	!==================================================================
	50	!
	51	! Purpose
	52	! -------
	53	! Run the physics package of the universal model in a 1D column.
	54	!
	55	! It can be compiled with a command like (e.g. for 25 layers):
	56	! "makegcm -p std -d 25 rcm1d"
	57	! It requires the files "callphys.def", "z2sig.def",
	58	! "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
	59	!
	60	! Authors
	61	! -------
	62	! Frederic Hourdin
	63	! R. Fournier
	64	! F. Forget
	65	! F. Montmessin (water ice added)
	66	! R. Wordsworth
[589]	67	! B. Charnay
	68	! A. Spiga
[728]	69	! J. Leconte (2012)
[253]	70	!
	71	!==================================================================
	72
	73	#include "dimensions.h"
[534]	74	#include "paramet.h"
[1308]	75	!include "dimphys.h"
[253]	76	#include "netcdf.inc"
[534]	77	#include "comgeom.h"
[253]	78
	79	c --------------------------------------------------------------
	80	c Declarations
	81	c --------------------------------------------------------------
	82	c
	83	INTEGER unitstart ! unite d'ecriture de "startfi"
	84	INTEGER nlayer,nlevel,nsoil,ndt
	85	INTEGER ilayer,ilevel,isoil,idt,iq
	86	LOGICAl firstcall,lastcall
[3562]	87	LOGICAL restart
	88	INTEGER nid_restart1D, nid_restartfi
	89	INTEGER controleid, uid, vid, tempid, tsurfid, did, tid
	90	INTEGER,ALLOCATABLE :: qid(:), qsurfid(:)
[253]	91	c
	92	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
[1539]	93	REAL day ! date durant le run
[253]	94	REAL time ! time (0<time<1 ; time=0.5 a midi)
[1539]	95	REAL play(llm) ! Pressure at the middle of the layers (Pa)
	96	REAL plev(llm+1) ! intermediate pressure levels (pa)
[1216]	97	REAL psurf,tsurf(1)
[1539]	98	REAL u(llm),v(llm) ! zonal, meridional wind
	99	REAL gru,grv ! prescribed "geostrophic" background wind
	100	REAL temp(llm) ! temperature at the middle of the layers
	101	REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
	102	REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
	103	REAL,ALLOCATABLE :: tsoil(:) ! subsurface soil temperature (K)
	104	! REAL co2ice ! co2ice layer (kg.m-2) !not used anymore
	105	integer :: i_co2_ice=0 ! tracer index of co2 ice
	106	integer :: i_h2o_ice=0 ! tracer index of h2o ice
	107	integer :: i_h2o_vap=0 ! tracer index of h2o vapor
[3335]	108	REAL emis(1) ! emissivity of surface
	109	real :: albedo(1,L_NSPECTV) ! surface albedo in various spectral bands
[1539]	110	REAL q2(llm+1) ! Turbulent Kinetic Energy
	111	REAL zlay(llm) ! altitude estimee dans les couches (km)
[253]	112
	113	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
[1308]	114	REAL du(llm),dv(llm),dtemp(llm)
	115	REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
[1549]	116	REAL dpsurf(1)
[1216]	117	REAL,ALLOCATABLE :: dq(:,:)
	118	REAL,ALLOCATABLE :: dqdyn(:,:)
[253]	119
	120	c Various intermediate variables
	121	! INTEGER thermo
	122	REAL zls
[1308]	123	REAL phi(llm),h(llm),s(llm)
	124	REAL pks, ptif, w(llm)
[253]	125	INTEGER ierr, aslun
[1308]	126	REAL tmp1(0:llm),tmp2(0:llm)
[787]	127	integer :: nq !=1 ! number of tracers
[253]	128
	129	character*2 str2
	130	character (len=7) :: str7
[1533]	131	character(len=44) :: txt
[2713]	132	character(len=44) :: tracer_profile_file_name
[253]	133
	134	logical oldcompare, earthhack,saveprofile
[3562]	135	real controle1D(length)
[253]	136
	137	! added by RW for zlay computation
	138	real Hscale, Hmax, rho, dz
	139
	140	! added by RW for autozlevs computation
[2354]	141	logical autozlevs
[253]	142	real nu, xx, pMIN, zlev, Htop
[1308]	143	real logplevs(llm)
[253]	144
	145	! added by BC
[1308]	146	REAL cloudfrac(1,llm)
[787]	147	REAL hice(1),totcloudfrac(1)
[253]	148
[1303]	149	! added by BC for ocean
	150	real rnat(1)
[3105]	151	REAL tslab(1,2),tsea_ice(1),sea_ice(1)
[3441]	152	real tice(1)
[1303]	153	real pctsrf_sic(1)
	154
[2621]	155	! added by BC to accelerate convergence
	156	logical accelerate_temperature
	157	real factor_acceleration
[1303]	158
[787]	159	! added by AS to avoid the use of adv trac common
[1984]	160	character*30,allocatable :: nametmp(:) !
[787]	161
[1543]	162	real :: latitude(1), longitude(1), cell_area(1)
[787]	163
[2436]	164	! added by JVO and YJ to read modern traceur.def
	165	character(len=500) :: line ! to store a line of text
	166	LOGICAL :: moderntracdef=.false. ! JVO, YJ : modern traceur.def
	167
[253]	168	c=======================================================================
	169	c INITIALISATION
	170	c=======================================================================
[3574]	171	if (command_argument_count() > 0) then ! Get the number of command-line arguments
[3576]	172	call get_command_argument(1,txt) ! Read the argument given to the program
	173	select case (trim(adjustl(txt)))
[3574]	174	case('version')
	175	call print_version_info()
	176	stop
	177	case default
	178	error stop 'The argument given to the program is '
	179	&//'unknown!'
	180	end select
	181	endif
	182
[2972]	183	! check if 'rcm1d.def' file is around
	184	open(90,file='rcm1d.def',status='old',form='formatted',
	185	& iostat=ierr)
	186	if (ierr.ne.0) then
	187	write(,) 'Cannot find required file "rcm1d.def"'
	188	write(,) 'which should contain some input parameters'
	189	write(,) ' ... might as well stop here ...'
	190	stop
	191	else
	192	close(90)
	193	endif
[253]	194
[2972]	195	! now, run.def is needed anyway. so we create a dummy temporary one
	196	! for ioipsl to work. if a run.def is already here, stop the
	197	! program and ask the user to do a bit of cleaning
	198	open(90,file='run.def',status='old',form='formatted',
	199	& iostat=ierr)
	200	if (ierr.eq.0) then
	201	close(90)
	202	write(,) 'There is already a run.def file.'
	203	write(,) 'This is not compatible with 1D runs.'
	204	write(,) 'Please remove the file and restart the run.'
	205	write(,) 'Runtime parameters are supposed to be in rcm1d.def'
	206	stop
	207	else
	208	call system('touch run.def')
	209	call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
	210	call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
	211	endif
	212
[3562]	213	! Check restart
	214	call getin("restart",restart)
	215	if (restart) then
	216	ierr = NF90_OPEN('start1D.nc', NF90_NOWRITE, nid_restart1D)
	217	if (ierr.NE.NF90_NOERR) then
	218	write(,)'readstart1D: problem opening 1D start file'
	219	write(,)trim(nf90_strerror(ierr))
	220	call abort_physic("readstart1D","Cannot open file",1)
	221	endif
	222	ierr = NF90_OPEN('startfi.nc', NF90_NOWRITE, nid_restartfi)
	223	if (ierr.NE.NF90_NOERR) then
	224	write(,)'readstart1D: problem opening startfi file'
	225	write(,)trim(nf90_strerror(ierr))
	226	call abort_physic("readstart1D","Cannot open file",1)
	227	endif
	228	endif ! restart
	229
[1883]	230	! read nq from traceur.def
	231	open(90,file='traceur.def',status='old',form='formatted',
	232	& iostat=ierr)
	233	if (ierr.eq.0) then
[2436]	234	! read number of tracers:
	235	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
	236	READ(90,'(A)') line
	237	IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
	238	READ(line,*,iostat=ierr) nq ! Try standard traceur.def
	239	ELSE
	240	moderntracdef = .true.
	241	DO
	242	READ(90,'(A)',iostat=ierr) line
	243	IF (ierr.eq.0) THEN
	244	IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
	245	READ(line,*,iostat=ierr) nq
	246	EXIT
	247	ENDIF
	248	ELSE ! If pb, or if reached EOF without having found nbtr
	249	write(,) "rcm1d: error reading number of tracers"
	250	write(,) " (first line of traceur.def) "
	251	stop
	252	ENDIF
	253	ENDDO
	254	ENDIF ! if modern or standard traceur.def
[1883]	255	else
	256	nq=0
	257	endif
	258	close(90)
	259
	260	! Initialize dimphy module
	261	call init_dimphy(1,llm)
[2972]	262
[1883]	263	! now initialize arrays using phys_state_var_init
[2972]	264	! but first initialise naerkind (from callphys.def)
	265	naerkind=0 !default
	266	call getin("naerkind",naerkind)
	267
[1883]	268	call phys_state_var_init(nq)
	269
[526]	270	saveprofile=.false.
[601]	271	saveprofile=.true.
[253]	272
[589]	273	c ----------------------------------------
	274	c Default values (corresponding to Mars)
	275	c ----------------------------------------
[253]	276
	277	pi=2.E+0*asin(1.E+0)
	278
	279	c Parametres Couche limite et Turbulence
	280	c --------------------------------------
	281	z0 = 1.e-2 ! surface roughness (m) ~0.01
	282
	283	c propriete optiques des calottes et emissivite du sol
	284	c ----------------------------------------------------
	285	emissiv= 0.95 ! Emissivite du sol martien ~.95
	286	emisice(1)=0.95 ! Emissivite calotte nord
	287	emisice(2)=0.95 ! Emissivite calotte sud
	288
	289	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
	290	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
	291	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
	292	dtemisice(2) = 2. ! time scale for snow metamorphism (south
	293	hybrid=.false.
	294
	295	c ------------------------------------------------------
[589]	296	c Load parameters from "run.def" and "gases.def"
[253]	297	c ------------------------------------------------------
	298
	299
	300	! check if we are going to run with or without tracers
	301	write(,) "Run with or without tracer transport ?"
	302	tracer=.false. ! default value
	303	call getin("tracer",tracer)
	304	write(,) " tracer = ",tracer
	305
[589]	306	! OK. now that run.def has been read once -- any variable is in memory.
	307	! so we can dump the dummy run.def
[1536]	308	! call system("rm -rf run.def") ! Ehouarn: delay this to after inifis
[589]	309
[253]	310	! while we're at it, check if there is a 'traceur.def' file
	311	! and preocess it, if necessary. Otherwise initialize tracer names
	312	if (tracer) then
	313	! load tracer names from file 'traceur.def'
	314	open(90,file='traceur.def',status='old',form='formatted',
	315	& iostat=ierr)
	316	if (ierr.eq.0) then
	317	write(,) "rcm1d: Reading file traceur.def"
	318	! read number of tracers:
[2436]	319	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
	320	READ(90,'(A)') line
	321	IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
	322	READ(line,*,iostat=ierr) nq ! Try standard traceur.def
	323	ELSE
	324	moderntracdef = .true.
	325	DO
	326	READ(90,'(A)',iostat=ierr) line
	327	IF (ierr.eq.0) THEN
	328	IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
	329	READ(line,*,iostat=ierr) nq
	330	EXIT
	331	ENDIF
	332	ELSE ! If pb, or if reached EOF without having found nbtr
	333	write(,) "rcm1d: error reading number of tracers"
	334	write(,) " (first line of traceur.def) "
	335	stop
	336	ENDIF
	337	ENDDO
	338	ENDIF ! if modern or standard traceur.def
[1216]	339	nqtot=nq ! set value of nqtot (in infotrac module) as nq
[253]	340	if (ierr.ne.0) then
	341	write(,) "rcm1d: error reading number of tracers"
	342	write(,) " (first line of traceur.def) "
	343	stop
[1216]	344	endif
	345	if (nq>0) then
	346	allocate(tname(nq))
[1802]	347	allocate(noms(nq))
[1308]	348	allocate(q(llm,nq))
[1216]	349	allocate(qsurf(nq))
[1308]	350	allocate(dq(llm,nq))
	351	allocate(dqdyn(llm,nq))
[253]	352	else
[1267]	353	write(,) "rcm1d: Error. You set tracer=.true."
	354	write(,) " but # of tracers in traceur.def is ",nq
	355	stop
[253]	356	endif
	357
	358	do iq=1,nq
	359	! minimal version, just read in the tracer names, 1 per line
[1216]	360	read(90,*,iostat=ierr) tname(iq)
[1802]	361	noms(iq)=tname(iq)
[253]	362	if (ierr.ne.0) then
	363	write(,) 'rcm1d: error reading tracer names...'
	364	stop
	365	endif
	366	enddo !of do iq=1,nq
	367	! check for co2_ice / h2o_ice tracers:
	368	i_co2_ice=0
	369	i_h2o_ice=0
[650]	370	i_h2o_vap=0
[253]	371	do iq=1,nq
[1216]	372	if (tname(iq)=="co2_ice") then
[253]	373	i_co2_ice=iq
[1216]	374	elseif (tname(iq)=="h2o_ice") then
[253]	375	i_h2o_ice=iq
[1216]	376	elseif (tname(iq)=="h2o_vap") then
[650]	377	i_h2o_vap=iq
[253]	378	endif
	379	enddo
	380	else
	381	write(,) 'Cannot find required file "traceur.def"'
	382	write(,) ' If you want to run with tracers, I need it'
	383	write(,) ' ... might as well stop here ...'
	384	stop
	385	endif
	386	close(90)
	387
[1802]	388
[1267]	389	else ! of if (tracer)
	390	nqtot=0
	391	nq=0
	392	! still, make allocations for 1 dummy tracer
	393	allocate(tname(1))
	394	allocate(qsurf(1))
[1308]	395	allocate(q(llm,1))
	396	allocate(dq(llm,1))
[1267]	397
[970]	398	! Check that tracer boolean is compliant with number of tracers
	399	! -- otherwise there is an error (and more generally we have to be consistent)
[1267]	400	if (nq .ge. 1) then
[970]	401	write(,) "------------------------------"
	402	write(,) "rcm1d: You set tracer=.false."
[1267]	403	write(,) " But set number of tracers to ",nq
[970]	404	write(,) " > If you want tracers, set tracer=.true."
	405	write(,) "------------------------------"
	406	stop
	407	endif
[253]	408	endif ! of if (tracer)
	409
[589]	410	!!! GEOGRAPHICAL INITIALIZATIONS
	411	!!! AREA. useless in 1D
[1542]	412	cell_area(1)=1.E+0 !JL+EM to have access to the area in the diagfi.nc files.
[650]	413	!!! GEOPOTENTIAL. useless in 1D because control by surface pressure
[589]	414	phisfi(1)=0.E+0
	415	!!! LATITUDE. can be set.
[1216]	416	latitude=0 ! default value for latitude
[589]	417	PRINT *,'latitude (in degrees) ?'
[1216]	418	call getin("latitude",latitude)
	419	write(,) " latitude = ",latitude
	420	latitude=latitude*pi/180.E+0
[589]	421	!!! LONGITUDE. useless in 1D.
[1216]	422	longitude=0.E+0
	423	longitude=longitude*pi/180.E+0
[589]	424
	425
	426	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	427	!!!! PLANETARY CONSTANTS !!!!
	428	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	429
	430	g = -99999.
	431	PRINT *,'GRAVITY in m s-2 ?'
	432	call getin("g",g)
	433	IF (g.eq.-99999.) THEN
[1275]	434	PRINT *,"STOP. I NEED g IN RCM1D.DEF."
[589]	435	STOP
	436	ELSE
	437	PRINT *,"--> g = ",g
	438	ENDIF
	439
[1275]	440	rad = -99999.
	441	PRINT *,'PLANETARY RADIUS in m ?'
	442	call getin("rad",rad)
	443	! Planetary radius is needed to compute shadow of the rings
[1371]	444	IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN
[1275]	445	PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
	446	STOP
	447	ELSE
	448	PRINT *,"--> rad = ",rad
	449	ENDIF
[589]	450
	451	daysec = -99999.
	452	PRINT *,'LENGTH OF A DAY in s ?'
	453	call getin("daysec",daysec)
	454	IF (daysec.eq.-99999.) THEN
[1275]	455	PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."
[589]	456	STOP
	457	ELSE
	458	PRINT *,"--> daysec = ",daysec
	459	ENDIF
	460	omeg=4.*asin(1.)/(daysec)
	461	PRINT *,"OK. FROM THIS I WORKED OUT:"
	462	PRINT *,"--> omeg = ",omeg
	463
	464	year_day = -99999.
	465	PRINT *,'LENGTH OF A YEAR in days ?'
	466	call getin("year_day",year_day)
	467	IF (year_day.eq.-99999.) THEN
[1275]	468	PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
[589]	469	STOP
	470	ELSE
	471	PRINT *,"--> year_day = ",year_day
	472	ENDIF
	473
	474	periastr = -99999.
	475	PRINT *,'MIN DIST STAR-PLANET in AU ?'
	476	call getin("periastr",periastr)
	477	IF (periastr.eq.-99999.) THEN
[1275]	478	PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."
[589]	479	STOP
	480	ELSE
	481	PRINT *,"--> periastr = ",periastr
	482	ENDIF
	483
	484	apoastr = -99999.
[591]	485	PRINT *,'MAX DIST STAR-PLANET in AU ?'
[589]	486	call getin("apoastr",apoastr)
	487	IF (apoastr.eq.-99999.) THEN
[1275]	488	PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."
[589]	489	STOP
	490	ELSE
	491	PRINT *,"--> apoastr = ",apoastr
	492	ENDIF
	493
	494	peri_day = -99999.
	495	PRINT *,'DATE OF PERIASTRON in days ?'
	496	call getin("peri_day",peri_day)
	497	IF (peri_day.eq.-99999.) THEN
[1275]	498	PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."
[589]	499	STOP
	500	ELSE IF (peri_day.gt.year_day) THEN
	501	PRINT *,"STOP. peri_day.gt.year_day"
	502	STOP
	503	ELSE
	504	PRINT *,"--> peri_day = ", peri_day
	505	ENDIF
	506
	507	obliquit = -99999.
	508	PRINT *,'OBLIQUITY in deg ?'
	509	call getin("obliquit",obliquit)
	510	IF (obliquit.eq.-99999.) THEN
[1275]	511	PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."
[589]	512	STOP
	513	ELSE
	514	PRINT *,"--> obliquit = ",obliquit
	515	ENDIF
	516
[3562]	517	if (restart) then
	518	ierr=NF90_INQ_VARID(nid_restart1D,'controle',controleid)
	519	if (ierr==NF90_NOERR) then
	520	ierr=NF90_GET_VAR(nid_restart1D,controleid,controle1D)
	521	if (ierr/=NF90_NOERR) then
	522	PRINT*, 'read restart 1D: Failed loading psurf'
	523	CALL abort_physic("readstart1D","Failed to read variable",1)
	524	endif
	525	endif
	526	psurf = controle1D(1)
	527	PRINT *,"In restart1D.nc, day0=",controle1D(2),
	528	& " and time=",controle1D(3)
	529	PRINT*,"(These values are NOT used, ",
	530	& "those given in rcm1d.def are used)"
	531	else ! restart
	532	psurf = -99999.
	533	PRINT *,'SURFACE PRESSURE in Pa ?'
	534	call getin("psurf",psurf)
	535	IF (psurf.eq.-99999.) THEN
	536	PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
	537	STOP
	538	ELSE
	539	PRINT *,"--> psurf = ",psurf
	540	ENDIF
	541	endif ! restart
[589]	542	!! we need reference pressures
	543	pa=psurf/30.
	544	preff=psurf ! these values are not needed in 1D (are you sure JL12)
	545
	546	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	547	!!!! END PLANETARY CONSTANTS !!!!
	548	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	549
[253]	550	c Date et heure locale du debut du run
	551	c ------------------------------------
	552	c Date (en sols depuis le solstice de printemps) du debut du run
[3562]	553	!if (restart) then
	554	! ! day
	555	! ierr=NF90_INQ_VARID(nid_restart1D,'day',did)
	556	! if (ierr==NF90_NOERR) then
	557	! ierr=NF90_GET_VAR(nid_restart1D,did,day)
	558	! if (ierr/=NF90_NOERR) then
	559	! PRINT*, 'read restart 1D: Failed loading day'
	560	! CALL abort_physic("readstart1D","Failed to read variable",1)
	561	! endif
	562	! endif
	563
	564	! ! time
	565	! ierr=NF90_INQ_VARID(nid_restart1D,'time',tid)
	566	! if (ierr==NF90_NOERR) then
	567	! ierr=NF90_GET_VAR(nid_restart1D,tid,time)
	568	! if (ierr/=NF90_NOERR) then
	569	! PRINT*, 'read restart 1D: Failed loading time'
	570	! CALL abort_physic("readstart1D","Failed to read variable",1)
	571	! endif
	572	! endif
	573
	574	!else
	575	day0 = 0 ! default value for day0
	576	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
	577	call getin("day0",day0)
	578	day=float(day0)
	579	write(,) " day0 = ",day0
[253]	580	c Heure de demarrage
[3562]	581	time=0 ! default value for time
	582	write(,)'Initial local time (in hours, between 0 and 24)?'
	583	call getin("time",time)
	584	write(,)" time = ",time
	585	time=time/24.E+0 ! convert time (hours) to fraction of sol
	586	!endif
[253]	587
	588	c Discretisation (Definition de la grille et des pas de temps)
	589	c --------------
	590	c
[1308]	591	nlayer=llm
[253]	592	nlevel=nlayer+1
	593
	594	day_step=48 ! default value for day_step
	595	PRINT *,'Number of time steps per sol ?'
	596	call getin("day_step",day_step)
	597	write(,) " day_step = ",day_step
	598
[1531]	599	iphysiq=1 ! in 1D model physics are called evry time step
[534]	600	ecritphy=day_step ! default value for ecritphy
	601	PRINT *,'Nunber of steps between writediagfi ?'
	602	call getin("ecritphy",ecritphy)
	603	write(,) " ecritphy = ",ecritphy
	604
[253]	605	ndt=10 ! default value for ndt
	606	PRINT *,'Number of sols to run ?'
	607	call getin("ndt",ndt)
	608	write(,) " ndt = ",ndt
[1525]	609	nday=ndt
[253]	610
	611	ndt=ndt*day_step
	612	dtphys=daysec/day_step
[589]	613	write(,)"-------------------------------------"
	614	write(,)"-------------------------------------"
	615	write(,)"--> Physical timestep is ",dtphys
	616	write(,)"-------------------------------------"
	617	write(,)"-------------------------------------"
[253]	618
[1531]	619	! initializations, as with iniphysiq.F90 for the 3D GCM
[1543]	620	call init_physics_distribution(regular_lonlat,4,
	621	& 1,1,1,nlayer,1)
	622	call init_interface_dyn_phys
[1531]	623	CALL init_regular_lonlat(1,1,longitude,latitude,
	624	& (/0.,0./),(/0.,0./))
[1543]	625	call init_geometry(1,longitude,latitude,
	626	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
	627	& cell_area)
[1621]	628	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
	629	! call init_vertical_layers(nlayer,preff,scaleheight,
	630	! & ap,bp,aps,bps,presnivs,pseudoalt)
[1883]	631	! call init_dimphy(1,nlayer) ! Initialize dimphy module
[1531]	632	call ini_planete_mod(nlayer,preff,ap,bp)
	633
[1524]	634	!!! CALL INIFIS
	635	!!! - read callphys.def
	636	!!! - calculate sine and cosine of longitude and latitude
	637	!!! - calculate mugaz and cp
	638	!!! NB: some operations are being done dummily in inifis in 1D
	639	!!! - physical constants: nevermind, things are done allright below
	640	!!! - physical frequency: nevermind, in inifis this is a simple print
[2076]	641	cpp=-9999. ! dummy init for inifis, will be rewrite later on
	642	r=-9999. ! dummy init for inifis, will be rewrite later on
[1525]	643	CALL inifis(1,llm,nq,day0,daysec,nday,dtphys,
[1542]	644	. latitude,longitude,cell_area,rad,g,r,cpp)
[1536]	645
[1539]	646	nsoil=nsoilmx
	647	allocate(tsoil(nsoilmx))
	648	!! those are defined in comsoil_h.F90
	649	IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx))
	650	IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1))
	651	IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx))
	652
[1536]	653	! At this point, both getin() and getin_p() functions have been used,
	654	! and the run.def file can be removed.
	655	call system("rm -rf run.def")
	656
[1524]	657	!!! We check everything is OK.
	658	PRINT *,"CHECK"
	659	PRINT *,"--> mugaz = ",mugaz
	660	PRINT *,"--> cpp = ",cpp
	661	r = 8.314511E+0 * 1000.E+0 / mugaz
	662	rcp = r / cpp
	663
[526]	664	c output spectrum?
	665	write(,)"Output spectral OLR?"
	666	specOLR=.false.
	667	call getin("specOLR",specOLR)
	668	write(,)" specOLR = ",specOLR
	669
[2618]	670	c option for temperature evolution?
	671	write(,)"accelerate temperature?"
	672	accelerate_temperature=.false.
	673	call getin("accelerate_temperature",accelerate_temperature)
	674	write(,)" accelerate_temperature = ",accelerate_temperature
	675
[2621]	676	write(,)"factor for temperature acceleration"
	677	factor_acceleration=1.0
	678	call getin("factor_acceleration",factor_acceleration)
	679	write(,)" factor_acceleration = ",factor_acceleration
[2618]	680
[2621]	681
	682
[253]	683	c
	684	c pour le schema d'ondes de gravite
	685	c ---------------------------------
	686	zmea(1)=0.E+0
	687	zstd(1)=0.E+0
	688	zsig(1)=0.E+0
	689	zgam(1)=0.E+0
	690	zthe(1)=0.E+0
	691
	692	c Initialisation des traceurs
	693	c ---------------------------
	694
[787]	695	DO iq=1,nq
[253]	696	DO ilayer=1,nlayer
	697	q(ilayer,iq) = 0.
	698	ENDDO
	699	ENDDO
	700
[787]	701	DO iq=1,nq
[253]	702	qsurf(iq) = 0.
	703	ENDDO
[1533]	704
	705	if (tracer) then ! Initialize tracers here.
	706
	707	write(,) "rcm1d : initializing tracers profiles"
[253]	708
[1533]	709	do iq=1,nq
[3562]	710
	711	if (restart) then
	712	if (iq.eq.1) then
	713	allocate(qid(nq))
	714	allocate(qsurfid(nq))
	715	endif
	716	! read q
	717	ierr=NF90_INQ_VARID(nid_restart1D,noms(iq),qid(iq))
	718	if (ierr==NF90_NOERR) then
	719	ierr=NF90_GET_VAR(nid_restart1D,qid(iq),q(:,iq),
	720	& count=(/1,nlayer,1,1/))
	721	if (ierr/=NF90_NOERR) then
	722	PRINT*, 'read restart 1D: Failed loading ', noms(iq)
	723	write(,)trim(nf90_strerror(ierr))
	724	CALL abort_physic("readstart1D",
	725	& "Failed to read variable",1)
	726	endif
	727	endif
	728	! read qsurf
	729	ierr=NF90_INQ_VARID(nid_restartfi,noms(iq),qsurfid(iq))
	730	if (ierr==NF90_NOERR) then
	731	ierr=NF90_GET_VAR(nid_restartfi,qsurfid(iq),qsurf(iq))
	732	if (ierr/=NF90_NOERR) then
	733	PRINT*, 'read restartfi: Failed loading ', noms(iq)
	734	write(,)trim(nf90_strerror(ierr))
	735	CALL abort_physic("readstartfi",
	736	& "Failed to read variable",1)
	737	endif
	738	endif
	739	else ! restart
[1533]	740
[3562]	741	txt=""
	742	write(txt,"(a)") tname(iq)
	743	write(,)" tracer:",trim(txt)
[1533]	744
[3562]	745	! CO2
	746	if (txt.eq."co2_ice") then
	747	q(:,iq)=0. ! kg/kg of atmosphere
	748	qsurf(iq)=0. ! kg/m2 at the surface
	749	! Look for a "profile_co2_ice" input file
	750	open(91,file='profile_co2_ice',status='old',
	751	& form='formatted',iostat=ierr)
	752	if (ierr.eq.0) then
	753	read(91,*) qsurf(iq)
	754	do ilayer=1,nlayer
	755	read(91,*) q(ilayer,iq)
	756	enddo
	757	else
	758	write(,) "No profile_co2_ice file!"
	759	endif
	760	close(91)
	761	endif ! of if (txt.eq."co2")
[1533]	762
[3562]	763	! WATER VAPOUR
	764	if (txt.eq."h2o_vap") then
	765	q(:,iq)=0. ! kg/kg of atmosphere
	766	qsurf(iq)=0. ! kg/m2 at the surface
	767	! Look for a "profile_h2o_vap" input file
	768	open(91,file='profile_h2o_vap',status='old',
	769	& form='formatted',iostat=ierr)
	770	if (ierr.eq.0) then
	771	read(91,*) qsurf(iq)
	772	do ilayer=1,nlayer
	773	read(91,*) q(ilayer,iq)
	774	enddo
	775	else
	776	write(,) "No profile_h2o_vap file!"
	777	endif
	778	close(91)
	779	endif ! of if (txt.eq."h2o_vap")
[1533]	780
[3562]	781	! WATER ICE
	782	if (txt.eq."h2o_ice") then
	783	q(:,iq)=0. ! kg/kg of atmosphere
	784	qsurf(iq)=0. ! kg/m2 at the surface
	785	! Look for a "profile_h2o_ice" input file
	786	open(91,file='profile_h2o_ice',status='old',
	787	& form='formatted',iostat=ierr)
	788	if (ierr.eq.0) then
	789	read(91,*) qsurf(iq)
	790	do ilayer=1,nlayer
	791	read(91,*) q(ilayer,iq)
	792	enddo
	793	else
	794	write(,) "No profile_h2o_ice file!"
	795	endif
	796	close(91)
	797	endif ! of if (txt.eq."h2o_ice")
[650]	798
[3562]	799	!_vap
	800	if((txt .ne. 'h2o_vap') .and.
	801	& (index(txt,'_vap' ) .ne. 0)) then
	802	q(:,iq)=0. !kg/kg of atmosphere
	803	qsurf(iq) = 0. !kg/kg of atmosphere
	804	! Look for a "profile_...._vap" input file
	805	tracer_profile_file_name=""
	806	tracer_profile_file_name='profile_'//txt
	807	open(91,file=tracer_profile_file_name,status='old',
	808	& form="formatted",iostat=ierr)
	809	if (ierr .eq. 0) then
	810	read(91,*)qsurf(iq)
	811	do ilayer=1,nlayer
	812	read(91,*)q(ilayer,iq)
	813	enddo
	814	else
	815	write(,) "No initial profile "
	816	write(,) " for this tracer :"
	817	write(,) txt
	818	endif
	819	close(91)
	820	endif ! (txt .eq. "_vap")
	821	!_ice
	822	if((txt.ne."h2o_ice") .and.
	823	& (index(txt,'_ice' ) /= 0)) then
	824	q(:,iq)=0. !kg/kg of atmosphere
	825	qsurf(iq) = 0. !kg/kg of atmosphere
	826	endif ! we only initialize the solid at 0
	827	endif ! restart
	828	enddo ! of do iq=1,nq
	829
[1533]	830	endif ! of tracer
	831
[253]	832	c Initialisation pour prendre en compte les vents en 1-D
	833	c ------------------------------------------------------
	834	ptif=2.E+0omegsinlat(1)
	835
	836
	837	c vent geostrophique
	838	gru=10. ! default value for gru
	839	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
	840	call getin("u",gru)
	841	write(,) " u = ",gru
	842	grv=0. !default value for grv
	843	PRINT *,'meridional northward component of the geostrophic',
	844	&' wind (m/s) ?'
	845	call getin("v",grv)
	846	write(,) " v = ",grv
	847
[1536]	848	! To be clean, also set vertical winds to zero
	849	w(1:nlayer)=0
	850
[3562]	851	if (restart) then
	852	ierr=NF90_INQ_VARID(nid_restart1D,'u',uid)
	853	if (ierr==NF90_NOERR) then
	854	ierr=NF90_GET_VAR(nid_restart1D,uid,u,
	855	& count=(/1,nlayer,1,1/))
	856	if (ierr/=NF90_NOERR) then
	857	PRINT*, 'read restart 1D: Failed loading u'
	858	write(,)trim(nf90_strerror(ierr))
	859	CALL abort_physic("readstart1D",
	860	& "Failed to read variable",1)
	861	endif
	862	endif
	863	ierr=NF90_INQ_VARID(nid_restart1D,'v',vid)
	864	if (ierr==NF90_NOERR) then
	865	ierr=NF90_GET_VAR(nid_restart1D,vid,v,
	866	& count=(/1,nlayer,1,1/))
	867	if (ierr/=NF90_NOERR) then
	868	PRINT*, 'read restart 1D: Failed loading v'
	869	write(,)trim(nf90_strerror(ierr))
	870	CALL abort_physic("readstart1D",
	871	& "Failed to read variable",1)
	872	endif
	873	endif
	874	else ! restart
[253]	875
[3562]	876	c Initialisation des vents au premier pas de temps
	877	DO ilayer=1,nlayer
	878	u(ilayer)=gru
	879	v(ilayer)=grv
	880	ENDDO
	881
	882	endif ! restart
	883
[253]	884	c energie cinetique turbulente
	885	DO ilevel=1,nlevel
	886	q2(ilevel)=0.E+0
	887	ENDDO
	888
	889	c emissivity / surface co2 ice ( + h2o ice??)
	890	c -------------------------------------------
[1216]	891	emis(1)=emissiv ! default value for emissivity
[253]	892	PRINT *,'Emissivity of bare ground ?'
[1216]	893	call getin("emis",emis(1))
	894	write(,) " emis = ",emis(1)
	895	emissiv=emis(1) ! we do this so that condense_co2 sets things to the right
[253]	896	! value if there is no snow
	897
	898	if(i_co2_ice.gt.0)then
	899	qsurf(i_co2_ice)=0 ! default value for co2ice
	900	print*,'Initial CO2 ice on the surface (kg.m-2)'
	901	call getin("co2ice",qsurf(i_co2_ice))
	902	write(,) " co2ice = ",qsurf(i_co2_ice)
	903	IF (qsurf(i_co2_ice).ge.1.E+0) THEN
	904	! if we have some CO2 ice on the surface, change emissivity
[1542]	905	if (latitude(1).ge.0) then ! northern hemisphere
[1216]	906	emis(1)=emisice(1)
[253]	907	else ! southern hemisphere
[1216]	908	emis(1)=emisice(2)
[253]	909	endif
	910	ENDIF
	911	endif
	912
	913	c calcul des pressions et altitudes en utilisant les niveaux sigma
	914	c ----------------------------------------------------------------
	915
	916	c Vertical Coordinates
	917	c """"""""""""""""""""
	918	hybrid=.true.
	919	PRINT *,'Hybrid coordinates ?'
	920	call getin("hybrid",hybrid)
	921	write(,) " hybrid = ", hybrid
	922
	923
	924	autozlevs=.false.
	925	PRINT *,'Auto-discretise vertical levels ?'
	926	call getin("autozlevs",autozlevs)
	927	write(,) " autozlevs = ", autozlevs
	928
[589]	929	pceil = psurf / 1000.0 ! Pascals
[253]	930	PRINT *,'Ceiling pressure (Pa) ?'
	931	call getin("pceil",pceil)
	932	write(,) " pceil = ", pceil
	933
	934	! Test of incompatibility:
	935	! if autozlevs used, cannot have hybrid too
	936
	937	if (autozlevs.and.hybrid) then
	938	print*,'Cannot use autozlevs and hybrid together!'
	939	call abort
	940	endif
	941
	942	if(autozlevs)then
	943
[305]	944	open(91,file="z2sig.def",form='formatted')
	945	read(91,*) Hscale
	946	DO ilayer=1,nlayer-2
	947	read(91,*) Hmax
	948	enddo
	949	close(91)
[253]	950
	951
	952	print*,'Hmax = ',Hmax,' km'
	953	print*,'Auto-shifting Hscale to:'
	954	! Hscale = Hmax / log(psurf/100.0)
	955	Hscale = Hmax / log(psurf/pceil)
	956	print*,'Hscale = ',Hscale,' km'
	957
	958	! none of this matters if we dont care about zlay
	959
	960	endif
	961
[2354]	962	call disvert_noterre
[1621]	963	! now that disvert has been called, initialize module vertical_layers_mod
	964	call init_vertical_layers(nlayer,preff,scaleheight,
	965	& ap,bp,aps,bps,presnivs,pseudoalt)
[253]	966
	967	if(.not.autozlevs)then
	968	! we want only the scale height from z2sig, in order to compute zlay
	969	open(91,file="z2sig.def",form='formatted')
	970	read(91,*) Hscale
	971	close(91)
	972	endif
	973
	974	! if(autozlevs)then
	975	! open(94,file="Hscale.temp",form='formatted')
	976	! read(94,*) Hscale
	977	! close(94)
	978	! endif
	979
	980	DO ilevel=1,nlevel
	981	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	982	ENDDO
	983
	984	DO ilayer=1,nlayer
	985	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	986	ENDDO
	987
	988
	989
	990	DO ilayer=1,nlayer
	991	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
	992	! & /g
	993	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
	994	ENDDO
	995
	996	! endif
	997
[3562]	998
[253]	999	c profil de temperature au premier appel
	1000	c --------------------------------------
	1001	pks=psurf**rcp
	1002
[3562]	1003	if (restart) then
	1004	! read temp
	1005	ierr=NF90_IN Q_VARID(nid_restart1D,'temp',tempid)
	1006	if (ierr==NF90_NOERR) then
	1007	ierr=NF90_GET_VAR(nid_restart1D,tempid,temp,
	1008	& count=(/1,nlayer,1,1/))
	1009	if (ierr/=NF90_NOERR) then
	1010	PRINT*, 'read restart 1D: Failed loading temp'
	1011	CALL abort_physic("readstart1D",
	1012	& "Failed to read variable",1)
	1013	endif
	1014	endif
	1015	! read tsurf
	1016	ierr=NF90_INQ_VARID(nid_restartfi,'tsurf',tsurfid)
	1017	if (ierr==NF90_NOERR) then
	1018	ierr=NF90_GET_VAR(nid_restartfi,tsurfid,tsurf)
	1019	if (ierr/=NF90_NOERR) then
	1020	PRINT*, 'read restartfi: Failed loading tsurf'
	1021	CALL abort_physic("readstartfi",
	1022	& "Failed to read variable",1)
	1023	endif
	1024	endif
	1025	else ! restart
	1026
[253]	1027	c altitude en km dans profile: on divise zlay par 1000
[3562]	1028	tmp1(0)=0.E+0
	1029	DO ilayer=1,nlayer
	1030	tmp1(ilayer)=zlay(ilayer)/1000.E+0
	1031	ENDDO
	1032	call profile(nlayer+1,tmp1,tmp2)
[253]	1033
[3562]	1034	tsurf(1)=tmp2(0)
	1035	DO ilayer=1,nlayer
	1036	temp(ilayer)=tmp2(ilayer)
	1037	ENDDO
	1038	print*,"check"
	1039	PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
	1040	PRINT*,"INPUT TEMPERATURE PROFILE",temp
[253]	1041
[3562]	1042	endif ! restart
	1043
[589]	1044	c Initialisation albedo / inertie du sol
	1045	c --------------------------------------
	1046	albedodat(1)=0.2 ! default value for albedodat
	1047	PRINT *,'Albedo of bare ground ?'
	1048	call getin("albedo",albedodat(1))
	1049	write(,) " albedo = ",albedodat(1)
[3335]	1050	! Initialize surface albedo to that of bare ground
	1051	albedo(1,:)=albedodat(1)
[253]	1052
[589]	1053	inertiedat(1,1)=400 ! default value for inertiedat
	1054	PRINT *,'Soil thermal inertia (SI) ?'
	1055	call getin("inertia",inertiedat(1,1))
	1056	write(,) " inertia = ",inertiedat(1,1)
	1057
[253]	1058	! Initialize soil properties and temperature
	1059	! ------------------------------------------
	1060	volcapa=1.e6 ! volumetric heat capacity
	1061	DO isoil=1,nsoil
	1062	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
[1216]	1063	tsoil(isoil)=tsurf(1) ! soil temperature
[253]	1064	ENDDO
	1065
	1066	! Initialize depths
	1067	! -----------------
	1068	do isoil=0,nsoil-1
[1802]	1069	mlayer(isoil)=3.e-2(2.*(isoil-0.5)) ! mid-layer depth
[253]	1070	enddo
	1071	do isoil=1,nsoil
[1802]	1072	layer(isoil)=3.e-2(2.*(isoil-1)) ! layer depth
[253]	1073	enddo
	1074
[1826]	1075	! Initialize cloud fraction and oceanic ice
	1076	! -----------------------------------------
	1077	hice=0.
	1078	do ilayer=1,llm
	1079	cloudfrac(1,ilayer)=0.
	1080	enddo
	1081	totcloudfrac=0.
[253]	1082
[1303]	1083	! Initialize slab ocean
	1084	! -----------------
	1085	rnat=1. ! default value for rnat
	1086	if(inertiedat(1,1).GE.10000.)then
	1087	rnat=0.
	1088	endif
	1089	if(ok_slab_ocean)then
	1090	rnat=0.
	1091	tslab(1,1)=tsurf(1)
	1092	tslab(1,2)=tsurf(1)
	1093	tsea_ice=tsurf
[3441]	1094	tice=0.
[1303]	1095	pctsrf_sic=0.
	1096	sea_ice=0.
	1097	endif
	1098
	1099
[1802]	1100	! Initialize chemical species
	1101	! -----------------
[2058]	1102	#ifndef MESOSCALE
[3562]	1103	if(tracer.and.photochem.and. .not.restart) then
	1104	print*, "Calling inichim_1D"
[2785]	1105	call initracer(1,nq)
[1802]	1106	allocate(nametmp(nq))
	1107	nametmp(1:nq)=tname(1:nq)
[2542]	1108	call inichim_1D(nlayer, nq, q, qsurf, play, 0, 0)
[1802]	1109	tname(1:nq)=nametmp(1:nq)
	1110	noms(1:nq)=nametmp(1:nq)
	1111	endif ! tracer and photochem
[2058]	1112	#endif
[1303]	1113
[3562]	1114	if (restart) then
	1115	ierr = NF90_CLOSE(nid_restart1D)
	1116	if (ierr/=NF90_NOERR) then
	1117	PRINT*, 'read restart1D: Failed closing restart1D.nc'
	1118	CALL abort_physic("readstart1D",
	1119	& "Failed closing file",1)
	1120	endif
	1121	ierr = NF90_CLOSE(nid_restartfi)
	1122	if (ierr/=NF90_NOERR) then
	1123	PRINT*, 'read restartfi: Failed closing restartfi.nc'
	1124	CALL abort_physic("readstartfi",
	1125	& "Failed closing file",1)
	1126	endif
	1127	endif
[1802]	1128
[1216]	1129	c Write a "startfi" file
[253]	1130	c --------------------
[1216]	1131	c This file will be read during the first call to "physiq".
	1132	c It is needed to transfert physics variables to "physiq"...
[253]	1133
[3562]	1134	if (.not. restart) then
	1135	call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
	1136	& dtphys,real(day0),time,cell_area,
	1137	& albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
	1138	call physdem1("startfi.nc",nsoilmx,1,llm,nq,
	1139	& dtphys,time,
	1140	& tsurf,tsoil,emis,albedo,q2,qsurf,
	1141	& cloudfrac,totcloudfrac,hice,
	1142	& rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice)
	1143	endif
[253]	1144	c=======================================================================
	1145	c BOUCLE TEMPORELLE DU MODELE 1D
	1146	c=======================================================================
	1147
	1148	firstcall=.true.
	1149	lastcall=.false.
	1150
	1151	DO idt=1,ndt
[486]	1152	IF (idt.eq.ndt) then !test
[253]	1153	lastcall=.true.
	1154	call stellarlong(day*1.0,zls)
[526]	1155	! write(103,) 'Ls=',zls180./pi
[1542]	1156	! write(103,) 'Lat=', latitude(1)180./pi
[526]	1157	! write(103,*) 'RunEnd - Atmos. Temp. File'
	1158	! write(103,*) 'RunEnd - Atmos. Temp. File'
	1159	! write(104,) 'Ls=',zls180./pi
[1542]	1160	! write(104,*) 'Lat=', latitude(1)
[526]	1161	! write(104,*) 'RunEnd - Atmos. Temp. File'
[253]	1162	ENDIF
	1163
	1164	c calcul du geopotentiel
	1165	c ~~~~~~~~~~~~~~~~~~~~~
	1166
	1167
[586]	1168	DO ilayer=1,nlayer
[253]	1169
	1170	! if(autozlevs)then
	1171	! s(ilayer)=(play(ilayer)/psurf)**rcp
	1172	! else
[586]	1173	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
[253]	1174	! endif
	1175	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
[586]	1176	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	1177	ENDDO
	1178
[253]	1179	! DO ilayer=1,nlayer
	1180	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	1181	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	1182	! ENDDO
	1183	phi(1)=pksh(1)(1.E+0-s(1))
	1184	DO ilayer=2,nlayer
	1185	phi(ilayer)=phi(ilayer-1)+
	1186	& pks(h(ilayer-1)+h(ilayer)).5E+0
	1187	& *(s(ilayer-1)-s(ilayer))
	1188
	1189	ENDDO
	1190
	1191	c appel de la physique
	1192	c --------------------
	1193
	1194
[787]	1195	CALL physiq (1,llm,nq,
[253]	1196	, firstcall,lastcall,
	1197	, day,time,dtphys,
	1198	, plev,play,phi,
	1199	, u, v,temp, q,
	1200	, w,
	1201	C - sorties
[1576]	1202	s du, dv, dtemp, dq,dpsurf)
[253]	1203
	1204
	1205	c evolution du vent : modele 1D
	1206	c -----------------------------
	1207
	1208	c la physique calcule les derivees temporelles de u et v.
	1209	c on y rajoute betement un effet Coriolis.
	1210	c
	1211	c DO ilayer=1,nlayer
	1212	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
	1213	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
	1214	c ENDDO
	1215
	1216	c Pour certain test : pas de coriolis a l'equateur
[1542]	1217	c if(latitude(1).eq.0.) then
[253]	1218	DO ilayer=1,nlayer
	1219	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
	1220	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
	1221	ENDDO
	1222	c end if
	1223	c
	1224	c
	1225	c Calcul du temps au pas de temps suivant
	1226	c ---------------------------------------
	1227	firstcall=.false.
	1228	time=time+dtphys/daysec
	1229	IF (time.gt.1.E+0) then
	1230	time=time-1.E+0
	1231	day=day+1
	1232	ENDIF
	1233
	1234	c calcul des vitesses et temperature au pas de temps suivant
	1235	c ----------------------------------------------------------
	1236
	1237	DO ilayer=1,nlayer
	1238	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
	1239	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
[2621]	1240	if(accelerate_temperature) then
	1241	temp(ilayer)=temp(ilayer)+dtphysdtemp(ilayer)
	1242	& max(1.0,play(ilayer)1e-5)*factor_acceleration
[2618]	1243	else
	1244	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
	1245	endif
[253]	1246	ENDDO
	1247
	1248	c calcul des pressions au pas de temps suivant
	1249	c ----------------------------------------------------------
	1250
[1549]	1251	psurf=psurf+dtphys*dpsurf(1) ! evolution de la pression de surface
[253]	1252	DO ilevel=1,nlevel
	1253	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	1254	ENDDO
	1255	DO ilayer=1,nlayer
	1256	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	1257	ENDDO
	1258
	1259	c calcul traceur au pas de temps suivant
	1260	c --------------------------------------
	1261
[787]	1262	DO iq = 1, nq
[253]	1263	DO ilayer=1,nlayer
	1264	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
	1265	ENDDO
	1266	END DO
	1267
	1268	c ========================================================
	1269	c GESTION DES SORTIE
	1270	c ========================================================
	1271	if(saveprofile)then
	1272	OPEN(12,file='profile.out',form='formatted')
[601]	1273	write(12,*) tsurf
[1308]	1274	DO ilayer=1,llm
[601]	1275	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
[253]	1276	ENDDO
	1277	CLOSE(12)
	1278	endif
	1279
[3562]	1280	c Produce the restart1D.nc file (MM)
	1281	if (lastcall) then
	1282	call writerestart1D('restart1D.nc',nlayer,nsoil,day,time,psurf,
	1283	& temp,tsoil,u,v,nq,q)
	1284	endif
[253]	1285
[601]	1286	ENDDO ! fin de la boucle temporelle
	1287
[1216]	1288	write(,) "rcm1d: Everything is cool."
	1289
[253]	1290	c ========================================================
	1291	end !rcm1d
	1292
	1293

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

Download in other formats:

Original Format