source: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F @ 2322

Last change on this file since 2322 was 2312, checked in by lrossi, 5 years ago

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