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

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

Bug correction in the initialization of specific gas constant in Saturn physics.
EM

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