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

Last change on this file since 310 was 298, checked in by milmd, 10 years ago

Less output messages are written. On 20000 cores it is better. In LMDZ, only master of MPI and OpenMP can write.

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