Context Navigation

newstart.F @ 1202

Last change on this file since 1202 was 1023, checked in by aslmd, 12 years ago
LMDZ.GENERIC. Removed useless inidissip in newstart.
File size: 52.5 KB

Rev	Line
[135]	1	C======================================================================
	2	PROGRAM newstart
	3	c=======================================================================
	4	c
	5	c
	6	c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick
	7	c ------
	8	c Derniere modif : 12/03
	9	c
	10	c
	11	c Objet: Create or modify the initial state for the LMD Mars GCM
	12	c ----- (fichiers NetCDF start et startfi)
	13	c
	14	c
	15	c=======================================================================
	16
[787]	17	USE tracer_h
	18	USE comsoil_h
	19	USE surfdat_h
	20	USE comgeomfi_h
[988]	21	use datafile_mod, only: datadir
	22	! to use 'getin'
	23	USE ioipsl_getincom, only: getin
[135]	24	implicit none
	25
	26	#include "dimensions.h"
	27	#include "dimphys.h"
	28	#include "planete.h"
	29	#include "paramet.h"
	30	#include "comconst.h"
	31	#include "comvert.h"
	32	#include "comgeom2.h"
	33	#include "control.h"
	34	#include "logic.h"
	35	#include "ener.h"
	36	#include "temps.h"
	37	#include "comdissnew.h"
	38	#include "serre.h"
	39	#include "netcdf.inc"
	40	#include "advtrac.h"
	41	c=======================================================================
	42	c Declarations
	43	c=======================================================================
	44
	45	c Variables dimension du fichier "start_archive"
	46	c------------------------------------
	47	CHARACTER relief*3
	48
	49
	50	c Variables pour les lectures NetCDF des fichiers "start_archive"
	51	c--------------------------------------------------
	52	INTEGER nid_dyn, nid_fi,nid,nvarid
	53	INTEGER length
	54	parameter (length = 100)
	55	INTEGER tab0
	56	INTEGER NB_ETATMAX
	57	parameter (NB_ETATMAX = 100)
	58
	59	REAL date
	60	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
	61
	62	c Variable histoire
	63	c------------------
	64	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	65	REAL phis(iip1,jjp1)
	66	REAL q(iip1,jjp1,llm,nqmx) ! champs advectes
	67
	68	c autre variables dynamique nouvelle grille
	69	c------------------------------------------
	70	REAL pks(iip1,jjp1)
	71	REAL w(iip1,jjp1,llm+1)
	72	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
	73	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	74	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
	75	REAL phi(iip1,jjp1,llm)
	76
	77	integer klatdat,klongdat
	78	PARAMETER (klatdat=180,klongdat=360)
	79
	80	c Physique sur grille scalaire
	81	c----------------------------
	82	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
	83	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
	84
	85	c variable physique
	86	c------------------
	87	REAL tsurf(ngridmx) ! surface temperature
	88	REAL tsoil(ngridmx,nsoilmx) ! soil temperature
	89	! REAL co2ice(ngridmx) ! CO2 ice layer
	90	REAL emis(ngridmx) ! surface emissivity
	91	real emisread ! added by RW
	92	REAL qsurf(ngridmx,nqmx)
	93	REAL q2(ngridmx,nlayermx+1)
	94	! REAL rnaturfi(ngridmx)
	95	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
	96	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
	97	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
	98	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
	99
	100	INTEGER i,j,l,isoil,ig,idum
	101	real mugaz ! molar mass of the atmosphere
	102
	103	integer ierr
	104
	105	c Variables on the new grid along scalar points
	106	c------------------------------------------------------
	107	! REAL p(iip1,jjp1)
	108	REAL t(iip1,jjp1,llm)
	109	REAL tset(iip1,jjp1,llm)
	110	real phisold_newgrid(iip1,jjp1)
	111	REAL :: teta(iip1, jjp1, llm)
	112	REAL :: pk(iip1,jjp1,llm)
	113	REAL :: pkf(iip1,jjp1,llm)
	114	REAL :: ps(iip1, jjp1)
	115	REAL :: masse(iip1,jjp1,llm)
	116	REAL :: xpn,xps,xppn(iim),xpps(iim)
	117	REAL :: p3d(iip1, jjp1, llm+1)
	118	REAL :: beta(iip1,jjp1,llm)
	119	! REAL dteta(ip1jmp1,llm)
	120
	121	c Variable de l'ancienne grille
	122	c------------------------------
	123	real time
	124	real tab_cntrl(100)
	125	real tab_cntrl_bis(100)
	126
	127	c variables diverses
	128	c-------------------
	129	real choix_1,pp
	130	character*80 fichnom
[988]	131	character*250 filestring
[135]	132	integer Lmodif,iq,thermo
	133	character modif*20
	134	real z_reel(iip1,jjp1)
	135	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
	136	real ptoto,pcap,patm,airetot,ptotn,patmn,psea
	137	! real ssum
	138	character*1 yes
	139	logical :: flagtset=.false. , flagps0=.false.
	140	real val, val2, val3 ! to store temporary variables
	141	real :: iceith=2000 ! thermal inertia of subterranean ice
	142	integer iref,jref
	143
	144	INTEGER :: itau
	145
	146	INTEGER :: nq,numvanle
	147	character(len=20) :: txt ! to store some text
[988]	148	character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
	149	character(len=150) :: longline
[135]	150	integer :: count
[699]	151	real :: profile(llm+1) ! to store an atmospheric profile + surface value
[135]	152
[253]	153	! added by RW for test
	154	real pmean, phi0
	155
	156	! added by BC for equilibrium temperature startup
	157	real teque
	158
	159	! added by BC for cloud fraction setup
	160	REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
	161	REAL totalfrac(ngridmx)
	162
	163	! added by RW for nuketharsis
	164	real fact1
	165	real fact2
	166
	167
[135]	168	c sortie visu pour les champs dynamiques
	169	c---------------------------------------
	170	! INTEGER :: visuid
	171	! real :: time_step,t_ops,t_wrt
	172	! CHARACTER*80 :: visu_file
	173
[535]	174	cpp = 0.
	175	preff = 0.
	176	pa = 0. ! to ensure disaster rather than minor error if we don`t
	177	! make deliberate choice of these values elsewhere.
[135]	178
	179	c=======================================================================
	180	c Choice of the start file(s) to use
	181	c=======================================================================
	182
	183	write(,) 'From which kind of files do you want to create new',
	184	. 'start and startfi files'
	185	write(,) ' 0 - from a file start_archive'
	186	write(,) ' 1 - from files start and startfi'
	187
	188	c-----------------------------------------------------------------------
	189	c Open file(s) to modify (start or start_archive)
	190	c-----------------------------------------------------------------------
	191
	192	DO
	193	read(,,iostat=ierr) choix_1
	194	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	195	ENDDO
	196
	197	c Open start_archive
	198	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	199	if (choix_1.eq.0) then
	200
	201	write(,) 'Creating start files from:'
	202	write(,) './start_archive.nc'
	203	write(,)
	204	fichnom = 'start_archive.nc'
	205	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	206	IF (ierr.NE.NF_NOERR) THEN
	207	write(6,*)' Problem opening file:',fichnom
	208	write(6,*)' ierr = ', ierr
	209	CALL ABORT
	210	ENDIF
	211	tab0 = 50
	212	Lmodif = 1
	213
	214	c OR open start and startfi files
	215	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	216	else
	217	write(,) 'Creating start files from:'
	218	write(,) './start.nc and ./startfi.nc'
	219	write(,)
	220	fichnom = 'start.nc'
	221	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	222	IF (ierr.NE.NF_NOERR) THEN
	223	write(6,*)' Problem opening file:',fichnom
	224	write(6,*)' ierr = ', ierr
	225	CALL ABORT
	226	ENDIF
	227
	228	fichnom = 'startfi.nc'
	229	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	230	IF (ierr.NE.NF_NOERR) THEN
	231	write(6,*)' Problem opening file:',fichnom
	232	write(6,*)' ierr = ', ierr
	233	CALL ABORT
	234	ENDIF
	235
	236	tab0 = 0
	237	Lmodif = 0
	238
	239	endif
	240
[649]	241
	242	c=======================================================================
	243	c INITIALISATIONS DIVERSES
	244	c=======================================================================
	245	! Load tracer names:
	246	call iniadvtrac(nq,numvanle)
	247	! tnom(:) now contains tracer names
	248	! JL12 we will need the tracer names to read start in dyneta0
	249
[787]	250	! Initialize global tracer indexes (stored in tracer.h)
	251	! ... this has to be done before phyetat0
	252	call initracer(ngridmx,nqmx,tnom)
	253
[837]	254	! Take care of arrays in common modules
	255	! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
	256	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
	257	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
	258	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
	259	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
	260	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
	261	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
	262	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
	263	! ALLOCATE ARRAYS in comsoil_h (if not already done)
	264	IF (.not.ALLOCATED(layer))
	265	. ALLOCATE(layer(nsoilmx))
	266	IF (.not.ALLOCATED(mlayer))
	267	. ALLOCATE(mlayer(0:nsoilmx-1))
	268	IF (.not.ALLOCATED(inertiedat))
	269	. ALLOCATE(inertiedat(ngridmx,nsoilmx))
	270	! ALLOCATE ARRAYS IN comgeomfi_h (done in inifis usually)
	271	IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngridmx))
	272	IF (.not. ALLOCATED(long)) ALLOCATE(long(ngridmx))
	273	IF (.not. ALLOCATED(area)) ALLOCATE(area(ngridmx))
	274	cursor = 1 ! added by AS in dimphys. 1 for sequential runs.
[787]	275
[135]	276	c-----------------------------------------------------------------------
	277	c Lecture du tableau des parametres du run (pour la dynamique)
	278	c-----------------------------------------------------------------------
	279
	280	if (choix_1.eq.0) then
	281
	282	write(,) 'reading tab_cntrl START_ARCHIVE'
	283	c
	284	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	285	#ifdef NC_DOUBLE
	286	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	287	#else
	288	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	289	#endif
	290	c
	291	else if (choix_1.eq.1) then
	292
	293	write(,) 'reading tab_cntrl START'
	294	c
	295	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	296	#ifdef NC_DOUBLE
	297	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	298	#else
	299	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	300	#endif
	301	c
	302	write(,) 'reading tab_cntrl STARTFI'
	303	c
	304	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	305	#ifdef NC_DOUBLE
	306	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	307	#else
	308	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	309	#endif
	310	c
	311	do i=1,50
	312	tab_cntrl(i+50)=tab_cntrl_bis(i)
	313	enddo
[649]	314	write(,) 'printing tab_cntrl', tab_cntrl
	315	do i=1,100
	316	write(,) i,tab_cntrl(i)
	317	enddo
	318
	319	! Lmodif set to 0 to disable modifications possibility in phyeta0
	320	write(,) 'Reading file START'
	321	fichnom = 'start.nc'
	322	CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
	323	. ps,phis,time)
	324
	325	write(,) 'Reading file STARTFI'
	326	fichnom = 'startfi.nc'
[787]	327	CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqmx,
[649]	328	. day_ini,time,
	329	. tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW
	330	. cloudfrac,totalfrac,hice)
	331
[905]	332	! copy albedo and soil thermal inertia on (local) physics grid
[649]	333	do i=1,ngridmx
	334	albfi(i) = albedodat(i)
	335	do j=1,nsoilmx
	336	ithfi(i,j) = inertiedat(i,j)
	337	enddo
	338	! build a surfithfi(:) using 1st layer of ithfi(:), which might
[905]	339	! be needed later on if reinitializing soil thermal inertia
[649]	340	surfithfi(i)=ithfi(i,1)
	341	enddo
[905]	342	! also copy albedo and soil thermal inertia on (local) dynamics grid
	343	! so that options below can manipulate either (but must then ensure
	344	! to correctly recast things on physics grid)
	345	call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
	346	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	347	call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
[135]	348
	349	endif
	350	c-----------------------------------------------------------------------
	351	c Initialisation des constantes dynamique
	352	c-----------------------------------------------------------------------
	353
	354	kappa = tab_cntrl(9)
	355	etot0 = tab_cntrl(12)
	356	ptot0 = tab_cntrl(13)
	357	ztot0 = tab_cntrl(14)
	358	stot0 = tab_cntrl(15)
	359	ang0 = tab_cntrl(16)
	360	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	361	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	362
	363	! for vertical coordinate
	364	preff=tab_cntrl(18) ! reference surface pressure
	365	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
	366	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
	367	write(,) "Newstart: preff=",preff," pa=",pa
	368	yes=' '
	369	do while ((yes.ne.'y').and.(yes.ne.'n'))
	370	write(,) "Change the values of preff and pa? (y/n)"
	371	read(*,fmt='(a)') yes
	372	end do
	373	if (yes.eq.'y') then
	374	write(,)"New value of reference surface pressure preff?"
	375	write(,)" (for Mars, typically preff=610)"
	376	read(,) preff
	377	write(,)"New value of reference pressure pa for purely"
	378	write(,)"pressure levels (for hybrid coordinates)?"
	379	write(,)" (for Mars, typically pa=20)"
	380	read(,) pa
	381	endif
	382	c-----------------------------------------------------------------------
	383	c Lecture du tab_cntrl et initialisation des constantes physiques
	384	c - pour start: Lmodif = 0 => pas de modifications possibles
	385	c (modif dans le tabfi de readfi + loin)
	386	c - pour start_archive: Lmodif = 1 => modifications possibles
	387	c-----------------------------------------------------------------------
	388	if (choix_1.eq.0) then
[787]	389	call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
[135]	390	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	391	else if (choix_1.eq.1) then
[649]	392	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
[787]	393	call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
[135]	394	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	395	endif
	396
	397	rad = p_rad
	398	omeg = p_omeg
	399	g = p_g
	400	cpp = p_cpp
	401	mugaz = p_mugaz
	402	daysec = p_daysec
	403
[535]	404	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
[135]	405
	406	c=======================================================================
	407	c INITIALISATIONS DIVERSES
	408	c=======================================================================
	409	! Load parameters from run.def file
	410	CALL defrun_new( 99, .TRUE. )
	411	CALL iniconst
	412	CALL inigeom
	413	idum=-1
	414	idum=0
	415
	416	c Initialisation coordonnees /aires
	417	c -------------------------------
	418	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
	419	! rlatu() and rlonv() are given in radians
	420	latfi(1)=rlatu(1)
	421	lonfi(1)=0.
	422	DO j=2,jjm
	423	DO i=1,iim
	424	latfi((j-2)*iim+1+i)=rlatu(j)
	425	lonfi((j-2)*iim+1+i)=rlonv(i)
	426	ENDDO
	427	ENDDO
	428	latfi(ngridmx)=rlatu(jjp1)
	429	lonfi(ngridmx)=0.
[699]	430
	431	! build airefi(), mesh area on physics grid
[135]	432	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
[699]	433	! Poles are single points on physics grid
	434	airefi(1)=airefi(1)*iim
	435	airefi(ngridmx)=airefi(ngridmx)*iim
[135]	436
	437	c=======================================================================
	438	c lecture topographie, albedo, inertie thermique, relief sous-maille
	439	c=======================================================================
	440
[988]	441	if (choix_1.eq.0) then ! for start_archive files,
	442	! where an external "surface.nc" file is needed
[135]	443
	444	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	445	c write(,)
	446	c write(,) 'choix du relief (mola,pla)'
	447	c write(,) '(Topographie MGS MOLA, plat)'
	448	c read(*,fmt='(a3)') relief
	449	relief="mola"
	450	c enddo
	451
[988]	452	! First get the correct value of datadir, if not already done:
	453	! default 'datadir' is set in "datafile_mod"
	454	call getin("datadir",datadir)
	455	write(,) 'Available surface data files are:'
	456	filestring='ls '//trim(datadir)//' \| grep .nc'
	457	call system(filestring)
[135]	458
[988]	459	write(,) ''
	460	write(,) 'Please choose the relevant file',
	461	& ' (e.g. "surface_mars.nc")'
	462	write(,) ' or "none" to not use any (and set planetary'
	463	write(,) ' albedo and surface thermal inertia)'
	464	read(*,fmt='(a50)') surfacefile
[135]	465
[988]	466	if (surfacefile.ne."none") then
[135]	467
[988]	468	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
	469	& zmeaS,zstdS,zsigS,zgamS,ztheS)
	470	else
	471	! specific case when not using a "surface.nc" file
	472	phis(:,:)=0
	473	zmeaS(:,:)=0
	474	zstdS(:,:)=0
	475	zsigS(:,:)=0
	476	zgamS(:,:)=0
	477	ztheS(:,:)=0
	478
	479	write(,) "Enter value of albedo of the bare ground:"
	480	read(,) alb(1,1)
	481	alb(:,:)=alb(1,1)
	482
	483	write(,) "Enter value of thermal inertia of soil:"
	484	read(,) surfith(1,1)
	485	surfith(:,:)=surfith(1,1)
	486
	487	endif ! of if (surfacefile.ne."none")
[135]	488
[988]	489	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	490	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	491	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	492	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	493	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	494	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	495	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	496	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
	497
	498	endif ! of if (choix_1.eq.0)
	499
	500
[135]	501	c=======================================================================
	502	c Lecture des fichiers (start ou start_archive)
	503	c=======================================================================
	504
	505	if (choix_1.eq.0) then
	506
	507	write(,) 'Reading file START_ARCHIVE'
	508	CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
	509	. t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
	510	& surfith,nid)
	511	write(,) "OK, read start_archive file"
	512	! copy soil thermal inertia
	513	ithfi(:,:)=inertiedat(:,:)
	514
	515	ierr= NF_CLOSE(nid)
	516
[649]	517	else if (choix_1.eq.1) then
	518	!do nothing, start and startfi have already been read
[135]	519	else
	520	CALL exit(1)
	521	endif
	522
	523	dtvr = daysec/FLOAT(day_step)
	524	dtphys = dtvr * FLOAT(iphysiq)
	525
	526	c=======================================================================
	527	c
	528	c=======================================================================
	529
	530	do ! infinite loop on list of changes
	531
	532	write(,)
	533	write(,)
	534	write(,) 'List of possible changes :'
	535	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	536	write(,)
	537	write(,) 'flat : no topography ("aquaplanet")'
[253]	538	write(,) 'nuketharsis : no Tharsis bulge'
[135]	539	write(,) 'bilball : uniform albedo and thermal inertia'
	540	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
	541	write(,) 'qname : change tracer name'
[996]	542	write(,) 't=profile : read temperature profile in profile.in'
[135]	543	write(,) 'q=0 : ALL tracer =zero'
	544	write(,) 'q=x : give a specific uniform value to one tracer'
[699]	545	write(,) 'q=profile : specify a profile for a tracer'
	546	! write(,) 'ini_q : tracers initialisation for chemistry, water an
	547	! $d ice '
	548	! write(,) 'ini_q-H2O : tracers initialisation for chemistry and
	549	! $ice '
	550	! write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
	551	! $ly '
[135]	552	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
	553	write(,) 'watercapn : H20 ice on permanent N polar cap '
	554	write(,) 'watercaps : H20 ice on permanent S polar cap '
[253]	555	write(,) 'noacglac : H2O ice across Noachis Terra'
[135]	556	write(,) 'oborealis : H2O ice across Vastitas Borealis'
[253]	557	write(,) 'iceball : Thick ice layer all over surface'
[135]	558	write(,) 'wetstart : start with a wet atmosphere'
[253]	559	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
	560	write(,) 'radequi : Earth-like radiative equilibrium temperature
	561	$ profile (lat-alt) and winds set to zero'
[135]	562	write(,) 'coldstart : Start X K above the CO2 frost point and
	563	$set wind to zero (assumes 100% CO2)'
	564	write(,) 'co2ice=0 : remove CO2 polar cap'
	565	write(,) 'ptot : change total pressure'
	566	write(,) 'emis : change surface emissivity'
[253]	567	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
[135]	568	&face values'
	569
	570	write(,)
	571	write(,) 'Change to perform ?'
	572	write(,) ' (enter keyword or return to end)'
	573	write(,)
	574
	575	read(*,fmt='(a20)') modif
	576	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
	577
	578	write(,)
	579	write(,) trim(modif) , ' : '
	580
	581	c 'flat : no topography ("aquaplanet")'
	582	c -------------------------------------
[699]	583	if (trim(modif) .eq. 'flat') then
[135]	584	c set topo to zero
[699]	585	z_reel(1:iip1,1:jjp1)=0
	586	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
[135]	587	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	588	write(,) 'topography set to zero.'
	589	write(,) 'WARNING : the subgrid topography parameters',
	590	& ' were not set to zero ! => set calllott to F'
	591
[253]	592	c Choice of surface pressure
[135]	593	yes=' '
	594	do while ((yes.ne.'y').and.(yes.ne.'n'))
	595	write(,) 'Do you wish to choose homogeneous surface',
	596	& 'pressure (y) or let newstart interpolate ',
	597	& ' the previous field (n)?'
	598	read(*,fmt='(a)') yes
	599	end do
	600	if (yes.eq.'y') then
	601	flagps0=.true.
	602	write(,) 'New value for ps (Pa) ?'
	603	201 read(,,iostat=ierr) patm
	604	if(ierr.ne.0) goto 201
	605	write(,)
	606	write(,) ' new ps everywhere (Pa) = ', patm
	607	write(,)
	608	do j=1,jjp1
	609	do i=1,iip1
	610	ps(i,j)=patm
	611	enddo
	612	enddo
	613	end if
	614
[253]	615	c 'nuketharsis : no tharsis bulge for Early Mars'
	616	c ---------------------------------------------
[699]	617	else if (trim(modif) .eq. 'nuketharsis') then
[253]	618
	619	DO j=1,jjp1
	620	DO i=1,iim
	621	ig=1+(j-2)*iim +i
	622	if(j.eq.1) ig=1
	623	if(j.eq.jjp1) ig=ngridmx
	624
	625	fact1=(((rlonv(i)180./pi)+100)*2 +
	626	& (rlatu(j)180./pi)2)/65*2
	627	fact2=exp( -fact1**2.5 )
	628
	629	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
	630
	631	! if(phis(i,j).gt.2500.*g)then
	632	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
	633	! phis(i,j)=2500.*g
	634	! endif
	635	! endif
	636
	637	ENDDO
	638	ENDDO
	639	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	640
	641
	642	c bilball : uniform albedo, thermal inertia
	643	c -----------------------------------------
[699]	644	else if (trim(modif) .eq. 'bilball') then
[135]	645	write(,) 'constante albedo and iner.therm:'
	646	write(,) 'New value for albedo (ex: 0.25) ?'
	647	101 read(,,iostat=ierr) alb_bb
	648	if(ierr.ne.0) goto 101
	649	write(,)
	650	write(,) ' uniform albedo (new value):',alb_bb
	651	write(,)
	652
	653	write(,) 'New value for thermal inertia (eg: 247) ?'
	654	102 read(,,iostat=ierr) ith_bb
	655	if(ierr.ne.0) goto 102
	656	write(,) 'uniform thermal inertia (new value):',ith_bb
	657	DO j=1,jjp1
	658	DO i=1,iip1
	659	alb(i,j) = alb_bb ! albedo
	660	do isoil=1,nsoilmx
	661	ith(i,j,isoil) = ith_bb ! thermal inertia
	662	enddo
	663	END DO
	664	END DO
	665	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	666	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	667	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	668
	669	c coldspole : sous-sol de la calotte sud toujours froid
	670	c -----------------------------------------------------
[699]	671	else if (trim(modif) .eq. 'coldspole') then
[135]	672	write(,)'new value for the subsurface temperature',
	673	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
	674	103 read(,,iostat=ierr) tsud
	675	if(ierr.ne.0) goto 103
	676	write(,)
	677	write(,) ' new value of the subsurface temperature:',tsud
	678	c nouvelle temperature sous la calotte permanente
	679	do l=2,nsoilmx
	680	tsoil(ngridmx,l) = tsud
	681	end do
	682
	683
	684	write(,)'new value for the albedo',
	685	& 'of the permanent southern polar cap ? (eg: 0.75)'
	686	104 read(,,iostat=ierr) albsud
	687	if(ierr.ne.0) goto 104
	688	write(,)
	689
	690	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	691	c Option 1: only the albedo of the pole is modified :
	692	albfi(ngridmx)=albsud
	693	write(,) 'ig=',ngridmx,' albedo perennial cap ',
	694	& albfi(ngridmx)
	695
	696	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	697	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
	698	c DO j=1,jjp1
	699	c DO i=1,iip1
	700	c ig=1+(j-2)*iim +i
	701	c if(j.eq.1) ig=1
	702	c if(j.eq.jjp1) ig=ngridmx
	703	c if ((rlatu(j)*180./pi.lt.-84.).and.
	704	c & (rlatu(j)*180./pi.gt.-91.).and.
	705	c & (rlonv(i)*180./pi.gt.-91.).and.
	706	c & (rlonv(i)*180./pi.lt.0.)) then
	707	cc albedo de la calotte permanente fixe a albsud
	708	c alb(i,j)=albsud
	709	c write(,) 'lat=',rlatu(j)*180./pi,
	710	c & ' lon=',rlonv(i)*180./pi
	711	cc fin de la condition sur les limites de la calotte permanente
	712	c end if
	713	c ENDDO
	714	c ENDDO
	715	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	716
	717	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	718
	719
	720	c ptot : Modification of the total pressure: ice + current atmosphere
	721	c -------------------------------------------------------------------
[699]	722	else if (trim(modif).eq.'ptot') then
[135]	723
	724	! check if we have a co2_ice surface tracer:
	725	if (igcm_co2_ice.eq.0) then
	726	write(,) " No surface CO2 ice !"
	727	write(,) " only atmospheric pressure will be considered!"
	728	endif
	729	c calcul de la pression totale glace + atm actuelle
	730	patm=0.
	731	airetot=0.
	732	pcap=0.
	733	DO j=1,jjp1
	734	DO i=1,iim
	735	ig=1+(j-2)*iim +i
	736	if(j.eq.1) ig=1
	737	if(j.eq.jjp1) ig=ngridmx
	738	patm = patm + ps(i,j)*aire(i,j)
	739	airetot= airetot + aire(i,j)
	740	if (igcm_co2_ice.ne.0) then
	741	!pcap = pcap + aire(i,j)co2ice(ig)g
	742	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2_ice)g
	743	endif
	744	ENDDO
	745	ENDDO
	746	ptoto = pcap + patm
	747
	748	print*,'Current total pressure at surface (co2 ice + atm) ',
	749	& ptoto/airetot
	750
	751	print*,'new value?'
	752	read(,) ptotn
	753	ptotn=ptotn*airetot
	754	patmn=ptotn-pcap
	755	print*,'ptoto,patm,ptotn,patmn'
	756	print*,ptoto,patm,ptotn,patmn
	757	print*,'Mult. factor for pressure (atm only)', patmn/patm
	758	do j=1,jjp1
	759	do i=1,iip1
	760	ps(i,j)=ps(i,j)*patmn/patm
	761	enddo
	762	enddo
	763
[253]	764
	765
[135]	766	c Correction pour la conservation des traceurs
	767	yes=' '
	768	do while ((yes.ne.'y').and.(yes.ne.'n'))
	769	write(,) 'Do you wish to conserve tracer total mass (y)',
	770	& ' or tracer mixing ratio (n) ?'
	771	read(*,fmt='(a)') yes
	772	end do
	773
	774	if (yes.eq.'y') then
	775	write(,) 'OK : conservation of tracer total mass'
	776	DO iq =1, nqmx
	777	DO l=1,llm
	778	DO j=1,jjp1
	779	DO i=1,iip1
	780	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	781	ENDDO
	782	ENDDO
	783	ENDDO
	784	ENDDO
	785	else
	786	write(,) 'OK : conservation of tracer mixing ratio'
	787	end if
	788
	789	c Correction pour la pression au niveau de la mer
	790	yes=' '
	791	do while ((yes.ne.'y').and.(yes.ne.'n'))
	792	write(,) 'Do you wish fix pressure at sea level (y)',
	793	& ' or not (n) ?'
	794	read(*,fmt='(a)') yes
	795	end do
	796
	797	if (yes.eq.'y') then
	798	write(,) 'Value?'
	799	read(,,iostat=ierr) psea
	800	DO i=1,iip1
	801	DO j=1,jjp1
	802	ps(i,j)=psea
	803
	804	ENDDO
	805	ENDDO
	806	write(,) 'psea=',psea
	807	else
	808	write(,) 'no'
	809	end if
	810
	811
	812	c emis : change surface emissivity (added by RW)
	813	c ----------------------------------------------
	814	else if (trim(modif).eq.'emis') then
	815
	816	print*,'new value?'
	817	read(,) emisread
	818
	819	do i=1,ngridmx
	820	emis(i)=emisread
	821	enddo
	822
	823	c qname : change tracer name
	824	c --------------------------
	825	else if (trim(modif).eq.'qname') then
	826	yes='y'
	827	do while (yes.eq.'y')
	828	write(,) 'Which tracer name do you want to change ?'
	829	do iq=1,nqmx
	830	write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
	831	enddo
	832	write(*,'(a35,i3)')
	833	& '(enter tracer number; between 1 and ',nqmx
	834	write(,)' or any other value to quit this option)'
	835	read(,) iq
	836	if ((iq.ge.1).and.(iq.le.nqmx)) then
	837	write(,)'Change tracer name ',trim(tnom(iq)),' to ?'
	838	read(,) txt
	839	tnom(iq)=txt
	840	write(,)'Do you want to change another tracer name (y/n)?'
	841	read(*,'(a)') yes
	842	else
	843	! inapropiate value of iq; quit this option
	844	yes='n'
	845	endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
	846	enddo ! of do while (yes.ne.'y')
	847
	848	c q=0 : set tracers to zero
	849	c -------------------------
[699]	850	else if (trim(modif).eq.'q=0') then
[135]	851	c mise a 0 des q (traceurs)
	852	write(,) 'Tracers set to 0 (1.E-30 in fact)'
	853	DO iq =1, nqmx
	854	DO l=1,llm
	855	DO j=1,jjp1
	856	DO i=1,iip1
	857	q(i,j,l,iq)=1.e-30
	858	ENDDO
	859	ENDDO
	860	ENDDO
	861	ENDDO
	862
	863	c set surface tracers to zero
	864	DO iq =1, nqmx
	865	DO ig=1,ngridmx
	866	qsurf(ig,iq)=0.
	867	ENDDO
	868	ENDDO
	869
	870	c q=x : initialise tracer manually
	871	c --------------------------------
[699]	872	else if (trim(modif).eq.'q=x') then
[135]	873	write(,) 'Which tracer do you want to modify ?'
	874	do iq=1,nqmx
	875	write(,)iq,' : ',trim(tnom(iq))
	876	enddo
	877	write(,) '(choose between 1 and ',nqmx,')'
	878	read(,) iq
	879	write(,)'mixing ratio of tracer ',trim(tnom(iq)),
	880	& ' ? (kg/kg)'
	881	read(,) val
	882	DO l=1,llm
	883	DO j=1,jjp1
	884	DO i=1,iip1
	885	q(i,j,l,iq)=val
	886	ENDDO
	887	ENDDO
	888	ENDDO
	889	write(,) 'SURFACE value of tracer ',trim(tnom(iq)),
	890	& ' ? (kg/m2)'
	891	read(,) val
	892	DO ig=1,ngridmx
	893	qsurf(ig,iq)=val
	894	ENDDO
	895
[996]	896	c t=profile : initialize temperature with a given profile
	897	c -------------------------------------------------------
	898	else if (trim(modif) .eq. 't=profile') then
	899	write(,) 'Temperature profile from ASCII file'
	900	write(,) "'profile.in' e.g. 1D output"
	901	write(,) "(one value per line in file; starting with"
	902	write(,) "surface value, the 1st atmospheric layer"
	903	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	904	txt="profile.in"
	905	open(33,file=trim(txt),status='old',form='formatted',
	906	& iostat=ierr)
	907	if (ierr.eq.0) then
	908	! OK, found file 'profile_...', load the profile
	909	do l=1,llm+1
	910	read(33,*,iostat=ierr) profile(l)
	911	write(,) profile(l)
	912	if (ierr.ne.0) then ! something went wrong
	913	exit ! quit loop
	914	endif
	915	enddo
	916	if (ierr.eq.0) then
	917	tsurf(1:ngridmx)=profile(1)
	918	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
	919	do l=1,llm
	920	Tset(1:iip1,1:jjp1,l)=profile(l+1)
	921	flagtset=.true.
	922	enddo
	923	ucov(1:iip1,1:jjp1,1:llm)=0.
	924	vcov(1:iip1,1:jjm,1:llm)=0.
	925	q2(1:ngridmx,1:nlayermx+1)=0.
	926	else
	927	write(,)'problem reading file ',trim(txt),' !'
	928	write(,)'No modifications to temperature'
	929	endif
	930	else
	931	write(,)'Could not find file ',trim(txt),' !'
	932	write(,)'No modifications to temperature'
	933	endif
	934
[699]	935	c q=profile : initialize tracer with a given profile
	936	c --------------------------------------------------
	937	else if (trim(modif) .eq. 'q=profile') then
	938	write(,) 'Tracer profile will be sought in ASCII file'
	939	write(,) "'profile_tracer' where 'tracer' is tracer name"
	940	write(,) "(one value per line in file; starting with"
	941	write(,) "surface value, the 1st atmospheric layer"
	942	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	943	write(,) 'Which tracer do you want to set?'
	944	do iq=1,nqmx
	945	write(,)iq,' : ',trim(tnom(iq))
	946	enddo
	947	write(,) '(choose between 1 and ',nqmx,')'
	948	read(,) iq
	949	if ((iq.lt.1).or.(iq.gt.nqmx)) then
	950	! wrong value for iq, go back to menu
	951	write(,) "wrong input value:",iq
	952	cycle
	953	endif
	954	! look for input file 'profile_tracer'
	955	txt="profile_"//trim(tnom(iq))
	956	open(41,file=trim(txt),status='old',form='formatted',
	957	& iostat=ierr)
	958	if (ierr.eq.0) then
	959	! OK, found file 'profile_...', load the profile
	960	do l=1,llm+1
	961	read(41,*,iostat=ierr) profile(l)
	962	if (ierr.ne.0) then ! something went wrong
	963	exit ! quit loop
	964	endif
	965	enddo
	966	if (ierr.eq.0) then
	967	! initialize tracer values
	968	qsurf(:,iq)=profile(1)
	969	do l=1,llm
	970	q(:,:,l,iq)=profile(l+1)
	971	enddo
	972	write(,)'OK, tracer ',trim(tnom(iq)),' initialized ',
	973	& 'using values from file ',trim(txt)
	974	else
	975	write(,)'problem reading file ',trim(txt),' !'
	976	write(,)'No modifications to tracer ',trim(tnom(iq))
	977	endif
	978	else
	979	write(,)'Could not find file ',trim(txt),' !'
	980	write(,)'No modifications to tracer ',trim(tnom(iq))
	981	endif
	982
[135]	983
	984	c wetstart : wet atmosphere with a north to south gradient
	985	c --------------------------------------------------------
[699]	986	else if (trim(modif) .eq. 'wetstart') then
[135]	987	! check that there is indeed a water vapor tracer
	988	if (igcm_h2o_vap.eq.0) then
	989	write(,) "No water vapour tracer! Can't use this option"
	990	stop
	991	endif
	992	DO l=1,llm
	993	DO j=1,jjp1
	994	DO i=1,iip1
	995	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
	996	ENDDO
	997	ENDDO
	998	ENDDO
	999
	1000	write(,) 'Water mass mixing ratio at north pole='
	1001	* ,q(1,1,1,igcm_h2o_vap)
	1002	write(,) '---------------------------south pole='
	1003	* ,q(1,jjp1,1,igcm_h2o_vap)
	1004
	1005	c noglacier : remove tropical water ice (to initialize high res sim)
	1006	c --------------------------------------------------
[699]	1007	else if (trim(modif) .eq. 'noglacier') then
[253]	1008	if (igcm_h2o_ice.eq.0) then
	1009	write(,) "No water ice tracer! Can't use this option"
	1010	stop
	1011	endif
[135]	1012	do ig=1,ngridmx
	1013	j=(ig-2)/iim +2
	1014	if(ig.eq.1) j=1
	1015	write(,) 'OK: remove surface ice for \|lat\|<45'
	1016	if (abs(rlatu(j)*180./pi).lt.45.) then
	1017	qsurf(ig,igcm_h2o_ice)=0.
	1018	end if
	1019	end do
	1020
	1021
	1022	c watercapn : H20 ice on permanent northern cap
	1023	c --------------------------------------------------
[699]	1024	else if (trim(modif) .eq. 'watercapn') then
[253]	1025	if (igcm_h2o_ice.eq.0) then
	1026	write(,) "No water ice tracer! Can't use this option"
	1027	stop
	1028	endif
[135]	1029
[253]	1030	DO j=1,jjp1
	1031	DO i=1,iim
	1032	ig=1+(j-2)*iim +i
	1033	if(j.eq.1) ig=1
	1034	if(j.eq.jjp1) ig=ngridmx
	1035
	1036	if (rlatu(j)*180./pi.gt.80.) then
[367]	1037	qsurf(ig,igcm_h2o_ice)=3.4e3
	1038	!do isoil=1,nsoilmx
	1039	! ith(i,j,isoil) = 18000. ! thermal inertia
	1040	!enddo
[253]	1041	write(,)' ==> Ice mesh North boundary (deg)= ',
	1042	& rlatv(j-1)*180./pi
	1043	end if
	1044	ENDDO
	1045	ENDDO
	1046	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1047
	1048	c$$$ do ig=1,ngridmx
	1049	c$$$ j=(ig-2)/iim +2
	1050	c$$$ if(ig.eq.1) j=1
	1051	c$$$ if (rlatu(j)*180./pi.gt.80.) then
	1052	c$$$
	1053	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
	1054	c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5
	1055	c$$$
	1056	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	1057	c$$$ & qsurf(ig,igcm_h2o_ice)
	1058	c$$$
	1059	c$$$ write(,)' ==> Ice mesh South boundary (deg)= ',
	1060	c$$$ & rlatv(j)*180./pi
	1061	c$$$ end if
	1062	c$$$ enddo
	1063
[135]	1064	c watercaps : H20 ice on permanent southern cap
	1065	c -------------------------------------------------
[699]	1066	else if (trim(modif) .eq. 'watercaps') then
[253]	1067	if (igcm_h2o_ice.eq.0) then
	1068	write(,) "No water ice tracer! Can't use this option"
	1069	stop
	1070	endif
[135]	1071
[253]	1072	DO j=1,jjp1
	1073	DO i=1,iim
	1074	ig=1+(j-2)*iim +i
	1075	if(j.eq.1) ig=1
	1076	if(j.eq.jjp1) ig=ngridmx
	1077
	1078	if (rlatu(j)*180./pi.lt.-80.) then
[367]	1079	qsurf(ig,igcm_h2o_ice)=3.4e3
	1080	!do isoil=1,nsoilmx
	1081	! ith(i,j,isoil) = 18000. ! thermal inertia
	1082	!enddo
[253]	1083	write(,)' ==> Ice mesh South boundary (deg)= ',
	1084	& rlatv(j-1)*180./pi
	1085	end if
	1086	ENDDO
	1087	ENDDO
	1088	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1089
	1090	c$$$ do ig=1,ngridmx
	1091	c$$$ j=(ig-2)/iim +2
	1092	c$$$ if(ig.eq.1) j=1
	1093	c$$$ if (rlatu(j)*180./pi.lt.-80.) then
	1094	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
	1095	c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
	1096	c$$$
	1097	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	1098	c$$$ & qsurf(ig,igcm_h2o_ice)
	1099	c$$$ write(,)' ==> Ice mesh North boundary (deg)= ',
	1100	c$$$ & rlatv(j-1)*180./pi
	1101	c$$$ end if
	1102	c$$$ enddo
	1103
	1104
	1105	c noacglac : H2O ice across highest terrain
[135]	1106	c --------------------------------------------
[699]	1107	else if (trim(modif) .eq. 'noacglac') then
[253]	1108	if (igcm_h2o_ice.eq.0) then
	1109	write(,) "No water ice tracer! Can't use this option"
	1110	stop
	1111	endif
	1112	DO j=1,jjp1
	1113	DO i=1,iim
	1114	ig=1+(j-2)*iim +i
	1115	if(j.eq.1) ig=1
	1116	if(j.eq.jjp1) ig=ngridmx
	1117
	1118	if(phis(i,j).gt.1000.*g)then
	1119	alb(i,j) = 0.5 ! snow value
	1120	do isoil=1,nsoilmx
	1121	ith(i,j,isoil) = 18000. ! thermal inertia
	1122	! this leads to rnat set to 'ocean' in physiq.F90
	1123	! actually no, because it is soil not surface
	1124	enddo
	1125	endif
	1126	ENDDO
	1127	ENDDO
	1128	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1129	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1130	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1131
	1132
	1133
	1134	c oborealis : H2O oceans across Vastitas Borealis
	1135	c -----------------------------------------------
[699]	1136	else if (trim(modif) .eq. 'oborealis') then
[253]	1137	if (igcm_h2o_ice.eq.0) then
	1138	write(,) "No water ice tracer! Can't use this option"
	1139	stop
	1140	endif
[135]	1141	DO j=1,jjp1
	1142	DO i=1,iim
	1143	ig=1+(j-2)*iim +i
	1144	if(j.eq.1) ig=1
	1145	if(j.eq.jjp1) ig=ngridmx
	1146
[253]	1147	if(phis(i,j).lt.-4000.*g)then
	1148	! if( (phis(i,j).lt.-4000.*g)
	1149	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
	1150	! phis(i,j)=-8000.0*g ! proper ocean
	1151	phis(i,j)=-4000.0*g ! scanty ocean
	1152
	1153	alb(i,j) = 0.07 ! oceanic value
	1154	do isoil=1,nsoilmx
	1155	ith(i,j,isoil) = 18000. ! thermal inertia
	1156	! this leads to rnat set to 'ocean' in physiq.F90
	1157	! actually no, because it is soil not surface
	1158	enddo
[135]	1159	endif
	1160	ENDDO
	1161	ENDDO
[253]	1162	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1163	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1164	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
[135]	1165
[253]	1166	c iborealis : H2O ice in Northern plains
	1167	c --------------------------------------
[699]	1168	else if (trim(modif) .eq. 'iborealis') then
[253]	1169	if (igcm_h2o_ice.eq.0) then
	1170	write(,) "No water ice tracer! Can't use this option"
	1171	stop
	1172	endif
	1173	DO j=1,jjp1
	1174	DO i=1,iim
	1175	ig=1+(j-2)*iim +i
	1176	if(j.eq.1) ig=1
	1177	if(j.eq.jjp1) ig=ngridmx
	1178
	1179	if(phis(i,j).lt.-4000.*g)then
[367]	1180	!qsurf(ig,igcm_h2o_ice)=1.e3
	1181	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
[253]	1182	endif
	1183	ENDDO
	1184	ENDDO
	1185	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1186	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1187	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1188
	1189
	1190	c oceanball : H2O liquid everywhere
	1191	c ----------------------------
[699]	1192	else if (trim(modif) .eq. 'oceanball') then
[253]	1193	if (igcm_h2o_ice.eq.0) then
	1194	write(,) "No water ice tracer! Can't use this option"
	1195	stop
	1196	endif
	1197	DO j=1,jjp1
	1198	DO i=1,iim
	1199	ig=1+(j-2)*iim +i
	1200	if(j.eq.1) ig=1
	1201	if(j.eq.jjp1) ig=ngridmx
	1202
	1203	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
	1204	qsurf(ig,igcm_h2o_ice)=0.0
	1205	alb(i,j) = 0.07 ! ocean value
	1206
	1207	do isoil=1,nsoilmx
	1208	ith(i,j,isoil) = 18000. ! thermal inertia
	1209	!ith(i,j,isoil) = 50. ! extremely low for test
	1210	! this leads to rnat set to 'ocean' in physiq.F90
	1211	enddo
	1212
	1213	ENDDO
	1214	ENDDO
	1215	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1216	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1217	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1218
	1219	c iceball : H2O ice everywhere
	1220	c ----------------------------
[699]	1221	else if (trim(modif) .eq. 'iceball') then
[253]	1222	if (igcm_h2o_ice.eq.0) then
	1223	write(,) "No water ice tracer! Can't use this option"
	1224	stop
	1225	endif
	1226	DO j=1,jjp1
	1227	DO i=1,iim
	1228	ig=1+(j-2)*iim +i
	1229	if(j.eq.1) ig=1
	1230	if(j.eq.jjp1) ig=ngridmx
	1231
	1232	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
	1233	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
	1234	alb(i,j) = 0.6 ! ice albedo value
	1235
	1236	do isoil=1,nsoilmx
	1237	!ith(i,j,isoil) = 18000. ! thermal inertia
	1238	! this leads to rnat set to 'ocean' in physiq.F90
	1239	enddo
	1240
	1241	ENDDO
	1242	ENDDO
	1243	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1244	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1245
[135]	1246	c isotherm : Isothermal temperatures and no winds
[253]	1247	c -----------------------------------------------
[699]	1248	else if (trim(modif) .eq. 'isotherm') then
[135]	1249
	1250	write(,)'Isothermal temperature of the atmosphere,
	1251	& surface and subsurface'
	1252	write(,) 'Value of this temperature ? :'
	1253	203 read(,,iostat=ierr) Tiso
	1254	if(ierr.ne.0) goto 203
	1255
[699]	1256	tsurf(1:ngridmx)=Tiso
	1257
	1258	tsoil(1:ngridmx,1:nsoilmx)=Tiso
	1259
	1260	Tset(1:iip1,1:jjp1,1:llm)=Tiso
[135]	1261	flagtset=.true.
[699]	1262
	1263	ucov(1:iip1,1:jjp1,1:llm)=0
	1264	vcov(1:iip1,1:jjm,1:llm)=0
	1265	q2(1:ngridmx,1:nlayermx+1)=0
[135]	1266
[253]	1267	c radequi : Radiative equilibrium profile of temperatures and no winds
	1268	c --------------------------------------------------------------------
[699]	1269	else if (trim(modif) .eq. 'radequi') then
[135]	1270
[253]	1271	write(,)'radiative equilibrium temperature profile'
	1272
	1273	do ig=1, ngridmx
	1274	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
	1275	& 10.0cos(latfi(ig))cos(latfi(ig))
	1276	tsurf(ig) = MAX(220.0,teque)
	1277	end do
	1278	do l=2,nsoilmx
	1279	do ig=1, ngridmx
	1280	tsoil(ig,l) = tsurf(ig)
	1281	end do
	1282	end do
	1283	DO j=1,jjp1
	1284	DO i=1,iim
	1285	Do l=1,llm
	1286	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
	1287	& 10.0cos(rlatu(j))cos(rlatu(j))
	1288	Tset(i,j,l)=MAX(220.0,teque)
	1289	end do
	1290	end do
	1291	end do
	1292	flagtset=.true.
[699]	1293	ucov(1:iip1,1:jjp1,1:llm)=0
	1294	vcov(1:iip1,1:jjm,1:llm)=0
	1295	q2(1:ngridmx,1:nlayermx+1)=0
[253]	1296
[135]	1297	c coldstart : T set 1K above CO2 frost point and no winds
	1298	c ------------------------------------------------
[699]	1299	else if (trim(modif) .eq. 'coldstart') then
[135]	1300
	1301	write(,)'set temperature of the atmosphere,'
	1302	&,'surface and subsurface how many degrees above CO2 frost point?'
	1303	204 read(,,iostat=ierr) Tabove
	1304	if(ierr.ne.0) goto 204
	1305
	1306	DO j=1,jjp1
	1307	DO i=1,iim
	1308	ig=1+(j-2)*iim +i
	1309	if(j.eq.1) ig=1
	1310	if(j.eq.jjp1) ig=ngridmx
	1311	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
	1312	END DO
	1313	END DO
	1314	do l=1,nsoilmx
	1315	do ig=1, ngridmx
	1316	tsoil(ig,l) = tsurf(ig)
	1317	end do
	1318	end do
	1319	DO j=1,jjp1
	1320	DO i=1,iim
	1321	Do l=1,llm
	1322	pp = aps(l) +bps(l)*ps(i,j)
	1323	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
	1324	end do
	1325	end do
	1326	end do
	1327
	1328	flagtset=.true.
[699]	1329	ucov(1:iip1,1:jjp1,1:llm)=0
	1330	vcov(1:iip1,1:jjm,1:llm)=0
	1331	q2(1:ngridmx,1:nlayermx+1)=0
[135]	1332
	1333
	1334	c co2ice=0 : remove CO2 polar ice caps'
	1335	c ------------------------------------------------
[699]	1336	else if (trim(modif) .eq. 'co2ice=0') then
[135]	1337	! check that there is indeed a co2_ice tracer ...
	1338	if (igcm_co2_ice.ne.0) then
	1339	do ig=1,ngridmx
	1340	!co2ice(ig)=0
	1341	qsurf(ig,igcm_co2_ice)=0
	1342	emis(ig)=emis(ngridmx/2)
	1343	end do
	1344	else
	1345	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
	1346	endif
	1347
	1348	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	1349	! ----------------------------------------------------------------------
	1350
[699]	1351	else if (trim(modif).eq.'therm_ini_s') then
[135]	1352	! write(,)"surfithfi(1):",surfithfi(1)
	1353	do isoil=1,nsoilmx
	1354	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1355	enddo
	1356	write(,)'OK: Soil thermal inertia has been reset to referenc
	1357	&e surface values'
	1358	! write(,)"inertiedat(1,1):",inertiedat(1,1)
	1359	ithfi(:,:)=inertiedat(:,:)
	1360	! recast ithfi() onto ith()
	1361	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1362	! Check:
	1363	! do i=1,iip1
	1364	! do j=1,jjp1
	1365	! do isoil=1,nsoilmx
	1366	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1367	! enddo
	1368	! enddo
	1369	! enddo
	1370
	1371
[699]	1372	else
	1373	write(,) ' Unknown (misspelled?) option!!!'
	1374	end if ! of if (trim(modif) .eq. '...') elseif ...
[253]	1375
[135]	1376
	1377
	1378
	1379
	1380
	1381
	1382
	1383
	1384
[253]	1385
[135]	1386	enddo ! of do ! infinite loop on liste of changes
	1387
	1388	999 continue
	1389
	1390
	1391	c=======================================================================
[253]	1392	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
	1393	c=======================================================================
	1394	DO ig=1, ngridmx
	1395	totalfrac(ig)=0.5
	1396	DO l=1,llm
	1397	cloudfrac(ig,l)=0.5
	1398	ENDDO
[588]	1399	! Ehouarn, march 2012: also add some initialisation for hice
	1400	hice(ig)=0.0
[253]	1401	ENDDO
	1402
	1403	c========================================================
	1404
	1405	! DO ig=1,ngridmx
	1406	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
	1407	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
	1408	! hice(ig)=min(hice(ig),1.0)
	1409	! ENDIF
	1410	! ENDDO
	1411
	1412	c=======================================================================
[135]	1413	c Correct pressure on the new grid (menu 0)
	1414	c=======================================================================
	1415
[588]	1416
[135]	1417	if ((choix_1.eq.0).and.(.not.flagps0)) then
	1418	r = 1000.*8.31/mugaz
	1419
	1420	do j=1,jjp1
	1421	do i=1,iip1
[588]	1422	ps(i,j) = ps(i,j) *
	1423	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1424	. (t(i,j,1) * r))
[253]	1425	end do
	1426	end do
	1427
[588]	1428	c periodicite de ps en longitude
[253]	1429	do j=1,jjp1
[588]	1430	ps(1,j) = ps(iip1,j)
[135]	1431	end do
	1432	end if
	1433
[588]	1434
[135]	1435	c=======================================================================
	1436	c=======================================================================
	1437
	1438	c=======================================================================
	1439	c Initialisation de la physique / ecriture de newstartfi :
	1440	c=======================================================================
	1441
[588]	1442
[135]	1443	CALL inifilr
	1444	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	1445
	1446	c-----------------------------------------------------------------------
	1447	c Initialisation pks:
	1448	c-----------------------------------------------------------------------
	1449
	1450	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
	1451	! Calcul de la temperature potentielle teta
	1452
	1453	if (flagtset) then
	1454	DO l=1,llm
	1455	DO j=1,jjp1
	1456	DO i=1,iim
	1457	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
	1458	ENDDO
	1459	teta (iip1,j,l)= teta (1,j,l)
	1460	ENDDO
	1461	ENDDO
	1462	else if (choix_1.eq.0) then
	1463	DO l=1,llm
	1464	DO j=1,jjp1
	1465	DO i=1,iim
	1466	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1467	ENDDO
	1468	teta (iip1,j,l)= teta (1,j,l)
	1469	ENDDO
	1470	ENDDO
	1471	endif
	1472
	1473	C Calcul intermediaire
	1474	c
	1475	if (choix_1.eq.0) then
	1476	CALL massdair( p3d, masse )
	1477	c
	1478	print *,' ALPHAX ',alphax
	1479
	1480	DO l = 1, llm
	1481	DO i = 1, iim
	1482	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1483	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1484	ENDDO
	1485	xpn = SUM(xppn)/apoln
	1486	xps = SUM(xpps)/apols
	1487	DO i = 1, iip1
	1488	masse( i , 1 , l ) = xpn
	1489	masse( i , jjp1 , l ) = xps
	1490	ENDDO
	1491	ENDDO
	1492	endif
	1493	phis(iip1,:) = phis(1,:)
	1494
	1495	itau=0
	1496	if (choix_1.eq.0) then
	1497	day_ini=int(date)
	1498	endif
	1499	c
	1500	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	1501
	1502	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	1503	* phi,w, pbaru,pbarv,day_ini+time )
	1504
[588]	1505
[135]	1506	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
	1507	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
	1508	C
	1509	C Ecriture etat initial physique
	1510	C
	1511
[253]	1512	! do ig=1,ngridmx
	1513	! print*,'surface ice in newstart=',qsurf(ig,igcm_h2o_ice)
	1514	! print*,'surface liq in newstart=',qsurf(ig,igcm_h2o_vap)
	1515	! enddo
	1516	! stop
	1517
	1518
[787]	1519	call physdem1(ngridmx,"restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
[135]	1520	. dtphys,float(day_ini),
	1521	. time,tsurf,tsoil,emis,q2,qsurf,
[253]	1522	. airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
[787]	1523	. cloudfrac,totalfrac,hice,tnom)
[135]	1524
	1525	c=======================================================================
	1526	c Formats
	1527	c=======================================================================
	1528
	1529	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	1530	*rrage est differente de la valeur parametree iim =',i4//)
	1531	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	1532	*rrage est differente de la valeur parametree jjm =',i4//)
	1533	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	1534	*rage est differente de la valeur parametree llm =',i4//)
	1535
	1536	end
	1537

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

Download in other formats:

Original Format