Context Navigation

newstart.F @ 1212

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

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

Download in other formats:

Original Format