source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/tabfi.F @ 270

Last change on this file since 270 was 270, checked in by millour, 10 years ago

Added features for the Saturn case:

  • Added possibility to run without startfi or restartfi.nc files
  • Added reference temperature "temp_profile.txt" profile to start from
  • More XIOS outputs, and put them on "presnivs (pressure) vertical coordinate
  • Added "-openmp-threadprivate compat" OpenMP option in Ada arch file

EM

File size: 23.0 KB
Line 
1c=======================================================================
2      SUBROUTINE tabfi(ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad,
3     .                 p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
4c=======================================================================
5c
6c   C. Hourdin 15/11/96
7c
8c   Object:        Lecture du tab_cntrl physique dans un fichier 
9c   ------            et initialisation des constantes physiques
10c
11c   Arguments:
12c   ----------
13c
14c     Inputs:
15c     ------
16c
17c      - nid:    unitne logique du fichier ou on va lire le tab_cntrl   
18c                      (ouvert dans le programme appellant) 
19c
20c                 si nid=0:
21c                       pas de lecture du tab_cntrl mais
22c                       Valeurs par default des constantes physiques
23c       
24c      - tab0:    Offset de tab_cntrl a partir duquel sont ranges 
25c                  les parametres physiques (50 pour start_archive)
26c
27c      - Lmodif:  si on souhaite modifier les constantes  Lmodif = 1 = TRUE
28c
29c
30c     Outputs:
31c     --------
32c
33c      - day_ini: tab_cntrl(tab0+3) (Dans les cas ou l'on souhaite
34c                              comparer avec le day_ini dynamique)
35c
36c      - lmax:    tab_cntrl(tab0+2) (pour test avec nlayer)
37c
38c      - p_rad
39c      - p_omeg   !
40c      - p_g      ! Constantes physiques ayant des
41c      - p_mugaz  ! homonymes dynamiques
42c      - p_daysec !
43c
44c=======================================================================
45! to use  'getin'
46      use ioipsl_getincom_p , only: getin_p
47
48      use surfdat_h, only: albedice, emisice, iceradius, dtemisice,
49     &                     emissiv
50      use comsoil_h, only: volcapa
51      use iostart, only: get_var
52      use mod_phys_lmdz_para, only: is_parallel
53      use planete_mod, only: year_day, periastr, apoastr, peri_day,
54     &                       obliquit, z0, lmixmin, emin_turb
55
56      implicit none
57 
58#include "comcstfi.h"
59!#include "planete.h"
60#include "netcdf.inc"
61#include "callkeys.h"
62
63c-----------------------------------------------------------------------
64c   Declarations
65c-----------------------------------------------------------------------
66
67c Arguments
68c ---------
69      INTEGER,INTENT(IN) :: ngrid,nid,tab0
70      INTEGER*4,INTENT(OUT) :: day_ini
71      INTEGER,INTENT(IN) :: Lmodif
72      INTEGER,INTENT(OUT) :: lmax
73      REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
74
75c Variables
76c ---------
77      INTEGER,PARAMETER :: length=100
78      REAL tab_cntrl(length) ! array in which are stored the run's parameters
79      INTEGER  ierr,nvarid
80      INTEGER size
81      CHARACTER modif*20
82      LOGICAL :: found
83     
84      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
85
86      IF (nid.eq.0) then
87c-----------------------------------------------------------------------
88c  Initialization of various physical constants to defaut values (nid = 0 case)
89c-----------------------------------------------------------------------
90        ! Ehouarn: Default Saturn values:
91        tab_cntrl(:)=0
92        write(*,*) "Using default Saturn values..."
93        ! these should be read in a def file I guess...
94        lmax=0 ! not used anyways
95        !day_ini=0
96        time=0
97        ! radius of the planet
98        rad=60268000
99        call getin_p("radius",rad)
100        ! Planetary rotation rate
101        omeg=0.00016512100410182
102        call getin_p("omega",omeg)
103        ! Gravity
104        g=10.44
105        call getin_p("g",g)
106        !mugaz=2.34 !EM: does not give cpp=11500
107        mugaz=2.53 ! with this value of mugaz, cpp=11500
108        call getin_p("mugaz",mugaz)
109        ! kappa
110        rcp=0.2857143
111        call getin_p("kappa",rcp)
112        cpp=(8.314511/(mugaz/1000.0))/rcp
113        call getin_p("cpp",cpp)
114!        write(*,*) "tabfi: cpp=",cpp
115        ! length (s) of a "standard" day
116        daysec=38052
117        call getin_p("day_length",daysec)
118        ! physics time step (s) ! not sure we need this here
119        dtphys=19026
120        ! length of year, in standard days
121        year_day=24430
122        ! Orbital parameters
123        periastr=9.02151966094971
124        apoastr=10.054479598999
125        peri_day=19280
126        obliquit=26.7299995422363
127        ! Other parameters some physical paréametrizations need
128        z0=1e-2
129        lmixmin=30
130        emin_turb=1.e-6
131        albedice(:)=0
132        emisice(:)=0
133        emissiv=0
134        iceradius(:)=1.e-6
135        dtemisice(:)=0
136        volcapa=1000000
137c-----------------------------------------------------------------------
138c       Save some constants for later use (as routine arguments)
139c-----------------------------------------------------------------------
140        p_omeg = omeg
141        p_g = g
142        p_cpp = cpp
143        p_mugaz = mugaz
144        p_daysec = daysec
145        p_rad=rad
146
147      ELSE
148c-----------------------------------------------------------------------
149c  Initialization of physical constants by reading array tab_cntrl(:)
150c               which contains these parameters (nid != 0 case)
151c-----------------------------------------------------------------------
152c Read 'controle' array
153c
154
155       call get_var("controle",tab_cntrl,found)
156       if (.not.found) then
157         write(*,*)"tabfi: Failed reading <controle> array"
158         call abort
159       else
160         write(*,*)'tabfi: tab_cntrl',tab_cntrl
161       endif
162c
163c  Initialization of some physical constants
164c informations on physics grid
165!      if(ngrid.ne.tab_cntrl(tab0+1)) then
166!         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
167!         print*,tab_cntrl(tab0+1),ngrid
168!      endif
169      lmax = nint(tab_cntrl(tab0+2))
170      day_ini = tab_cntrl(tab0+3)
171      time = tab_cntrl(tab0+4)
172      write (*,*) 'IN tabfi day_ini=',day_ini
173c Informations about planet for dynamics and physics
174      rad = tab_cntrl(tab0+5)
175      omeg = tab_cntrl(tab0+6)
176      g = tab_cntrl(tab0+7)
177      mugaz = tab_cntrl(tab0+8)
178      rcp = tab_cntrl(tab0+9)
179      cpp=(8.314511/(mugaz/1000.0))/rcp
180      daysec = tab_cntrl(tab0+10)
181      dtphys = tab_cntrl(tab0+11)
182c Informations about planet for the physics only
183      year_day = tab_cntrl(tab0+14)
184      periastr = tab_cntrl(tab0+15)
185      apoastr = tab_cntrl(tab0+16)
186      peri_day = tab_cntrl(tab0+17)
187      obliquit = tab_cntrl(tab0+18)
188c boundary layer and turbeulence
189      z0 = tab_cntrl(tab0+19)
190      lmixmin = tab_cntrl(tab0+20)
191      emin_turb = tab_cntrl(tab0+21)
192c optical properties of polar caps and ground emissivity
193      albedice(1)= tab_cntrl(tab0+22)
194      albedice(2)= tab_cntrl(tab0+23)
195      emisice(1) = tab_cntrl(tab0+24)
196      emisice(2) = tab_cntrl(tab0+25)
197      emissiv    = tab_cntrl(tab0+26)
198      iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north)
199      iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south)
200      dtemisice(1)= tab_cntrl(tab0+33) !time scale for snow metamorphism (north)
201      dtemisice(2)= tab_cntrl(tab0+34) !time scale for snow metamorphism (south)
202c soil properties
203      volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
204c-----------------------------------------------------------------------
205c       Save some constants for later use (as routine arguments)
206c-----------------------------------------------------------------------
207      p_omeg = omeg
208      p_g = g
209      p_cpp = cpp
210      p_mugaz = mugaz
211      p_daysec = daysec
212      p_rad=rad
213
214      ENDIF    ! end of (nid = 0)
215
216c-----------------------------------------------------------------------
217c       Write physical constants to output before modifying them
218c-----------------------------------------------------------------------
219 
220   6  FORMAT(a20,e15.6,e15.6)
221   5  FORMAT(a20,f12.2,f12.2)
222 
223      write(*,*) '*****************************************************'
224      write(*,*) 'Reading tab_cntrl when calling tabfi before changes'
225      write(*,*) '*****************************************************'
226      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
227      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
228      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
229      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
230      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
231      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
232      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
233      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
234      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
235      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
236
237      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
238      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
239      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
240      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
241      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
242
243      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
244      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
245      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
246
247      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
248      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
249      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
250      write(*,5) '(22)    albedice(1)',tab_cntrl(tab0+22),albedice(1)
251      write(*,5) '(23)    albedice(2)',tab_cntrl(tab0+23),albedice(2)
252      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
253      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
254      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
255      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
256
257      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
258
259      write(*,*)
260      write(*,*) 'Lmodif in tabfi!!!!!!!',Lmodif
261
262c-----------------------------------------------------------------------
263c        Modifications...
264! NB: Modifying controls should only be done by newstart, and in seq mode
265      if ((Lmodif.eq.1).and.is_parallel) then
266        write(*,*) "tabfi: Error modifying tab_control should",
267     &             " only happen in serial mode (eg: by newstart)"
268        stop
269      endif
270c-----------------------------------------------------------------------
271
272      IF(Lmodif.eq.1) then
273
274      write(*,*)
275      write(*,*) 'Change values in tab_cntrl ? :'
276      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
277      write(*,*) '(Current values given above)'
278      write(*,*)
279      write(*,*) '(3)          day_ini : Initial day (=0 at Ls=0)'
280      write(*,*) '(19)              z0 :  surface roughness (m)'
281      write(*,*) '(21)       emin_turb :  minimal energy (PBL)'
282      write(*,*) '(20)         lmixmin : mixing length (PBL)'
283      write(*,*) '(26)         emissiv : ground emissivity'
284      write(*,*) '(24 et 25)   emisice : CO2 ice max emissivity '
285      write(*,*) '(22 et 23)  albedice : CO2 ice cap albedos'
286      write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow'
287      write(*,*) '(33 et 34) dtemisice : time scale for snow',
288     &           'metamorphism'
289      write(*,*) '(35)      volcapa : soil volumetric heat capacity'
290      write(*,*) '(18)     obliquit : planet obliquity (deg)'
291      write(*,*) '(17)     peri_day : periastron date (sols since Ls=0)'
292      write(*,*) '(15)     periastr : min. star-planet dist (UA)'
293      write(*,*) '(16)     apoastr  : max. star-planet (UA)'
294      write(*,*) '(14)     year_day : length of year (in sols)'
295      write(*,*) '(5)      rad      : radius of the planet (m)'
296      write(*,*) '(6)      omeg     : planet rotation rate (rad/s)'
297      write(*,*) '(7)      g        : gravity (m/s2)'
298      write(*,*) '(8)      mugaz    : molecular mass '
299      write(*,*) '                       of the atmosphere (g/mol)'
300      write(*,*) '(9)      rcp      : r/Cp'
301      write(*,*) '(8)+(9)  calc_cpp_mugaz : r/Cp and mugaz '
302      write(*,*) '                 computed from gases.def'
303      write(*,*) '(10)     daysec   : length of a sol (s)'
304      write(*,*)
305 
306 
307      do while(modif(1:1).ne.'hello')
308        write(*,*)
309        write(*,*)
310        write(*,*) 'Changes to perform ?'
311        write(*,*) '   (enter keyword or return )'
312        write(*,*)
313        read(*,fmt='(a20)') modif
314        if (modif(1:1) .eq. ' ') goto 999
315 
316        write(*,*)
317        write(*,*) modif(1:len_trim(modif)) , ' : '
318
319        if (modif(1:len_trim(modif)) .eq. 'day_ini') then
320          write(*,*) 'current value:',day_ini
321          write(*,*) 'enter new value:'
322 101      read(*,*,iostat=ierr) day_ini
323          if(ierr.ne.0) goto 101
324          write(*,*) ' '
325          write(*,*) 'day_ini (new value):',day_ini
326
327        else if (modif(1:len_trim(modif)) .eq. 'z0') then
328          write(*,*) 'current value:',z0
329          write(*,*) 'enter new value:'
330 102      read(*,*,iostat=ierr) z0
331          if(ierr.ne.0) goto 102
332          write(*,*) ' '
333          write(*,*) ' z0 (new value):',z0
334
335        else if (modif(1:len_trim(modif)) .eq. 'emin_turb') then
336          write(*,*) 'current value:',emin_turb
337          write(*,*) 'enter new value:'
338 103      read(*,*,iostat=ierr) emin_turb
339          if(ierr.ne.0) goto 103
340          write(*,*) ' '
341          write(*,*) ' emin_turb (new value):',emin_turb
342
343        else if (modif(1:len_trim(modif)) .eq. 'lmixmin') then
344          write(*,*) 'current value:',lmixmin
345          write(*,*) 'enter new value:'
346 104      read(*,*,iostat=ierr) lmixmin
347          if(ierr.ne.0) goto 104
348          write(*,*) ' '
349          write(*,*) ' lmixmin (new value):',lmixmin
350
351        else if (modif(1:len_trim(modif)) .eq. 'emissiv') then
352          write(*,*) 'current value:',emissiv
353          write(*,*) 'enter new value:'
354 105      read(*,*,iostat=ierr) emissiv
355          if(ierr.ne.0) goto 105
356          write(*,*) ' '
357          write(*,*) ' emissiv (new value):',emissiv
358
359        else if (modif(1:len_trim(modif)) .eq. 'emisice') then
360          write(*,*) 'current value emisice(1) North:',emisice(1)
361          write(*,*) 'enter new value:'
362 106      read(*,*,iostat=ierr) emisice(1)
363          if(ierr.ne.0) goto 106
364          write(*,*) 
365          write(*,*) ' emisice(1) (new value):',emisice(1)
366          write(*,*)
367
368          write(*,*) 'current value emisice(2) South:',emisice(2)
369          write(*,*) 'enter new value:'
370 107      read(*,*,iostat=ierr) emisice(2)
371          if(ierr.ne.0) goto 107
372          write(*,*) 
373          write(*,*) ' emisice(2) (new value):',emisice(2)
374
375        else if (modif(1:len_trim(modif)) .eq. 'albedice') then
376          write(*,*) 'current value albedice(1) North:',albedice(1)
377          write(*,*) 'enter new value:'
378 108      read(*,*,iostat=ierr) albedice(1)
379          if(ierr.ne.0) goto 108
380          write(*,*) 
381          write(*,*) ' albedice(1) (new value):',albedice(1)
382          write(*,*)
383
384          write(*,*) 'current value albedice(2) South:',albedice(2)
385          write(*,*) 'enter new value:'
386 109      read(*,*,iostat=ierr) albedice(2)
387          if(ierr.ne.0) goto 109
388          write(*,*) 
389          write(*,*) ' albedice(2) (new value):',albedice(2)
390
391        else if (modif(1:len_trim(modif)) .eq. 'iceradius') then
392          write(*,*) 'current value iceradius(1) North:',iceradius(1)
393          write(*,*) 'enter new value:'
394 110      read(*,*,iostat=ierr) iceradius(1)
395          if(ierr.ne.0) goto 110
396          write(*,*) 
397          write(*,*) ' iceradius(1) (new value):',iceradius(1)
398          write(*,*)
399
400          write(*,*) 'current value iceradius(2) South:',iceradius(2)
401          write(*,*) 'enter new value:'
402 111      read(*,*,iostat=ierr) iceradius(2)
403          if(ierr.ne.0) goto 111
404          write(*,*) 
405          write(*,*) ' iceradius(2) (new value):',iceradius(2)
406
407        else if (modif(1:len_trim(modif)) .eq. 'dtemisice') then
408          write(*,*) 'current value dtemisice(1) North:',dtemisice(1)
409          write(*,*) 'enter new value:'
410 112      read(*,*,iostat=ierr) dtemisice(1)
411          if(ierr.ne.0) goto 112
412          write(*,*) 
413          write(*,*) ' dtemisice(1) (new value):',dtemisice(1)
414          write(*,*)
415
416          write(*,*) 'current value dtemisice(2) South:',dtemisice(2)
417          write(*,*) 'enter new value:'
418 113      read(*,*,iostat=ierr) dtemisice(2)
419          if(ierr.ne.0) goto 113
420          write(*,*) 
421          write(*,*) ' dtemisice(2) (new value):',dtemisice(2)
422
423        else if (modif(1:len_trim(modif)) .eq. 'obliquit') then
424          write(*,*) 'current value:',obliquit
425          write(*,*) 'obliquit should be 25.19 on current Mars'
426          write(*,*) 'enter new value:'
427 115      read(*,*,iostat=ierr) obliquit
428          if(ierr.ne.0) goto 115
429          write(*,*) 
430          write(*,*) ' obliquit (new value):',obliquit
431
432        else if (modif(1:len_trim(modif)) .eq. 'peri_day') then
433          write(*,*) 'current value:',peri_day
434          write(*,*) 'peri_day should be 485 on current Mars'
435          write(*,*) 'enter new value:'
436 116      read(*,*,iostat=ierr) peri_day
437          if(ierr.ne.0) goto 116
438          write(*,*) 
439          write(*,*) ' peri_day (new value):',peri_day
440
441        else if (modif(1:len_trim(modif)) .eq. 'periastr') then
442          write(*,*) 'current value:',periastr
443          write(*,*) 'periastr should be 206.66 on present-day Mars'
444          write(*,*) 'enter new value:'
445 117      read(*,*,iostat=ierr) periastr
446          if(ierr.ne.0) goto 117
447          write(*,*) 
448          write(*,*) ' periastr (new value):',periastr
449 
450        else if (modif(1:len_trim(modif)) .eq. 'apoastr') then
451          write(*,*) 'current value:',apoastr
452          write(*,*) 'apoastr should be 249.22 on present-day Mars'
453          write(*,*) 'enter new value:'
454 118      read(*,*,iostat=ierr) apoastr
455          if(ierr.ne.0) goto 118
456          write(*,*) 
457          write(*,*) ' apoastr (new value):',apoastr
458 
459        else if (modif(1:len_trim(modif)) .eq. 'volcapa') then
460          write(*,*) 'current value:',volcapa
461          write(*,*) 'enter new value:'
462 119      read(*,*,iostat=ierr) volcapa
463          if(ierr.ne.0) goto 119
464          write(*,*) 
465          write(*,*) ' volcapa (new value):',volcapa
466       
467        else if (modif(1:len_trim(modif)).eq.'rad') then
468          write(*,*) 'current value:',rad
469          write(*,*) 'enter new value:'
470 120      read(*,*,iostat=ierr) rad
471          if(ierr.ne.0) goto 120
472          write(*,*) 
473          write(*,*) ' rad (new value):',rad
474
475        else if (modif(1:len_trim(modif)).eq.'omeg') then
476          write(*,*) 'current value:',omeg
477          write(*,*) 'enter new value:'
478 121      read(*,*,iostat=ierr) omeg
479          if(ierr.ne.0) goto 121
480          write(*,*) 
481          write(*,*) ' omeg (new value):',omeg
482       
483        else if (modif(1:len_trim(modif)).eq.'g') then
484          write(*,*) 'current value:',g
485          write(*,*) 'enter new value:'
486 122      read(*,*,iostat=ierr) g
487          if(ierr.ne.0) goto 122
488          write(*,*) 
489          write(*,*) ' g (new value):',g
490
491        else if (modif(1:len_trim(modif)).eq.'mugaz') then
492          write(*,*) 'current value:',mugaz
493          write(*,*) 'enter new value:'
494 123      read(*,*,iostat=ierr) mugaz
495          if(ierr.ne.0) goto 123
496          write(*,*) 
497          write(*,*) ' mugaz (new value):',mugaz
498          r=8.314511/(mugaz/1000.0)
499          write(*,*) ' R (new value):',r
500
501        else if (modif(1:len_trim(modif)).eq.'rcp') then
502          write(*,*) 'current value:',rcp
503          write(*,*) 'enter new value:'
504 124      read(*,*,iostat=ierr) rcp
505          if(ierr.ne.0) goto 124
506          write(*,*) 
507          write(*,*) ' rcp (new value):',rcp
508          r=8.314511/(mugaz/1000.0)
509          cpp=r/rcp
510          write(*,*) ' cpp (new value):',cpp
511
512        else if (modif(1:len_trim(modif)).eq.'calc_cpp_mugaz') then
513          write(*,*) 'current value rcp, mugaz:',rcp,mugaz
514          check_cpp_match=.false.
515          force_cpp=.false.
516          call su_gases
517          call calc_cpp_mugaz
518          write(*,*) 
519          write(*,*) ' cpp (new value):',cpp
520          write(*,*) ' mugaz (new value):',mugaz
521          r=8.314511/(mugaz/1000.0)
522          rcp=r/cpp
523          write(*,*) ' rcp (new value):',rcp
524         
525        else if (modif(1:len_trim(modif)).eq.'daysec') then
526          write(*,*) 'current value:',daysec
527          write(*,*) 'enter new value:'
528 125      read(*,*,iostat=ierr) daysec
529          if(ierr.ne.0) goto 125
530          write(*,*) 
531          write(*,*) ' daysec (new value):',daysec
532
533!         added by RW!
534        else if (modif(1:len_trim(modif)).eq.'year_day') then
535          write(*,*) 'current value:',year_day
536          write(*,*) 'enter new value:' 
537 126      read(*,*,iostat=ierr) year_day
538          if(ierr.ne.0) goto 126
539          write(*,*)
540          write(*,*) ' year_day (new value):',year_day
541
542        endif
543      enddo ! of do while(modif(1:1).ne.'hello')
544
545 999  continue
546
547c-----------------------------------------------------------------------
548c       Write values of physical constants after modifications
549c-----------------------------------------------------------------------
550 
551      write(*,*) '*****************************************************'
552      write(*,*) 'Reading tab_cntrl when calling tabfi AFTER changes'
553      write(*,*) '*****************************************************'
554      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
555      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
556      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
557      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
558      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
559      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
560      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
561      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
562      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
563      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
564 
565      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
566      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
567      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
568      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
569      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
570 
571      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
572      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
573      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
574 
575      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
576      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
577      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
578      write(*,5) '(22)    albedice(1)',tab_cntrl(tab0+22),albedice(1)
579      write(*,5) '(23)    albedice(2)',tab_cntrl(tab0+23),albedice(2)
580      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
581      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
582      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
583      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
584 
585      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
586
587      write(*,*) 
588      write(*,*) 
589
590      ENDIF                     !       of if (Lmodif == 1)
591
592c-----------------------------------------------------------------------
593c       Save some constants for later use (as routine arguments)
594c-----------------------------------------------------------------------
595      p_omeg = omeg
596      p_g = g
597      p_cpp = cpp
598      p_mugaz = mugaz
599      p_daysec = daysec
600      p_rad=rad
601
602
603      end
Note: See TracBrowser for help on using the repository browser.