Context Navigation

newstart.F @ 1130

Last change on this file since 1130 was 1130, checked in by emillour, 11 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
dyn3d/control.h is now module control_mod.F90
rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

File size: 56.7 KB

Rev	Line
[38]	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
[164]	17	! to use 'getin'
[1036]	18	use ioipsl_getincom, only: getin
[1130]	19	use infotrac, only: iniadvtrac, nqtot, tname
	20	use tracer_mod, only: noms, igcm_dust_number, igcm_dust_mass,
	21	& igcm_ccn_number, igcm_ccn_mass,
	22	& igcm_h2o_vap, igcm_h2o_ice
[1047]	23	use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
	24	& albedodat, z0_default
	25	use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx
[1130]	26	use control_mod, only: day_step, iphysiq, anneeref
	27	use phyredem, only: physdem0, physdem1
	28	use iostart, only: open_startphy
	29	use comgeomphy, only: initcomgeomphy
[38]	30	implicit none
	31
	32	#include "dimensions.h"
[1047]	33	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
[38]	34	#include "dimphys.h"
[1047]	35	!#include "surfdat.h"
	36	!#include "comsoil.h"
	37	!#include "dimradmars.h"
	38	!#include "yomaer.h"
[38]	39	#include "planete.h"
	40	#include "paramet.h"
	41	#include "comconst.h"
	42	#include "comvert.h"
	43	#include "comgeom2.h"
[1130]	44	!#include "control.h"
[38]	45	#include "logic.h"
	46	#include "description.h"
	47	#include "ener.h"
	48	#include "temps.h"
[1130]	49	!#include "lmdstd.h"
[38]	50	#include "comdissnew.h"
	51	#include "clesph0.h"
	52	#include "serre.h"
	53	#include "netcdf.inc"
[1036]	54	!#include"advtrac.h"
	55	!#include"tracer.h"
[164]	56	#include "datafile.h"
[38]	57	c=======================================================================
	58	c Declarations
	59	c=======================================================================
	60
	61	c Variables dimension du fichier "start_archive"
	62	c------------------------------------
	63	CHARACTER relief*3
	64
	65	c et autres:
	66	c----------
	67
	68	c Variables pour les lectures NetCDF des fichiers "start_archive"
	69	c--------------------------------------------------
	70	INTEGER nid_dyn, nid_fi,nid,nvarid
	71	! INTEGER size
	72	INTEGER length
	73	parameter (length = 100)
	74	INTEGER tab0
	75	INTEGER NB_ETATMAX
	76	parameter (NB_ETATMAX = 100)
	77
	78	REAL date
	79	REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
	80
	81	c Variable histoire
	82	c------------------
	83	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	84	REAL phis(iip1,jjp1)
[1036]	85	REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes
[38]	86
	87	c autre variables dynamique nouvelle grille
	88	c------------------------------------------
	89	REAL pks(iip1,jjp1)
	90	REAL w(iip1,jjp1,llm+1)
	91	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
	92	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	93	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
	94	REAL phi(iip1,jjp1,llm)
	95
	96	integer klatdat,klongdat
	97	PARAMETER (klatdat=180,klongdat=360)
	98
	99	c Physique sur grille scalaire
	100	c----------------------------
	101	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
	102	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
[224]	103	real z0S(iip1,jjp1)
[38]	104
	105	c variable physique
	106	c------------------
	107	REAL tsurf(ngridmx) ! surface temperature
	108	REAL tsoil(ngridmx,nsoilmx) ! soil temperature
	109	REAL co2ice(ngridmx) ! CO2 ice layer
	110	REAL emis(ngridmx) ! surface emissivity
[1036]	111	REAL,ALLOCATABLE :: qsurf(:,:)
[38]	112	REAL q2(ngridmx,nlayermx+1)
	113	! REAL rnaturfi(ngridmx)
	114	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
	115	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
	116	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
	117	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
	118
	119	INTEGER i,j,l,isoil,ig,idum
	120	real mugaz ! molar mass of the atmosphere
	121
	122	EXTERNAL iniconst,geopot,inigeom
	123	integer ierr !, nbetat
	124	integer ismin
	125	external ismin
	126
	127	c Variables on the new grid along scalar points
	128	c------------------------------------------------------
	129	! REAL p(iip1,jjp1)
	130	REAL t(iip1,jjp1,llm)
	131	real phisold_newgrid(iip1,jjp1)
	132	REAL :: teta(iip1, jjp1, llm)
	133	REAL :: pk(iip1,jjp1,llm)
	134	REAL :: pkf(iip1,jjp1,llm)
	135	REAL :: ps(iip1, jjp1)
	136	REAL :: masse(iip1,jjp1,llm)
	137	REAL :: xpn,xps,xppn(iim),xpps(iim)
	138	REAL :: p3d(iip1, jjp1, llm+1)
	139	REAL :: beta(iip1,jjp1,llm)
	140	! REAL dteta(ip1jmp1,llm)
	141
	142	c Variable de l'ancienne grille
	143	c------------------------------
	144	real time
	145	real tab_cntrl(100)
	146	real tab_cntrl_bis(100)
	147
	148	c variables diverses
	149	c-------------------
[1130]	150	real choix_1 ! ==0 : read start_archive file ; ==1: read start files
[38]	151	character*80 fichnom
[618]	152	integer Lmodif,iq
	153	integer flagthermo, flagh2o
[38]	154	character modif*20
	155	real z_reel(iip1,jjp1)
	156	real tsud,albsud,alb_bb,ith_bb,Tiso
	157	real ptoto,pcap,patm,airetot,ptotn,patmn
	158	! real ssum
	159	character*1 yes
	160	logical :: flagiso=.false. , flagps0=.false.
	161	real val, val2, val3 ! to store temporary variables
	162	real :: iceith=2000 ! thermal inertia of subterranean ice
	163	real :: iceithN,iceithS ! values of thermal inertias in N & S hemispheres
	164	integer iref,jref
	165
	166	INTEGER :: itau
	167
	168	INTEGER :: nq,numvanle
[563]	169	character(len=50) :: txt ! to store some text
[38]	170	integer :: count
[563]	171	real :: profile(llm+1) ! to store an atmospheric profile + surface value
[38]	172
	173	! MONS data:
	174	real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
	175	real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
	176	! coefficient to apply to convert d21 to 'true' depth (m)
	177	real :: MONS_coeff
	178	real :: MONS_coeffS ! coeff for southern hemisphere
	179	real :: MONS_coeffN ! coeff for northern hemisphere
	180	! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
	181
	182	c sortie visu pour les champs dynamiques
	183	c---------------------------------------
	184	! INTEGER :: visuid
	185	! real :: time_step,t_ops,t_wrt
	186	! CHARACTER*80 :: visu_file
	187
	188	cpp = 744.499 ! for Mars, instead of 1004.70885 (Earth)
	189	preff = 610. ! for Mars, instead of 101325. (Earth)
	190	pa= 20 ! for Mars, instead of 500 (Earth)
	191
[1130]	192	! initialize "serial/parallel" related stuff
	193	CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
	194	call initcomgeomphy
	195
[1036]	196	! Load tracer number and names:
	197	call iniadvtrac(nqtot,numvanle)
	198	! allocate arrays
	199	allocate(q(iip1,jjp1,llm,nqtot))
	200	allocate(qsurf(ngridmx,nqtot))
	201
[38]	202	c=======================================================================
	203	c Choice of the start file(s) to use
	204	c=======================================================================
	205
	206	write(,) 'From which kind of files do you want to create new',
	207	. 'start and startfi files'
	208	write(,) ' 0 - from a file start_archive'
	209	write(,) ' 1 - from files start and startfi'
	210
	211	c-----------------------------------------------------------------------
	212	c Open file(s) to modify (start or start_archive)
	213	c-----------------------------------------------------------------------
	214
	215	DO
	216	read(,,iostat=ierr) choix_1
	217	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	218	ENDDO
	219
	220	c Open start_archive
	221	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	222	if (choix_1.eq.0) then
	223
	224	write(,) 'Creating start files from:'
	225	write(,) './start_archive.nc'
	226	write(,)
	227	fichnom = 'start_archive.nc'
	228	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	229	IF (ierr.NE.NF_NOERR) THEN
	230	write(6,*)' Problem opening file:',fichnom
	231	write(6,*)' ierr = ', ierr
	232	CALL ABORT
	233	ENDIF
	234	tab0 = 50
	235	Lmodif = 1
	236
	237	c OR open start and startfi files
	238	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	239	else
	240	write(,) 'Creating start files from:'
	241	write(,) './start.nc and ./startfi.nc'
	242	write(,)
	243	fichnom = 'start.nc'
	244	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	245	IF (ierr.NE.NF_NOERR) THEN
	246	write(6,*)' Problem opening file:',fichnom
	247	write(6,*)' ierr = ', ierr
	248	CALL ABORT
	249	ENDIF
	250
	251	fichnom = 'startfi.nc'
	252	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	253	IF (ierr.NE.NF_NOERR) THEN
	254	write(6,*)' Problem opening file:',fichnom
	255	write(6,*)' ierr = ', ierr
	256	CALL ABORT
	257	ENDIF
	258
	259	tab0 = 0
	260	Lmodif = 0
	261
	262	endif
	263
	264	c-----------------------------------------------------------------------
	265	c Lecture du tableau des parametres du run (pour la dynamique)
	266	c-----------------------------------------------------------------------
	267
	268	if (choix_1.eq.0) then
	269
	270	write(,) 'reading tab_cntrl START_ARCHIVE'
	271	c
	272	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	273	#ifdef NC_DOUBLE
	274	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	275	#else
	276	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	277	#endif
	278	c
	279	else if (choix_1.eq.1) then
	280
	281	write(,) 'reading tab_cntrl START'
	282	c
	283	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	284	#ifdef NC_DOUBLE
	285	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	286	#else
	287	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	288	#endif
	289	c
	290	write(,) 'reading tab_cntrl STARTFI'
	291	c
	292	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	293	#ifdef NC_DOUBLE
	294	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	295	#else
	296	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	297	#endif
	298	c
	299	do i=1,50
	300	tab_cntrl(i+50)=tab_cntrl_bis(i)
	301	enddo
	302	write(,) 'printing tab_cntrl', tab_cntrl
	303	do i=1,100
	304	write(,) i,tab_cntrl(i)
	305	enddo
	306
	307	endif
	308	c-----------------------------------------------------------------------
	309	c Initialisation des constantes dynamique
	310	c-----------------------------------------------------------------------
	311
	312	kappa = tab_cntrl(9)
	313	etot0 = tab_cntrl(12)
	314	ptot0 = tab_cntrl(13)
	315	ztot0 = tab_cntrl(14)
	316	stot0 = tab_cntrl(15)
	317	ang0 = tab_cntrl(16)
	318	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	319	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	320
	321	c-----------------------------------------------------------------------
	322	c Lecture du tab_cntrl et initialisation des constantes physiques
	323	c - pour start: Lmodif = 0 => pas de modifications possibles
	324	c (modif dans le tabfi de readfi + loin)
	325	c - pour start_archive: Lmodif = 1 => modifications possibles
	326	c-----------------------------------------------------------------------
	327	if (choix_1.eq.0) then
[1130]	328	! tabfi requires that input file be first opened by open_startphy(fichnom)
	329	fichnom = 'start_archive.nc'
	330	call open_startphy(fichnom)
[38]	331	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
	332	. p_omeg,p_g,p_mugaz,p_daysec,time)
	333	else if (choix_1.eq.1) then
[1130]	334	fichnom = 'startfi.nc'
	335	call open_startphy(fichnom)
[38]	336	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
	337	. p_omeg,p_g,p_mugaz,p_daysec,time)
	338	endif
	339
	340	rad = p_rad
	341	omeg = p_omeg
	342	g = p_g
	343	mugaz = p_mugaz
	344	daysec = p_daysec
	345	! write(,) 'aire',aire
	346
	347
	348	c=======================================================================
	349	c INITIALISATIONS DIVERSES
	350	c=======================================================================
	351
	352	day_step=180 !?! Note: day_step is a common in "control.h"
	353	CALL defrun_new( 99, .TRUE. )
	354	CALL iniconst
	355	CALL inigeom
	356	idum=-1
	357	idum=0
	358
	359	c Initialisation coordonnees /aires
	360	c -------------------------------
	361	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
	362	! rlatu() and rlonv() are given in radians
	363	latfi(1)=rlatu(1)
	364	lonfi(1)=0.
	365	DO j=2,jjm
	366	DO i=1,iim
	367	latfi((j-2)*iim+1+i)=rlatu(j)
	368	lonfi((j-2)*iim+1+i)=rlonv(i)
	369	ENDDO
	370	ENDDO
	371	latfi(ngridmx)=rlatu(jjp1)
	372	lonfi(ngridmx)=0.
[697]	373
	374	! build airefi(), mesh area on physics grid
[38]	375	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
[697]	376	! Poles are single points on physics grid
	377	airefi(1)=airefi(1)*iim
	378	airefi(ngridmx)=airefi(ngridmx)*iim
[38]	379
[321]	380	! also initialize various physics flags/settings which might be needed
	381	! (for instance initracer needs to know about some flags, and/or
	382	! 'datafile' path may be changed by user)
[1036]	383	call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
[321]	384	& latfi,lonfi,airefi,rad,g,r,cpp)
	385
[38]	386	c=======================================================================
	387	c lecture topographie, albedo, inertie thermique, relief sous-maille
	388	c=======================================================================
	389
	390	if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier
	391	! surface.dat dans le cas des start
	392
	393	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	394	c write(,)
	395	c write(,) 'choix du relief (mola,pla)'
	396	c write(,) '(Topographie MGS MOLA, plat)'
	397	c read(*,fmt='(a3)') relief
	398	relief="mola"
	399	c enddo
	400
[224]	401	CALL datareadnc(relief,phis,alb,surfith,z0S,
	402	& zmeaS,zstdS,zsigS,zgamS,ztheS)
[38]	403
	404	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	405	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	406	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	407	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
[224]	408	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,z0S,z0)
[38]	409	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	410	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	411	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	412	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	413	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
	414
	415	endif ! of if (choix_1.ne.1)
	416
	417
	418	c=======================================================================
	419	c Lecture des fichiers (start ou start_archive)
	420	c=======================================================================
	421
	422	if (choix_1.eq.0) then
	423
	424	write(,) 'Reading file START_ARCHIVE'
[1047]	425	CALL lect_start_archive(ngridmx,llm,nqtot,
	426	& date,tsurf,tsoil,emis,q2,
	427	& t,ucov,vcov,ps,co2ice,teta,phisold_newgrid,q,qsurf,
[38]	428	& surfith,nid)
	429	write(,) "OK, read start_archive file"
	430	! copy soil thermal inertia
	431	ithfi(:,:)=inertiedat(:,:)
	432
	433	ierr= NF_CLOSE(nid)
	434
	435	else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui
	436	! permet de changer les valeurs du
	437	! tab_cntrl Lmodif=1
	438	tab0=0
	439	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
	440	write(,) 'Reading file START'
	441	fichnom = 'start.nc'
[1036]	442	CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
[38]	443	. ps,phis,time)
	444
	445	write(,) 'Reading file STARTFI'
	446	fichnom = 'startfi.nc'
[1047]	447	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
[38]	448	. day_ini,time,
	449	. tsurf,tsoil,emis,q2,qsurf,co2ice)
	450
	451	! copy albedo and soil thermal inertia
	452	do i=1,ngridmx
	453	albfi(i) = albedodat(i)
	454	do j=1,nsoilmx
	455	ithfi(i,j) = inertiedat(i,j)
	456	enddo
	457	! build a surfithfi(:) using 1st layer of ithfi(:), which might
	458	! be neede later on if reinitializing soil thermal inertia
	459	surfithfi(i)=ithfi(i,1)
	460	enddo
	461
	462	else
	463	CALL exit(1)
	464	endif
	465
	466	dtvr = daysec/REAL(day_step)
	467	dtphys = dtvr * REAL(iphysiq)
	468
	469	c=======================================================================
	470	c
	471	c=======================================================================
	472	! If tracer names follow 'old' convention (q01, q02, ...) then
	473	! rename them
	474	count=0
[1036]	475	do iq=1,nqtot
[38]	476	txt=" "
	477	write(txt,'(a1,i2.2)') 'q',iq
[1130]	478	if (txt.eq.tname(iq)) then
[38]	479	count=count+1
	480	endif
[1036]	481	enddo ! of do iq=1,nqtot
[38]	482
[618]	483	! initialize tracer names noms(:) and indexes (igcm_co2, igcm_h2o_vap, ...)
[1036]	484	call initracer(ngridmx,nqtot,qsurf,co2ice)
[618]	485
[1036]	486	if (count.eq.nqtot) then
[38]	487	write(,) 'Newstart: updating tracer names'
[1130]	488	! copy noms(:) to tname(:) to have matching tracer names in physics
[618]	489	! and dynamics
[1130]	490	tname(1:nqtot)=noms(1:nqtot)
[38]	491	endif
	492
	493	c=======================================================================
	494	c
	495	c=======================================================================
	496
	497	do ! infinite loop on list of changes
	498
	499	write(,)
	500	write(,)
	501	write(,) 'List of possible changes :'
[618]	502	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~'
[38]	503	write(,)
[618]	504	write(,) 'flat : no topography ("aquaplanet")'
	505	write(,) 'bilball : uniform albedo and thermal inertia'
	506	write(,) 'z0 : set a uniform surface roughness length'
	507	write(,) 'coldspole : cold subsurface and high albedo at
	508	$ S.Pole'
	509	write(,) 'qname : change tracer name'
	510	write(,) 'q=0 : ALL tracer =zero'
	511	write(,) 'q=x : give a specific uniform value to one
	512	$ tracer'
	513	write(,) 'q=profile : specify a profile for a tracer'
[1130]	514	write(,) 'freedust : rescale dust to a true value'
[618]	515	write(,) 'ini_q : tracers initialization for chemistry
	516	$ and water vapour'
	517	write(,) 'ini_q-h2o : tracers initialization for chemistry
	518	$ only'
	519	write(,) 'ini_h2osurf : reinitialize surface water ice '
	520	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
	521	write(,) 'watercapn : H20 ice on permanent N polar cap '
	522	write(,) 'watercaps : H20 ice on permanent S polar cap '
	523	write(,) 'wetstart : start with a wet atmosphere'
	524	write(,) 'isotherm : Isothermal Temperatures, wind set to
	525	$ zero'
	526	write(,) 'co2ice=0 : remove CO2 polar cap'
	527	write(,) 'ptot : change total pressure'
	528	write(,) 'therm_ini_s : set soil thermal inertia to reference
	529	$ surface values'
	530	write(,) 'subsoilice_n : put deep underground ice layer in
	531	$ northern hemisphere'
	532	write(,) 'subsoilice_s : put deep underground ice layer in
	533	$ southern hemisphere'
	534	write(,) 'mons_ice : put underground ice layer according
	535	$ to MONS derived data'
[38]	536
	537	write(,)
	538	write(,) 'Change to perform ?'
	539	write(,) ' (enter keyword or return to end)'
	540	write(,)
	541
	542	read(*,fmt='(a20)') modif
	543	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
	544
	545	write(,)
[164]	546	write(,) trim(modif) , ' : '
[38]	547
	548	c 'flat : no topography ("aquaplanet")'
	549	c -------------------------------------
[164]	550	if (trim(modif) .eq. 'flat') then
[38]	551	c set topo to zero
	552	CALL initial0(ip1jmp1,z_reel)
	553	CALL multscal(ip1jmp1,z_reel,g,phis)
	554	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	555	write(,) 'topography set to zero.'
	556	write(,) 'WARNING : the subgrid topography parameters',
	557	& ' were not set to zero ! => set calllott to F'
	558
	559	c Choice for surface pressure
	560	yes=' '
	561	do while ((yes.ne.'y').and.(yes.ne.'n'))
	562	write(,) 'Do you wish to choose homogeneous surface',
	563	& 'pressure (y) or let newstart interpolate ',
	564	& ' the previous field (n)?'
	565	read(*,fmt='(a)') yes
	566	end do
	567	if (yes.eq.'y') then
	568	flagps0=.true.
	569	write(,) 'New value for ps (Pa) ?'
	570	201 read(,,iostat=ierr) patm
	571	if(ierr.ne.0) goto 201
	572	write(,)
	573	write(,) ' new ps everywhere (Pa) = ', patm
	574	write(,)
	575	do j=1,jjp1
	576	do i=1,iip1
	577	ps(i,j)=patm
	578	enddo
	579	enddo
	580	end if
	581
	582	c bilball : albedo, inertie thermique uniforme
	583	c --------------------------------------------
[164]	584	else if (trim(modif) .eq. 'bilball') then
[38]	585	write(,) 'constante albedo and iner.therm:'
	586	write(,) 'New value for albedo (ex: 0.25) ?'
	587	101 read(,,iostat=ierr) alb_bb
	588	if(ierr.ne.0) goto 101
	589	write(,)
	590	write(,) ' uniform albedo (new value):',alb_bb
	591	write(,)
	592
	593	write(,) 'New value for thermal inertia (eg: 247) ?'
	594	102 read(,,iostat=ierr) ith_bb
	595	if(ierr.ne.0) goto 102
	596	write(,) 'uniform thermal inertia (new value):',ith_bb
	597	DO j=1,jjp1
	598	DO i=1,iip1
	599	alb(i,j) = alb_bb ! albedo
	600	do isoil=1,nsoilmx
	601	ith(i,j,isoil) = ith_bb ! thermal inertia
	602	enddo
	603	END DO
	604	END DO
	605	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	606	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	607	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
[224]	608
	609	! also reset surface roughness length to default value
	610	write(,) 'surface roughness length set to:',z0_default,' m'
	611	z0(:)=z0_default
[38]	612
[224]	613	! z0 : set surface roughness length to a constant value
	614	! -----------------------------------------------------
	615	else if (trim(modif) .eq. 'z0') then
	616	write(,) 'set a uniform surface roughness length'
	617	write(,) ' value for z0_default (ex: ',z0_default,')?'
	618	ierr=1
	619	do while (ierr.ne.0)
	620	read(,,iostat=ierr) z0_default
	621	enddo
	622	z0(:)=z0_default
	623
[38]	624	c coldspole : sous-sol de la calotte sud toujours froid
	625	c -----------------------------------------------------
[164]	626	else if (trim(modif) .eq. 'coldspole') then
[38]	627	write(,)'new value for the subsurface temperature',
	628	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
	629	103 read(,,iostat=ierr) tsud
	630	if(ierr.ne.0) goto 103
	631	write(,)
	632	write(,) ' new value of the subsurface temperature:',tsud
	633	c nouvelle temperature sous la calotte permanente
	634	do l=2,nsoilmx
	635	tsoil(ngridmx,l) = tsud
	636	end do
	637
	638
	639	write(,)'new value for the albedo',
	640	& 'of the permanent southern polar cap ? (eg: 0.75)'
	641	104 read(,,iostat=ierr) albsud
	642	if(ierr.ne.0) goto 104
	643	write(,)
	644
	645	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	646	c Option 1: only the albedo of the pole is modified :
	647	albfi(ngridmx)=albsud
	648	write(,) 'ig=',ngridmx,' albedo perennial cap ',
	649	& albfi(ngridmx)
	650
	651	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	652	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
	653	c DO j=1,jjp1
	654	c DO i=1,iip1
	655	c ig=1+(j-2)*iim +i
	656	c if(j.eq.1) ig=1
	657	c if(j.eq.jjp1) ig=ngridmx
	658	c if ((rlatu(j)*180./pi.lt.-84.).and.
	659	c & (rlatu(j)*180./pi.gt.-91.).and.
	660	c & (rlonv(i)*180./pi.gt.-91.).and.
	661	c & (rlonv(i)*180./pi.lt.0.)) then
	662	cc albedo de la calotte permanente fixe a albsud
	663	c alb(i,j)=albsud
	664	c write(,) 'lat=',rlatu(j)*180./pi,
	665	c & ' lon=',rlonv(i)*180./pi
	666	cc fin de la condition sur les limites de la calotte permanente
	667	c end if
	668	c ENDDO
	669	c ENDDO
	670	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	671
	672	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	673
	674
	675	c ptot : Modification of the total pressure: ice + current atmosphere
	676	c -------------------------------------------------------------------
[164]	677	else if (trim(modif) .eq. 'ptot') then
[38]	678
	679	c calcul de la pression totale glace + atm actuelle
	680	patm=0.
	681	airetot=0.
	682	pcap=0.
	683	DO j=1,jjp1
	684	DO i=1,iim
	685	ig=1+(j-2)*iim +i
	686	if(j.eq.1) ig=1
	687	if(j.eq.jjp1) ig=ngridmx
	688	patm = patm + ps(i,j)*aire(i,j)
	689	airetot= airetot + aire(i,j)
	690	pcap = pcap + aire(i,j)co2ice(ig)g
	691	ENDDO
	692	ENDDO
	693	ptoto = pcap + patm
	694
	695	print*,'Current total pressure at surface (co2 ice + atm) ',
	696	& ptoto/airetot
	697
	698	print*,'new value?'
	699	read(,) ptotn
	700	ptotn=ptotn*airetot
	701	patmn=ptotn-pcap
	702	print*,'ptoto,patm,ptotn,patmn'
	703	print*,ptoto,patm,ptotn,patmn
	704	print*,'Mult. factor for pressure (atm only)', patmn/patm
	705	do j=1,jjp1
	706	do i=1,iip1
	707	ps(i,j)=ps(i,j)*patmn/patm
	708	enddo
	709	enddo
	710
	711	c Correction pour la conservation des traceurs
	712	yes=' '
	713	do while ((yes.ne.'y').and.(yes.ne.'n'))
	714	write(,) 'Do you wish to conserve tracer total mass (y)',
	715	& ' or tracer mixing ratio (n) ?'
	716	read(*,fmt='(a)') yes
	717	end do
	718
	719	if (yes.eq.'y') then
	720	write(,) 'OK : conservation of tracer total mass'
[1036]	721	DO iq =1, nqtot
[38]	722	DO l=1,llm
	723	DO j=1,jjp1
	724	DO i=1,iip1
	725	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	726	ENDDO
	727	ENDDO
	728	ENDDO
	729	ENDDO
	730	else
	731	write(,) 'OK : conservation of tracer mixing ratio'
	732	end if
	733
	734	c qname : change tracer name
	735	c --------------------------
	736	else if (trim(modif).eq.'qname') then
	737	yes='y'
	738	do while (yes.eq.'y')
	739	write(,) 'Which tracer name do you want to change ?'
[1036]	740	do iq=1,nqtot
[1130]	741	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
[38]	742	enddo
	743	write(*,'(a35,i3)')
[1036]	744	& '(enter tracer number; between 1 and ',nqtot
[38]	745	write(,)' or any other value to quit this option)'
	746	read(,) iq
[1036]	747	if ((iq.ge.1).and.(iq.le.nqtot)) then
[1130]	748	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
[38]	749	read(,) txt
[1130]	750	tname(iq)=txt
[38]	751	write(,)'Do you want to change another tracer name (y/n)?'
	752	read(*,'(a)') yes
	753	else
	754	! inapropiate value of iq; quit this option
	755	yes='n'
[1036]	756	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
[38]	757	enddo ! of do while (yes.ne.'y')
	758
	759	c q=0 : set tracers to zero
	760	c -------------------------
[164]	761	else if (trim(modif) .eq. 'q=0') then
[38]	762	c mise a 0 des q (traceurs)
	763	write(,) 'Tracers set to 0 (1.E-30 in fact)'
[1036]	764	DO iq =1, nqtot
[38]	765	DO l=1,llm
	766	DO j=1,jjp1
	767	DO i=1,iip1
	768	q(i,j,l,iq)=1.e-30
	769	ENDDO
	770	ENDDO
	771	ENDDO
	772	ENDDO
	773
	774	c set surface tracers to zero
[1036]	775	DO iq =1, nqtot
[38]	776	DO ig=1,ngridmx
	777	qsurf(ig,iq)=0.
	778	ENDDO
	779	ENDDO
	780
	781	c q=x : initialise tracer manually
	782	c --------------------------------
[164]	783	else if (trim(modif) .eq. 'q=x') then
[38]	784	write(,) 'Which tracer do you want to modify ?'
[1036]	785	do iq=1,nqtot
[1130]	786	write(,)iq,' : ',trim(tname(iq))
[38]	787	enddo
[1036]	788	write(,) '(choose between 1 and ',nqtot,')'
[38]	789	read(,) iq
[1036]	790	if ((iq.lt.1).or.(iq.gt.nqtot)) then
[171]	791	! wrong value for iq, go back to menu
	792	write(,) "wrong input value:",iq
	793	cycle
	794	endif
[1130]	795	write(,)'mixing ratio of tracer ',trim(tname(iq)),
[38]	796	& ' ? (kg/kg)'
	797	read(,) val
	798	DO l=1,llm
	799	DO j=1,jjp1
	800	DO i=1,iip1
	801	q(i,j,l,iq)=val
	802	ENDDO
	803	ENDDO
	804	ENDDO
[1130]	805	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
[38]	806	& ' ? (kg/m2)'
	807	read(,) val
	808	DO ig=1,ngridmx
	809	qsurf(ig,iq)=val
	810	ENDDO
	811
[563]	812	c q=profile : initialize tracer with a given profile
	813	c --------------------------------------------------
	814	else if (trim(modif) .eq. 'q=profile') then
	815	write(,) 'Tracer profile will be sought in ASCII file'
	816	write(,) "'profile_tracer' where 'tracer' is tracer name"
	817	write(,) "(one value per line in file; starting with"
	818	write(,) "surface value, the 1st atmospheric layer"
	819	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	820	write(,) 'Which tracer do you want to set?'
[1036]	821	do iq=1,nqtot
[1130]	822	write(,)iq,' : ',trim(tname(iq))
[563]	823	enddo
[1036]	824	write(,) '(choose between 1 and ',nqtot,')'
[563]	825	read(,) iq
[1036]	826	if ((iq.lt.1).or.(iq.gt.nqtot)) then
[563]	827	! wrong value for iq, go back to menu
	828	write(,) "wrong input value:",iq
	829	cycle
	830	endif
	831	! look for input file 'profile_tracer'
[1130]	832	txt="profile_"//trim(tname(iq))
[563]	833	open(41,file=trim(txt),status='old',form='formatted',
	834	& iostat=ierr)
	835	if (ierr.eq.0) then
	836	! OK, found file 'profile_...', load the profile
	837	do l=1,llm+1
	838	read(41,*,iostat=ierr) profile(l)
	839	if (ierr.ne.0) then ! something went wrong
	840	exit ! quit loop
	841	endif
	842	enddo
	843	if (ierr.eq.0) then
	844	! initialize tracer values
	845	qsurf(:,iq)=profile(1)
	846	do l=1,llm
	847	q(:,:,l,iq)=profile(l+1)
	848	enddo
[1130]	849	write(,)'OK, tracer ',trim(tname(iq)),
	850	& ' initialized ','using values from file ',trim(txt)
[563]	851	else
	852	write(,)'problem reading file ',trim(txt),' !'
[1130]	853	write(,)'No modifications to tracer ',trim(tname(iq))
[563]	854	endif
	855	else
	856	write(,)'Could not find file ',trim(txt),' !'
[1130]	857	write(,)'No modifications to tracer ',trim(tname(iq))
[563]	858	endif
	859
[1130]	860	c q=profile : initialize tracer with a given profile
	861	c --------------------------------------------------
	862	else if (trim(modif) .eq. 'freedust') then
	863	do l=1,llm
	864	do j=1,jjp1
	865	do i=1,iip1
	866	if (igcm_dust_number .ne. 0)
	867	& q(i,j,l,igcm_dust_number)=
	868	& q(i,j,l,igcm_dust_number) * 1e-3 ! grosso modo
	869	if (igcm_dust_mass .ne. 0)
	870	& q(i,j,l,igcm_dust_mass)=
	871	& q(i,j,l,igcm_dust_mass) * 1e-3 ! grosso modo
	872	if (igcm_ccn_number .ne. 0)
	873	& q(i,j,l,igcm_ccn_number)=
	874	& q(i,j,l,igcm_ccn_number) * 1e-3 ! grosso modo
	875	if (igcm_ccn_mass .ne. 0)
	876	& q(i,j,l,igcm_ccn_mass)=
	877	& q(i,j,l,igcm_ccn_mass) * 1e-3 ! grosso modo
	878	end do
	879	end do
	880	end do
[563]	881
[1130]	882	! We want to have the very same value at lon -180 and lon 180
	883	do l = 1,llm
	884	do j = 1,jjp1
	885	do iq = 1,nqtot
	886	q(iip1,j,l,iq) = q(1,j,l,iq)
	887	end do
	888	end do
	889	end do
	890
[38]	891	c ini_q : Initialize tracers for chemistry
	892	c -----------------------------------------------
[164]	893	else if (trim(modif) .eq. 'ini_q') then
[618]	894	flagh2o = 1
	895	flagthermo = 0
	896	yes=' '
[38]	897	c For more than 32 layers, possible to initiate thermosphere only
	898	if (llm.gt.32) then
	899	do while ((yes.ne.'y').and.(yes.ne.'n'))
	900	write(,)'',
	901	& 'initialisation for thermosphere only? (y/n)'
	902	read(*,fmt='(a)') yes
	903	if (yes.eq.'y') then
[618]	904	flagthermo=1
[38]	905	else
[618]	906	flagthermo=0
[38]	907	endif
	908	enddo
	909	endif
	910
[1047]	911	call inichim_newstart(ngridmx, nq, q, qsurf, ps,
	912	& flagh2o, flagthermo)
[664]	913
	914	! We want to have the very same value at lon -180 and lon 180
	915	do l = 1,llm
	916	do j = 1,jjp1
[1036]	917	do iq = 1,nqtot
[664]	918	q(iip1,j,l,iq) = q(1,j,l,iq)
	919	end do
	920	end do
	921	end do
	922
[618]	923	write(,) 'inichim_newstart: chemical species and
	924	$ water vapour initialised'
[38]	925
[664]	926	c ini_q-h2o : as above except for the water vapour tracer
[38]	927	c ------------------------------------------------------
[618]	928	else if (trim(modif) .eq. 'ini_q-h2o') then
	929	flagh2o = 0
	930	flagthermo = 0
	931	yes=' '
[38]	932	! for more than 32 layers, possible to initiate thermosphere only
	933	if(llm.gt.32) then
	934	do while ((yes.ne.'y').and.(yes.ne.'n'))
	935	write(,)'',
	936	& 'initialisation for thermosphere only? (y/n)'
	937	read(*,fmt='(a)') yes
	938	if (yes.eq.'y') then
[618]	939	flagthermo=1
[38]	940	else
[618]	941	flagthermo=0
[38]	942	endif
	943	enddo
	944	endif
[664]	945
[618]	946	call inichim_newstart(q, qsurf, ps, flagh2o, flagthermo)
[664]	947
	948	! We want to have the very same value at lon -180 and lon 180
	949	do l = 1,llm
	950	do j = 1,jjp1
[1036]	951	do iq = 1,nqtot
[664]	952	q(iip1,j,l,iq) = q(1,j,l,iq)
	953	end do
	954	end do
	955	end do
	956
[618]	957	write(,) 'inichim_newstart: chemical species initialised
	958	$ (except water vapour)'
[38]	959
	960	c wetstart : wet atmosphere with a north to south gradient
	961	c --------------------------------------------------------
[164]	962	else if (trim(modif) .eq. 'wetstart') then
[38]	963	! check that there is indeed a water vapor tracer
	964	if (igcm_h2o_vap.eq.0) then
	965	write(,) "No water vapour tracer! Can't use this option"
	966	stop
	967	endif
	968	DO l=1,llm
	969	DO j=1,jjp1
[321]	970	DO i=1,iip1-1
[38]	971	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
	972	ENDDO
[321]	973	! We want to have the very same value at lon -180 and lon 180
	974	q(iip1,j,l,igcm_h2o_vap) = q(1,j,l,igcm_h2o_vap)
[38]	975	ENDDO
	976	ENDDO
	977
	978	write(,) 'Water mass mixing ratio at north pole='
	979	* ,q(1,1,1,igcm_h2o_vap)
	980	write(,) '---------------------------south pole='
	981	* ,q(1,jjp1,1,igcm_h2o_vap)
	982
[480]	983	c ini_h2osurf : reinitialize surface water ice
	984	c --------------------------------------------------
	985	else if (trim(modif) .eq. 'ini_h2osurf') then
	986	write(,)'max surface ice left?(e.g. 0.2 kg/m2=200microns)'
	987	207 read(,,iostat=ierr) val
	988	if(ierr.ne.0) goto 207
	989	write(,)'also set negative values of surf ice to 0'
	990	do ig=1,ngridmx
	991	qsurf(ig,igcm_h2o_ice)=min(val,qsurf(ig,igcm_h2o_ice))
	992	qsurf(ig,igcm_h2o_ice)=max(0.,qsurf(ig,igcm_h2o_ice))
	993	end do
	994
[38]	995	c noglacier : remove tropical water ice (to initialize high res sim)
	996	c --------------------------------------------------
[164]	997	else if (trim(modif) .eq. 'noglacier') then
[38]	998	do ig=1,ngridmx
	999	j=(ig-2)/iim +2
	1000	if(ig.eq.1) j=1
	1001	write(,) 'OK: remove surface ice for \|lat\|<45'
	1002	if (abs(rlatu(j)*180./pi).lt.45.) then
	1003	qsurf(ig,igcm_h2o_ice)=0.
	1004	end if
	1005	end do
	1006
	1007
	1008	c watercapn : H20 ice on permanent northern cap
	1009	c --------------------------------------------------
[164]	1010	else if (trim(modif) .eq. 'watercapn') then
[38]	1011	do ig=1,ngridmx
	1012	j=(ig-2)/iim +2
	1013	if(ig.eq.1) j=1
	1014	if (rlatu(j)*180./pi.gt.80.) then
	1015	qsurf(ig,igcm_h2o_ice)=1.e5
	1016	write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	1017	& qsurf(ig,igcm_h2o_ice)
	1018	write(,)' ==> Ice mesh South boundary (deg)= ',
	1019	& rlatv(j)*180./pi
	1020	end if
	1021	enddo
	1022
	1023	c watercaps : H20 ice on permanent southern cap
	1024	c -------------------------------------------------
[164]	1025	else if (trim(modif) .eq. 'watercaps') then
[38]	1026	do ig=1,ngridmx
	1027	j=(ig-2)/iim +2
	1028	if(ig.eq.1) j=1
	1029	if (rlatu(j)*180./pi.lt.-80.) then
	1030	qsurf(ig,igcm_h2o_ice)=1.e5
	1031	write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	1032	& qsurf(ig,igcm_h2o_ice)
	1033	write(,)' ==> Ice mesh North boundary (deg)= ',
	1034	& rlatv(j-1)*180./pi
	1035	end if
	1036	enddo
	1037
	1038	c isotherm : Isothermal temperatures and no winds
	1039	c ------------------------------------------------
[164]	1040	else if (trim(modif) .eq. 'isotherm') then
[38]	1041
	1042	write(,)'Isothermal temperature of the atmosphere,
	1043	& surface and subsurface'
	1044	write(,) 'Value of this temperature ? :'
	1045	203 read(,,iostat=ierr) Tiso
	1046	if(ierr.ne.0) goto 203
	1047
	1048	do ig=1, ngridmx
	1049	tsurf(ig) = Tiso
	1050	end do
	1051	do l=2,nsoilmx
	1052	do ig=1, ngridmx
	1053	tsoil(ig,l) = Tiso
	1054	end do
	1055	end do
	1056	flagiso=.true.
	1057	call initial0(llm*ip1jmp1,ucov)
	1058	call initial0(llm*ip1jm,vcov)
	1059	call initial0(ngridmx*(llm+1),q2)
	1060
	1061	c co2ice=0 : remove CO2 polar ice caps'
	1062	c ------------------------------------------------
[164]	1063	else if (trim(modif) .eq. 'co2ice=0') then
[38]	1064	do ig=1,ngridmx
	1065	co2ice(ig)=0
	1066	emis(ig)=emis(ngridmx/2)
	1067	end do
	1068
	1069	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	1070	! ----------------------------------------------------------------------
	1071
[164]	1072	else if (trim(modif).eq.'therm_ini_s') then
[38]	1073	! write(,)"surfithfi(1):",surfithfi(1)
	1074	do isoil=1,nsoilmx
	1075	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1076	enddo
	1077	write(,)'OK: Soil thermal inertia has been reset to referenc
	1078	&e surface values'
	1079	! write(,)"inertiedat(1,1):",inertiedat(1,1)
	1080	ithfi(:,:)=inertiedat(:,:)
	1081	! recast ithfi() onto ith()
	1082	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1083	! Check:
	1084	! do i=1,iip1
	1085	! do j=1,jjp1
	1086	! do isoil=1,nsoilmx
	1087	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1088	! enddo
	1089	! enddo
	1090	! enddo
	1091
	1092	! subsoilice_n: Put deep ice layer in northern hemisphere soil
	1093	! ------------------------------------------------------------
	1094
[164]	1095	else if (trim(modif).eq.'subsoilice_n') then
[38]	1096
	1097	write(,)'From which latitude (in deg.), up to the north pole,
	1098	&should we put subterranean ice?'
	1099	ierr=1
	1100	do while (ierr.ne.0)
	1101	read(,,iostat=ierr) val
	1102	if (ierr.eq.0) then ! got a value
	1103	! do a sanity check
	1104	if((val.lt.0.).or.(val.gt.90)) then
	1105	write(,)'Latitude should be between 0 and 90 deg. !!!'
	1106	ierr=1
	1107	else ! find corresponding jref (nearest latitude)
	1108	! note: rlatu(:) contains decreasing values of latitude
	1109	! starting from PI/2 to -PI/2
	1110	do j=1,jjp1
	1111	if ((rlatu(j)*180./pi.ge.val).and.
	1112	& (rlatu(j+1)*180./pi.le.val)) then
	1113	! find which grid point is nearest to val:
	1114	if (abs(rlatu(j)*180./pi-val).le.
	1115	& abs((rlatu(j+1)*180./pi-val))) then
	1116	jref=j
	1117	else
	1118	jref=j+1
	1119	endif
	1120
	1121	write(,)'Will use nearest grid latitude which is:',
	1122	& rlatu(jref)*180./pi
	1123	endif
	1124	enddo ! of do j=1,jjp1
	1125	endif ! of if((val.lt.0.).or.(val.gt.90))
	1126	endif !of if (ierr.eq.0)
	1127	enddo ! of do while
	1128
	1129	! Build layers() (as in soil_settings.F)
	1130	val2=sqrt(mlayer(0)*mlayer(1))
	1131	val3=mlayer(1)/mlayer(0)
	1132	do isoil=1,nsoilmx
	1133	layer(isoil)=val2(val3*(isoil-1))
	1134	enddo
	1135
	1136	write(,)'At which depth (in m.) does the ice layer begin?'
	1137	write(,)'(currently, the deepest soil layer extends down to:'
	1138	& ,layer(nsoilmx),')'
	1139	ierr=1
	1140	do while (ierr.ne.0)
	1141	read(,,iostat=ierr) val2
	1142	! write(,)'val2:',val2,'ierr=',ierr
	1143	if (ierr.eq.0) then ! got a value, but do a sanity check
	1144	if(val2.gt.layer(nsoilmx)) then
	1145	write(,)'Depth should be less than ',layer(nsoilmx)
	1146	ierr=1
	1147	endif
	1148	if(val2.lt.layer(1)) then
	1149	write(,)'Depth should be more than ',layer(1)
	1150	ierr=1
	1151	endif
	1152	endif
	1153	enddo ! of do while
	1154
	1155	! find the reference index iref the depth corresponds to
	1156	! if (val2.lt.layer(1)) then
	1157	! iref=1
	1158	! else
	1159	do isoil=1,nsoilmx-1
	1160	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	1161	& then
	1162	iref=isoil
	1163	exit
	1164	endif
	1165	enddo
	1166	! endif
	1167
	1168	! write(,)'iref:',iref,' jref:',jref
	1169	! write(,)'layer',layer
	1170	! write(,)'mlayer',mlayer
	1171
	1172	! thermal inertia of the ice:
	1173	ierr=1
	1174	do while (ierr.ne.0)
	1175	write(,)'What is the value of subterranean ice thermal inert
	1176	&ia? (e.g.: 2000)'
	1177	read(,,iostat=ierr)iceith
	1178	enddo ! of do while
	1179
	1180	! recast ithfi() onto ith()
	1181	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1182
	1183	do j=1,jref
	1184	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	1185	do i=1,iip1 ! loop on longitudes
	1186	! Build "equivalent" thermal inertia for the mixed layer
	1187	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1188	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	1189	& ((layer(iref+1)-val2)/(iceith)**2)))
	1190	! Set thermal inertia of lower layers
	1191	do isoil=iref+2,nsoilmx
	1192	ith(i,j,isoil)=iceith ! ice
	1193	enddo
	1194	enddo ! of do i=1,iip1
	1195	enddo ! of do j=1,jjp1
	1196
	1197
	1198	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1199
	1200	! do i=1,nsoilmx
	1201	! write(,)'i:',i,'ithfi(1,i):',ithfi(1,i)
	1202	! enddo
	1203
	1204
	1205	! subsoilice_s: Put deep ice layer in southern hemisphere soil
	1206	! ------------------------------------------------------------
	1207
[164]	1208	else if (trim(modif).eq.'subsoilice_s') then
[38]	1209
	1210	write(,)'From which latitude (in deg.), down to the south pol
	1211	&e, should we put subterranean ice?'
	1212	ierr=1
	1213	do while (ierr.ne.0)
	1214	read(,,iostat=ierr) val
	1215	if (ierr.eq.0) then ! got a value
	1216	! do a sanity check
	1217	if((val.gt.0.).or.(val.lt.-90)) then
	1218	write(,)'Latitude should be between 0 and -90 deg. !!!'
	1219	ierr=1
	1220	else ! find corresponding jref (nearest latitude)
	1221	! note: rlatu(:) contains decreasing values of latitude
	1222	! starting from PI/2 to -PI/2
	1223	do j=1,jjp1
	1224	if ((rlatu(j)*180./pi.ge.val).and.
	1225	& (rlatu(j+1)*180./pi.le.val)) then
	1226	! find which grid point is nearest to val:
	1227	if (abs(rlatu(j)*180./pi-val).le.
	1228	& abs((rlatu(j+1)*180./pi-val))) then
	1229	jref=j
	1230	else
	1231	jref=j+1
	1232	endif
	1233
	1234	write(,)'Will use nearest grid latitude which is:',
	1235	& rlatu(jref)*180./pi
	1236	endif
	1237	enddo ! of do j=1,jjp1
	1238	endif ! of if((val.lt.0.).or.(val.gt.90))
	1239	endif !of if (ierr.eq.0)
	1240	enddo ! of do while
	1241
	1242	! Build layers() (as in soil_settings.F)
	1243	val2=sqrt(mlayer(0)*mlayer(1))
	1244	val3=mlayer(1)/mlayer(0)
	1245	do isoil=1,nsoilmx
	1246	layer(isoil)=val2(val3*(isoil-1))
	1247	enddo
	1248
	1249	write(,)'At which depth (in m.) does the ice layer begin?'
	1250	write(,)'(currently, the deepest soil layer extends down to:'
	1251	& ,layer(nsoilmx),')'
	1252	ierr=1
	1253	do while (ierr.ne.0)
	1254	read(,,iostat=ierr) val2
	1255	! write(,)'val2:',val2,'ierr=',ierr
	1256	if (ierr.eq.0) then ! got a value, but do a sanity check
	1257	if(val2.gt.layer(nsoilmx)) then
	1258	write(,)'Depth should be less than ',layer(nsoilmx)
	1259	ierr=1
	1260	endif
	1261	if(val2.lt.layer(1)) then
	1262	write(,)'Depth should be more than ',layer(1)
	1263	ierr=1
	1264	endif
	1265	endif
	1266	enddo ! of do while
	1267
	1268	! find the reference index iref the depth corresponds to
	1269	do isoil=1,nsoilmx-1
	1270	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	1271	& then
	1272	iref=isoil
	1273	exit
	1274	endif
	1275	enddo
	1276
	1277	! write(,)'iref:',iref,' jref:',jref
	1278
	1279	! thermal inertia of the ice:
	1280	ierr=1
	1281	do while (ierr.ne.0)
	1282	write(,)'What is the value of subterranean ice thermal inert
	1283	&ia? (e.g.: 2000)'
	1284	read(,,iostat=ierr)iceith
	1285	enddo ! of do while
	1286
	1287	! recast ithfi() onto ith()
	1288	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1289
	1290	do j=jref,jjp1
	1291	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	1292	do i=1,iip1 ! loop on longitudes
	1293	! Build "equivalent" thermal inertia for the mixed layer
	1294	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1295	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	1296	& ((layer(iref+1)-val2)/(iceith)**2)))
	1297	! Set thermal inertia of lower layers
	1298	do isoil=iref+2,nsoilmx
	1299	ith(i,j,isoil)=iceith ! ice
	1300	enddo
	1301	enddo ! of do i=1,iip1
	1302	enddo ! of do j=jref,jjp1
	1303
	1304
	1305	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1306
	1307	c 'mons_ice' : use MONS data to build subsurface ice table
	1308	c --------------------------------------------------------
[164]	1309	else if (trim(modif).eq.'mons_ice') then
[38]	1310
	1311	! 1. Load MONS data
	1312	call load_MONS_data(MONS_Hdn,MONS_d21)
	1313
	1314	! 2. Get parameters from user
	1315	ierr=1
	1316	do while (ierr.ne.0)
	1317	write(,) "Coefficient to apply to MONS 'depth' in Northern",
	1318	& " Hemisphere?"
	1319	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
	1320	read(,,iostat=ierr) MONS_coeffN
	1321	enddo
	1322	ierr=1
	1323	do while (ierr.ne.0)
	1324	write(,) "Coefficient to apply to MONS 'depth' in Southern",
	1325	& " Hemisphere?"
	1326	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
	1327	read(,,iostat=ierr) MONS_coeffS
	1328	enddo
	1329	ierr=1
	1330	do while (ierr.ne.0)
	1331	write(,) "Value of subterranean ice thermal inertia ",
	1332	& " in Northern hemisphere?"
	1333	write(,) " (e.g.: 2000, or perhaps 2290)"
	1334	! read(,,iostat=ierr) iceith
	1335	read(,,iostat=ierr) iceithN
	1336	enddo
	1337	ierr=1
	1338	do while (ierr.ne.0)
	1339	write(,) "Value of subterranean ice thermal inertia ",
	1340	& " in Southern hemisphere?"
	1341	write(,) " (e.g.: 2000, or perhaps 2290)"
	1342	! read(,,iostat=ierr) iceith
	1343	read(,,iostat=ierr) iceithS
	1344	enddo
	1345
	1346	! 3. Build subterranean thermal inertia
	1347
	1348	! initialise subsurface inertia with reference surface values
	1349	do isoil=1,nsoilmx
	1350	ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1351	enddo
	1352	! recast ithfi() onto ith()
	1353	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1354
	1355	do i=1,iip1 ! loop on longitudes
	1356	do j=1,jjp1 ! loop on latitudes
	1357	! set MONS_coeff
	1358	if (rlatu(j).ge.0) then ! northern hemisphere
	1359	! N.B: rlatu(:) contains decreasing values of latitude
	1360	! starting from PI/2 to -PI/2
	1361	MONS_coeff=MONS_coeffN
	1362	iceith=iceithN
	1363	else ! southern hemisphere
	1364	MONS_coeff=MONS_coeffS
	1365	iceith=iceithS
	1366	endif
	1367	! check if we should put subterranean ice
	1368	if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14%
	1369	! compute depth at which ice lies:
	1370	val=MONS_d21(i,j)*MONS_coeff
	1371	! compute val2= the diurnal skin depth of surface inertia
	1372	! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1
	1373	val2=ith(i,j,1)1.e-6sqrt(88775./3.14159)
	1374	if (val.lt.val2) then
	1375	! ice must be below the (surface inertia) diurnal skin depth
	1376	val=val2
	1377	endif
	1378	if (val.lt.layer(nsoilmx)) then ! subterranean ice
	1379	! find the reference index iref that depth corresponds to
	1380	iref=0
	1381	do isoil=1,nsoilmx-1
	1382	if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1)))
	1383	& then
	1384	iref=isoil
	1385	exit
	1386	endif
	1387	enddo
	1388	! Build "equivalent" thermal inertia for the mixed layer
	1389	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1390	& (((val-layer(iref))/(ith(i,j,iref+1)**2))+
	1391	& ((layer(iref+1)-val)/(iceith)**2)))
	1392	! Set thermal inertia of lower layers
	1393	do isoil=iref+2,nsoilmx
	1394	ith(i,j,isoil)=iceith
	1395	enddo
	1396	endif ! of if (val.lt.layer(nsoilmx))
	1397	endif ! of if (MONS_Hdn(i,j).lt.14.0)
	1398	enddo ! do j=1,jjp1
	1399	enddo ! do i=1,iip1
	1400
	1401	! Check:
	1402	! do i=1,iip1
	1403	! do j=1,jjp1
	1404	! do isoil=1,nsoilmx
	1405	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1406	! enddo
	1407	! enddo
	1408	! enddo
	1409
	1410	! recast ith() into ithfi()
	1411	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1412
	1413	else
	1414	write(,) ' Unknown (misspelled?) option!!!'
[164]	1415	end if ! of if (trim(modif) .eq. '...') elseif ...
[38]	1416
	1417	enddo ! of do ! infinite loop on liste of changes
	1418
	1419	999 continue
	1420
	1421
	1422	c=======================================================================
	1423	c Correct pressure on the new grid (menu 0)
	1424	c=======================================================================
	1425
	1426	if (choix_1.eq.0) then
	1427	r = 1000.*8.31/mugaz
	1428
	1429	do j=1,jjp1
	1430	do i=1,iip1
	1431	ps(i,j) = ps(i,j) *
	1432	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1433	. (t(i,j,1) * r))
	1434	end do
	1435	end do
	1436
	1437	c periodicity of surface ps in longitude
	1438	do j=1,jjp1
	1439	ps(1,j) = ps(iip1,j)
	1440	end do
	1441	end if
	1442
	1443	c=======================================================================
	1444	c=======================================================================
	1445
	1446	c=======================================================================
	1447	c Initialisation de la physique / ecriture de newstartfi :
	1448	c=======================================================================
	1449
	1450
	1451	CALL inifilr
	1452	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	1453
	1454	c-----------------------------------------------------------------------
	1455	c Initialisation pks:
	1456	c-----------------------------------------------------------------------
	1457
	1458	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
	1459	! Calcul de la temperature potentielle teta
	1460
	1461	if (flagiso) then
	1462	DO l=1,llm
	1463	DO j=1,jjp1
	1464	DO i=1,iim
	1465	teta(i,j,l) = Tiso * cpp/pk(i,j,l)
	1466	ENDDO
	1467	teta (iip1,j,l)= teta (1,j,l)
	1468	ENDDO
	1469	ENDDO
	1470	else if (choix_1.eq.0) then
	1471	DO l=1,llm
	1472	DO j=1,jjp1
	1473	DO i=1,iim
	1474	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1475	ENDDO
	1476	teta (iip1,j,l)= teta (1,j,l)
	1477	ENDDO
	1478	ENDDO
	1479	endif
	1480
	1481	C Calcul intermediaire
	1482	c
	1483	if (choix_1.eq.0) then
	1484	CALL massdair( p3d, masse )
	1485	c
	1486	print *,' ALPHAX ',alphax
	1487
	1488	DO l = 1, llm
	1489	DO i = 1, iim
	1490	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1491	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1492	ENDDO
	1493	xpn = SUM(xppn)/apoln
	1494	xps = SUM(xpps)/apols
	1495	DO i = 1, iip1
	1496	masse( i , 1 , l ) = xpn
	1497	masse( i , jjp1 , l ) = xps
	1498	ENDDO
	1499	ENDDO
	1500	endif
	1501	phis(iip1,:) = phis(1,:)
	1502
	1503	CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
	1504	* tetagdiv, tetagrot , tetatemp )
	1505	itau=0
	1506	if (choix_1.eq.0) then
	1507	day_ini=int(date)
[999]	1508	hour_ini=date-int(date)
[38]	1509	endif
	1510	c
	1511	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	1512
	1513	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	1514	* phi,w, pbaru,pbarv,day_ini+time )
	1515	c CALL caldyn
	1516	c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	1517	c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
	1518
[1036]	1519	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqtot)
[999]	1520	CALL dynredem1("restart.nc",hour_ini,vcov,ucov,teta,q,
[1036]	1521	. nqtot,masse,ps)
[38]	1522	C
	1523	C Ecriture etat initial physique
	1524	C
	1525
[1047]	1526	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
	1527	. nqtot,dtphys,real(day_ini),0.0,
[38]	1528	. airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
[1047]	1529	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
[999]	1530	. dtphys,hour_ini,
	1531	. tsurf,tsoil,co2ice,emis,q2,qsurf)
[38]	1532
	1533	c=======================================================================
	1534	c Formats
	1535	c=======================================================================
	1536
	1537	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	1538	*rrage est differente de la valeur parametree iim =',i4//)
	1539	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	1540	*rrage est differente de la valeur parametree jjm =',i4//)
	1541	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	1542	*rage est differente de la valeur parametree llm =',i4//)
	1543
[171]	1544	write(,) "newstart: All is well that ends well."
	1545
[38]	1546	end
	1547
	1548	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1549	subroutine load_MONS_data(MONS_Hdn,MONS_d21)
	1550	implicit none
	1551	! routine to load Benedicte Diez MONS dataset, fill in date in southern
	1552	! polar region, and interpolate the result onto the GCM grid
	1553	#include"dimensions.h"
	1554	#include"paramet.h"
	1555	#include"datafile.h"
	1556	#include"comgeom.h"
	1557	! arguments:
	1558	real,intent(out) :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
	1559	real,intent(out) :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
	1560	! N.B MONS datasets should be of dimension (iip1,jjp1)
	1561	! local variables:
	1562	character(len=88) :: filename="results_MONS_lat_lon_H_depth.txt"
	1563	character(len=88) :: txt ! to store some text
	1564	integer :: ierr,i,j
	1565	integer,parameter :: nblon=180 ! number of longitudes of MONS datasets
	1566	integer,parameter :: nblat=90 ! number of latitudes of MONS datasets
	1567	real :: pi
	1568	real :: longitudes(nblon) ! MONS dataset longitudes
	1569	real :: latitudes(nblat) ! MONS dataset latitudes
	1570	! MONS dataset: mass fraction of H2O where H is assumed to be in H2O
	1571	real :: Hdn(nblon,nblat)
	1572	real :: d21(nblon,nblat)! MONS dataset "depth" (g/cm2)
	1573
	1574	! Extended MONS dataset (for interp_horiz)
	1575	real :: Hdnx(nblon+1,nblat)
	1576	real :: d21x(nblon+1,nblat)
	1577	real :: lon_bound(nblon+1) ! longitude boundaries
	1578	real :: lat_bound(nblat-1) ! latitude boundaries
	1579
	1580	! 1. Initializations:
	1581
	1582	write(,) "Loading MONS data"
	1583
	1584	! Open MONS datafile:
	1585	open(42,file=trim(datafile)//"/"//trim(filename),
	1586	& status="old",iostat=ierr)
	1587	if (ierr/=0) then
	1588	write(,) "Error in load_MONS_data:"
	1589	write(,) "Failed opening file ",
	1590	& trim(datafile)//"/"//trim(filename)
	1591	write(,)'1) You can change the path to the file in '
	1592	write(,)' file phymars/datafile.h'
	1593	write(,)'2) If necessary ',trim(filename),
	1594	& ' (and other datafiles)'
	1595	write(,)' can be obtained online at:'
	1596	write(,)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
	1597	CALL ABORT
	1598	else ! skip first line of file (dummy read)
	1599	read(42,*) txt
	1600	endif
	1601
	1602	pi=2.*asin(1.)
	1603
	1604	!2. Load MONS data (on MONS grid)
	1605	do j=1,nblat
	1606	do i=1,nblon
	1607	! swap latitude index so latitudes go from north pole to south pole:
	1608	read(42,*) latitudes(nblat-j+1),longitudes(i),
	1609	& Hdn(i,nblat-j+1),d21(i,nblat-j+1)
	1610	! multiply d21 by 10 to convert from g/cm2 to kg/m2
	1611	d21(i,nblat-j+1)=d21(i,nblat-j+1)*10.0
	1612	enddo
	1613	enddo
	1614	close(42)
	1615
	1616	! there is unfortunately no d21 data for latitudes -77 to -90
	1617	! so we build some by linear interpolation between values at -75
	1618	! and assuming d21=0 at the pole
	1619	do j=84,90 ! latitudes(84)=-77 ; latitudes(83)=-75
	1620	do i=1,nblon
	1621	d21(i,j)=d21(i,83)*((latitudes(j)+90)/15.0)
	1622	enddo
	1623	enddo
	1624
	1625	! 3. Build extended MONS dataset & boundaries (for interp_horiz)
	1626	! longitude boundaries (in radians):
	1627	do i=1,nblon
	1628	! NB: MONS data is every 2 degrees in longitude
	1629	lon_bound(i)=(longitudes(i)+1.0)*pi/180.0
	1630	enddo
	1631	! extra 'modulo' value
	1632	lon_bound(nblon+1)=lon_bound(1)+2.0*pi
	1633
	1634	! latitude boundaries (in radians):
	1635	do j=1,nblat-1
	1636	! NB: Mons data is every 2 degrees in latitude
	1637	lat_bound(j)=(latitudes(j)-1.0)*pi/180.0
	1638	enddo
	1639
	1640	! MONS datasets:
	1641	do j=1,nblat
	1642	Hdnx(1:nblon,j)=Hdn(1:nblon,j)
	1643	Hdnx(nblon+1,j)=Hdnx(1,j)
	1644	d21x(1:nblon,j)=d21(1:nblon,j)
	1645	d21x(nblon+1,j)=d21x(1,j)
	1646	enddo
	1647
	1648	! Interpolate onto GCM grid
	1649	call interp_horiz(Hdnx,MONS_Hdn,nblon,nblat-1,iim,jjm,1,
	1650	& lon_bound,lat_bound,rlonu,rlatv)
	1651	call interp_horiz(d21x,MONS_d21,nblon,nblat-1,iim,jjm,1,
	1652	& lon_bound,lat_bound,rlonu,rlatv)
	1653
	1654	end subroutine

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

Download in other formats:

Original Format