source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/newstart.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 125.4 KB
Line 
1C======================================================================
2      PROGRAM newstart
3c=======================================================================
4c
5c
6c   Auteur:   Christophe Hourdin/Francois Forget/Yann Wanherdrick
7c   ------
8c             Derniere modif : 12/03
9c
10c
11c   Objet:  Create or modify the initial state for the LMD Mars GCM
12c   -----           (fichiers NetCDF start et startfi)
13c
14c
15c=======================================================================
16
17      use infotrac, only: iniadvtrac, nqtot, tname
18      USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
19      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
20      USE surfdat_h, ONLY: phisfi, albedodat,
21     &                     zmea, zstd, zsig, zgam, zthe
22      USE comgeomfi_h, ONLY: lati, long, area
23      use datafile_mod, only: datadir
24! to use  'getin'
25      USE ioipsl_getincom, only: getin
26      use control_mod, only: day_step, iphysiq, anneeref
27      use phyredem, only: physdem0, physdem1
28      use iostart, only: open_startphy
29      use comgeomphy, only: initcomgeomphy
30      use slab_ice_h, only:noceanmx
31      implicit none
32
33!-----------------------------------------------------------------------
34!   INCLUDE 'dimensions.h'
35!
36!   dimensions.h contient les dimensions du modele
37!   ndm est tel que iim=2**ndm
38!-----------------------------------------------------------------------
39
40      INTEGER iim,jjm,llm,ndm
41
42      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
43
44!-----------------------------------------------------------------------
45!-----------------------------------------------------------------------
46!   INCLUDE 'dimphys.h'
47
48! ngridmx : number of horizontal grid points
49! note: the -1/jjm term will be 0; unless jj=1
50      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)   
51! nlayermx : number of atmospheric layers
52      integer, parameter :: nlayermx = llm 
53! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
54      !integer, parameter :: nsoilmx = 4 ! for a test
55      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
56      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
57!-----------------------------------------------------------------------
58
59!-----------------------------------------------------------------------
60! INCLUDE planet.h
61
62      COMMON/planet/apoastr,periastr,year_day,peri_day,                 &
63     &       obliquit,nres,                                             &
64     &       z0,lmixmin,emin_turb,coefvis,coefir,                       &
65     &       timeperi,e_elips,p_elips
66
67      real apoastr,periastr,year_day,peri_day,                          &
68     &     obliquit,nres,                                               &
69     &     z0,lmixmin,emin_turb,coefvis,coefir,                         &
70     &       timeperi,e_elips,p_elips
71     
72
73!-----------------------------------------------------------------------
74!
75! $Header$
76!
77!
78!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
79!                 veillez  n'utiliser que des ! pour les commentaires
80!                 et  bien positionner les & des lignes de continuation
81!                 (les placer en colonne 6 et en colonne 73)
82!
83!
84!-----------------------------------------------------------------------
85!   INCLUDE 'paramet.h'
86
87      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
88      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
89      INTEGER  ijmllm,mvar
90      INTEGER jcfil,jcfllm
91
92      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
93     &    ,jjp1=jjm+1-1/jjm)
94      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
95      PARAMETER( kftd  = iim/2 -ndm )
96      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
97      PARAMETER( ip1jmi1= ip1jm - iip1 )
98      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
99      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
100      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
101
102!-----------------------------------------------------------------------
103!
104! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
105!
106!-----------------------------------------------------------------------
107! INCLUDE comconst.h
108
109      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
110     &                 iflag_top_bound,mode_top_bound
111      COMMON/comconstr/dtvr,daysec,                                     &
112     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
113     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
114     & ,dissip_pupstart  ,tau_top_bound,                                &
115     & daylen,molmass, ihf
116      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
117
118      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
119      REAL dtvr ! dynamical time step (in s)
120      REAL daysec !length (in s) of a standard day
121      REAL pi    ! something like 3.14159....
122      REAL dtphys ! (s) time step for the physics
123      REAL dtdiss ! (s) time step for the dissipation
124      REAL rad ! (m) radius of the planet
125      REAL r ! Reduced Gas constant r=R/mu
126             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
127      REAL cpp   ! Cp
128      REAL kappa ! kappa=R/Cp
129      REAL cotot
130      REAL unsim ! = 1./iim
131      REAL g ! (m/s2) gravity
132      REAL omeg ! (rad/s) rotation rate of the planet
133! Dissipation factors, for Earth model:
134      REAL dissip_factz,dissip_zref !dissip_deltaz
135! Dissipation factors, for other planets:
136      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
137      REAL dissip_pupstart
138      INTEGER iflag_top_bound,mode_top_bound
139      REAL tau_top_bound
140      REAL daylen ! length of solar day, in 'standard' day length
141      REAL molmass ! (g/mol) molar mass of the atmosphere
142
143      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
144      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
145
146
147!-----------------------------------------------------------------------
148!
149! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
150!
151!-----------------------------------------------------------------------
152!   INCLUDE 'comvert.h'
153
154      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
155     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
156     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
157
158      common/comverti/disvert_type, pressure_exner
159
160      real ap     ! hybrid pressure contribution at interlayers
161      real bp     ! hybrid sigma contribution at interlayer
162      real presnivs ! (reference) pressure at mid-layers
163      real dpres
164      real pa     ! reference pressure (Pa) at which hybrid coordinates
165                  ! become purely pressure
166      real preff  ! reference surface pressure (Pa)
167      real nivsigs
168      real nivsig
169      real aps    ! hybrid pressure contribution at mid-layers
170      real bps    ! hybrid sigma contribution at mid-layers
171      real scaleheight ! atmospheric (reference) scale height (km)
172      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
173                     ! preff and scaleheight
174
175      integer disvert_type ! type of vertical discretization:
176                           ! 1: Earth (default for planet_type==earth),
177                           !     automatic generation
178                           ! 2: Planets (default for planet_type!=earth),
179                           !     using 'z2sig.def' (or 'esasig.def) file
180
181      logical pressure_exner
182!     compute pressure inside layers using Exner function, else use mean
183!     of pressure values at interfaces
184
185 !-----------------------------------------------------------------------
186!
187! $Header$
188!
189!CDK comgeom2
190      COMMON/comgeom/                                                   &
191     & cu(iip1,jjp1),cv(iip1,jjm),unscu2(iip1,jjp1),unscv2(iip1,jjm)  , &
192     & aire(iip1,jjp1),airesurg(iip1,jjp1),aireu(iip1,jjp1)           , &
193     & airev(iip1,jjm),unsaire(iip1,jjp1),apoln,apols                 , &
194     & unsairez(iip1,jjm),airuscv2(iip1,jjm),airvscu2(iip1,jjm)       , &
195     & aireij1(iip1,jjp1),aireij2(iip1,jjp1),aireij3(iip1,jjp1)       , &
196     & aireij4(iip1,jjp1),alpha1(iip1,jjp1),alpha2(iip1,jjp1)         , &
197     & alpha3(iip1,jjp1),alpha4(iip1,jjp1),alpha1p2(iip1,jjp1)        , &
198     & alpha1p4(iip1,jjp1),alpha2p3(iip1,jjp1),alpha3p4(iip1,jjp1)    , &
199     & fext(iip1,jjm),constang(iip1,jjp1), rlatu(jjp1),rlatv(jjm),      &
200     & rlonu(iip1),rlonv(iip1),cuvsurcv(iip1,jjm),cvsurcuv(iip1,jjm)  , &
201     & cvusurcu(iip1,jjp1),cusurcvu(iip1,jjp1)                        , &
202     & cuvscvgam1(iip1,jjm),cuvscvgam2(iip1,jjm),cvuscugam1(iip1,jjp1), &
203     & cvuscugam2(iip1,jjp1),cvscuvgam(iip1,jjm),cuscvugam(iip1,jjp1) , &
204     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2                , &
205     & unsair_gam1(iip1,jjp1),unsair_gam2(iip1,jjp1)                  , &
206     & unsairz_gam(iip1,jjm),aivscu2gam(iip1,jjm),aiuscv2gam(iip1,jjm)  &
207     & , xprimu(iip1),xprimv(iip1)
208
209
210      REAL                                                               &
211     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,apoln,apols,unsaire &
212     & ,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4     , &
213     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 , &
214     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     , &
215     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1           , &
216     & unsapolnga2,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2     , &
217     & unsairz_gam,aivscu2gam,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu    , &
218     & cusurcvu,xprimu,xprimv
219!#include "control.h"
220!
221! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
222!
223!
224! NB: keep items of different kinds in seperate common blocs to avoid
225!     "misaligned commons" issues
226!-----------------------------------------------------------------------
227! INCLUDE 'logic.h'
228
229      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
230     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
231     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
232     &  ,ok_limit,ok_etat0,hybrid                                       &
233     &  ,moyzon_mu,moyzon_ch
234
235      COMMON/logici/ iflag_phys,iflag_trac
236     
237      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
238     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
239     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
240     &  ,ok_limit,ok_etat0
241      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
242                     ! (only used if disvert_type==2)
243      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
244
245      integer iflag_phys,iflag_trac
246!$OMP THREADPRIVATE(/logicl/)
247!$OMP THREADPRIVATE(/logici/)
248!-----------------------------------------------------------------------
249!
250! $Id: ener.h 1447 2010-10-22 16:18:27Z jghattas $
251!
252!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
253!                 veillez à n'utiliser que des ! pour les commentaires
254!                 et à bien positionner les & des lignes de continuation
255!                 (les placer en colonne 6 et en colonne 73)
256!
257! INCLUDE 'ener.h'
258
259      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,                         &
260     &            ang,etot,ptot,ztot,stot,rmsdpdt ,                     &
261     &            rmsv,gtot(llmm1)
262
263      REAL ang0,etot0,ptot0,ztot0,stot0,                                &
264     &     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
265
266!-----------------------------------------------------------------------
267!
268! $Id: temps.h 1577 2011-10-20 15:06:47Z fairhead $
269!
270!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
271!                 veillez  n'utiliser que des ! pour les commentaires
272!                 et  bien positionner les & des lignes de continuation
273!                 (les placer en colonne 6 et en colonne 73)
274!
275!
276! jD_ref = jour julien de la date de reference (lancement de l'experience)
277! hD_ref = "heure" julienne de la date de reference
278!-----------------------------------------------------------------------
279! INCLUDE 'temps.h'
280
281      COMMON/temps_r/dt,jD_ref,jH_ref,start_time,hour_ini
282      COMMON/temps_i/day_ini,day_end,annee_ref,day_ref,                 &
283     &             itau_dyn,itau_phy,itaufin
284      COMMON/temps_c/calend
285
286
287      INTEGER   itaufin ! total number of dynamical steps for the run
288      INTEGER   itau_dyn, itau_phy
289      INTEGER   day_ini ! initial day # of simulation sequence
290      INTEGER   day_end ! final day # ; i.e. day # when this simulation ends
291      INTEGER   annee_ref
292      INTEGER   day_ref
293      REAL      dt ! (dynamics) time step (changes if doing Matsuno or LF step)
294      REAL      jD_ref, jH_ref, start_time
295      CHARACTER (len=10) :: calend
296
297      ! Additionnal Mars stuff:
298      real hour_ini ! initial fraction of day of simulation sequence (0=<hour_ini<1)
299
300!-----------------------------------------------------------------------
301!
302! $Id: comdissnew.h 1319 2010-02-23 21:29:54Z fairhead $
303!
304!
305!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
306!                 veillez à n'utiliser que des ! pour les commentaires
307!                 et à bien positionner les & des lignes de continuation
308!                 (les placer en colonne 6 et en colonne 73)
309!
310!-----------------------------------------------------------------------
311! INCLUDE 'comdissnew.h'
312
313      COMMON/comdissnew/ lstardis,nitergdiv,nitergrot,niterh,tetagdiv,  &
314     &                   tetagrot,tetatemp,coefdis, vert_prof_dissip
315
316      LOGICAL lstardis
317      INTEGER nitergdiv, nitergrot, niterh
318
319! For the Earth model:
320      integer vert_prof_dissip ! vertical profile of horizontal dissipation
321!     Allowed values:
322!     0: rational fraction, function of pressure
323!     1: tanh of altitude
324
325      REAL     tetagdiv, tetagrot,  tetatemp, coefdis
326
327!
328! ... Les parametres de ce common comdissnew sont  lues par defrun_new
329!              sur le fichier  run.def    ....
330!
331!-----------------------------------------------------------------------
332!
333! $Header$
334!
335!c
336!c
337!c..include serre.h
338!c
339       REAL clon,clat,transx,transy,alphax,alphay,pxo,pyo,              &
340     &  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
341       COMMON/serre/clon,clat,transx,transy,alphax,alphay,pxo,pyo ,     &
342     &  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
343!     NetCDF-3.
344!
345! netcdf version 3 fortran interface:
346!
347
348!
349! external netcdf data types:
350!
351      integer nf_byte
352      integer nf_int1
353      integer nf_char
354      integer nf_short
355      integer nf_int2
356      integer nf_int
357      integer nf_float
358      integer nf_real
359      integer nf_double
360
361      parameter (nf_byte = 1)
362      parameter (nf_int1 = nf_byte)
363      parameter (nf_char = 2)
364      parameter (nf_short = 3)
365      parameter (nf_int2 = nf_short)
366      parameter (nf_int = 4)
367      parameter (nf_float = 5)
368      parameter (nf_real = nf_float)
369      parameter (nf_double = 6)
370
371!
372! default fill values:
373!
374      integer           nf_fill_byte
375      integer           nf_fill_int1
376      integer           nf_fill_char
377      integer           nf_fill_short
378      integer           nf_fill_int2
379      integer           nf_fill_int
380      real              nf_fill_float
381      real              nf_fill_real
382      doubleprecision   nf_fill_double
383
384      parameter (nf_fill_byte = -127)
385      parameter (nf_fill_int1 = nf_fill_byte)
386      parameter (nf_fill_char = 0)
387      parameter (nf_fill_short = -32767)
388      parameter (nf_fill_int2 = nf_fill_short)
389      parameter (nf_fill_int = -2147483647)
390      parameter (nf_fill_float = 9.9692099683868690e+36)
391      parameter (nf_fill_real = nf_fill_float)
392      parameter (nf_fill_double = 9.9692099683868690d+36)
393
394!
395! mode flags for opening and creating a netcdf dataset:
396!
397      integer nf_nowrite
398      integer nf_write
399      integer nf_clobber
400      integer nf_noclobber
401      integer nf_fill
402      integer nf_nofill
403      integer nf_lock
404      integer nf_share
405      integer nf_64bit_offset
406      integer nf_sizehint_default
407      integer nf_align_chunk
408      integer nf_format_classic
409      integer nf_format_64bit
410
411      parameter (nf_nowrite = 0)
412      parameter (nf_write = 1)
413      parameter (nf_clobber = 0)
414      parameter (nf_noclobber = 4)
415      parameter (nf_fill = 0)
416      parameter (nf_nofill = 256)
417      parameter (nf_lock = 1024)
418      parameter (nf_share = 2048)
419      parameter (nf_64bit_offset = 512)
420      parameter (nf_sizehint_default = 0)
421      parameter (nf_align_chunk = -1)
422      parameter (nf_format_classic = 1)
423      parameter (nf_format_64bit = 2)
424
425!
426! size argument for defining an unlimited dimension:
427!
428      integer nf_unlimited
429      parameter (nf_unlimited = 0)
430
431!
432! global attribute id:
433!
434      integer nf_global
435      parameter (nf_global = 0)
436
437!
438! implementation limits:
439!
440      integer nf_max_dims
441      integer nf_max_attrs
442      integer nf_max_vars
443      integer nf_max_name
444      integer nf_max_var_dims
445
446      parameter (nf_max_dims = 1024)
447      parameter (nf_max_attrs = 8192)
448      parameter (nf_max_vars = 8192)
449      parameter (nf_max_name = 256)
450      parameter (nf_max_var_dims = nf_max_dims)
451
452!
453! error codes:
454!
455      integer nf_noerr
456      integer nf_ebadid
457      integer nf_eexist
458      integer nf_einval
459      integer nf_eperm
460      integer nf_enotindefine
461      integer nf_eindefine
462      integer nf_einvalcoords
463      integer nf_emaxdims
464      integer nf_enameinuse
465      integer nf_enotatt
466      integer nf_emaxatts
467      integer nf_ebadtype
468      integer nf_ebaddim
469      integer nf_eunlimpos
470      integer nf_emaxvars
471      integer nf_enotvar
472      integer nf_eglobal
473      integer nf_enotnc
474      integer nf_ests
475      integer nf_emaxname
476      integer nf_eunlimit
477      integer nf_enorecvars
478      integer nf_echar
479      integer nf_eedge
480      integer nf_estride
481      integer nf_ebadname
482      integer nf_erange
483      integer nf_enomem
484      integer nf_evarsize
485      integer nf_edimsize
486      integer nf_etrunc
487
488      parameter (nf_noerr = 0)
489      parameter (nf_ebadid = -33)
490      parameter (nf_eexist = -35)
491      parameter (nf_einval = -36)
492      parameter (nf_eperm = -37)
493      parameter (nf_enotindefine = -38)
494      parameter (nf_eindefine = -39)
495      parameter (nf_einvalcoords = -40)
496      parameter (nf_emaxdims = -41)
497      parameter (nf_enameinuse = -42)
498      parameter (nf_enotatt = -43)
499      parameter (nf_emaxatts = -44)
500      parameter (nf_ebadtype = -45)
501      parameter (nf_ebaddim = -46)
502      parameter (nf_eunlimpos = -47)
503      parameter (nf_emaxvars = -48)
504      parameter (nf_enotvar = -49)
505      parameter (nf_eglobal = -50)
506      parameter (nf_enotnc = -51)
507      parameter (nf_ests = -52)
508      parameter (nf_emaxname = -53)
509      parameter (nf_eunlimit = -54)
510      parameter (nf_enorecvars = -55)
511      parameter (nf_echar = -56)
512      parameter (nf_eedge = -57)
513      parameter (nf_estride = -58)
514      parameter (nf_ebadname = -59)
515      parameter (nf_erange = -60)
516      parameter (nf_enomem = -61)
517      parameter (nf_evarsize = -62)
518      parameter (nf_edimsize = -63)
519      parameter (nf_etrunc = -64)
520!
521! error handling modes:
522!
523      integer  nf_fatal
524      integer nf_verbose
525
526      parameter (nf_fatal = 1)
527      parameter (nf_verbose = 2)
528
529!
530! miscellaneous routines:
531!
532      character*80   nf_inq_libvers
533      external       nf_inq_libvers
534
535      character*80   nf_strerror
536!                         (integer             ncerr)
537      external       nf_strerror
538
539      logical        nf_issyserr
540!                         (integer             ncerr)
541      external       nf_issyserr
542
543!
544! control routines:
545!
546      integer         nf_inq_base_pe
547!                         (integer             ncid,
548!                          integer             pe)
549      external        nf_inq_base_pe
550
551      integer         nf_set_base_pe
552!                         (integer             ncid,
553!                          integer             pe)
554      external        nf_set_base_pe
555
556      integer         nf_create
557!                         (character*(*)       path,
558!                          integer             cmode,
559!                          integer             ncid)
560      external        nf_create
561
562      integer         nf__create
563!                         (character*(*)       path,
564!                          integer             cmode,
565!                          integer             initialsz,
566!                          integer             chunksizehint,
567!                          integer             ncid)
568      external        nf__create
569
570      integer         nf__create_mp
571!                         (character*(*)       path,
572!                          integer             cmode,
573!                          integer             initialsz,
574!                          integer             basepe,
575!                          integer             chunksizehint,
576!                          integer             ncid)
577      external        nf__create_mp
578
579      integer         nf_open
580!                         (character*(*)       path,
581!                          integer             mode,
582!                          integer             ncid)
583      external        nf_open
584
585      integer         nf__open
586!                         (character*(*)       path,
587!                          integer             mode,
588!                          integer             chunksizehint,
589!                          integer             ncid)
590      external        nf__open
591
592      integer         nf__open_mp
593!                         (character*(*)       path,
594!                          integer             mode,
595!                          integer             basepe,
596!                          integer             chunksizehint,
597!                          integer             ncid)
598      external        nf__open_mp
599
600      integer         nf_set_fill
601!                         (integer             ncid,
602!                          integer             fillmode,
603!                          integer             old_mode)
604      external        nf_set_fill
605
606      integer         nf_set_default_format
607!                          (integer             format,
608!                          integer             old_format)
609      external        nf_set_default_format
610
611      integer         nf_redef
612!                         (integer             ncid)
613      external        nf_redef
614
615      integer         nf_enddef
616!                         (integer             ncid)
617      external        nf_enddef
618
619      integer         nf__enddef
620!                         (integer             ncid,
621!                          integer             h_minfree,
622!                          integer             v_align,
623!                          integer             v_minfree,
624!                          integer             r_align)
625      external        nf__enddef
626
627      integer         nf_sync
628!                         (integer             ncid)
629      external        nf_sync
630
631      integer         nf_abort
632!                         (integer             ncid)
633      external        nf_abort
634
635      integer         nf_close
636!                         (integer             ncid)
637      external        nf_close
638
639      integer         nf_delete
640!                         (character*(*)       ncid)
641      external        nf_delete
642
643!
644! general inquiry routines:
645!
646
647      integer         nf_inq
648!                         (integer             ncid,
649!                          integer             ndims,
650!                          integer             nvars,
651!                          integer             ngatts,
652!                          integer             unlimdimid)
653      external        nf_inq
654
655      integer         nf_inq_ndims
656!                         (integer             ncid,
657!                          integer             ndims)
658      external        nf_inq_ndims
659
660      integer         nf_inq_nvars
661!                         (integer             ncid,
662!                          integer             nvars)
663      external        nf_inq_nvars
664
665      integer         nf_inq_natts
666!                         (integer             ncid,
667!                          integer             ngatts)
668      external        nf_inq_natts
669
670      integer         nf_inq_unlimdim
671!                         (integer             ncid,
672!                          integer             unlimdimid)
673      external        nf_inq_unlimdim
674
675      integer         nf_inq_format
676!                         (integer             ncid,
677!                          integer             format)
678      external        nf_inq_format
679
680!
681! dimension routines:
682!
683
684      integer         nf_def_dim
685!                         (integer             ncid,
686!                          character(*)        name,
687!                          integer             len,
688!                          integer             dimid)
689      external        nf_def_dim
690
691      integer         nf_inq_dimid
692!                         (integer             ncid,
693!                          character(*)        name,
694!                          integer             dimid)
695      external        nf_inq_dimid
696
697      integer         nf_inq_dim
698!                         (integer             ncid,
699!                          integer             dimid,
700!                          character(*)        name,
701!                          integer             len)
702      external        nf_inq_dim
703
704      integer         nf_inq_dimname
705!                         (integer             ncid,
706!                          integer             dimid,
707!                          character(*)        name)
708      external        nf_inq_dimname
709
710      integer         nf_inq_dimlen
711!                         (integer             ncid,
712!                          integer             dimid,
713!                          integer             len)
714      external        nf_inq_dimlen
715
716      integer         nf_rename_dim
717!                         (integer             ncid,
718!                          integer             dimid,
719!                          character(*)        name)
720      external        nf_rename_dim
721
722!
723! general attribute routines:
724!
725
726      integer         nf_inq_att
727!                         (integer             ncid,
728!                          integer             varid,
729!                          character(*)        name,
730!                          integer             xtype,
731!                          integer             len)
732      external        nf_inq_att
733
734      integer         nf_inq_attid
735!                         (integer             ncid,
736!                          integer             varid,
737!                          character(*)        name,
738!                          integer             attnum)
739      external        nf_inq_attid
740
741      integer         nf_inq_atttype
742!                         (integer             ncid,
743!                          integer             varid,
744!                          character(*)        name,
745!                          integer             xtype)
746      external        nf_inq_atttype
747
748      integer         nf_inq_attlen
749!                         (integer             ncid,
750!                          integer             varid,
751!                          character(*)        name,
752!                          integer             len)
753      external        nf_inq_attlen
754
755      integer         nf_inq_attname
756!                         (integer             ncid,
757!                          integer             varid,
758!                          integer             attnum,
759!                          character(*)        name)
760      external        nf_inq_attname
761
762      integer         nf_copy_att
763!                         (integer             ncid_in,
764!                          integer             varid_in,
765!                          character(*)        name,
766!                          integer             ncid_out,
767!                          integer             varid_out)
768      external        nf_copy_att
769
770      integer         nf_rename_att
771!                         (integer             ncid,
772!                          integer             varid,
773!                          character(*)        curname,
774!                          character(*)        newname)
775      external        nf_rename_att
776
777      integer         nf_del_att
778!                         (integer             ncid,
779!                          integer             varid,
780!                          character(*)        name)
781      external        nf_del_att
782
783!
784! attribute put/get routines:
785!
786
787      integer         nf_put_att_text
788!                         (integer             ncid,
789!                          integer             varid,
790!                          character(*)        name,
791!                          integer             len,
792!                          character(*)        text)
793      external        nf_put_att_text
794
795      integer         nf_get_att_text
796!                         (integer             ncid,
797!                          integer             varid,
798!                          character(*)        name,
799!                          character(*)        text)
800      external        nf_get_att_text
801
802      integer         nf_put_att_int1
803!                         (integer             ncid,
804!                          integer             varid,
805!                          character(*)        name,
806!                          integer             xtype,
807!                          integer             len,
808!                          nf_int1_t           i1vals(1))
809      external        nf_put_att_int1
810
811      integer         nf_get_att_int1
812!                         (integer             ncid,
813!                          integer             varid,
814!                          character(*)        name,
815!                          nf_int1_t           i1vals(1))
816      external        nf_get_att_int1
817
818      integer         nf_put_att_int2
819!                         (integer             ncid,
820!                          integer             varid,
821!                          character(*)        name,
822!                          integer             xtype,
823!                          integer             len,
824!                          nf_int2_t           i2vals(1))
825      external        nf_put_att_int2
826
827      integer         nf_get_att_int2
828!                         (integer             ncid,
829!                          integer             varid,
830!                          character(*)        name,
831!                          nf_int2_t           i2vals(1))
832      external        nf_get_att_int2
833
834      integer         nf_put_att_int
835!                         (integer             ncid,
836!                          integer             varid,
837!                          character(*)        name,
838!                          integer             xtype,
839!                          integer             len,
840!                          integer             ivals(1))
841      external        nf_put_att_int
842
843      integer         nf_get_att_int
844!                         (integer             ncid,
845!                          integer             varid,
846!                          character(*)        name,
847!                          integer             ivals(1))
848      external        nf_get_att_int
849
850      integer         nf_put_att_real
851!                         (integer             ncid,
852!                          integer             varid,
853!                          character(*)        name,
854!                          integer             xtype,
855!                          integer             len,
856!                          real                rvals(1))
857      external        nf_put_att_real
858
859      integer         nf_get_att_real
860!                         (integer             ncid,
861!                          integer             varid,
862!                          character(*)        name,
863!                          real                rvals(1))
864      external        nf_get_att_real
865
866      integer         nf_put_att_double
867!                         (integer             ncid,
868!                          integer             varid,
869!                          character(*)        name,
870!                          integer             xtype,
871!                          integer             len,
872!                          double              dvals(1))
873      external        nf_put_att_double
874
875      integer         nf_get_att_double
876!                         (integer             ncid,
877!                          integer             varid,
878!                          character(*)        name,
879!                          double              dvals(1))
880      external        nf_get_att_double
881
882!
883! general variable routines:
884!
885
886      integer         nf_def_var
887!                         (integer             ncid,
888!                          character(*)        name,
889!                          integer             datatype,
890!                          integer             ndims,
891!                          integer             dimids(1),
892!                          integer             varid)
893      external        nf_def_var
894
895      integer         nf_inq_var
896!                         (integer             ncid,
897!                          integer             varid,
898!                          character(*)        name,
899!                          integer             datatype,
900!                          integer             ndims,
901!                          integer             dimids(1),
902!                          integer             natts)
903      external        nf_inq_var
904
905      integer         nf_inq_varid
906!                         (integer             ncid,
907!                          character(*)        name,
908!                          integer             varid)
909      external        nf_inq_varid
910
911      integer         nf_inq_varname
912!                         (integer             ncid,
913!                          integer             varid,
914!                          character(*)        name)
915      external        nf_inq_varname
916
917      integer         nf_inq_vartype
918!                         (integer             ncid,
919!                          integer             varid,
920!                          integer             xtype)
921      external        nf_inq_vartype
922
923      integer         nf_inq_varndims
924!                         (integer             ncid,
925!                          integer             varid,
926!                          integer             ndims)
927      external        nf_inq_varndims
928
929      integer         nf_inq_vardimid
930!                         (integer             ncid,
931!                          integer             varid,
932!                          integer             dimids(1))
933      external        nf_inq_vardimid
934
935      integer         nf_inq_varnatts
936!                         (integer             ncid,
937!                          integer             varid,
938!                          integer             natts)
939      external        nf_inq_varnatts
940
941      integer         nf_rename_var
942!                         (integer             ncid,
943!                          integer             varid,
944!                          character(*)        name)
945      external        nf_rename_var
946
947      integer         nf_copy_var
948!                         (integer             ncid_in,
949!                          integer             varid,
950!                          integer             ncid_out)
951      external        nf_copy_var
952
953!
954! entire variable put/get routines:
955!
956
957      integer         nf_put_var_text
958!                         (integer             ncid,
959!                          integer             varid,
960!                          character(*)        text)
961      external        nf_put_var_text
962
963      integer         nf_get_var_text
964!                         (integer             ncid,
965!                          integer             varid,
966!                          character(*)        text)
967      external        nf_get_var_text
968
969      integer         nf_put_var_int1
970!                         (integer             ncid,
971!                          integer             varid,
972!                          nf_int1_t           i1vals(1))
973      external        nf_put_var_int1
974
975      integer         nf_get_var_int1
976!                         (integer             ncid,
977!                          integer             varid,
978!                          nf_int1_t           i1vals(1))
979      external        nf_get_var_int1
980
981      integer         nf_put_var_int2
982!                         (integer             ncid,
983!                          integer             varid,
984!                          nf_int2_t           i2vals(1))
985      external        nf_put_var_int2
986
987      integer         nf_get_var_int2
988!                         (integer             ncid,
989!                          integer             varid,
990!                          nf_int2_t           i2vals(1))
991      external        nf_get_var_int2
992
993      integer         nf_put_var_int
994!                         (integer             ncid,
995!                          integer             varid,
996!                          integer             ivals(1))
997      external        nf_put_var_int
998
999      integer         nf_get_var_int
1000!                         (integer             ncid,
1001!                          integer             varid,
1002!                          integer             ivals(1))
1003      external        nf_get_var_int
1004
1005      integer         nf_put_var_real
1006!                         (integer             ncid,
1007!                          integer             varid,
1008!                          real                rvals(1))
1009      external        nf_put_var_real
1010
1011      integer         nf_get_var_real
1012!                         (integer             ncid,
1013!                          integer             varid,
1014!                          real                rvals(1))
1015      external        nf_get_var_real
1016
1017      integer         nf_put_var_double
1018!                         (integer             ncid,
1019!                          integer             varid,
1020!                          doubleprecision     dvals(1))
1021      external        nf_put_var_double
1022
1023      integer         nf_get_var_double
1024!                         (integer             ncid,
1025!                          integer             varid,
1026!                          doubleprecision     dvals(1))
1027      external        nf_get_var_double
1028
1029!
1030! single variable put/get routines:
1031!
1032
1033      integer         nf_put_var1_text
1034!                         (integer             ncid,
1035!                          integer             varid,
1036!                          integer             index(1),
1037!                          character*1         text)
1038      external        nf_put_var1_text
1039
1040      integer         nf_get_var1_text
1041!                         (integer             ncid,
1042!                          integer             varid,
1043!                          integer             index(1),
1044!                          character*1         text)
1045      external        nf_get_var1_text
1046
1047      integer         nf_put_var1_int1
1048!                         (integer             ncid,
1049!                          integer             varid,
1050!                          integer             index(1),
1051!                          nf_int1_t           i1val)
1052      external        nf_put_var1_int1
1053
1054      integer         nf_get_var1_int1
1055!                         (integer             ncid,
1056!                          integer             varid,
1057!                          integer             index(1),
1058!                          nf_int1_t           i1val)
1059      external        nf_get_var1_int1
1060
1061      integer         nf_put_var1_int2
1062!                         (integer             ncid,
1063!                          integer             varid,
1064!                          integer             index(1),
1065!                          nf_int2_t           i2val)
1066      external        nf_put_var1_int2
1067
1068      integer         nf_get_var1_int2
1069!                         (integer             ncid,
1070!                          integer             varid,
1071!                          integer             index(1),
1072!                          nf_int2_t           i2val)
1073      external        nf_get_var1_int2
1074
1075      integer         nf_put_var1_int
1076!                         (integer             ncid,
1077!                          integer             varid,
1078!                          integer             index(1),
1079!                          integer             ival)
1080      external        nf_put_var1_int
1081
1082      integer         nf_get_var1_int
1083!                         (integer             ncid,
1084!                          integer             varid,
1085!                          integer             index(1),
1086!                          integer             ival)
1087      external        nf_get_var1_int
1088
1089      integer         nf_put_var1_real
1090!                         (integer             ncid,
1091!                          integer             varid,
1092!                          integer             index(1),
1093!                          real                rval)
1094      external        nf_put_var1_real
1095
1096      integer         nf_get_var1_real
1097!                         (integer             ncid,
1098!                          integer             varid,
1099!                          integer             index(1),
1100!                          real                rval)
1101      external        nf_get_var1_real
1102
1103      integer         nf_put_var1_double
1104!                         (integer             ncid,
1105!                          integer             varid,
1106!                          integer             index(1),
1107!                          doubleprecision     dval)
1108      external        nf_put_var1_double
1109
1110      integer         nf_get_var1_double
1111!                         (integer             ncid,
1112!                          integer             varid,
1113!                          integer             index(1),
1114!                          doubleprecision     dval)
1115      external        nf_get_var1_double
1116
1117!
1118! variable array put/get routines:
1119!
1120
1121      integer         nf_put_vara_text
1122!                         (integer             ncid,
1123!                          integer             varid,
1124!                          integer             start(1),
1125!                          integer             count(1),
1126!                          character(*)        text)
1127      external        nf_put_vara_text
1128
1129      integer         nf_get_vara_text
1130!                         (integer             ncid,
1131!                          integer             varid,
1132!                          integer             start(1),
1133!                          integer             count(1),
1134!                          character(*)        text)
1135      external        nf_get_vara_text
1136
1137      integer         nf_put_vara_int1
1138!                         (integer             ncid,
1139!                          integer             varid,
1140!                          integer             start(1),
1141!                          integer             count(1),
1142!                          nf_int1_t           i1vals(1))
1143      external        nf_put_vara_int1
1144
1145      integer         nf_get_vara_int1
1146!                         (integer             ncid,
1147!                          integer             varid,
1148!                          integer             start(1),
1149!                          integer             count(1),
1150!                          nf_int1_t           i1vals(1))
1151      external        nf_get_vara_int1
1152
1153      integer         nf_put_vara_int2
1154!                         (integer             ncid,
1155!                          integer             varid,
1156!                          integer             start(1),
1157!                          integer             count(1),
1158!                          nf_int2_t           i2vals(1))
1159      external        nf_put_vara_int2
1160
1161      integer         nf_get_vara_int2
1162!                         (integer             ncid,
1163!                          integer             varid,
1164!                          integer             start(1),
1165!                          integer             count(1),
1166!                          nf_int2_t           i2vals(1))
1167      external        nf_get_vara_int2
1168
1169      integer         nf_put_vara_int
1170!                         (integer             ncid,
1171!                          integer             varid,
1172!                          integer             start(1),
1173!                          integer             count(1),
1174!                          integer             ivals(1))
1175      external        nf_put_vara_int
1176
1177      integer         nf_get_vara_int
1178!                         (integer             ncid,
1179!                          integer             varid,
1180!                          integer             start(1),
1181!                          integer             count(1),
1182!                          integer             ivals(1))
1183      external        nf_get_vara_int
1184
1185      integer         nf_put_vara_real
1186!                         (integer             ncid,
1187!                          integer             varid,
1188!                          integer             start(1),
1189!                          integer             count(1),
1190!                          real                rvals(1))
1191      external        nf_put_vara_real
1192
1193      integer         nf_get_vara_real
1194!                         (integer             ncid,
1195!                          integer             varid,
1196!                          integer             start(1),
1197!                          integer             count(1),
1198!                          real                rvals(1))
1199      external        nf_get_vara_real
1200
1201      integer         nf_put_vara_double
1202!                         (integer             ncid,
1203!                          integer             varid,
1204!                          integer             start(1),
1205!                          integer             count(1),
1206!                          doubleprecision     dvals(1))
1207      external        nf_put_vara_double
1208
1209      integer         nf_get_vara_double
1210!                         (integer             ncid,
1211!                          integer             varid,
1212!                          integer             start(1),
1213!                          integer             count(1),
1214!                          doubleprecision     dvals(1))
1215      external        nf_get_vara_double
1216
1217!
1218! strided variable put/get routines:
1219!
1220
1221      integer         nf_put_vars_text
1222!                         (integer             ncid,
1223!                          integer             varid,
1224!                          integer             start(1),
1225!                          integer             count(1),
1226!                          integer             stride(1),
1227!                          character(*)        text)
1228      external        nf_put_vars_text
1229
1230      integer         nf_get_vars_text
1231!                         (integer             ncid,
1232!                          integer             varid,
1233!                          integer             start(1),
1234!                          integer             count(1),
1235!                          integer             stride(1),
1236!                          character(*)        text)
1237      external        nf_get_vars_text
1238
1239      integer         nf_put_vars_int1
1240!                         (integer             ncid,
1241!                          integer             varid,
1242!                          integer             start(1),
1243!                          integer             count(1),
1244!                          integer             stride(1),
1245!                          nf_int1_t           i1vals(1))
1246      external        nf_put_vars_int1
1247
1248      integer         nf_get_vars_int1
1249!                         (integer             ncid,
1250!                          integer             varid,
1251!                          integer             start(1),
1252!                          integer             count(1),
1253!                          integer             stride(1),
1254!                          nf_int1_t           i1vals(1))
1255      external        nf_get_vars_int1
1256
1257      integer         nf_put_vars_int2
1258!                         (integer             ncid,
1259!                          integer             varid,
1260!                          integer             start(1),
1261!                          integer             count(1),
1262!                          integer             stride(1),
1263!                          nf_int2_t           i2vals(1))
1264      external        nf_put_vars_int2
1265
1266      integer         nf_get_vars_int2
1267!                         (integer             ncid,
1268!                          integer             varid,
1269!                          integer             start(1),
1270!                          integer             count(1),
1271!                          integer             stride(1),
1272!                          nf_int2_t           i2vals(1))
1273      external        nf_get_vars_int2
1274
1275      integer         nf_put_vars_int
1276!                         (integer             ncid,
1277!                          integer             varid,
1278!                          integer             start(1),
1279!                          integer             count(1),
1280!                          integer             stride(1),
1281!                          integer             ivals(1))
1282      external        nf_put_vars_int
1283
1284      integer         nf_get_vars_int
1285!                         (integer             ncid,
1286!                          integer             varid,
1287!                          integer             start(1),
1288!                          integer             count(1),
1289!                          integer             stride(1),
1290!                          integer             ivals(1))
1291      external        nf_get_vars_int
1292
1293      integer         nf_put_vars_real
1294!                         (integer             ncid,
1295!                          integer             varid,
1296!                          integer             start(1),
1297!                          integer             count(1),
1298!                          integer             stride(1),
1299!                          real                rvals(1))
1300      external        nf_put_vars_real
1301
1302      integer         nf_get_vars_real
1303!                         (integer             ncid,
1304!                          integer             varid,
1305!                          integer             start(1),
1306!                          integer             count(1),
1307!                          integer             stride(1),
1308!                          real                rvals(1))
1309      external        nf_get_vars_real
1310
1311      integer         nf_put_vars_double
1312!                         (integer             ncid,
1313!                          integer             varid,
1314!                          integer             start(1),
1315!                          integer             count(1),
1316!                          integer             stride(1),
1317!                          doubleprecision     dvals(1))
1318      external        nf_put_vars_double
1319
1320      integer         nf_get_vars_double
1321!                         (integer             ncid,
1322!                          integer             varid,
1323!                          integer             start(1),
1324!                          integer             count(1),
1325!                          integer             stride(1),
1326!                          doubleprecision     dvals(1))
1327      external        nf_get_vars_double
1328
1329!
1330! mapped variable put/get routines:
1331!
1332
1333      integer         nf_put_varm_text
1334!                         (integer             ncid,
1335!                          integer             varid,
1336!                          integer             start(1),
1337!                          integer             count(1),
1338!                          integer             stride(1),
1339!                          integer             imap(1),
1340!                          character(*)        text)
1341      external        nf_put_varm_text
1342
1343      integer         nf_get_varm_text
1344!                         (integer             ncid,
1345!                          integer             varid,
1346!                          integer             start(1),
1347!                          integer             count(1),
1348!                          integer             stride(1),
1349!                          integer             imap(1),
1350!                          character(*)        text)
1351      external        nf_get_varm_text
1352
1353      integer         nf_put_varm_int1
1354!                         (integer             ncid,
1355!                          integer             varid,
1356!                          integer             start(1),
1357!                          integer             count(1),
1358!                          integer             stride(1),
1359!                          integer             imap(1),
1360!                          nf_int1_t           i1vals(1))
1361      external        nf_put_varm_int1
1362
1363      integer         nf_get_varm_int1
1364!                         (integer             ncid,
1365!                          integer             varid,
1366!                          integer             start(1),
1367!                          integer             count(1),
1368!                          integer             stride(1),
1369!                          integer             imap(1),
1370!                          nf_int1_t           i1vals(1))
1371      external        nf_get_varm_int1
1372
1373      integer         nf_put_varm_int2
1374!                         (integer             ncid,
1375!                          integer             varid,
1376!                          integer             start(1),
1377!                          integer             count(1),
1378!                          integer             stride(1),
1379!                          integer             imap(1),
1380!                          nf_int2_t           i2vals(1))
1381      external        nf_put_varm_int2
1382
1383      integer         nf_get_varm_int2
1384!                         (integer             ncid,
1385!                          integer             varid,
1386!                          integer             start(1),
1387!                          integer             count(1),
1388!                          integer             stride(1),
1389!                          integer             imap(1),
1390!                          nf_int2_t           i2vals(1))
1391      external        nf_get_varm_int2
1392
1393      integer         nf_put_varm_int
1394!                         (integer             ncid,
1395!                          integer             varid,
1396!                          integer             start(1),
1397!                          integer             count(1),
1398!                          integer             stride(1),
1399!                          integer             imap(1),
1400!                          integer             ivals(1))
1401      external        nf_put_varm_int
1402
1403      integer         nf_get_varm_int
1404!                         (integer             ncid,
1405!                          integer             varid,
1406!                          integer             start(1),
1407!                          integer             count(1),
1408!                          integer             stride(1),
1409!                          integer             imap(1),
1410!                          integer             ivals(1))
1411      external        nf_get_varm_int
1412
1413      integer         nf_put_varm_real
1414!                         (integer             ncid,
1415!                          integer             varid,
1416!                          integer             start(1),
1417!                          integer             count(1),
1418!                          integer             stride(1),
1419!                          integer             imap(1),
1420!                          real                rvals(1))
1421      external        nf_put_varm_real
1422
1423      integer         nf_get_varm_real
1424!                         (integer             ncid,
1425!                          integer             varid,
1426!                          integer             start(1),
1427!                          integer             count(1),
1428!                          integer             stride(1),
1429!                          integer             imap(1),
1430!                          real                rvals(1))
1431      external        nf_get_varm_real
1432
1433      integer         nf_put_varm_double
1434!                         (integer             ncid,
1435!                          integer             varid,
1436!                          integer             start(1),
1437!                          integer             count(1),
1438!                          integer             stride(1),
1439!                          integer             imap(1),
1440!                          doubleprecision     dvals(1))
1441      external        nf_put_varm_double
1442
1443      integer         nf_get_varm_double
1444!                         (integer             ncid,
1445!                          integer             varid,
1446!                          integer             start(1),
1447!                          integer             count(1),
1448!                          integer             stride(1),
1449!                          integer             imap(1),
1450!                          doubleprecision     dvals(1))
1451      external        nf_get_varm_double
1452
1453
1454!     NetCDF-2.
1455!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1456! begin netcdf 2.4 backward compatibility:
1457!
1458
1459!     
1460! functions in the fortran interface
1461!
1462      integer nccre
1463      integer ncopn
1464      integer ncddef
1465      integer ncdid
1466      integer ncvdef
1467      integer ncvid
1468      integer nctlen
1469      integer ncsfil
1470
1471      external nccre
1472      external ncopn
1473      external ncddef
1474      external ncdid
1475      external ncvdef
1476      external ncvid
1477      external nctlen
1478      external ncsfil
1479
1480
1481      integer ncrdwr
1482      integer nccreat
1483      integer ncexcl
1484      integer ncindef
1485      integer ncnsync
1486      integer nchsync
1487      integer ncndirty
1488      integer nchdirty
1489      integer nclink
1490      integer ncnowrit
1491      integer ncwrite
1492      integer ncclob
1493      integer ncnoclob
1494      integer ncglobal
1495      integer ncfill
1496      integer ncnofill
1497      integer maxncop
1498      integer maxncdim
1499      integer maxncatt
1500      integer maxncvar
1501      integer maxncnam
1502      integer maxvdims
1503      integer ncnoerr
1504      integer ncebadid
1505      integer ncenfile
1506      integer nceexist
1507      integer nceinval
1508      integer nceperm
1509      integer ncenotin
1510      integer nceindef
1511      integer ncecoord
1512      integer ncemaxds
1513      integer ncename
1514      integer ncenoatt
1515      integer ncemaxat
1516      integer ncebadty
1517      integer ncebadd
1518      integer ncests
1519      integer nceunlim
1520      integer ncemaxvs
1521      integer ncenotvr
1522      integer nceglob
1523      integer ncenotnc
1524      integer ncfoobar
1525      integer ncsyserr
1526      integer ncfatal
1527      integer ncverbos
1528      integer ncentool
1529
1530
1531!
1532! netcdf data types:
1533!
1534      integer ncbyte
1535      integer ncchar
1536      integer ncshort
1537      integer nclong
1538      integer ncfloat
1539      integer ncdouble
1540
1541      parameter(ncbyte = 1)
1542      parameter(ncchar = 2)
1543      parameter(ncshort = 3)
1544      parameter(nclong = 4)
1545      parameter(ncfloat = 5)
1546      parameter(ncdouble = 6)
1547
1548!     
1549!     masks for the struct nc flag field; passed in as 'mode' arg to
1550!     nccreate and ncopen.
1551!     
1552
1553!     read/write, 0 => readonly
1554      parameter(ncrdwr = 1)
1555!     in create phase, cleared by ncendef
1556      parameter(nccreat = 2)
1557!     on create destroy existing file
1558      parameter(ncexcl = 4)
1559!     in define mode, cleared by ncendef
1560      parameter(ncindef = 8)
1561!     synchronise numrecs on change (x'10')
1562      parameter(ncnsync = 16)
1563!     synchronise whole header on change (x'20')
1564      parameter(nchsync = 32)
1565!     numrecs has changed (x'40')
1566      parameter(ncndirty = 64) 
1567!     header info has changed (x'80')
1568      parameter(nchdirty = 128)
1569!     prefill vars on endef and increase of record, the default behavior
1570      parameter(ncfill = 0)
1571!     do not fill vars on endef and increase of record (x'100')
1572      parameter(ncnofill = 256)
1573!     isa link (x'8000')
1574      parameter(nclink = 32768)
1575
1576!     
1577!     'mode' arguments for nccreate and ncopen
1578!     
1579      parameter(ncnowrit = 0)
1580      parameter(ncwrite = ncrdwr)
1581      parameter(ncclob = nf_clobber)
1582      parameter(ncnoclob = nf_noclobber)
1583
1584!     
1585!     'size' argument to ncdimdef for an unlimited dimension
1586!     
1587      integer ncunlim
1588      parameter(ncunlim = 0)
1589
1590!     
1591!     attribute id to put/get a global attribute
1592!     
1593      parameter(ncglobal  = 0)
1594
1595!     
1596!     advisory maximums:
1597!     
1598      parameter(maxncop = 64)
1599      parameter(maxncdim = 1024)
1600      parameter(maxncatt = 8192)
1601      parameter(maxncvar = 8192)
1602!     not enforced
1603      parameter(maxncnam = 256)
1604      parameter(maxvdims = maxncdim)
1605
1606!     
1607!     global netcdf error status variable
1608!     initialized in error.c
1609!     
1610
1611!     no error
1612      parameter(ncnoerr = nf_noerr)
1613!     not a netcdf id
1614      parameter(ncebadid = nf_ebadid)
1615!     too many netcdfs open
1616      parameter(ncenfile = -31)   ! nc_syserr
1617!     netcdf file exists && ncnoclob
1618      parameter(nceexist = nf_eexist)
1619!     invalid argument
1620      parameter(nceinval = nf_einval)
1621!     write to read only
1622      parameter(nceperm = nf_eperm)
1623!     operation not allowed in data mode
1624      parameter(ncenotin = nf_enotindefine )   
1625!     operation not allowed in define mode
1626      parameter(nceindef = nf_eindefine)   
1627!     coordinates out of domain
1628      parameter(ncecoord = nf_einvalcoords)
1629!     maxncdims exceeded
1630      parameter(ncemaxds = nf_emaxdims)
1631!     string match to name in use
1632      parameter(ncename = nf_enameinuse)   
1633!     attribute not found
1634      parameter(ncenoatt = nf_enotatt)
1635!     maxncattrs exceeded
1636      parameter(ncemaxat = nf_emaxatts)
1637!     not a netcdf data type
1638      parameter(ncebadty = nf_ebadtype)
1639!     invalid dimension id
1640      parameter(ncebadd = nf_ebaddim)
1641!     ncunlimited in the wrong index
1642      parameter(nceunlim = nf_eunlimpos)
1643!     maxncvars exceeded
1644      parameter(ncemaxvs = nf_emaxvars)
1645!     variable not found
1646      parameter(ncenotvr = nf_enotvar)
1647!     action prohibited on ncglobal varid
1648      parameter(nceglob = nf_eglobal)
1649!     not a netcdf file
1650      parameter(ncenotnc = nf_enotnc)
1651      parameter(ncests = nf_ests)
1652      parameter (ncentool = nf_emaxname) 
1653      parameter(ncfoobar = 32)
1654      parameter(ncsyserr = -31)
1655
1656!     
1657!     global options variable. used to determine behavior of error handler.
1658!     initialized in lerror.c
1659!     
1660      parameter(ncfatal = 1)
1661      parameter(ncverbos = 2)
1662
1663!
1664!     default fill values.  these must be the same as in the c interface.
1665!
1666      integer filbyte
1667      integer filchar
1668      integer filshort
1669      integer fillong
1670      real filfloat
1671      doubleprecision fildoub
1672
1673      parameter (filbyte = -127)
1674      parameter (filchar = 0)
1675      parameter (filshort = -32767)
1676      parameter (fillong = -2147483647)
1677      parameter (filfloat = 9.9692099683868690e+36)
1678      parameter (fildoub = 9.9692099683868690d+36)
1679
1680!     NetCDF-4.
1681!     This is part of netCDF-4. Copyright 2006, UCAR, See COPYRIGHT
1682!     file for distribution information.
1683
1684!     Netcdf version 4 fortran interface.
1685
1686!     $Id: netcdf4.inc,v 1.28 2010/05/25 13:53:02 ed Exp $
1687
1688!     New netCDF-4 types.
1689      integer nf_ubyte
1690      integer nf_ushort
1691      integer nf_uint
1692      integer nf_int64
1693      integer nf_uint64
1694      integer nf_string
1695      integer nf_vlen
1696      integer nf_opaque
1697      integer nf_enum
1698      integer nf_compound
1699
1700      parameter (nf_ubyte = 7)
1701      parameter (nf_ushort = 8)
1702      parameter (nf_uint = 9)
1703      parameter (nf_int64 = 10)
1704      parameter (nf_uint64 = 11)
1705      parameter (nf_string = 12)
1706      parameter (nf_vlen = 13)
1707      parameter (nf_opaque = 14)
1708      parameter (nf_enum = 15)
1709      parameter (nf_compound = 16)
1710
1711!     New netCDF-4 fill values.
1712      integer           nf_fill_ubyte
1713      integer           nf_fill_ushort
1714!      real              nf_fill_uint
1715!      real              nf_fill_int64
1716!      real              nf_fill_uint64
1717      parameter (nf_fill_ubyte = 255)
1718      parameter (nf_fill_ushort = 65535)
1719
1720!     New constants.
1721      integer nf_format_netcdf4
1722      parameter (nf_format_netcdf4 = 3)
1723
1724      integer nf_format_netcdf4_classic
1725      parameter (nf_format_netcdf4_classic = 4)
1726
1727      integer nf_netcdf4
1728      parameter (nf_netcdf4 = 4096)
1729
1730      integer nf_classic_model
1731      parameter (nf_classic_model = 256)
1732
1733      integer nf_chunk_seq
1734      parameter (nf_chunk_seq = 0)
1735      integer nf_chunk_sub
1736      parameter (nf_chunk_sub = 1)
1737      integer nf_chunk_sizes
1738      parameter (nf_chunk_sizes = 2)
1739
1740      integer nf_endian_native
1741      parameter (nf_endian_native = 0)
1742      integer nf_endian_little
1743      parameter (nf_endian_little = 1)
1744      integer nf_endian_big
1745      parameter (nf_endian_big = 2)
1746
1747!     For NF_DEF_VAR_CHUNKING
1748      integer nf_chunked
1749      parameter (nf_chunked = 0)
1750      integer nf_contiguous
1751      parameter (nf_contiguous = 1)
1752
1753!     For NF_DEF_VAR_FLETCHER32
1754      integer nf_nochecksum
1755      parameter (nf_nochecksum = 0)
1756      integer nf_fletcher32
1757      parameter (nf_fletcher32 = 1)
1758
1759!     For NF_DEF_VAR_DEFLATE
1760      integer nf_noshuffle
1761      parameter (nf_noshuffle = 0)
1762      integer nf_shuffle
1763      parameter (nf_shuffle = 1)
1764
1765!     For NF_DEF_VAR_SZIP
1766      integer nf_szip_ec_option_mask
1767      parameter (nf_szip_ec_option_mask = 4)
1768      integer nf_szip_nn_option_mask
1769      parameter (nf_szip_nn_option_mask = 32)
1770
1771!     For parallel I/O.
1772      integer nf_mpiio     
1773      parameter (nf_mpiio = 8192)
1774      integer nf_mpiposix
1775      parameter (nf_mpiposix = 16384)
1776      integer nf_pnetcdf
1777      parameter (nf_pnetcdf = 32768)
1778
1779!     For NF_VAR_PAR_ACCESS.
1780      integer nf_independent
1781      parameter (nf_independent = 0)
1782      integer nf_collective
1783      parameter (nf_collective = 1)
1784
1785!     New error codes.
1786      integer nf_ehdferr        ! Error at HDF5 layer.
1787      parameter (nf_ehdferr = -101)
1788      integer nf_ecantread      ! Can't read.
1789      parameter (nf_ecantread = -102)
1790      integer nf_ecantwrite     ! Can't write.
1791      parameter (nf_ecantwrite = -103)
1792      integer nf_ecantcreate    ! Can't create.
1793      parameter (nf_ecantcreate = -104)
1794      integer nf_efilemeta      ! Problem with file metadata.
1795      parameter (nf_efilemeta = -105)
1796      integer nf_edimmeta       ! Problem with dimension metadata.
1797      parameter (nf_edimmeta = -106)
1798      integer nf_eattmeta       ! Problem with attribute metadata.
1799      parameter (nf_eattmeta = -107)
1800      integer nf_evarmeta       ! Problem with variable metadata.
1801      parameter (nf_evarmeta = -108)
1802      integer nf_enocompound    ! Not a compound type.
1803      parameter (nf_enocompound = -109)
1804      integer nf_eattexists     ! Attribute already exists.
1805      parameter (nf_eattexists = -110)
1806      integer nf_enotnc4        ! Attempting netcdf-4 operation on netcdf-3 file.   
1807      parameter (nf_enotnc4 = -111)
1808      integer nf_estrictnc3     ! Attempting netcdf-4 operation on strict nc3 netcdf-4 file.   
1809      parameter (nf_estrictnc3 = -112)
1810      integer nf_enotnc3        ! Attempting netcdf-3 operation on netcdf-4 file.   
1811      parameter (nf_enotnc3 = -113)
1812      integer nf_enopar         ! Parallel operation on file opened for non-parallel access.   
1813      parameter (nf_enopar = -114)
1814      integer nf_eparinit       ! Error initializing for parallel access.   
1815      parameter (nf_eparinit = -115)
1816      integer nf_ebadgrpid      ! Bad group ID.   
1817      parameter (nf_ebadgrpid = -116)
1818      integer nf_ebadtypid      ! Bad type ID.   
1819      parameter (nf_ebadtypid = -117)
1820      integer nf_etypdefined    ! Type has already been defined and may not be edited.
1821      parameter (nf_etypdefined = -118)
1822      integer nf_ebadfield      ! Bad field ID.   
1823      parameter (nf_ebadfield = -119)
1824      integer nf_ebadclass      ! Bad class.   
1825      parameter (nf_ebadclass = -120)
1826      integer nf_emaptype       ! Mapped access for atomic types only.   
1827      parameter (nf_emaptype = -121)
1828      integer nf_elatefill      ! Attempt to define fill value when data already exists.
1829      parameter (nf_elatefill = -122)
1830      integer nf_elatedef       ! Attempt to define var properties, like deflate, after enddef.
1831      parameter (nf_elatedef = -123)
1832      integer nf_edimscale      ! Probem with HDF5 dimscales.
1833      parameter (nf_edimscale = -124)
1834      integer nf_enogrp       ! No group found.
1835      parameter (nf_enogrp = -125)
1836
1837
1838!     New functions.
1839
1840!     Parallel I/O.
1841      integer nf_create_par
1842      external nf_create_par
1843
1844      integer nf_open_par
1845      external nf_open_par
1846
1847      integer nf_var_par_access
1848      external nf_var_par_access
1849
1850!     Functions to handle groups.
1851      integer nf_inq_ncid
1852      external nf_inq_ncid
1853
1854      integer nf_inq_grps
1855      external nf_inq_grps
1856
1857      integer nf_inq_grpname
1858      external nf_inq_grpname
1859
1860      integer nf_inq_grpname_full
1861      external nf_inq_grpname_full
1862
1863      integer nf_inq_grpname_len
1864      external nf_inq_grpname_len
1865
1866      integer nf_inq_grp_parent
1867      external nf_inq_grp_parent
1868
1869      integer nf_inq_grp_ncid
1870      external nf_inq_grp_ncid
1871
1872      integer nf_inq_grp_full_ncid
1873      external nf_inq_grp_full_ncid
1874
1875      integer nf_inq_varids
1876      external nf_inq_varids
1877
1878      integer nf_inq_dimids
1879      external nf_inq_dimids
1880
1881      integer nf_def_grp
1882      external nf_def_grp
1883
1884!     New options for netCDF variables.
1885      integer nf_def_var_deflate
1886      external nf_def_var_deflate
1887
1888      integer nf_inq_var_deflate
1889      external nf_inq_var_deflate
1890
1891      integer nf_def_var_fletcher32
1892      external nf_def_var_fletcher32
1893
1894      integer nf_inq_var_fletcher32
1895      external nf_inq_var_fletcher32
1896
1897      integer nf_def_var_chunking
1898      external nf_def_var_chunking
1899
1900      integer nf_inq_var_chunking
1901      external nf_inq_var_chunking
1902
1903      integer nf_def_var_fill
1904      external nf_def_var_fill
1905
1906      integer nf_inq_var_fill
1907      external nf_inq_var_fill
1908
1909      integer nf_def_var_endian
1910      external nf_def_var_endian
1911
1912      integer nf_inq_var_endian
1913      external nf_inq_var_endian
1914
1915!     User defined types.
1916      integer nf_inq_typeids
1917      external nf_inq_typeids
1918
1919      integer nf_inq_typeid
1920      external nf_inq_typeid
1921
1922      integer nf_inq_type
1923      external nf_inq_type
1924
1925      integer nf_inq_user_type
1926      external nf_inq_user_type
1927
1928!     User defined types - compound types.
1929      integer nf_def_compound
1930      external nf_def_compound
1931
1932      integer nf_insert_compound
1933      external nf_insert_compound
1934
1935      integer nf_insert_array_compound
1936      external nf_insert_array_compound
1937
1938      integer nf_inq_compound
1939      external nf_inq_compound
1940
1941      integer nf_inq_compound_name
1942      external nf_inq_compound_name
1943
1944      integer nf_inq_compound_size
1945      external nf_inq_compound_size
1946
1947      integer nf_inq_compound_nfields
1948      external nf_inq_compound_nfields
1949
1950      integer nf_inq_compound_field
1951      external nf_inq_compound_field
1952
1953      integer nf_inq_compound_fieldname
1954      external nf_inq_compound_fieldname
1955
1956      integer nf_inq_compound_fieldindex
1957      external nf_inq_compound_fieldindex
1958
1959      integer nf_inq_compound_fieldoffset
1960      external nf_inq_compound_fieldoffset
1961
1962      integer nf_inq_compound_fieldtype
1963      external nf_inq_compound_fieldtype
1964
1965      integer nf_inq_compound_fieldndims
1966      external nf_inq_compound_fieldndims
1967
1968      integer nf_inq_compound_fielddim_sizes
1969      external nf_inq_compound_fielddim_sizes
1970
1971!     User defined types - variable length arrays.
1972      integer nf_def_vlen
1973      external nf_def_vlen
1974
1975      integer nf_inq_vlen
1976      external nf_inq_vlen
1977
1978      integer nf_free_vlen
1979      external nf_free_vlen
1980
1981!     User defined types - enums.
1982      integer nf_def_enum
1983      external nf_def_enum
1984
1985      integer nf_insert_enum
1986      external nf_insert_enum
1987
1988      integer nf_inq_enum
1989      external nf_inq_enum
1990
1991      integer nf_inq_enum_member
1992      external nf_inq_enum_member
1993
1994      integer nf_inq_enum_ident
1995      external nf_inq_enum_ident
1996
1997!     User defined types - opaque.
1998      integer nf_def_opaque
1999      external nf_def_opaque
2000
2001      integer nf_inq_opaque
2002      external nf_inq_opaque
2003
2004!     Write and read attributes of any type, including user defined
2005!     types.
2006      integer nf_put_att
2007      external nf_put_att
2008      integer nf_get_att
2009      external nf_get_att
2010
2011!     Write and read variables of any type, including user defined
2012!     types.
2013      integer nf_put_var
2014      external nf_put_var
2015      integer nf_put_var1
2016      external nf_put_var1
2017      integer nf_put_vara
2018      external nf_put_vara
2019      integer nf_put_vars
2020      external nf_put_vars
2021      integer nf_get_var
2022      external nf_get_var
2023      integer nf_get_var1
2024      external nf_get_var1
2025      integer nf_get_vara
2026      external nf_get_vara
2027      integer nf_get_vars
2028      external nf_get_vars
2029
2030!     64-bit int functions.
2031      integer nf_put_var1_int64
2032      external nf_put_var1_int64
2033      integer nf_put_vara_int64
2034      external nf_put_vara_int64
2035      integer nf_put_vars_int64
2036      external nf_put_vars_int64
2037      integer nf_put_varm_int64
2038      external nf_put_varm_int64
2039      integer nf_put_var_int64
2040      external nf_put_var_int64
2041      integer nf_get_var1_int64
2042      external nf_get_var1_int64
2043      integer nf_get_vara_int64
2044      external nf_get_vara_int64
2045      integer nf_get_vars_int64
2046      external nf_get_vars_int64
2047      integer nf_get_varm_int64
2048      external nf_get_varm_int64
2049      integer nf_get_var_int64
2050      external nf_get_var_int64
2051
2052!     For helping F77 users with VLENs.
2053      integer nf_get_vlen_element
2054      external nf_get_vlen_element
2055      integer nf_put_vlen_element
2056      external nf_put_vlen_element
2057
2058!     For dealing with file level chunk cache.
2059      integer nf_set_chunk_cache
2060      external nf_set_chunk_cache
2061      integer nf_get_chunk_cache
2062      external nf_get_chunk_cache
2063
2064!     For dealing with per variable chunk cache.
2065      integer nf_set_var_chunk_cache
2066      external nf_set_var_chunk_cache
2067      integer nf_get_var_chunk_cache
2068      external nf_get_var_chunk_cache
2069!#include "advtrac.h"
2070c=======================================================================
2071c   Declarations
2072c=======================================================================
2073
2074c Variables dimension du fichier "start_archive"
2075c------------------------------------
2076      CHARACTER relief*3
2077
2078
2079c Variables pour les lectures NetCDF des fichiers "start_archive" 
2080c--------------------------------------------------
2081      INTEGER nid_dyn, nid_fi,nid,nvarid
2082      INTEGER length
2083      parameter (length = 100)
2084      INTEGER tab0
2085      INTEGER   NB_ETATMAX
2086      parameter (NB_ETATMAX = 100)
2087
2088      REAL  date
2089      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
2090
2091c Variable histoire 
2092c------------------
2093      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
2094      REAL phis(iip1,jjp1)
2095      REAL,ALLOCATABLE :: q(:,:,:,:)               ! champs advectes
2096
2097c autre variables dynamique nouvelle grille
2098c------------------------------------------
2099      REAL pks(iip1,jjp1)
2100      REAL w(iip1,jjp1,llm+1)
2101      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
2102!      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
2103!      REAL dh(ip1jmp1,llm),dp(ip1jmp1)
2104      REAL phi(iip1,jjp1,llm)
2105
2106      integer klatdat,klongdat
2107      PARAMETER (klatdat=180,klongdat=360)
2108
2109c Physique sur grille scalaire 
2110c----------------------------
2111      real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
2112      real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
2113
2114c variable physique
2115c------------------
2116      REAL tsurf(ngridmx)       ! surface temperature
2117      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
2118!      REAL co2ice(ngridmx)     ! CO2 ice layer
2119      REAL emis(ngridmx)        ! surface emissivity
2120      real emisread             ! added by RW
2121      REAL,ALLOCATABLE :: qsurf(:,:)
2122      REAL q2(ngridmx,nlayermx+1)
2123!      REAL rnaturfi(ngridmx)
2124      real alb(iip1,jjp1),albfi(ngridmx) ! albedos
2125      real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
2126      real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
2127      REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
2128
2129      INTEGER i,j,l,isoil,ig,idum
2130      real mugaz ! molar mass of the atmosphere
2131
2132      integer ierr 
2133
2134      REAL rnat(ngridmx)
2135      REAL tslab(ngridmx,nsoilmx) ! slab ocean temperature
2136      REAL pctsrf_sic(ngridmx) ! sea ice cover
2137      REAL tsea_ice(ngridmx) ! temperature sea_ice
2138      REAL sea_ice(ngridmx) ! mass of sea ice (kg/m2)
2139
2140c Variables on the new grid along scalar points 
2141c------------------------------------------------------
2142!      REAL p(iip1,jjp1)
2143      REAL t(iip1,jjp1,llm)
2144      REAL tset(iip1,jjp1,llm)
2145      real phisold_newgrid(iip1,jjp1)
2146      REAL :: teta(iip1, jjp1, llm)
2147      REAL :: pk(iip1,jjp1,llm)
2148      REAL :: pkf(iip1,jjp1,llm)
2149      REAL :: ps(iip1, jjp1)
2150      REAL :: masse(iip1,jjp1,llm)
2151      REAL :: xpn,xps,xppn(iim),xpps(iim)
2152      REAL :: p3d(iip1, jjp1, llm+1)
2153      REAL :: beta(iip1,jjp1,llm)
2154!      REAL dteta(ip1jmp1,llm)
2155
2156c Variable de l'ancienne grille
2157c------------------------------
2158      real time
2159      real tab_cntrl(100)
2160      real tab_cntrl_bis(100)
2161
2162c variables diverses
2163c-------------------
2164      real choix_1,pp
2165      character*80      fichnom
2166      character*250  filestring
2167      integer Lmodif,iq,thermo
2168      character modif*20
2169      real z_reel(iip1,jjp1)
2170      real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
2171      real ptoto,pcap,patm,airetot,ptotn,patmn,psea
2172!      real ssum
2173      character*1 yes
2174      logical :: flagtset=.false. ,  flagps0=.false.
2175      real val, val2, val3 ! to store temporary variables
2176      real :: iceith=2000 ! thermal inertia of subterranean ice
2177      integer iref,jref
2178
2179      INTEGER :: itau
2180     
2181      INTEGER :: nq,numvanle
2182      character(len=20) :: txt ! to store some text
2183      character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
2184      character(len=150) :: longline
2185      integer :: count
2186      real :: profile(llm+1) ! to store an atmospheric profile + surface value
2187
2188!     added by RW for test
2189      real pmean, phi0
2190
2191!     added by BC for equilibrium temperature startup
2192      real teque
2193
2194!     added by BC for cloud fraction setup
2195      REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
2196      REAL totalfrac(ngridmx)
2197
2198
2199!     added by RW for nuketharsis
2200      real fact1
2201      real fact2
2202
2203
2204c sortie visu pour les champs dynamiques
2205c---------------------------------------
2206!      INTEGER :: visuid
2207!      real :: time_step,t_ops,t_wrt
2208!      CHARACTER*80 :: visu_file
2209
2210      cpp    = 0.
2211      preff  = 0.
2212      pa     = 0. ! to ensure disaster rather than minor error if we don`t
2213                  ! make deliberate choice of these values elsewhere.
2214
2215! initialize "serial/parallel" related stuff
2216      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
2217      call initcomgeomphy
2218! Load tracer number and names:
2219      call iniadvtrac(nqtot,numvanle)
2220! allocate arrays
2221      allocate(q(iip1,jjp1,llm,nqtot))
2222      allocate(qsurf(ngridmx,nqtot))
2223
2224c=======================================================================
2225c   Choice of the start file(s) to use
2226c=======================================================================
2227      write(*,*) 'From which kind of files do you want to create new',
2228     .  'start and startfi files'
2229      write(*,*) '    0 - from a file start_archive'
2230      write(*,*) '    1 - from files start and startfi'
2231 
2232c-----------------------------------------------------------------------
2233c   Open file(s) to modify (start or start_archive)
2234c-----------------------------------------------------------------------
2235
2236      DO
2237         read(*,*,iostat=ierr) choix_1
2238         if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
2239      ENDDO
2240
2241c     Open start_archive
2242c     ~~~~~~~~~~~~~~~~~~~~~~~~~~
2243      if (choix_1.eq.0) then
2244
2245        write(*,*) 'Creating start files from:'
2246        write(*,*) './start_archive.nc'
2247        write(*,*)
2248        fichnom = 'start_archive.nc'
2249        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
2250        IF (ierr.NE.NF_NOERR) THEN
2251          write(6,*)' Problem opening file:',fichnom
2252          write(6,*)' ierr = ', ierr
2253          CALL ABORT
2254        ENDIF
2255        tab0 = 50
2256        Lmodif = 1
2257
2258c     OR open start and startfi files
2259c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2260      else
2261        write(*,*) 'Creating start files from:'
2262        write(*,*) './start.nc and ./startfi.nc'
2263        write(*,*)
2264        fichnom = 'start.nc'
2265        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
2266        IF (ierr.NE.NF_NOERR) THEN
2267          write(6,*)' Problem opening file:',fichnom
2268          write(6,*)' ierr = ', ierr
2269          CALL ABORT
2270        ENDIF
2271 
2272        fichnom = 'startfi.nc'
2273        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
2274        IF (ierr.NE.NF_NOERR) THEN
2275          write(6,*)' Problem opening file:',fichnom
2276          write(6,*)' ierr = ', ierr
2277          CALL ABORT
2278        ENDIF
2279
2280        tab0 = 0
2281        Lmodif = 0
2282
2283      endif
2284
2285
2286c=======================================================================
2287c  INITIALISATIONS DIVERSES
2288c=======================================================================
2289
2290! Initialize global tracer indexes (stored in tracer.h)
2291! ... this has to be done before phyetat0
2292      call initracer(ngridmx,nqtot,tname)
2293
2294! Take care of arrays in common modules
2295      ! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
2296      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
2297      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
2298      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
2299      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
2300      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
2301      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
2302      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
2303      ! ALLOCATE ARRAYS in comsoil_h (if not already done)
2304      IF (.not.ALLOCATED(layer))
2305     . ALLOCATE(layer(nsoilmx))
2306      IF (.not.ALLOCATED(mlayer))
2307     . ALLOCATE(mlayer(0:nsoilmx-1))
2308      IF (.not.ALLOCATED(inertiedat))
2309     . ALLOCATE(inertiedat(ngridmx,nsoilmx))
2310      ! ALLOCATE ARRAYS IN comgeomfi_h (done in inifis usually)
2311      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngridmx))
2312      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngridmx))
2313      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngridmx))
2314
2315c-----------------------------------------------------------------------
2316c Lecture du tableau des parametres du run (pour la dynamique)
2317c-----------------------------------------------------------------------
2318      if (choix_1.eq.0) then
2319
2320        write(*,*) 'reading tab_cntrl START_ARCHIVE'
2321c
2322        ierr = NF_INQ_VARID (nid, "controle", nvarid)
2323        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
2324c
2325      else if (choix_1.eq.1) then
2326
2327        write(*,*) 'reading tab_cntrl START'
2328c
2329        ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
2330        ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
2331c
2332        write(*,*) 'reading tab_cntrl STARTFI'
2333c
2334        ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
2335        ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
2336c
2337        do i=1,50
2338          tab_cntrl(i+50)=tab_cntrl_bis(i)
2339        enddo
2340        write(*,*) 'printing tab_cntrl', tab_cntrl
2341        do i=1,100
2342          write(*,*) i,tab_cntrl(i)
2343        enddo
2344       
2345        ! Lmodif set to 0 to disable modifications possibility in phyeta0                           
2346        write(*,*) 'Reading file START'
2347        fichnom = 'start.nc'
2348        CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
2349     .       ps,phis,time)
2350
2351        write(*,*) 'Reading file STARTFI'
2352        fichnom = 'startfi.nc'
2353        CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqtot,
2354     .        day_ini,time,
2355     .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
2356     .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
2357     .        sea_ice)
2358
2359        ! copy albedo and soil thermal inertia on (local) physics grid
2360        do i=1,ngridmx
2361          albfi(i) = albedodat(i)
2362          do j=1,nsoilmx
2363           ithfi(i,j) = inertiedat(i,j)
2364          enddo
2365        ! build a surfithfi(:) using 1st layer of ithfi(:), which might
2366        ! be needed later on if reinitializing soil thermal inertia
2367          surfithfi(i)=ithfi(i,1)
2368        enddo
2369        ! also copy albedo and soil thermal inertia on (local) dynamics grid
2370        ! so that options below can manipulate either (but must then ensure
2371        ! to correctly recast things on physics grid)
2372        call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
2373        call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
2374        call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
2375     
2376      endif
2377c-----------------------------------------------------------------------
2378c               Initialisation des constantes dynamique
2379c-----------------------------------------------------------------------
2380
2381      kappa = tab_cntrl(9)
2382      etot0 = tab_cntrl(12)
2383      ptot0 = tab_cntrl(13)
2384      ztot0 = tab_cntrl(14)
2385      stot0 = tab_cntrl(15)
2386      ang0 = tab_cntrl(16)
2387      write(*,*) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
2388      write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0
2389
2390      ! for vertical coordinate
2391      preff=tab_cntrl(18) ! reference surface pressure
2392      pa=tab_cntrl(17)  ! reference pressure at which coord is purely pressure
2393      !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
2394      write(*,*) "Newstart: preff=",preff," pa=",pa
2395      yes=' '
2396      do while ((yes.ne.'y').and.(yes.ne.'n'))
2397        write(*,*) "Change the values of preff and pa? (y/n)"
2398        read(*,fmt='(a)') yes
2399      end do
2400      if (yes.eq.'y') then
2401        write(*,*)"New value of reference surface pressure preff?"
2402        write(*,*)"   (for Mars, typically preff=610)"
2403        read(*,*) preff
2404        write(*,*)"New value of reference pressure pa for purely"
2405        write(*,*)"pressure levels (for hybrid coordinates)?"
2406        write(*,*)"   (for Mars, typically pa=20)"
2407        read(*,*) pa
2408      endif
2409c-----------------------------------------------------------------------
2410c   Lecture du tab_cntrl et initialisation des constantes physiques
2411c  - pour start:  Lmodif = 0 => pas de modifications possibles
2412c                  (modif dans le tabfi de readfi + loin)
2413c  - pour start_archive:  Lmodif = 1 => modifications possibles
2414c-----------------------------------------------------------------------
2415      if (choix_1.eq.0) then
2416         ! tabfi requires that input file be first opened by open_startphy(fichnom)
2417         fichnom = 'start_archive.nc'
2418         call open_startphy(fichnom)
2419         call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
2420     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
2421      else if (choix_1.eq.1) then
2422         fichnom = 'startfi.nc'
2423         call open_startphy(fichnom)
2424         Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
2425         call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
2426     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
2427      endif
2428
2429      rad = p_rad
2430      omeg = p_omeg
2431      g = p_g
2432      cpp = p_cpp
2433      mugaz = p_mugaz
2434      daysec = p_daysec
2435
2436      kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW
2437
2438c=======================================================================
2439c  INITIALISATIONS DIVERSES
2440c=======================================================================
2441! Load parameters from run.def file
2442      CALL defrun_new( 99, .TRUE. )
2443      CALL iniconst
2444      CALL inigeom
2445      idum=-1
2446      idum=0
2447
2448c Initialisation coordonnees /aires
2449c -------------------------------
2450! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
2451!       rlatu() and rlonv() are given in radians
2452      latfi(1)=rlatu(1)
2453      lonfi(1)=0.
2454      DO j=2,jjm
2455         DO i=1,iim
2456            latfi((j-2)*iim+1+i)=rlatu(j)
2457            lonfi((j-2)*iim+1+i)=rlonv(i)
2458         ENDDO
2459      ENDDO
2460      latfi(ngridmx)=rlatu(jjp1)
2461      lonfi(ngridmx)=0.
2462       
2463      ! build airefi(), mesh area on physics grid
2464      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
2465      ! Poles are single points on physics grid
2466      airefi(1)=airefi(1)*iim
2467      airefi(ngridmx)=airefi(ngridmx)*iim
2468
2469! also initialize various physics flags/settings which might be needed
2470!    (for instance initracer needs to know about some flags, and/or
2471!      'datafile' path may be changed by user)
2472      call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
2473     &                latfi,lonfi,airefi,rad,g,r,cpp)
2474
2475c=======================================================================
2476c   lecture topographie, albedo, inertie thermique, relief sous-maille
2477c=======================================================================
2478
2479      if (choix_1.eq.0) then  ! for start_archive files,
2480                              ! where an external "surface.nc" file is needed
2481
2482c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
2483c       write(*,*)
2484c       write(*,*) 'choix du relief (mola,pla)'
2485c       write(*,*) '(Topographie MGS MOLA, plat)'
2486c       read(*,fmt='(a3)') relief
2487        relief="mola"
2488c     enddo
2489
2490!    First get the correct value of datadir, if not already done:
2491        ! default 'datadir' is set in "datafile_mod"
2492        call getin("datadir",datadir)
2493        write(*,*) 'Available surface data files are:'
2494        filestring='ls '//trim(datadir)//' | grep .nc'
2495        call system(filestring)
2496
2497        write(*,*) ''
2498        write(*,*) 'Please choose the relevant file',
2499     &  ' (e.g. "surface_mars.nc")'
2500        write(*,*) ' or "none" to not use any (and set planetary'
2501        write(*,*) '  albedo and surface thermal inertia)'
2502        read(*,fmt='(a50)') surfacefile
2503
2504        if (surfacefile.ne."none") then
2505
2506          CALL datareadnc(relief,surfacefile,phis,alb,surfith,
2507     &          zmeaS,zstdS,zsigS,zgamS,ztheS)
2508        else
2509        ! specific case when not using a "surface.nc" file
2510          phis(:,:)=0
2511          zmeaS(:,:)=0
2512          zstdS(:,:)=0
2513          zsigS(:,:)=0
2514          zgamS(:,:)=0
2515          ztheS(:,:)=0
2516         
2517          write(*,*) "Enter value of albedo of the bare ground:"
2518          read(*,*) alb(1,1)
2519          alb(:,:)=alb(1,1)
2520         
2521          write(*,*) "Enter value of thermal inertia of soil:"
2522          read(*,*) surfith(1,1)
2523          surfith(:,:)=surfith(1,1)
2524       
2525        endif ! of if (surfacefile.ne."none")
2526
2527        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
2528        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
2529        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
2530        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
2531        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
2532        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
2533        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
2534        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
2535
2536      endif ! of if (choix_1.eq.0)
2537
2538
2539c=======================================================================
2540c  Lecture des fichiers (start ou start_archive)
2541c=======================================================================
2542
2543      if (choix_1.eq.0) then
2544
2545        write(*,*) 'Reading file START_ARCHIVE'
2546        CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
2547     &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
2548     &   surfith,nid,
2549     &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
2550        write(*,*) "OK, read start_archive file"
2551        ! copy soil thermal inertia
2552        ithfi(:,:)=inertiedat(:,:)
2553       
2554        ierr= NF_CLOSE(nid)
2555
2556      else if (choix_1.eq.1) then
2557         !do nothing, start and startfi have already been read
2558      else
2559        CALL exit(1)
2560      endif
2561
2562      dtvr   = daysec/FLOAT(day_step)
2563      dtphys   = dtvr * FLOAT(iphysiq)
2564
2565c=======================================================================
2566c
2567c=======================================================================
2568
2569      do ! infinite loop on list of changes
2570
2571      write(*,*)
2572      write(*,*)
2573      write(*,*) 'List of possible changes :'
2574      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
2575      write(*,*)
2576      write(*,*) 'flat : no topography ("aquaplanet")'
2577      write(*,*) 'nuketharsis : no Tharsis bulge'
2578      write(*,*) 'bilball : uniform albedo and thermal inertia'
2579      write(*,*) 'coldspole : cold subsurface and high albedo at S.pole'
2580      write(*,*) 'qname : change tracer name'
2581      write(*,*) 't=profile  : read temperature profile in profile.in'
2582      write(*,*) 'q=0 : ALL tracer =zero'
2583      write(*,*) 'q=x : give a specific uniform value to one tracer'
2584      write(*,*) 'q=profile    : specify a profile for a tracer'
2585!      write(*,*) 'ini_q : tracers initialisation for chemistry, water an
2586!     $d ice   '
2587!      write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and
2588!     $ice '
2589!      write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on
2590!     $ly '
2591      write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45'
2592      write(*,*) 'watercapn : H20 ice on permanent N polar cap '
2593      write(*,*) 'watercaps : H20 ice on permanent S polar cap '
2594      write(*,*) 'noacglac  : H2O ice across Noachis Terra'
2595      write(*,*) 'oborealis : H2O ice across Vastitas Borealis'
2596      write(*,*) 'iceball   : Thick ice layer all over surface'
2597      write(*,*) 'wetstart  : start with a wet atmosphere'
2598      write(*,*) 'isotherm  : Isothermal Temperatures, wind set to zero'
2599      write(*,*) 'radequi   : Earth-like radiative equilibrium temperature
2600     $ profile (lat-alt) and winds set to zero'
2601      write(*,*) 'coldstart : Start X K above the CO2 frost point and
2602     $set wind to zero (assumes 100% CO2)'
2603      write(*,*) 'co2ice=0 : remove CO2 polar cap'
2604      write(*,*) 'ptot : change total pressure'
2605      write(*,*) 'emis : change surface emissivity'
2606      write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur
2607     &face values'
2608      write(*,*) 'slab_ocean_0 : initialisation of slab
2609     $ocean variables'
2610
2611        write(*,*)
2612        write(*,*) 'Change to perform ?'
2613        write(*,*) '   (enter keyword or return to end)'
2614        write(*,*)
2615
2616        read(*,fmt='(a20)') modif
2617        if (modif(1:1) .eq. ' ')then
2618         exit ! exit loop on changes
2619        endif
2620
2621        write(*,*)
2622        write(*,*) trim(modif) , ' : '
2623
2624c       'flat : no topography ("aquaplanet")'
2625c       -------------------------------------
2626        if (trim(modif) .eq. 'flat') then
2627c         set topo to zero 
2628          z_reel(1:iip1,1:jjp1)=0
2629          phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
2630          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
2631          write(*,*) 'topography set to zero.'
2632          write(*,*) 'WARNING : the subgrid topography parameters',
2633     &    ' were not set to zero ! => set calllott to F'                   
2634
2635c        Choice of surface pressure
2636         yes=' '
2637         do while ((yes.ne.'y').and.(yes.ne.'n'))
2638            write(*,*) 'Do you wish to choose homogeneous surface',
2639     &                 'pressure (y) or let newstart interpolate ',
2640     &                 ' the previous field  (n)?'
2641             read(*,fmt='(a)') yes
2642         end do
2643         if (yes.eq.'y') then
2644           flagps0=.true.
2645           write(*,*) 'New value for ps (Pa) ?'
2646 201       read(*,*,iostat=ierr) patm
2647            if(ierr.ne.0) goto 201
2648             write(*,*)
2649             write(*,*) ' new ps everywhere (Pa) = ', patm
2650             write(*,*)
2651             do j=1,jjp1
2652               do i=1,iip1
2653                 ps(i,j)=patm
2654               enddo
2655             enddo
2656         end if
2657
2658c       'nuketharsis : no tharsis bulge for Early Mars'
2659c       ---------------------------------------------
2660        else if (trim(modif) .eq. 'nuketharsis') then
2661
2662           DO j=1,jjp1       
2663              DO i=1,iim
2664                 ig=1+(j-2)*iim +i
2665                 if(j.eq.1) ig=1
2666                 if(j.eq.jjp1) ig=ngridmx
2667
2668                 fact1=(((rlonv(i)*180./pi)+100)**2 + 
2669     &                (rlatu(j)*180./pi)**2)/65**2 
2670                 fact2=exp( -fact1**2.5 )
2671
2672                 phis(i,j) = phis(i,j) - (phis(i,j)+4000.*g)*fact2
2673
2674!                 if(phis(i,j).gt.2500.*g)then
2675!                    if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
2676!                       phis(i,j)=2500.*g
2677!                    endif
2678!                 endif
2679
2680              ENDDO
2681           ENDDO
2682          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
2683
2684
2685c       bilball : uniform albedo, thermal inertia
2686c       -----------------------------------------
2687        else if (trim(modif) .eq. 'bilball') then
2688          write(*,*) 'constante albedo and iner.therm:'
2689          write(*,*) 'New value for albedo (ex: 0.25) ?'
2690 101      read(*,*,iostat=ierr) alb_bb
2691          if(ierr.ne.0) goto 101
2692          write(*,*)
2693          write(*,*) ' uniform albedo (new value):',alb_bb
2694          write(*,*)
2695
2696          write(*,*) 'New value for thermal inertia (eg: 247) ?'
2697 102      read(*,*,iostat=ierr) ith_bb
2698          if(ierr.ne.0) goto 102
2699          write(*,*) 'uniform thermal inertia (new value):',ith_bb
2700          DO j=1,jjp1
2701             DO i=1,iip1
2702                alb(i,j) = alb_bb       ! albedo
2703                do isoil=1,nsoilmx
2704                  ith(i,j,isoil) = ith_bb       ! thermal inertia
2705                enddo
2706             END DO
2707          END DO
2708!          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
2709          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
2710          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
2711
2712c       coldspole : sous-sol de la calotte sud toujours froid
2713c       -----------------------------------------------------
2714        else if (trim(modif) .eq. 'coldspole') then
2715          write(*,*)'new value for the subsurface temperature',
2716     &   ' beneath the permanent southern polar cap ? (eg: 141 K)'
2717 103      read(*,*,iostat=ierr) tsud
2718          if(ierr.ne.0) goto 103
2719          write(*,*)
2720          write(*,*) ' new value of the subsurface temperature:',tsud
2721c         nouvelle temperature sous la calotte permanente
2722          do l=2,nsoilmx
2723               tsoil(ngridmx,l) =  tsud
2724          end do
2725
2726
2727          write(*,*)'new value for the albedo',
2728     &   'of the permanent southern polar cap ? (eg: 0.75)'
2729 104      read(*,*,iostat=ierr) albsud
2730          if(ierr.ne.0) goto 104
2731          write(*,*)
2732
2733c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2734c         Option 1:  only the albedo of the pole is modified :   
2735          albfi(ngridmx)=albsud
2736          write(*,*) 'ig=',ngridmx,'   albedo perennial cap ',
2737     &    albfi(ngridmx)
2738
2739c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
2740c          Option 2 A haute resolution : coordonnee de la vrai calotte ~   
2741c           DO j=1,jjp1
2742c             DO i=1,iip1
2743c                ig=1+(j-2)*iim +i
2744c                if(j.eq.1) ig=1
2745c                if(j.eq.jjp1) ig=ngridmx
2746c                if ((rlatu(j)*180./pi.lt.-84.).and.
2747c     &            (rlatu(j)*180./pi.gt.-91.).and.
2748c     &            (rlonv(i)*180./pi.gt.-91.).and.
2749c     &            (rlonv(i)*180./pi.lt.0.))         then
2750cc    albedo de la calotte permanente fixe a albsud
2751c                   alb(i,j)=albsud
2752c                   write(*,*) 'lat=',rlatu(j)*180./pi,
2753c     &                      ' lon=',rlonv(i)*180./pi
2754cc     fin de la condition sur les limites de la calotte permanente
2755c                end if
2756c             ENDDO
2757c          ENDDO
2758c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2759
2760c         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
2761
2762
2763c       ptot : Modification of the total pressure: ice + current atmosphere 
2764c       -------------------------------------------------------------------
2765        else if (trim(modif).eq.'ptot') then
2766
2767          ! check if we have a co2_ice surface tracer:
2768          if (igcm_co2_ice.eq.0) then
2769            write(*,*) " No surface CO2 ice !"
2770            write(*,*) " only atmospheric pressure will be considered!"
2771          endif
2772c         calcul de la pression totale glace + atm actuelle
2773          patm=0.
2774          airetot=0.
2775          pcap=0.
2776          DO j=1,jjp1
2777             DO i=1,iim
2778                ig=1+(j-2)*iim +i
2779                if(j.eq.1) ig=1
2780                if(j.eq.jjp1) ig=ngridmx
2781                patm = patm + ps(i,j)*aire(i,j)
2782                airetot= airetot + aire(i,j)
2783                if (igcm_co2_ice.ne.0) then
2784                  !pcap = pcap + aire(i,j)*co2ice(ig)*g
2785                  pcap = pcap + aire(i,j)*qsurf(ig,igcm_co2_ice)*g
2786                endif
2787             ENDDO
2788          ENDDO
2789          ptoto = pcap + patm
2790
2791          print*,'Current total pressure at surface (co2 ice + atm) ',
2792     &       ptoto/airetot
2793
2794          print*,'new value?'
2795          read(*,*) ptotn
2796          ptotn=ptotn*airetot
2797          patmn=ptotn-pcap
2798          print*,'ptoto,patm,ptotn,patmn'
2799          print*,ptoto,patm,ptotn,patmn
2800          print*,'Mult. factor for pressure (atm only)', patmn/patm
2801          do j=1,jjp1
2802             do i=1,iip1
2803                ps(i,j)=ps(i,j)*patmn/patm
2804             enddo
2805          enddo
2806
2807
2808
2809c        Correction pour la conservation des traceurs
2810         yes=' '
2811         do while ((yes.ne.'y').and.(yes.ne.'n'))
2812            write(*,*) 'Do you wish to conserve tracer total mass (y)',
2813     &              ' or tracer mixing ratio (n) ?'
2814             read(*,fmt='(a)') yes
2815         end do
2816
2817         if (yes.eq.'y') then
2818           write(*,*) 'OK : conservation of tracer total mass'
2819           DO iq =1, nqtot
2820             DO l=1,llm
2821               DO j=1,jjp1
2822                  DO i=1,iip1
2823                    q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
2824                  ENDDO
2825               ENDDO
2826             ENDDO
2827           ENDDO
2828          else
2829            write(*,*) 'OK : conservation of tracer mixing ratio'
2830          end if
2831
2832c        Correction pour la pression au niveau de la mer
2833         yes=' '
2834         do while ((yes.ne.'y').and.(yes.ne.'n'))
2835            write(*,*) 'Do you wish fix pressure at sea level (y)',
2836     &              ' or not (n) ?'
2837             read(*,fmt='(a)') yes
2838         end do
2839
2840         if (yes.eq.'y') then
2841           write(*,*) 'Value?'
2842                read(*,*,iostat=ierr) psea
2843             DO i=1,iip1
2844               DO j=1,jjp1
2845                    ps(i,j)=psea
2846
2847                  ENDDO
2848               ENDDO
2849                write(*,*) 'psea=',psea
2850          else
2851            write(*,*) 'no'
2852          end if
2853
2854
2855c       emis : change surface emissivity (added by RW)
2856c       ----------------------------------------------
2857        else if (trim(modif).eq.'emis') then
2858
2859           print*,'new value?'
2860           read(*,*) emisread
2861
2862           do i=1,ngridmx
2863              emis(i)=emisread
2864           enddo
2865
2866c       qname : change tracer name
2867c       --------------------------
2868        else if (trim(modif).eq.'qname') then
2869         yes='y'
2870         do while (yes.eq.'y')
2871          write(*,*) 'Which tracer name do you want to change ?'
2872          do iq=1,nqtot
2873            write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
2874          enddo
2875          write(*,'(a35,i3)')
2876     &            '(enter tracer number; between 1 and ',nqtot
2877          write(*,*)' or any other value to quit this option)'
2878          read(*,*) iq
2879          if ((iq.ge.1).and.(iq.le.nqtot)) then
2880            write(*,*)'Change tracer name ',trim(tname(iq)),' to ?'
2881            read(*,*) txt
2882            tname(iq)=txt
2883            write(*,*)'Do you want to change another tracer name (y/n)?'
2884            read(*,'(a)') yes 
2885          else
2886! inapropiate value of iq; quit this option
2887            yes='n'
2888          endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
2889         enddo ! of do while (yes.ne.'y')
2890
2891c       q=0 : set tracers to zero
2892c       -------------------------
2893        else if (trim(modif).eq.'q=0') then
2894c          mise a 0 des q (traceurs)
2895          write(*,*) 'Tracers set to 0 (1.E-30 in fact)'
2896           DO iq =1, nqtot
2897             DO l=1,llm
2898               DO j=1,jjp1
2899                  DO i=1,iip1
2900                    q(i,j,l,iq)=1.e-30
2901                  ENDDO
2902               ENDDO
2903             ENDDO
2904           ENDDO
2905
2906c          set surface tracers to zero
2907           DO iq =1, nqtot
2908             DO ig=1,ngridmx
2909                 qsurf(ig,iq)=0.
2910             ENDDO
2911           ENDDO
2912
2913c       q=x : initialise tracer manually 
2914c       --------------------------------
2915        else if (trim(modif).eq.'q=x') then
2916             write(*,*) 'Which tracer do you want to modify ?'
2917             do iq=1,nqtot
2918               write(*,*)iq,' : ',trim(tname(iq))
2919             enddo
2920             write(*,*) '(choose between 1 and ',nqtot,')'
2921             read(*,*) iq 
2922             write(*,*)'mixing ratio of tracer ',trim(tname(iq)),
2923     &                 ' ? (kg/kg)'
2924             read(*,*) val
2925             DO l=1,llm
2926               DO j=1,jjp1
2927                  DO i=1,iip1
2928                    q(i,j,l,iq)=val
2929                  ENDDO
2930               ENDDO
2931             ENDDO
2932             write(*,*) 'SURFACE value of tracer ',trim(tname(iq)),
2933     &                   ' ? (kg/m2)'
2934             read(*,*) val
2935             DO ig=1,ngridmx
2936                 qsurf(ig,iq)=val
2937             ENDDO
2938
2939c       t=profile : initialize temperature with a given profile
2940c       -------------------------------------------------------
2941        else if (trim(modif) .eq. 't=profile') then
2942             write(*,*) 'Temperature profile from ASCII file'
2943             write(*,*) "'profile.in' e.g. 1D output"
2944             write(*,*) "(one value per line in file; starting with"
2945             write(*,*) "surface value, the 1st atmospheric layer"
2946             write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
2947             txt="profile.in"
2948             open(33,file=trim(txt),status='old',form='formatted',
2949     &            iostat=ierr)
2950             if (ierr.eq.0) then
2951               ! OK, found file 'profile_...', load the profile
2952               do l=1,llm+1
2953                 read(33,*,iostat=ierr) profile(l)
2954                 write(*,*) profile(l)
2955                 if (ierr.ne.0) then ! something went wrong
2956                   exit ! quit loop
2957                 endif
2958               enddo
2959               if (ierr.eq.0) then
2960                 tsurf(1:ngridmx)=profile(1)
2961                 tsoil(1:ngridmx,1:nsoilmx)=profile(1)
2962                 do l=1,llm
2963                   Tset(1:iip1,1:jjp1,l)=profile(l+1)
2964                   flagtset=.true.
2965                 enddo
2966                 ucov(1:iip1,1:jjp1,1:llm)=0.
2967                 vcov(1:iip1,1:jjm,1:llm)=0.
2968                 q2(1:ngridmx,1:nlayermx+1)=0.
2969               else
2970                 write(*,*)'problem reading file ',trim(txt),' !'
2971                 write(*,*)'No modifications to temperature'
2972               endif
2973             else
2974               write(*,*)'Could not find file ',trim(txt),' !'
2975               write(*,*)'No modifications to temperature'
2976             endif
2977
2978c       q=profile : initialize tracer with a given profile
2979c       --------------------------------------------------
2980        else if (trim(modif) .eq. 'q=profile') then
2981             write(*,*) 'Tracer profile will be sought in ASCII file'
2982             write(*,*) "'profile_tracer' where 'tracer' is tracer name"
2983             write(*,*) "(one value per line in file; starting with"
2984             write(*,*) "surface value, the 1st atmospheric layer"
2985             write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
2986             write(*,*) 'Which tracer do you want to set?'
2987             do iq=1,nqtot
2988               write(*,*)iq,' : ',trim(tname(iq))
2989             enddo
2990             write(*,*) '(choose between 1 and ',nqtot,')'
2991             read(*,*) iq 
2992             if ((iq.lt.1).or.(iq.gt.nqtot)) then
2993               ! wrong value for iq, go back to menu
2994               write(*,*) "wrong input value:",iq
2995               cycle
2996             endif
2997             ! look for input file 'profile_tracer'
2998             txt="profile_"//trim(tname(iq))
2999             open(41,file=trim(txt),status='old',form='formatted',
3000     &            iostat=ierr)
3001             if (ierr.eq.0) then
3002               ! OK, found file 'profile_...', load the profile
3003               do l=1,llm+1
3004                 read(41,*,iostat=ierr) profile(l)
3005                 if (ierr.ne.0) then ! something went wrong
3006                   exit ! quit loop
3007                 endif
3008               enddo
3009               if (ierr.eq.0) then
3010                 ! initialize tracer values
3011                 qsurf(:,iq)=profile(1)
3012                 do l=1,llm
3013                   q(:,:,l,iq)=profile(l+1)
3014                 enddo
3015                 write(*,*)'OK, tracer ',trim(tname(iq)),' initialized '
3016     &                    ,'using values from file ',trim(txt)
3017               else
3018                 write(*,*)'problem reading file ',trim(txt),' !'
3019                 write(*,*)'No modifications to tracer ',trim(tname(iq))
3020               endif
3021             else
3022               write(*,*)'Could not find file ',trim(txt),' !'
3023               write(*,*)'No modifications to tracer ',trim(tname(iq))
3024             endif
3025             
3026
3027c      wetstart : wet atmosphere with a north to south gradient
3028c      --------------------------------------------------------
3029       else if (trim(modif) .eq. 'wetstart') then
3030        ! check that there is indeed a water vapor tracer
3031        if (igcm_h2o_vap.eq.0) then
3032          write(*,*) "No water vapour tracer! Can't use this option"
3033          stop
3034        endif
3035          DO l=1,llm
3036            DO j=1,jjp1
3037              DO i=1,iip1
3038                q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
3039              ENDDO
3040            ENDDO
3041          ENDDO
3042
3043         write(*,*) 'Water mass mixing ratio at north pole='
3044     *               ,q(1,1,1,igcm_h2o_vap)
3045         write(*,*) '---------------------------south pole='
3046     *               ,q(1,jjp1,1,igcm_h2o_vap)
3047
3048c      noglacier : remove tropical water ice (to initialize high res sim)
3049c      --------------------------------------------------
3050        else if (trim(modif) .eq. 'noglacier') then
3051           if (igcm_h2o_ice.eq.0) then
3052             write(*,*) "No water ice tracer! Can't use this option"
3053             stop
3054           endif
3055           do ig=1,ngridmx
3056             j=(ig-2)/iim +2
3057              if(ig.eq.1) j=1
3058              write(*,*) 'OK: remove surface ice for |lat|<45'
3059              if (abs(rlatu(j)*180./pi).lt.45.) then
3060                   qsurf(ig,igcm_h2o_ice)=0.
3061              end if
3062           end do
3063
3064
3065c      watercapn : H20 ice on permanent northern cap
3066c      --------------------------------------------------
3067        else if (trim(modif) .eq. 'watercapn') then
3068           if (igcm_h2o_ice.eq.0) then
3069             write(*,*) "No water ice tracer! Can't use this option"
3070             stop
3071           endif
3072
3073           DO j=1,jjp1       
3074              DO i=1,iim
3075                 ig=1+(j-2)*iim +i
3076                 if(j.eq.1) ig=1
3077                 if(j.eq.jjp1) ig=ngridmx
3078
3079                 if (rlatu(j)*180./pi.gt.80.) then
3080                    qsurf(ig,igcm_h2o_ice)=3.4e3
3081                    !do isoil=1,nsoilmx
3082                    !   ith(i,j,isoil) = 18000. ! thermal inertia
3083                    !enddo
3084                   write(*,*)'     ==> Ice mesh North boundary (deg)= ',
3085     &                   rlatv(j-1)*180./pi
3086                 end if
3087              ENDDO
3088           ENDDO
3089           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
3090
3091c$$$           do ig=1,ngridmx
3092c$$$             j=(ig-2)/iim +2
3093c$$$              if(ig.eq.1) j=1
3094c$$$              if (rlatu(j)*180./pi.gt.80.) then
3095c$$$
3096c$$$                   qsurf(ig,igcm_h2o_ice)=1.e5
3097c$$$                   qsurf(ig,igcm_h2o_vap)=0.0!1.e5
3098c$$$
3099c$$$                   write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
3100c$$$     &              qsurf(ig,igcm_h2o_ice)
3101c$$$
3102c$$$                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
3103c$$$     &              rlatv(j)*180./pi
3104c$$$                end if
3105c$$$           enddo
3106
3107c      watercaps : H20 ice on permanent southern cap
3108c      -------------------------------------------------
3109        else if (trim(modif) .eq. 'watercaps') then
3110           if (igcm_h2o_ice.eq.0) then
3111              write(*,*) "No water ice tracer! Can't use this option"
3112              stop
3113           endif
3114
3115           DO j=1,jjp1       
3116              DO i=1,iim
3117                 ig=1+(j-2)*iim +i
3118                 if(j.eq.1) ig=1
3119                 if(j.eq.jjp1) ig=ngridmx
3120
3121                 if (rlatu(j)*180./pi.lt.-80.) then
3122                    qsurf(ig,igcm_h2o_ice)=3.4e3
3123                    !do isoil=1,nsoilmx
3124                    !   ith(i,j,isoil) = 18000. ! thermal inertia
3125                    !enddo
3126                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
3127     &                   rlatv(j-1)*180./pi
3128                 end if
3129              ENDDO
3130           ENDDO
3131           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
3132
3133c$$$           do ig=1,ngridmx
3134c$$$               j=(ig-2)/iim +2
3135c$$$               if(ig.eq.1) j=1
3136c$$$               if (rlatu(j)*180./pi.lt.-80.) then
3137c$$$                  qsurf(ig,igcm_h2o_ice)=1.e5
3138c$$$                  qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
3139c$$$
3140c$$$                  write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
3141c$$$     &                 qsurf(ig,igcm_h2o_ice)
3142c$$$                  write(*,*)'     ==> Ice mesh North boundary (deg)= ',
3143c$$$     &                 rlatv(j-1)*180./pi
3144c$$$               end if
3145c$$$           enddo
3146
3147
3148c       noacglac : H2O ice across highest terrain
3149c       --------------------------------------------
3150        else if (trim(modif) .eq. 'noacglac') then
3151           if (igcm_h2o_ice.eq.0) then
3152             write(*,*) "No water ice tracer! Can't use this option"
3153             stop
3154           endif
3155          DO j=1,jjp1       
3156             DO i=1,iim
3157                ig=1+(j-2)*iim +i
3158                if(j.eq.1) ig=1
3159                if(j.eq.jjp1) ig=ngridmx
3160
3161                if(phis(i,j).gt.1000.*g)then
3162                    alb(i,j) = 0.5 ! snow value
3163                    do isoil=1,nsoilmx
3164                       ith(i,j,isoil) = 18000. ! thermal inertia
3165                       ! this leads to rnat set to 'ocean' in physiq.F90
3166                       ! actually no, because it is soil not surface
3167                    enddo
3168                endif
3169             ENDDO
3170          ENDDO
3171          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
3172          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
3173          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
3174
3175
3176
3177c       oborealis : H2O oceans across Vastitas Borealis
3178c       -----------------------------------------------
3179        else if (trim(modif) .eq. 'oborealis') then
3180           if (igcm_h2o_ice.eq.0) then
3181             write(*,*) "No water ice tracer! Can't use this option"
3182             stop
3183           endif
3184          DO j=1,jjp1       
3185             DO i=1,iim
3186                ig=1+(j-2)*iim +i
3187                if(j.eq.1) ig=1
3188                if(j.eq.jjp1) ig=ngridmx
3189
3190                if(phis(i,j).lt.-4000.*g)then
3191!                if( (phis(i,j).lt.-4000.*g)
3192!     &               .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
3193!                    phis(i,j)=-8000.0*g ! proper ocean
3194                    phis(i,j)=-4000.0*g ! scanty ocean
3195
3196                    alb(i,j) = 0.07 ! oceanic value
3197                    do isoil=1,nsoilmx
3198                       ith(i,j,isoil) = 18000. ! thermal inertia
3199                       ! this leads to rnat set to 'ocean' in physiq.F90
3200                       ! actually no, because it is soil not surface
3201                    enddo
3202                endif
3203             ENDDO
3204          ENDDO
3205          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
3206          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
3207          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
3208
3209c       iborealis : H2O ice in Northern plains
3210c       --------------------------------------
3211        else if (trim(modif) .eq. 'iborealis') then
3212           if (igcm_h2o_ice.eq.0) then
3213             write(*,*) "No water ice tracer! Can't use this option"
3214             stop
3215           endif
3216          DO j=1,jjp1       
3217             DO i=1,iim
3218                ig=1+(j-2)*iim +i
3219                if(j.eq.1) ig=1
3220                if(j.eq.jjp1) ig=ngridmx
3221
3222                if(phis(i,j).lt.-4000.*g)then
3223                   !qsurf(ig,igcm_h2o_ice)=1.e3
3224                   qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
3225                endif
3226             ENDDO
3227          ENDDO
3228          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
3229          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
3230          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
3231
3232
3233c       oceanball : H2O liquid everywhere
3234c       ----------------------------
3235        else if (trim(modif) .eq. 'oceanball') then
3236           if (igcm_h2o_ice.eq.0) then
3237             write(*,*) "No water ice tracer! Can't use this option"
3238             stop
3239           endif
3240          DO j=1,jjp1       
3241             DO i=1,iim
3242                ig=1+(j-2)*iim +i
3243                if(j.eq.1) ig=1
3244                if(j.eq.jjp1) ig=ngridmx
3245
3246                qsurf(ig,igcm_h2o_vap)=0.0    ! for ocean, this is infinite source
3247                qsurf(ig,igcm_h2o_ice)=0.0
3248                alb(i,j) = 0.07 ! ocean value
3249
3250                do isoil=1,nsoilmx
3251                   ith(i,j,isoil) = 18000. ! thermal inertia
3252                   !ith(i,j,isoil) = 50. ! extremely low for test
3253                   ! this leads to rnat set to 'ocean' in physiq.F90
3254                enddo
3255
3256             ENDDO
3257          ENDDO
3258          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
3259          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
3260          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
3261
3262c       iceball : H2O ice everywhere
3263c       ----------------------------
3264        else if (trim(modif) .eq. 'iceball') then
3265           if (igcm_h2o_ice.eq.0) then
3266             write(*,*) "No water ice tracer! Can't use this option"
3267             stop
3268           endif
3269          DO j=1,jjp1       
3270             DO i=1,iim
3271                ig=1+(j-2)*iim +i
3272                if(j.eq.1) ig=1
3273                if(j.eq.jjp1) ig=ngridmx
3274
3275                qsurf(ig,igcm_h2o_vap)=-50.    ! for ocean, this is infinite source
3276                qsurf(ig,igcm_h2o_ice)=50.     ! == to 0.05 m of oceanic ice
3277                alb(i,j) = 0.6 ! ice albedo value
3278
3279                do isoil=1,nsoilmx
3280                   !ith(i,j,isoil) = 18000. ! thermal inertia
3281                   ! this leads to rnat set to 'ocean' in physiq.F90
3282                enddo
3283
3284             ENDDO
3285          ENDDO
3286          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
3287          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
3288
3289c       isotherm : Isothermal temperatures and no winds
3290c       -----------------------------------------------
3291        else if (trim(modif) .eq. 'isotherm') then
3292
3293          write(*,*)'Isothermal temperature of the atmosphere,
3294     &           surface and subsurface'
3295          write(*,*) 'Value of this temperature ? :'
3296 203      read(*,*,iostat=ierr) Tiso
3297          if(ierr.ne.0) goto 203
3298
3299          tsurf(1:ngridmx)=Tiso
3300         
3301          tsoil(1:ngridmx,1:nsoilmx)=Tiso
3302         
3303          Tset(1:iip1,1:jjp1,1:llm)=Tiso
3304          flagtset=.true.
3305         
3306          ucov(1:iip1,1:jjp1,1:llm)=0
3307          vcov(1:iip1,1:jjm,1:llm)=0
3308          q2(1:ngridmx,1:nlayermx+1)=0
3309
3310c       radequi : Radiative equilibrium profile of temperatures and no winds
3311c       --------------------------------------------------------------------
3312        else if (trim(modif) .eq. 'radequi') then
3313
3314          write(*,*)'radiative equilibrium temperature profile'       
3315
3316          do ig=1, ngridmx
3317             teque= 335.0-60.0*sin(latfi(ig))*sin(latfi(ig))-
3318     &            10.0*cos(latfi(ig))*cos(latfi(ig))
3319             tsurf(ig) = MAX(220.0,teque)
3320          end do
3321          do l=2,nsoilmx
3322             do ig=1, ngridmx
3323                tsoil(ig,l) = tsurf(ig)
3324             end do
3325          end do
3326          DO j=1,jjp1
3327             DO i=1,iim
3328                Do l=1,llm
3329                   teque=335.-60.0*sin(rlatu(j))*sin(rlatu(j))-
3330     &                  10.0*cos(rlatu(j))*cos(rlatu(j))
3331                   Tset(i,j,l)=MAX(220.0,teque)
3332                end do
3333             end do
3334          end do
3335          flagtset=.true.
3336          ucov(1:iip1,1:jjp1,1:llm)=0
3337          vcov(1:iip1,1:jjm,1:llm)=0
3338          q2(1:ngridmx,1:nlayermx+1)=0
3339
3340c       coldstart : T set 1K above CO2 frost point and no winds
3341c       ------------------------------------------------
3342        else if (trim(modif) .eq. 'coldstart') then
3343
3344          write(*,*)'set temperature of the atmosphere,' 
3345     &,'surface and subsurface how many degrees above CO2 frost point?'
3346 204      read(*,*,iostat=ierr) Tabove
3347          if(ierr.ne.0) goto 204
3348
3349            DO j=1,jjp1
3350             DO i=1,iim
3351                ig=1+(j-2)*iim +i
3352                if(j.eq.1) ig=1
3353                if(j.eq.jjp1) ig=ngridmx
3354            tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
3355             END DO
3356            END DO
3357          do l=1,nsoilmx
3358            do ig=1, ngridmx
3359              tsoil(ig,l) = tsurf(ig)
3360            end do
3361          end do
3362          DO j=1,jjp1
3363           DO i=1,iim
3364            Do l=1,llm
3365               pp = aps(l) +bps(l)*ps(i,j) 
3366               Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
3367            end do
3368           end do
3369          end do
3370
3371          flagtset=.true.
3372          ucov(1:iip1,1:jjp1,1:llm)=0
3373          vcov(1:iip1,1:jjm,1:llm)=0
3374          q2(1:ngridmx,1:nlayermx+1)=0
3375
3376
3377c       co2ice=0 : remove CO2 polar ice caps'
3378c       ------------------------------------------------
3379        else if (trim(modif) .eq. 'co2ice=0') then
3380         ! check that there is indeed a co2_ice tracer ...
3381          if (igcm_co2_ice.ne.0) then
3382           do ig=1,ngridmx
3383              !co2ice(ig)=0
3384              qsurf(ig,igcm_co2_ice)=0
3385              emis(ig)=emis(ngridmx/2)
3386           end do
3387          else
3388            write(*,*) "Can't remove CO2 ice!! (no co2_ice tracer)"
3389          endif
3390       
3391!       therm_ini_s: (re)-set soil thermal inertia to reference surface values
3392!       ----------------------------------------------------------------------
3393
3394        else if (trim(modif) .eq. 'therm_ini_s') then
3395!          write(*,*)"surfithfi(1):",surfithfi(1)
3396          do isoil=1,nsoilmx
3397            inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
3398          enddo
3399          write(*,*)'OK: Soil thermal inertia has been reset to referenc
3400     &e surface values'
3401!         write(*,*)"inertiedat(1,1):",inertiedat(1,1)
3402          ithfi(:,:)=inertiedat(:,:)
3403         ! recast ithfi() onto ith()
3404         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
3405! Check:
3406!         do i=1,iip1
3407!           do j=1,jjp1
3408!             do isoil=1,nsoilmx
3409!               write(77,*) i,j,isoil,"  ",ith(i,j,isoil)
3410!             enddo
3411!           enddo
3412!        enddo
3413
3414
3415
3416c       slab_ocean_initialisation
3417c       ------------------------------------------------
3418        else if (trim(modif) .eq. 'slab_ocean_0') then
3419        write(*,*)'OK: initialisation of slab ocean' 
3420
3421      DO ig=1, ngridmx
3422         rnat(ig)=1.
3423         tslab(ig,1)=0.
3424         tslab(ig,2)=0.
3425         tsea_ice(ig)=0.
3426         sea_ice(ig)=0.
3427         pctsrf_sic(ig)=0.
3428         
3429         if(ithfi(ig,1).GT.10000.)then
3430           rnat(ig)=0.
3431           phisfi(ig)=0.
3432           tsea_ice(ig)=273.15-1.8
3433           tslab(ig,1)=tsurf(ig)
3434           tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
3435          endif
3436
3437      ENDDO
3438          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
3439
3440
3441
3442        else
3443          write(*,*) '       Unknown (misspelled?) option!!!'
3444        end if ! of if (trim(modif) .eq. '...') elseif ...
3445
3446
3447
3448       enddo ! of do ! infinite loop on liste of changes
3449
3450 999  continue
3451
3452 
3453c=======================================================================
3454c   Initialisation for cloud fraction and oceanic ice (added by BC 2010)
3455c=======================================================================
3456      DO ig=1, ngridmx
3457         totalfrac(ig)=0.5
3458         DO l=1,llm
3459            cloudfrac(ig,l)=0.5
3460         ENDDO
3461! Ehouarn, march 2012: also add some initialisation for hice
3462         hice(ig)=0.0
3463      ENDDO
3464     
3465c========================================================
3466
3467!      DO ig=1,ngridmx
3468!         IF(tsurf(ig) .LT. 273.16-1.8) THEN
3469!            hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
3470!            hice(ig)=min(hice(ig),1.0)
3471!         ENDIF
3472!      ENDDO
3473
3474
3475
3476
3477c=======================================================================
3478c   Correct pressure on the new grid (menu 0)
3479c=======================================================================
3480
3481
3482      if ((choix_1.eq.0).and.(.not.flagps0)) then
3483        r = 1000.*8.31/mugaz
3484
3485        do j=1,jjp1
3486          do i=1,iip1
3487             ps(i,j) = ps(i,j) *
3488     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
3489     .                                  (t(i,j,1) * r))
3490          end do
3491        end do
3492
3493c periodicite de ps en longitude
3494        do j=1,jjp1
3495          ps(1,j) = ps(iip1,j)
3496        end do
3497      end if
3498
3499         
3500c=======================================================================
3501c=======================================================================
3502
3503c=======================================================================
3504c    Initialisation de la physique / ecriture de newstartfi :
3505c=======================================================================
3506
3507
3508      CALL inifilr 
3509      CALL pression(ip1jmp1, ap, bp, ps, p3d)
3510
3511c-----------------------------------------------------------------------
3512c   Initialisation  pks:
3513c-----------------------------------------------------------------------
3514
3515      CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
3516! Calcul de la temperature potentielle teta
3517
3518      if (flagtset) then
3519          DO l=1,llm
3520             DO j=1,jjp1
3521                DO i=1,iim
3522                   teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
3523                ENDDO
3524                teta (iip1,j,l)= teta (1,j,l)
3525             ENDDO
3526          ENDDO
3527      else if (choix_1.eq.0) then
3528         DO l=1,llm
3529            DO j=1,jjp1
3530               DO i=1,iim
3531                  teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
3532               ENDDO
3533               teta (iip1,j,l)= teta (1,j,l)
3534            ENDDO
3535         ENDDO
3536      endif
3537
3538C Calcul intermediaire
3539c
3540      if (choix_1.eq.0) then
3541         CALL massdair( p3d, masse  )
3542c
3543         print *,' ALPHAX ',alphax
3544
3545         DO  l = 1, llm
3546           DO  i    = 1, iim
3547             xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
3548             xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
3549           ENDDO
3550             xpn      = SUM(xppn)/apoln
3551             xps      = SUM(xpps)/apols
3552           DO i   = 1, iip1
3553             masse(   i   ,   1     ,  l )   = xpn
3554             masse(   i   ,   jjp1  ,  l )   = xps
3555           ENDDO
3556         ENDDO
3557      endif
3558      phis(iip1,:) = phis(1,:)
3559
3560      itau=0
3561      if (choix_1.eq.0) then
3562         day_ini=int(date)
3563      endif
3564c
3565      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
3566
3567      CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
3568     *                phi,w, pbaru,pbarv,day_ini+time )
3569
3570         
3571      CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqtot)
3572      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqtot,masse,ps) 
3573C
3574C Ecriture etat initial physique
3575C
3576
3577      call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
3578     &              nqtot,dtphys,real(day_ini),0.0,
3579     &              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
3580      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nq,
3581     &                dtphys,real(day_ini),
3582     &                tsurf,tsoil,emis,q2,qsurf,
3583     &                cloudfrac,totalfrac,hice,
3584     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
3585
3586c=======================================================================
3587c       Formats 
3588c=======================================================================
3589
3590   1  FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
3591     *rrage est differente de la valeur parametree iim =',i4//)
3592   2  FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
3593     *rrage est differente de la valeur parametree jjm =',i4//)
3594   3  FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
3595     *rage est differente de la valeur parametree llm =',i4//)
3596
3597      write(*,*) "newstart: All is well that ends well."
3598
3599      end
3600
Note: See TracBrowser for help on using the repository browser.