source: trunk/LMDZ.GENERIC/libf/phystd/tabfi.F @ 1145

Last change on this file since 1145 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
File size: 20.5 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 nlayermx)
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
47
48      USE surfdat_h
49      use comsoil_h
50      USE comgeomfi_h
51
52      implicit none
53 
54#include "dimensions.h"
55#include "dimphys.h"
56#include "comcstfi.h"
57#include "planete.h"
58#include "netcdf.inc"
59#include "callkeys.h"
60
61c-----------------------------------------------------------------------
62c   Declarations
63c-----------------------------------------------------------------------
64
65      INTEGER ngrid
66
67c Arguments
68c ---------
69      INTEGER nid,nvarid,tab0
70      INTEGER*4 day_ini
71      INTEGER Lmodif
72      INTEGER lmax
73      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
74
75c Variables
76c ---------
77      INTEGER length
78      parameter (length = 100)
79      REAL tab_cntrl(length) ! array in which are stored the run's parameters
80      INTEGER  ierr
81      INTEGER size
82      CHARACTER modif*20
83
84c-----------------------------------------------------------------------
85c  Initialization of physical constants by reading array tab_cntrl(:)
86c               which contains these parameters (nid != 0 case)
87c-----------------------------------------------------------------------
88c Read 'controle' array
89c
90      ierr = NF_INQ_VARID (nid, "controle", nvarid)
91      IF (ierr .NE. NF_NOERR) THEN
92         PRINT*, "Tabfi: Could not find <controle> data"
93         CALL abort
94      ENDIF
95#ifdef NC_DOUBLE
96      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
97#else
98      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
99#endif
100      IF (ierr .NE. NF_NOERR) THEN
101         PRINT*, "Tabfi: Failed reading <controle> array"
102         CALL abort
103      ENDIF
104
105       print*,'tabfi: tab_cntrl',tab_cntrl
106c
107c  Initialization of some physical constants
108c informations on physics grid
109      if(ngrid.ne.tab_cntrl(tab0+1)) then
110         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
111         print*,tab_cntrl(tab0+1),ngrid
112      endif
113      lmax = nint(tab_cntrl(tab0+2))
114      day_ini = tab_cntrl(tab0+3)
115      time = tab_cntrl(tab0+4)
116      write (*,*) 'IN tabfi day_ini=',day_ini
117c Informations about planet for dynamics and physics
118      rad = tab_cntrl(tab0+5)
119      omeg = tab_cntrl(tab0+6)
120      g = tab_cntrl(tab0+7)
121      mugaz = tab_cntrl(tab0+8)
122      rcp = tab_cntrl(tab0+9)
123      cpp=(8.314511/(mugaz/1000.0))/rcp
124      daysec = tab_cntrl(tab0+10)
125      dtphys = tab_cntrl(tab0+11)
126c Informations about planet for the physics only
127      year_day = tab_cntrl(tab0+14)
128      periastr = tab_cntrl(tab0+15)
129      apoastr = tab_cntrl(tab0+16)
130      peri_day = tab_cntrl(tab0+17)
131      obliquit = tab_cntrl(tab0+18)
132c boundary layer and turbeulence
133      z0 = tab_cntrl(tab0+19)
134      lmixmin = tab_cntrl(tab0+20)
135      emin_turb = tab_cntrl(tab0+21)
136c optical properties of polar caps and ground emissivity
137      albedice(1)= tab_cntrl(tab0+22)
138      albedice(2)= tab_cntrl(tab0+23)
139      emisice(1) = tab_cntrl(tab0+24)
140      emisice(2) = tab_cntrl(tab0+25)
141      emissiv    = tab_cntrl(tab0+26)
142      iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north)
143      iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south)
144      dtemisice(1)= tab_cntrl(tab0+33) !time scale for snow metamorphism (north)
145      dtemisice(2)= tab_cntrl(tab0+34) !time scale for snow metamorphism (south)
146c soil properties
147      volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
148
149
150
151
152c-----------------------------------------------------------------------
153c       Save some constants for later use (as routine arguments)
154c-----------------------------------------------------------------------
155      p_omeg = omeg
156      p_g = g
157      p_cpp = cpp
158      p_mugaz = mugaz
159      p_daysec = daysec
160      p_rad=rad
161
162
163c-----------------------------------------------------------------------
164c       Write physical constants to output before modifying them
165c-----------------------------------------------------------------------
166 
167   6  FORMAT(a20,e15.6,e15.6)
168   5  FORMAT(a20,f12.2,f12.2)
169 
170      write(*,*) '*****************************************************'
171      write(*,*) 'Reading tab_cntrl when calling tabfi before changes'
172      write(*,*) '*****************************************************'
173      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
174      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
175      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
176      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
177      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
178      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
179      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
180      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
181      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
182      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
183
184      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
185      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
186      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
187      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
188      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
189
190      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
191      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
192      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
193
194      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
195      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
196      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
197      write(*,5) '(22)    albedice(1)',tab_cntrl(tab0+22),albedice(1)
198      write(*,5) '(23)    albedice(2)',tab_cntrl(tab0+23),albedice(2)
199      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
200      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
201      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
202      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
203
204      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
205
206      write(*,*)
207      write(*,*) 'Lmodif in tabfi!!!!!!!',Lmodif
208
209c-----------------------------------------------------------------------
210c        Modifications...
211c-----------------------------------------------------------------------
212
213      IF(Lmodif.eq.1) then
214
215      write(*,*)
216      write(*,*) 'Change values in tab_cntrl ? :'
217      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
218      write(*,*) '(Current values given above)'
219      write(*,*)
220      write(*,*) '(3)          day_ini : Initial day (=0 at Ls=0)'
221      write(*,*) '(19)              z0 :  surface roughness (m)'
222      write(*,*) '(21)       emin_turb :  minimal energy (PBL)'
223      write(*,*) '(20)         lmixmin : mixing length (PBL)'
224      write(*,*) '(26)         emissiv : ground emissivity'
225      write(*,*) '(24 et 25)   emisice : CO2 ice max emissivity '
226      write(*,*) '(22 et 23)  albedice : CO2 ice cap albedos'
227      write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow'
228      write(*,*) '(33 et 34) dtemisice : time scale for snow',
229     &           'metamorphism'
230      write(*,*) '(35)      volcapa : soil volumetric heat capacity'
231      write(*,*) '(18)     obliquit : planet obliquity (deg)'
232      write(*,*) '(17)     peri_day : periastron date (sols since Ls=0)'
233      write(*,*) '(15)     periastr : min. star-planet dist (UA)'
234      write(*,*) '(16)     apoastr  : max. star-planet (UA)'
235      write(*,*) '(14)     year_day : length of year (in sols)'
236      write(*,*) '(5)      rad      : radius of the planet (m)'
237      write(*,*) '(6)      omeg     : planet rotation rate (rad/s)'
238      write(*,*) '(7)      g        : gravity (m/s2)'
239      write(*,*) '(8)      mugaz    : molecular mass '
240      write(*,*) '                       of the atmosphere (g/mol)'
241      write(*,*) '(9)      rcp      : r/Cp'
242      write(*,*) '(8)+(9)  calc_cpp_mugaz : r/Cp and mugaz '
243      write(*,*) '                 computed from gases.def'
244      write(*,*) '(10)     daysec   : length of a sol (s)'
245      write(*,*)
246 
247 
248      do while(modif(1:1).ne.'hello')
249        write(*,*)
250        write(*,*)
251        write(*,*) 'Changes to perform ?'
252        write(*,*) '   (enter keyword or return )'
253        write(*,*)
254        read(*,fmt='(a20)') modif
255        if (modif(1:1) .eq. ' ') goto 999
256 
257        write(*,*)
258        write(*,*) modif(1:len_trim(modif)) , ' : '
259
260        if (modif(1:len_trim(modif)) .eq. 'day_ini') then
261          write(*,*) 'current value:',day_ini
262          write(*,*) 'enter new value:'
263 101      read(*,*,iostat=ierr) day_ini
264          if(ierr.ne.0) goto 101
265          write(*,*) ' '
266          write(*,*) 'day_ini (new value):',day_ini
267
268        else if (modif(1:len_trim(modif)) .eq. 'z0') then
269          write(*,*) 'current value:',z0
270          write(*,*) 'enter new value:'
271 102      read(*,*,iostat=ierr) z0
272          if(ierr.ne.0) goto 102
273          write(*,*) ' '
274          write(*,*) ' z0 (new value):',z0
275
276        else if (modif(1:len_trim(modif)) .eq. 'emin_turb') then
277          write(*,*) 'current value:',emin_turb
278          write(*,*) 'enter new value:'
279 103      read(*,*,iostat=ierr) emin_turb
280          if(ierr.ne.0) goto 103
281          write(*,*) ' '
282          write(*,*) ' emin_turb (new value):',emin_turb
283
284        else if (modif(1:len_trim(modif)) .eq. 'lmixmin') then
285          write(*,*) 'current value:',lmixmin
286          write(*,*) 'enter new value:'
287 104      read(*,*,iostat=ierr) lmixmin
288          if(ierr.ne.0) goto 104
289          write(*,*) ' '
290          write(*,*) ' lmixmin (new value):',lmixmin
291
292        else if (modif(1:len_trim(modif)) .eq. 'emissiv') then
293          write(*,*) 'current value:',emissiv
294          write(*,*) 'enter new value:'
295 105      read(*,*,iostat=ierr) emissiv
296          if(ierr.ne.0) goto 105
297          write(*,*) ' '
298          write(*,*) ' emissiv (new value):',emissiv
299
300        else if (modif(1:len_trim(modif)) .eq. 'emisice') then
301          write(*,*) 'current value emisice(1) North:',emisice(1)
302          write(*,*) 'enter new value:'
303 106      read(*,*,iostat=ierr) emisice(1)
304          if(ierr.ne.0) goto 106
305          write(*,*)
306          write(*,*) ' emisice(1) (new value):',emisice(1)
307          write(*,*)
308
309          write(*,*) 'current value emisice(2) South:',emisice(2)
310          write(*,*) 'enter new value:'
311 107      read(*,*,iostat=ierr) emisice(2)
312          if(ierr.ne.0) goto 107
313          write(*,*)
314          write(*,*) ' emisice(2) (new value):',emisice(2)
315
316        else if (modif(1:len_trim(modif)) .eq. 'albedice') then
317          write(*,*) 'current value albedice(1) North:',albedice(1)
318          write(*,*) 'enter new value:'
319 108      read(*,*,iostat=ierr) albedice(1)
320          if(ierr.ne.0) goto 108
321          write(*,*)
322          write(*,*) ' albedice(1) (new value):',albedice(1)
323          write(*,*)
324
325          write(*,*) 'current value albedice(2) South:',albedice(2)
326          write(*,*) 'enter new value:'
327 109      read(*,*,iostat=ierr) albedice(2)
328          if(ierr.ne.0) goto 109
329          write(*,*)
330          write(*,*) ' albedice(2) (new value):',albedice(2)
331
332        else if (modif(1:len_trim(modif)) .eq. 'iceradius') then
333          write(*,*) 'current value iceradius(1) North:',iceradius(1)
334          write(*,*) 'enter new value:'
335 110      read(*,*,iostat=ierr) iceradius(1)
336          if(ierr.ne.0) goto 110
337          write(*,*)
338          write(*,*) ' iceradius(1) (new value):',iceradius(1)
339          write(*,*)
340
341          write(*,*) 'current value iceradius(2) South:',iceradius(2)
342          write(*,*) 'enter new value:'
343 111      read(*,*,iostat=ierr) iceradius(2)
344          if(ierr.ne.0) goto 111
345          write(*,*)
346          write(*,*) ' iceradius(2) (new value):',iceradius(2)
347
348        else if (modif(1:len_trim(modif)) .eq. 'dtemisice') then
349          write(*,*) 'current value dtemisice(1) North:',dtemisice(1)
350          write(*,*) 'enter new value:'
351 112      read(*,*,iostat=ierr) dtemisice(1)
352          if(ierr.ne.0) goto 112
353          write(*,*)
354          write(*,*) ' dtemisice(1) (new value):',dtemisice(1)
355          write(*,*)
356
357          write(*,*) 'current value dtemisice(2) South:',dtemisice(2)
358          write(*,*) 'enter new value:'
359 113      read(*,*,iostat=ierr) dtemisice(2)
360          if(ierr.ne.0) goto 113
361          write(*,*)
362          write(*,*) ' dtemisice(2) (new value):',dtemisice(2)
363
364        else if (modif(1:len_trim(modif)) .eq. 'obliquit') then
365          write(*,*) 'current value:',obliquit
366          write(*,*) 'obliquit should be 25.19 on current Mars'
367          write(*,*) 'enter new value:'
368 115      read(*,*,iostat=ierr) obliquit
369          if(ierr.ne.0) goto 115
370          write(*,*)
371          write(*,*) ' obliquit (new value):',obliquit
372
373        else if (modif(1:len_trim(modif)) .eq. 'peri_day') then
374          write(*,*) 'current value:',peri_day
375          write(*,*) 'peri_day should be 485 on current Mars'
376          write(*,*) 'enter new value:'
377 116      read(*,*,iostat=ierr) peri_day
378          if(ierr.ne.0) goto 116
379          write(*,*)
380          write(*,*) ' peri_day (new value):',peri_day
381
382        else if (modif(1:len_trim(modif)) .eq. 'periastr') then
383          write(*,*) 'current value:',periastr
384          write(*,*) 'periastr should be 206.66 on present-day Mars'
385          write(*,*) 'enter new value:'
386 117      read(*,*,iostat=ierr) periastr
387          if(ierr.ne.0) goto 117
388          write(*,*)
389          write(*,*) ' periastr (new value):',periastr
390 
391        else if (modif(1:len_trim(modif)) .eq. 'apoastr') then
392          write(*,*) 'current value:',apoastr
393          write(*,*) 'apoastr should be 249.22 on present-day Mars'
394          write(*,*) 'enter new value:'
395 118      read(*,*,iostat=ierr) apoastr
396          if(ierr.ne.0) goto 118
397          write(*,*)
398          write(*,*) ' apoastr (new value):',apoastr
399 
400        else if (modif(1:len_trim(modif)) .eq. 'volcapa') then
401          write(*,*) 'current value:',volcapa
402          write(*,*) 'enter new value:'
403 119      read(*,*,iostat=ierr) volcapa
404          if(ierr.ne.0) goto 119
405          write(*,*)
406          write(*,*) ' volcapa (new value):',volcapa
407       
408        else if (modif(1:len_trim(modif)).eq.'rad') then
409          write(*,*) 'current value:',rad
410          write(*,*) 'enter new value:'
411 120      read(*,*,iostat=ierr) rad
412          if(ierr.ne.0) goto 120
413          write(*,*)
414          write(*,*) ' rad (new value):',rad
415
416        else if (modif(1:len_trim(modif)).eq.'omeg') then
417          write(*,*) 'current value:',omeg
418          write(*,*) 'enter new value:'
419 121      read(*,*,iostat=ierr) omeg
420          if(ierr.ne.0) goto 121
421          write(*,*)
422          write(*,*) ' omeg (new value):',omeg
423       
424        else if (modif(1:len_trim(modif)).eq.'g') then
425          write(*,*) 'current value:',g
426          write(*,*) 'enter new value:'
427 122      read(*,*,iostat=ierr) g
428          if(ierr.ne.0) goto 122
429          write(*,*)
430          write(*,*) ' g (new value):',g
431
432        else if (modif(1:len_trim(modif)).eq.'mugaz') then
433          write(*,*) 'current value:',mugaz
434          write(*,*) 'enter new value:'
435 123      read(*,*,iostat=ierr) mugaz
436          if(ierr.ne.0) goto 123
437          write(*,*)
438          write(*,*) ' mugaz (new value):',mugaz
439          r=8.314511/(mugaz/1000.0)
440          write(*,*) ' R (new value):',r
441
442        else if (modif(1:len_trim(modif)).eq.'rcp') then
443          write(*,*) 'current value:',rcp
444          write(*,*) 'enter new value:'
445 124      read(*,*,iostat=ierr) rcp
446          if(ierr.ne.0) goto 124
447          write(*,*)
448          write(*,*) ' rcp (new value):',rcp
449          r=8.314511/(mugaz/1000.0)
450          cpp=r/rcp
451          write(*,*) ' cpp (new value):',cpp
452
453        else if (modif(1:len_trim(modif)).eq.'calc_cpp_mugaz') then
454          write(*,*) 'current value rcp, mugaz:',rcp,mugaz
455          check_cpp_match=.false.
456          force_cpp=.false.
457          call su_gases
458          call calc_cpp_mugaz
459          write(*,*)
460          write(*,*) ' cpp (new value):',cpp
461          write(*,*) ' mugaz (new value):',mugaz
462          r=8.314511/(mugaz/1000.0)
463          rcp=r/cpp
464          write(*,*) ' rcp (new value):',rcp
465         
466        else if (modif(1:len_trim(modif)).eq.'daysec') then
467          write(*,*) 'current value:',daysec
468          write(*,*) 'enter new value:'
469 125      read(*,*,iostat=ierr) daysec
470          if(ierr.ne.0) goto 125
471          write(*,*)
472          write(*,*) ' daysec (new value):',daysec
473
474!         added by RW!
475        else if (modif(1:len_trim(modif)).eq.'year_day') then
476          write(*,*) 'current value:',year_day
477          write(*,*) 'enter new value:'
478 126      read(*,*,iostat=ierr) year_day
479          if(ierr.ne.0) goto 126
480          write(*,*)
481          write(*,*) ' year_day (new value):',year_day
482
483        endif
484      enddo ! of do while(modif(1:1).ne.'hello')
485
486 999  continue
487
488c-----------------------------------------------------------------------
489c       Write values of physical constants after modifications
490c-----------------------------------------------------------------------
491 
492      write(*,*) '*****************************************************'
493      write(*,*) 'Reading tab_cntrl when calling tabfi AFTER changes'
494      write(*,*) '*****************************************************'
495      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
496      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
497      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
498      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
499      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
500      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
501      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
502      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
503      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
504      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
505 
506      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
507      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
508      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
509      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
510      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
511 
512      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
513      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
514      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
515 
516      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
517      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
518      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
519      write(*,5) '(22)    albedice(1)',tab_cntrl(tab0+22),albedice(1)
520      write(*,5) '(23)    albedice(2)',tab_cntrl(tab0+23),albedice(2)
521      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
522      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
523      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
524      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
525 
526      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
527
528      write(*,*) 
529      write(*,*)
530
531      ENDIF                     !       of if (Lmodif == 1)
532
533c-----------------------------------------------------------------------
534c       Save some constants for later use (as routine arguments)
535c-----------------------------------------------------------------------
536      p_omeg = omeg
537      p_g = g
538      p_cpp = cpp
539      p_mugaz = mugaz
540      p_daysec = daysec
541      p_rad=rad
542
543
544      end
Note: See TracBrowser for help on using the repository browser.