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

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