Context Navigation

newstart.F @ 830

Last change on this file since 830 was 787, checked in by aslmd, 13 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

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

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

Download in other formats:

Original Format