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

Last change on this file since 313 was 313, checked in by ymipsl, 10 years ago
  • implement splitting of XIOS file for lmdz physics
  • Termination is done properly in parallel by calling MPI_ABORT instead of abort or stop

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