source: trunk/LMDZ.GENERIC/libf/dyn3d/newstart.F @ 889

Last change on this file since 889 was 837, checked in by aslmd, 13 years ago

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