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

Last change on this file since 1057 was 313, checked in by ymipsl, 10 years ago
  • implement splitting of XIOS file for lmdz physics
  • Termination is done properly in parallel by calling MPI_ABORT instead of abort or stop

YM

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_physiq
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        CALL abort_physiq
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.