source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/newstart.F @ 222

Last change on this file since 222 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

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