source: trunk/mars/libf/dyn3d/newstart.F @ 77

Last change on this file since 77 was 38, checked in by emillour, 14 years ago

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