source: trunk/LMDZ.MARS/libf/dyn3d/newstart.F @ 480

Last change on this file since 480 was 480, checked in by tnavarro, 13 years ago

new option in newstart: ini_h2osurf, water ice surface initialisation to keep seasonal frost

but get rid of polar cap positive or negative values.

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