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

Last change on this file since 2609 was 2565, checked in by emillour, 4 years ago

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