source: trunk/MESOSCALE/LMDZ.MARS/libf_gcm/dyn3d/newstart.F @ 3586

Last change on this file since 3586 was 1858, checked in by aslmd, 7 years ago

added gfortran support for old GCM. commented part of ch.f which pose problems (not used anyway in low-atmosphere GCM for mesoscale applications. added an adapted makegcm_gnu. corrected a problem of condition in newstart picked by picky gfortran.
File size: 30.9 KB
RevLine 
[57]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      implicit none
18
19#include "dimensions.h"
20#include "dimphys.h"
21#include "surfdat.h"
22#include "dimradmars.h"
23#include "yomaer.h"
24#include "planete.h"
25#include "paramet.h"
26#include "comconst.h"
27#include "comvert.h"
28#include "comgeom2.h"
29#include "control.h"
30#include "logic.h"
31#include "description.h"
32#include "ener.h"
33#include "temps.h"
34#include "lmdstd.h"
35#include "comdissnew.h"
36#include "clesph0.h"
37#include "serre.h"
38#include "netcdf.inc"
39
40c=======================================================================
41c   Declarations
42c=======================================================================
43
44c Variables dimension du fichier "start_archive"
45c------------------------------------
46      CHARACTER relief*3
47
48c et autres:
49c----------
50      INTEGER lnblnk
51      EXTERNAL lnblnk
52
53c Variables pour les lectures NetCDF des fichiers "start_archive"
54c--------------------------------------------------
55      INTEGER nid_dyn, nid_fi,nid,nvarid
56      INTEGER size
57      INTEGER length
58      parameter (length = 100)
59      INTEGER tab0
60      INTEGER   NB_ETATMAX
61      parameter (NB_ETATMAX = 100)
62
63      REAL  date
64      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
65
66c Variable histoire
67c------------------
68      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
69      REAL phis(iip1,jjp1)
70      REAL q(iip1,jjp1,llm,nqmx)               ! champs advectes
71
72c autre variables dynamique nouvelle grille
73c------------------------------------------
74      REAL pks(iip1,jjp1)
75      REAL w(iip1,jjp1,llm+1)
76      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
77      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
78      REAL dh(ip1jmp1,llm),dp(ip1jmp1)
79      REAL phi(iip1,jjp1,llm)
80
81      integer klatdat,klongdat
82      PARAMETER (klatdat=180,klongdat=360)
83
84c Physique sur grille scalaire
85c----------------------------
86      real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
87      real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
88
89c variable physique
90c------------------
91      REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx),co2ice(ngridmx)
92      REAL emis(ngridmx),qsurf(ngridmx,nqmx)
93      REAL q2(ngridmx,nlayermx+1)
94      REAL rnaturfi(ngridmx)
95      real alb(iip1,jjp1),albfi(ngridmx)
96      real ith(iip1,jjp1),ithfi(ngridmx)
97      REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
98
99      INTEGER i,j,l,idum,ig
100      real mugaz ! masse molaire de l''atmosphere
101
102      EXTERNAL iniconst,geopot,inigeom
103      integer ierr, nbetat
104      integer ismin
105      external ismin
106
107c Variable nouvelle grille naturelle au point scalaire
108c------------------------------------------------------
109      REAL p(iip1,jjp1)
110      REAL t(iip1,jjp1,llm)
111      real phisold_newgrid(iip1,jjp1)
112      REAL :: teta(iip1, jjp1, llm)
113      REAL :: pk(iip1,jjp1,llm)
114      REAL :: pkf(iip1,jjp1,llm)
115      REAL :: ps(iip1, jjp1)
116      REAL :: masse(iip1,jjp1,llm)
117      REAL :: xpn,xps,xppn(iim),xpps(iim)
118      REAL :: p3d(iip1, jjp1, llm+1)
119      REAL :: beta(iip1,jjp1,llm)
120      REAL dteta(ip1jmp1,llm)
121
122c Variable de l'ancienne grille
123c------------------------------
124      real time
125      real tab_cntrl(100)
126      real tab_cntrl_bis(100)
127
128c variables diverses
129c-------------------
130      real choix_1
131      character*80      fichnom
132      integer Lmodif,iq,thermo
133      character modif*20
134      real z_reel(iip1,jjp1)
135      real tsud,albsud,alb_bb,ith_bb,Tiso
136      real ptoto,pcap,patm,airetot,ptotn,patmn
137      real ssum
138      character*1 yes
139      logical :: flagiso=.false. ,  flagps0=.false.
140      real val
141
142      INTEGER :: itau
143
144c sortie visu pour les champs dynamiques
145c---------------------------------------
146      INTEGER :: visuid
147      real :: time_step,t_ops,t_wrt
148      CHARACTER*80 :: visu_file
149
150      cpp    = 744.499 ! au lieu de 1004.70885 (Terre)
151      preff  = 610.    ! au lieu de 101325. (Terre)
152      pa= 20           ! au lieu de 500 (Terre)
153
154c=======================================================================
155c   Choix du fichier de demarrage a modifier
156c=======================================================================
157
158      write(*,*) 'From which kind of files do you want to create new',
159     .  'start and startfi files'
160      write(*,*) '    0 - from a file start_archive'
161      write(*,*) '    1 - from files start and startfi'
162 
163c-----------------------------------------------------------------------
164c   Ouverture des fichiers a modifier (start ou start_archive)
165c-----------------------------------------------------------------------
166
167      DO
168         read(*,*,iostat=ierr) choix_1
169         if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
170      ENDDO
171
172c     Ouverture de start_archive
173c     ~~~~~~~~~~~~~~~~~~~~~~~~~~
174      if (choix_1.eq.0) then
175
176        write(*,*) 'Creating start files from:'
177        write(*,*) './start_archive.dat and ./start_archive.dic'
178        write(*,*)
179        fichnom = 'start_archive.nc'
180        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
181        IF (ierr.NE.NF_NOERR) THEN
182          write(6,*)' Problem opening: ',fichnom
183          write(6,*)' ierr = ', ierr
184          CALL ABORT
185        ENDIF
186        tab0 = 50
187        Lmodif = 1
188
189c     OU Ouverture de start
190c     ~~~~~~~~~~~~~~~~~~~~~
191      else
192        write(*,*) 'Creating start files from:'
193        write(*,*) './start.nc and ./startfi.nc'
194        write(*,*)
195        fichnom = 'start.nc'
196        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
197        IF (ierr.NE.NF_NOERR) THEN
198          write(6,*)' Problem opening:  ',fichnom
199          write(6,*)' ierr = ', ierr
200          CALL ABORT
201        ENDIF
202 
203        fichnom = 'startfi.nc'
204        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
205        IF (ierr.NE.NF_NOERR) THEN
206          write(6,*)' Problem opening:  ',fichnom
207          write(6,*)' ierr = ', ierr
208          CALL ABORT
209        ENDIF
210
211        tab0 = 0
212        Lmodif = 0
213
214      endif
215
216c-----------------------------------------------------------------------
217c Lecture du tableau des parametres du run (pour la dynamique)
218c-----------------------------------------------------------------------
219
220      if (choix_1.eq.0) then
221
222        write(*,*) 'reading tab_cntrl START_ARCHIVE'
223c
224        ierr = NF_INQ_VARID (nid, "controle", nvarid)
225#ifdef NC_DOUBLE
226        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
227#else
228        ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
229#endif
230c
231      else if (choix_1.eq.1) then
232
233        write(*,*) 'reading tab_cntrl START'
234c
235        ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
236#ifdef NC_DOUBLE
237        ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
238#else
239        ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
240#endif
241c
242        write(*,*) 'reading tab_cntrl STARTFI'
243c
244        ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
245#ifdef NC_DOUBLE
246        ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
247#else
248        ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
249#endif
250c
251        do i=1,50
252          tab_cntrl(i+50)=tab_cntrl_bis(i)
253        enddo
254      write(*,*) 'printing tab_cntrl', tab_cntrl
255      do i=1,100
256        write(*,*) i,tab_cntrl(i)
257      enddo
258     
259      endif
260c-----------------------------------------------------------------------
261c               Initialisation des constantes dynamique
262c-----------------------------------------------------------------------
263
264      kappa = tab_cntrl(9)
265      etot0 = tab_cntrl(12)
266      ptot0 = tab_cntrl(13)
267      ztot0 = tab_cntrl(14)
268      stot0 = tab_cntrl(15)
269      ang0 = tab_cntrl(16)
270      write(*,*) "kappa,etot0,ptot0,ztot0,stot0,ang0"
271      write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0
272
273c-----------------------------------------------------------------------
274c   Lecture du tab_cntrl et initialisation des constantes physiques
275c  - pour start:  Lmodif = 0 => pas de modifications possibles
276c                  (modif dans le tabfi de readfi + loin)
277c  - pour start_archive:  Lmodif = 1 => modifications possibles
278c-----------------------------------------------------------------------
279      if (choix_1.eq.0) then
280         call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
281     .            p_omeg,p_g,p_mugaz,p_daysec,time)
282      else if (choix_1.eq.1) then
283         call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
284     .            p_omeg,p_g,p_mugaz,p_daysec,time)
285      endif
286
287      rad = p_rad
288      omeg = p_omeg
289      g = p_g
290      mugaz = p_mugaz
291      daysec = p_daysec
292      write(*,*) 'aire',aire
293
294
295c=======================================================================
296c  INITIALISATIONS DIVERSES
297c=======================================================================
298
299      day_step=180
300      CALL defrun_new( 99, .TRUE. )
301      CALL iniconst
302      CALL inigeom
303      idum=-1
304      idum=0
305
306c Initialisation coordonnees /aires
307c -------------------------------
308
309      latfi(1)=rlatu(1)
310      lonfi(1)=0.
311      DO j=2,jjm
312         DO i=1,iim
313            latfi((j-2)*iim+1+i)=rlatu(j)
314            lonfi((j-2)*iim+1+i)=rlonv(i)
315         ENDDO
316      ENDDO
317      latfi(ngridmx)=rlatu(jjp1)
318      lonfi(ngridmx)=0.
319      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
320
321
322c=======================================================================
323c  Lecture des fichiers (start ou start_archive)
324c=======================================================================
325
326      if (choix_1.eq.0) then
327
328        write(*,*) 'Reading file START_ARCHIVE'
329        CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
330     .   t,ucov,vcov,ps,co2ice,teta,phisold_newgrid,q,qsurf,nid)
331       
332        ierr= NF_CLOSE(nid)
333
334      else if (choix_1.eq.1) then !  c'est l'appel a tabfi de phyeta0 qui
335                                  !  permet de changer les valeurs du
336                                  !  tab_cntrl Lmodif=1
337        tab0=0
338        Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
339        write(*,*) 'Reading file START'
340        fichnom = 'start.nc'
341        CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
342     .       ps,phis,time)
343
344        write(*,*) 'Reading file STARTFI'
345        fichnom = 'startfi.nc'
346        CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
347     .        day_ini,time,
348     .        tsurf,tsoil,emis,q2,qsurf,co2ice)
349       
350        do i=1,ngridmx
351          albfi(i) = albedodat(i)
352          ithfi(i) = inertiedat(i)
353        enddo
354c       save old topography
355          do j=1,jjp1
356            do i=1,iip1
357               phisold_newgrid(i,j)=phis(i,j)         
358            enddo
359          enddo
360      else
361        CALL exit(1)
362      endif
363
364      dtvr   = daysec/FLOAT(day_step)
365      dtphys   = dtvr * FLOAT(iphysiq)
366
367c=======================================================================
368c   lecture topographie, albedo, inertie thermique, relief sous-maille
369c=======================================================================
370
371      if (choix_1.ne.1) then  ! pour ne pas avoir besoin du fichier
372                              ! surface.dat dans le cas des start
373
374c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
375c       write(*,*)
376c       write(*,*) 'choix du relief (mola,pla)'
377c       write(*,*) '(Topographie MGS MOLA, plat)'
378c       read(*,fmt='(a3)') relief
379        relief="mola"
380c     enddo
381
382      CALL datareadnc(relief,phis,alb,ith,zmeaS,zstdS,zsigS,zgamS,
383     .          ztheS)
384
385      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
386      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
387      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
388      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
389      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
390      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
391      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
392      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
393
394      endif
395
396c=======================================================================
397c
398c=======================================================================
399 
400 888  continue
401
402      write(*,*)
403      write(*,*)
404      write(*,*) 'Other possible changes :'
405      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
406      write(*,*)
407      write(*,*) 'flat : no topography ("aquaplanet")'
408      write(*,*) 'bilball : uniform albedo and thermal inertia '
409      write(*,*) 'coldspole : cold subsurface and high albedo at S.pole'
410      write(*,*) 'q=0 : ALL tracer =zero'
411      write(*,*) 'q=x : give a specific uniform value to one tracer'
412      write(*,*) 'ini_q : tracers initialisation for chemistry, water an
413     $d ice   '
414      write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and
415     $ice '
416      write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on
417     $ly '
418      write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45'
419      write(*,*) 'watercapn : H20 ice on permanent N polar cap '
420      write(*,*) 'watercaps : H20 ice on permanent S polar cap '
421      write(*,*) 'wetstart  : start with a wet atmosphere'
422      write(*,*) 'isotherm : Isothermal Temperatures, wind set to zero'
423      write(*,*) 'co2ice=0 : remove CO2 polar cap'
424      write(*,*) 'ptot : change total pressure'
425
426
427      do while(modif(1:1).ne.'hello')
428        write(*,*)
429        write(*,*) 'Changes to perform ?'
430        write(*,*) '   (enter keyword or return to end)'
431        write(*,*)
432
433        read(*,fmt='(a20)') modif
434        if (modif(1:1) .eq. ' ') goto 999
435        write(*,*)
436        write(*,*) modif(1:lnblnk(modif)) , ' : '
437
438c       'flat : no topography ("aquaplanet")'
439c       -------------------------------------
440        if (modif(1:lnblnk(modif)) .eq. 'flat') then
441c         set topo to zero
442          CALL initial0(ip1jmp1,z_reel)
443          CALL multscal(ip1jmp1,z_reel,g,phis)
444          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
445          write(*,*) 'topography set to zero.'
446          write(*,*) 'WARNING : the subgrid topography parameters',
447     &    ' were not set to zero ! => set calllott to F'                   
448
449c        Choice for surface pressure
450         yes=' '
451         do while ((yes.ne.'y').and.(yes.ne.'n'))
452            write(*,*) 'Do you wish to choose homogeneous surface',
453     &                 'pressure (y) or let newstart interpolate ',
454     &                 ' the previous field  (n)?'
455             read(*,fmt='(a)') yes
456         end do
457         if (yes.eq.'y') then
458           flagps0=.true.
459           write(*,*) 'New value for ps (Pa) ?'
460 201       read(*,*,iostat=ierr) patm
461            if(ierr.ne.0) goto 201
462             write(*,*)
463             write(*,*) ' new ps everywhere (Pa) = ', patm
464             write(*,*)
465             do j=1,jjp1
466               do i=1,iip1
467                 ps(i,j)=patm
468               enddo
469             enddo
470         end if
471
472c       bilball : albedo, inertie thermique uniforme
473c       --------------------------------------------
474        else if (modif(1:lnblnk(modif)) .eq. 'bilball') then
475          write(*,*) 'constante albedo and iner.therm:'
476          write(*,*) 'New value for albedo (ex: 0.25) ?'
477 101      read(*,*,iostat=ierr) alb_bb
478          if(ierr.ne.0) goto 101
479          write(*,*)
480          write(*,*) ' albedo uniform (new value):',alb_bb
481          write(*,*)
482
483          write(*,*) 'New value for iner.therm (ex: 247) ?'
484 102      read(*,*,iostat=ierr) ith_bb
485          if(ierr.ne.0) goto 102
486          write(*,*) 'iner.therm uniform (new value):',ith_bb
487          DO j=1,jjp1
488             DO i=1,iip1
489                alb(i,j) = alb_bb
490                ith(i,j) = ith_bb
491             END DO
492          END DO
493          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
494          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
495
496c       coldspole : sous-sol de la calotte sud toujours froid
497c       -----------------------------------------------------
498        else if (modif(1:lnblnk(modif)) .eq. 'coldspole') then
499          write(*,*)'nouvelle valeur de la temperature du',
500     &   'sous sol de la calotte permanente sud ? (ex: 141 K)'
501 103      read(*,*,iostat=ierr) tsud
502          if(ierr.ne.0) goto 103
503          write(*,*)
504          write(*,*) ' nouvelle valeur de la temperature:',tsud
505c         nouvelle temperature sous la calotte permanente
506          do l=2,nsoilmx
507               tsoil(ngridmx,l) =  tsud
508          end do
509
510
511          write(*,*)'nouvelle valeur de l albedo',
512     &   'de la calotte permanente sud ? (ex: 0.75 K)'
513 104      read(*,*,iostat=ierr) albsud
514          if(ierr.ne.0) goto 104
515          write(*,*)
516
517c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
518c         Option 1:  Seul l'albedo du pole est modifié :   
519          albfi(ngridmx)=albsud
520          write(*,*) 'ig=',ngridmx,'   albedo perennial cap ',
521     &    albfi(ngridmx)
522
523c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~é
524c          Option 2 A haute resolution : coordonnée de la vrai calotte ~   
525c           DO j=1,jjp1
526c             DO i=1,iip1
527c                ig=1+(j-2)*iim +i
528c                if(j.eq.1) ig=1
529c                if(j.eq.jjp1) ig=ngridmx
530c                if ((rlatu(j)*180./pi.lt.-84.).and.
531c     &            (rlatu(j)*180./pi.gt.-91.).and.
532c     &            (rlonv(i)*180./pi.gt.-91.).and.
533c     &            (rlonv(i)*180./pi.lt.0.))         then
534cc    albedo de la calotte permanente fixe a albsud
535c                   alb(i,j)=albsud
536c                   write(*,*) 'lat=',rlatu(j)*180./pi,
537c     &                      ' lon=',rlonv(i)*180./pi
538cc     fin de la condition sur les limites de la calotte permanente
539c                end if
540c             ENDDO
541c          ENDDO
542c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
543
544c         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
545
546
547c       ptot : Modification de la pression totale glace + atm actuelle
548c       --------------------------------------------------------------
549        else if (modif(1:lnblnk(modif)) .eq. 'ptot') then
550
551c         calcul de la pression totale glace + atm actuelle
552          patm=0.
553          airetot=0.
554          pcap=0.
555          DO j=1,jjp1
556             DO i=1,iim
557                ig=1+(j-2)*iim +i
558                if(j.eq.1) ig=1
559                if(j.eq.jjp1) ig=ngridmx
560                patm = patm + ps(i,j)*aire(i,j)
561                airetot= airetot + aire(i,j)
562                pcap = pcap + aire(i,j)*co2ice(ig)*g
563             ENDDO
564          ENDDO
565          ptoto = pcap + patm
566
567          print*,'Pression totale au sol actuelle (co2 ice + atm) ',
568     &       ptoto/airetot
569
570          print*,'nouvelle valeur?'
571          read(*,*) ptotn
572          ptotn=ptotn*airetot
573          patmn=ptotn-pcap
574          print*,'ptoto,patm,ptotn,patmn'
575          print*,ptoto,patm,ptotn,patmn
576          print*,'Facteur mult de la pression (atm only)', patmn/patm
577          do j=1,jjp1
578             do i=1,iip1
579                ps(i,j)=ps(i,j)*patmn/patm
580             enddo
581          enddo
582
583c        Correction pour la conservation des traceurs
584         yes=' '
585         do while ((yes.ne.'y').and.(yes.ne.'n'))
586            write(*,*) 'Do you wish to conserve tracer total mass (y)',
587     &              ' or tracer mixing ratio (n) ?'
588             read(*,fmt='(a)') yes
589         end do
590
591         if (yes.eq.'y') then
592           write(*,*) 'OK : conservation of tracer total mass'
593           DO iq =1, nqmx
594             DO l=1,llm
595               DO j=1,jjp1
596                  DO i=1,iip1
597                    q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
598                  ENDDO
599               ENDDO
600             ENDDO
601           ENDDO
602          else
603            write(*,*) 'OK : conservation of tracer mixing ratio'
604          end if
605
606c       q=0 : traceurs a zero
607c       ---------------------
608        else if (modif(1:lnblnk(modif)) .eq. 'q=0') then
609c          mise a 0 des q (traceurs)
610          write(*,*) 'Traceurs mis a 0 (1.E-30 en fait)'
611           DO iq =1, nqmx
612             DO l=1,llm
613               DO j=1,jjp1
614                  DO i=1,iip1
615                    q(i,j,l,iq)=1.e-30
616                  ENDDO
617               ENDDO
618             ENDDO
619           ENDDO
620
621c          mise a 0 des qsurf (traceurs a la surface)
622           DO iq =1, nqmx
623             DO ig=1,ngridmx   
624                 qsurf(ig,iq)=0.
625             ENDDO
626           ENDDO
627
628c       q=x : initialise tracer manually
629c       --------------------------------
630        else if (modif(1:lnblnk(modif)) .eq. 'q=x') then
631             write(*,*) 'Which tracer do you want to modify ?'
632             write(*,*) '(choose between 1 and ',nqmx,')'
633             read(*,*) iq
634             write(*,*)'mixing ratio of tracer #',iq, ' ? (kg/kg)'
635             read(*,*) val
636             DO l=1,llm
637               DO j=1,jjp1
638                  DO i=1,iip1
639                    q(i,j,l,iq)=val
640                  ENDDO
641               ENDDO
642             ENDDO
643             write(*,*) 'SURFACE value of tracer #', iq , ' ? (kg/m2)'
644             read(*,*) val
645             DO ig=1,ngridmx
646                 qsurf(ig,iq)=val
647             ENDDO
648
649c       ini_q : traceurs initialises pour la chimie
650c       -----------------------------------------------
651        else if (modif(1:lnblnk(modif)) .eq. 'ini_q') then
652c         For more than 32 layers, possible to initiate thermosphere only     
653          thermo=0
654          yes=' '
655          if (llm.gt.32) then
656            do while ((yes.ne.'y').and.(yes.ne.'n'))
657            write(*,*)'',
658     &     'initialisation for thermosphere only? (y/n)'
659            read(*,fmt='(a)') yes
660            if (yes.eq.'y') then
661            thermo=1
662            else
663            thermo=0
664            endif
665            enddo 
666          endif
667         
668              call inichim_newstart(q,ps,sig,nqmx,latfi,lonfi,airefi,
669     $   thermo)
670          write(*,*) 'Especes chimiques initialisees'
671
672        if (thermo.eq.0) then
673c          mise a 0 des qsurf (traceurs a la surface)
674           DO iq =1, nqmx
675             DO ig=1,ngridmx
676                 qsurf(ig,iq)=0.
677             ENDDO
678           ENDDO
679        endif   
680
681c       ini_q-H2O : idem sauf le dernier traceur (H2O)
682c       -----------------------------------------------
683        else if (modif(1:lnblnk(modif)) .eq. 'ini_q-H2O') then
684          ! for more than 32 layers, possible to initiate thermosphere only     
685          thermo=0
686          yes=' '
687          if(llm.gt.32) then
688            do while ((yes.ne.'y').and.(yes.ne.'n'))
689            write(*,*)'',
690     &      'initialisation for thermosphere only? (y/n)'
691            read(*,fmt='(a)') yes
692            if (yes.eq.'y') then
693            thermo=1
694            else
695            thermo=0
696            endif
697            enddo
698          endif
699              call inichim_newstart(q,ps,sig,nqmx-1,latfi,lonfi,airefi,
700     $   thermo)
701          write(*,*) 'Esp. chim. initialisees sauf derniere (H2O)'
702
703        if (thermo.eq.0) then
704c          mise a 0 des qsurf (traceurs a la surface)
705           DO iq =1, nqmx-1
706             DO ig=1,ngridmx
707                 qsurf(ig,iq)=0.
708             ENDDO
709           ENDDO
710        endif
711
712c       ini_q-iceH2O : idem sauf ice et H2O
713c       -----------------------------------------------
714        else if (modif(1:lnblnk(modif)) .eq. 'ini_q-iceH2O') then
715c         For more than 32 layers, possible to initiate thermosphere only     
716          thermo=0
717          yes=' '
718          if(llm.gt.32) then
719            do while ((yes.ne.'y').and.(yes.ne.'n'))
720            write(*,*)'',
721     &      'initialisation for thermosphere only? (y/n)'
722            read(*,fmt='(a)') yes
723            if (yes.eq.'y') then
724            thermo=1
725            else
726            thermo=0
727            endif
728            enddo
729          endif
730     
731         call inichim_newstart(q,ps,sig,nqmx-2,latfi,lonfi,airefi,
732     $   thermo)
733          write(*,*) 'Esp. chim. initialisees sauf ice et H2O'
734
735        if (thermo.eq.0) then
736c          mise a 0 des qsurf (traceurs a la surface)
737           DO iq =1, nqmx-2
738             DO ig=1,ngridmx
739                 qsurf(ig,iq)=0.
740             ENDDO
741           ENDDO
742        endif
743c      wetstart : wet atmosphere with a north to south gradient
744c      --------------------------------------------------------
745       else if (modif(1:lnblnk(modif)) .eq. 'wetstart') then
746          DO l=1,llm
747            DO j=1,jjp1
748              DO i=1,iip1
749                q(i,j,l,nqmx)=150.e-6 * (rlatu(j)+pi/2.) / pi
750              ENDDO
751            ENDDO
752          ENDDO
753
754         write(*,*) 'Water mass mixing ratio at north pole='
755     *               ,q(1,1,1,nqmx)
756         write(*,*) '---------------------------south pole='
757     *               ,q(1,jjp1,1,nqmx)
758
759c      noglacier : remove tropical water ice (to initialize high res sim)
760c      --------------------------------------------------
761        else if (modif(1:lnblnk(modif)) .eq. 'noglacier') then
762           do ig=1,ngridmx
763             j=(ig-2)/iim +2
764              if(ig.eq.1) j=1
765              write(*,*) 'OK: remove surface ice for |lat|<45'
766              if (abs(rlatu(j)*180./pi).lt.45.) then
767                   qsurf(ig,nqmx)=0.
768              end if
769           end do
770c      watercapn : H20 ice sur la calotte permanente nord
771c      --------------------------------------------------
772        else if (modif(1:lnblnk(modif)) .eq. 'watercapn') then
773           do ig=1,ngridmx
774             j=(ig-2)/iim +2
775              if(ig.eq.1) j=1
776              if (rlatu(j)*180./pi.gt.80.) then
777                   qsurf(ig,nqmx)=1.e5
778                   write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
779     &              qsurf(ig,nqmx)
780                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
781     &              rlatv(j)*180./pi
782                end if
783           enddo
784
785c      watercaps : H20 ice sur la calotte permanente sud
786c      -------------------------------------------------
787        else if (modif(1:lnblnk(modif)) .eq. 'watercaps') then
788           do ig=1,ngridmx
789               j=(ig-2)/iim +2
790               if(ig.eq.1) j=1
791               if (rlatu(j)*180./pi.lt.-80.) then
792                   qsurf(ig,nqmx)=1.e5
793                   write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
794     &              qsurf(ig,nqmx)
795                   write(*,*)'     ==> Ice mesh North boundary (deg)= ',
796     &              rlatv(j-1)*180./pi
797               end if
798           enddo
799
800c       isotherm : Temperatures isothermes et vents nuls
801c       ------------------------------------------------
802        else if (modif(1:lnblnk(modif)) .eq. 'isotherm') then
803
804          write(*,*)'Isothermal temperature of the atmosphere,
805     &           surface and subsurface'
806          write(*,*) 'Value of THIS temperature ? :'
807 203      read(*,*,iostat=ierr) Tiso
808          if(ierr.ne.0) goto 203
809
810          do ig=1, ngridmx
811            tsurf(ig) = Tiso
812          end do
813          do l=2,nsoilmx
814            do ig=1, ngridmx
815              tsoil(ig,l) = Tiso
816            end do
817          end do
818          flagiso=.true.
819          call initial0(llm*ip1jmp1,ucov)
820          call initial0(llm*ip1jm,vcov)
821          call initial0(ngridmx*(llm+1),q2)
822
823c       co2ice=0 : ellimination des calottes polaires de CO2'
824c       ------------------------------------------------
825        else if (modif(1:lnblnk(modif)) .eq. 'co2ice=0') then
826           do ig=1,ngridmx
827              co2ice(ig)=0
828              emis(ig)=emis(ngridmx/2)
829           end do
830        else
831          goto 888
832        end if
833      end do
834
835 999  continue
836
837 
838c=======================================================================
839c   Correction de la pression pour nouvelle grille (menu 0)
840c=======================================================================
841
[1858]842      if (flagps0.eqv..false.) then
[57]843        r = 1000.*8.31/mugaz
844
845        do j=1,jjp1
846          do i=1,iip1
847             ps(i,j) = ps(i,j) *
848     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
849     .                                  (t(i,j,1) * r))
850          end do
851        end do
852 
853c periodicite de ps en longitude
854        do j=1,jjp1
855          ps(1,j) = ps(iip1,j)
856        end do
857      end if
858
859c=======================================================================
860c=======================================================================
861
862c=======================================================================
863c    Initialisation de la physique / ecriture de newstartfi :
864c=======================================================================
865
866
867      CALL inifilr
868      CALL pression(ip1jmp1, ap, bp, ps, p3d)
869
870c-----------------------------------------------------------------------
871c   Initialisation  pks:
872c-----------------------------------------------------------------------
873
874      CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
875! Calcul de la temperature potentielle teta
876
877      if (flagiso) then
878          DO l=1,llm
879             DO j=1,jjp1
880                DO i=1,iim
881                   teta(i,j,l) = Tiso * cpp/pk(i,j,l)
882                ENDDO
883                teta (iip1,j,l)= teta (1,j,l)
884             ENDDO
885          ENDDO
886      else if (choix_1.eq.0) then
887         DO l=1,llm
888            DO j=1,jjp1
889               DO i=1,iim
890                  teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
891               ENDDO
892               teta (iip1,j,l)= teta (1,j,l)
893            ENDDO
894         ENDDO
895      endif
896
897C Calcul intermediaire
898c
899      if (choix_1.eq.0) then
900         CALL massdair( p3d, masse  )
901c
902         print *,' ALPHAX ',alphax
903
904         DO  l = 1, llm
905           DO  i    = 1, iim
906             xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
907             xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
908           ENDDO
909             xpn      = SUM(xppn)/apoln
910             xps      = SUM(xpps)/apols
911           DO i   = 1, iip1
912             masse(   i   ,   1     ,  l )   = xpn
913             masse(   i   ,   jjp1  ,  l )   = xps
914           ENDDO
915         ENDDO
916      endif
917      phis(iip1,:) = phis(1,:)
918
919      CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
920     *                tetagdiv, tetagrot , tetatemp  )
921      itau=0
922      if (choix_1.eq.0) then
923         day_ini=int(date)
924      endif
925c
926      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
927
928      CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
929     *                phi,w, pbaru,pbarv,day_ini+time )
930c     CALL caldyn
931c    $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
932c    $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
933
934      CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
935      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
936C
937C Ecriture etat initial physique
938C
939      call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
940     .              dtphys,float(day_ini),
941     .              time,tsurf,tsoil,co2ice,emis,q2,qsurf,
942     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
943
944c=======================================================================
945c       Formats
946c=======================================================================
947
948   1  FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
949     *rrage est differente de la valeur parametree iim =',i4//)
950   2  FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
951     *rrage est differente de la valeur parametree jjm =',i4//)
952   3  FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
953     *rage est differente de la valeur parametree llm =',i4//)
954
955      end
Note: See TracBrowser for help on using the repository browser.