c======================================================================= SUBROUTINE tabfi(nid,Lmodif,tab0,day_ini,lmax,p_rad, . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) c======================================================================= c c C. Hourdin 15/11/96 c c Object: Lecture du tab_cntrl physique dans un fichier c ------ et initialisation des constantes physiques c c Arguments: c ---------- c c Inputs: c ------ c c - nid: unitne logique du fichier ou on va lire le tab_cntrl c (ouvert dans le programme appellant) c c si nid=0: c pas de lecture du tab_cntrl mais c Valeurs par default des constantes physiques c c - tab0: Offset de tab_cntrl a partir duquel sont ranges c les parametres physiques (50 pour start_archive) c c - Lmodif: si on souhaite modifier les constantes Lmodif = 1 = TRUE c c c Outputs: c -------- c c - day_ini: tab_cntrl(tab0+3) (Dans les cas ou l'on souhaite c comparer avec le day_ini dynamique) c c - lmax: tab_cntrl(tab0+2) (pour test avec nlayermx) c c - p_rad c - p_omeg ! c - p_g ! Constantes physiques ayant des c - p_mugaz ! homonymes dynamiques c - p_daysec ! c c======================================================================= ! to use 'getin' use ioipsl_getincom use planet_h use comgeomfi_h use comsoil_h implicit none #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "surfdat.h" #include "netcdf.inc" #include "callkeys.h" c----------------------------------------------------------------------- c Declarations c----------------------------------------------------------------------- c Arguments c --------- INTEGER nid,nvarid,tab0 INTEGER*4 day_ini INTEGER Lmodif INTEGER lmax REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time REAL rm !TB c Variables c --------- INTEGER length parameter (length = 100) REAL tab_cntrl(length) ! array in which are stored the run's parameters INTEGER ierr INTEGER size CHARACTER modif*20,choice*1 REAL eccpal,sma ! paleo c----------------------------------------------------------------------- c Initialization of various physical constants to defaut values (nid = 0 case) c----------------------------------------------------------------------- c TB : modif : suppression de la boucle nid c IF (nid.eq.0) then c Reference pressure c------------------------------------- c pressrf = 1. ! Pression de reference (Pa) ~1 c Default (Pluton) parameters for the dynamics and physics c---------------------------------------------------------- rad=1187000 !old= 1173000. ! radius of Pluton (m) daysec=551854.08 ! length of a sol (s) ~551837 s omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1) g=0.617189 ! gravity (m.s-2) GM/r2 old: 0.65 mugaz=28. ! Molar mass of the atmosphere (g.mol-1) ~28 rcp=.2853 ! = r/cp ~0.2853 c Default (Pluton) parameters for physics only c---------------------------------------------- year_day = 14178.30343 ! length of year (sols) ~14182.245 periheli = 4436.82 ! min. Star-Planet distance (Mkm) ~4436 aphelie = 7375.93 ! max. Star-Planet distance (Mkm) ~7375 peri_day = 87.362 ! date of perihelion (sols since N. spring) obliquit=119.591 ! Obliquity of the planet (deg) old ~119.6 tpal=0. ! time (Myr) to start paleoclimate run (e.g -10) adjust=0. ! N2 ice albedo adjustment c Boundary layer and turbulence c---------------------------- z0 = 1.e-2 ! surface roughness (m) ~0.01 c soil properties volcapa = 1.e6 ! soil volumetric heat capacity (in comsoil.h) c----------------------------------------------------------------------- c Initialization of physical constants by reading array tab_cntrl(:) c which contains these parameters (nid != 0 case) c----------------------------------------------------------------------- c Read 'controle' array c ierr = NF_INQ_VARID (nid, "controle", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Tabfi: Could not find data" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) #else ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "Tabfi: Failed reading array" CALL abort ENDIF c Write physical constants to output before modifying them c----------------------------------------------------------------------- 6 FORMAT(a20,e15.6,e15.6) 5 FORMAT(a20,f12.2,f12.2) write(*,*) '*****************************************************' write(*,*) 'Reading tab_cntrl when calling tabfi before changes' write(*,*) '*****************************************************' write(*,*) ' VALUE IN NC DEFAULT VALUE' write(*,5) '(1) = ngridmx?',tab_cntrl(tab0+1),float(ngridmx) write(*,5) '(2) lmax',tab_cntrl(tab0+2),float(lmax) write(*,5) '(3) day_ini',tab_cntrl(tab0+3),float(day_ini) write(*,5) '(5) rad',tab_cntrl(tab0+5),rad write(*,5) '(10) daysec',tab_cntrl(tab0+10),daysec write(*,6) '(6) omeg',tab_cntrl(tab0+6),omeg write(*,5) '(7) g',tab_cntrl(tab0+7),g write(*,5) '(8) mugaz',tab_cntrl(tab0+8),mugaz write(*,5) '(9) rcp',tab_cntrl(tab0+9),rcp write(*,6) '(11) dtphys?',tab_cntrl(tab0+11),dtphys write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day write(*,5) '(15) periheli',tab_cntrl(tab0+15),periheli write(*,5) '(16) aphelie',tab_cntrl(tab0+16),aphelie write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit write(*,6) '(19) z0',tab_cntrl(tab0+19),z0 write(*,5) '(20) tpal',tab_cntrl(tab0+20),tpal write(*,5) '(21) adjust',tab_cntrl(tab0+21),adjust write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa write(*,*) c c Initialization of some physical constants c informations on physics grid if(ngridmx.ne.tab_cntrl(tab0+1)) then print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngridmx' print*,tab_cntrl(tab0+1),ngridmx endif lmax = nint(tab_cntrl(tab0+2)) day_ini = tab_cntrl(tab0+3) time = tab_cntrl(tab0+4) c Informations about planet for dynamics and physics rad = tab_cntrl(tab0+5) omeg = tab_cntrl(tab0+6) g = tab_cntrl(tab0+7) mugaz = tab_cntrl(tab0+8) rcp = tab_cntrl(tab0+9) daysec = tab_cntrl(tab0+10) dtphys = tab_cntrl(tab0+11) c Informations about planet for the physics only year_day = tab_cntrl(tab0+14) periheli = tab_cntrl(tab0+15) aphelie = tab_cntrl(tab0+16) peri_day = tab_cntrl(tab0+17) obliquit = tab_cntrl(tab0+18) c boundary layer and turbeulence z0 = tab_cntrl(tab0+19) c for paleoclimate tpal = tab_cntrl(tab0+20) adjust = tab_cntrl(tab0+21) ! for Triton albedo adjustment c soil properties volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity c----------------------------------------------------------------------- c Save some constants for later use (as routine arguments) c----------------------------------------------------------------------- p_omeg = omeg p_g = g rm= 8.314511E+0 *1000.E+0/mugaz p_cpp = rm/rcp p_mugaz = mugaz p_daysec = daysec p_rad=rad write(*,*),'In particular in the nc file: rad, g = ',p_rad,p_g write(*,*),'cpp calculated from rcp is : ',p_cpp c ENDIF ! end of (nid = 0) c----------------------------------------------------------------------- c Modifications... c----------------------------------------------------------------------- IF(Lmodif.eq.1) then write(*,*) write(*,*) 'Change values in tab_cntrl ? :' write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(*,*) '(Current values given above)' write(*,*) write(*,*) '(2) lmax' write(*,*) '(3) day_ini : Initial day (=0 at Ls=0)' write(*,*) '(5) rad : radius of the planet (m)' write(*,*) '(6) omeg : planet rotation rate (rad/s)' write(*,*) '(7) g : gravity (m/s2)' write(*,*) '(8) mugaz : molecular mass of atm (g/mol)' write(*,*) '(9) rcp : r/Cp' write(*,*) '(10) daysec : length of a sol (s)' write(*,*) '(14) year_day : length of year (in sols)' write(*,*) '(15) periheli : min. sun-pluto dist (Mkm)' write(*,*) '(16) aphelie : max. sun-pluto dist (Mkm)' write(*,*) '(17) peri_day : perihelion date (sol since Ls=0)' write(*,*) '(18) obliquit : planet obliquity (deg)' write(*,*) '(19) z0 : surface roughness (m)' write(*,*) '(20) tpal: paleoclimatic date in million years' write(*,*) '(21) adjust: N2 ice albedo adjustment' write(*,*) '(35) volcapa : soil volumetric heat capacity' write(*,*) do while(modif(1:1).ne.'hello') write(*,*) write(*,*) write(*,*) 'Changes to perform ?' write(*,*) ' (enter keyword or return )' write(*,*) read(*,fmt='(a20)') modif if (modif(1:1) .eq. ' ') goto 999 write(*,*) write(*,*) modif(1:len_trim(modif)) , ' : ' if (modif(1:len_trim(modif)) .eq. 'day_ini') then write(*,*) 'current value:',day_ini write(*,*) 'enter new value:' 101 read(*,*,iostat=ierr) day_ini if(ierr.ne.0) goto 101 write(*,*) ' ' write(*,*) 'day_ini (new value):',day_ini else if (modif(1:len_trim(modif)) .eq. 'z0') then write(*,*) 'current value:',z0 write(*,*) 'enter new value:' 102 read(*,*,iostat=ierr) z0 if(ierr.ne.0) goto 102 write(*,*) ' ' write(*,*) ' z0 (new value):',z0 else if (modif(1:len_trim(modif)) .eq. 'obliquit') then write(*,*) 'current value:',obliquit write(*,*) 'obliquit should be ... on current Pluto' write(*,*) 'enter new value:' 115 read(*,*,iostat=ierr) obliquit if(ierr.ne.0) goto 115 write(*,*) write(*,*) ' obliquit (new value):',obliquit else if (modif(1:len_trim(modif)) .eq. 'peri_day') then write(*,*) 'current value:',peri_day write(*,*) 'peri_day should be 72 on current Pluto' write(*,*) 'enter new value:' 116 read(*,*,iostat=ierr) peri_day if(ierr.ne.0) goto 116 write(*,*) write(*,*) ' peri_day (new value):',peri_day else if (modif(1:len_trim(modif)) .eq. 'periheli') then write(*,*) 'current value:',periheli,' e6 km' write(*,*) 'perihelion should be 4436 on current Pluto' write(*,*) 'enter new value:' 117 read(*,*,iostat=ierr) periheli if(ierr.ne.0) goto 117 write(*,*) write(*,*) ' periheli (new value):',periheli else if (modif(1:len_trim(modif)) .eq. 'aphelie') then write(*,*) 'current value:',aphelie,' e6 km' write(*,*) 'aphelion should be 7375 on current Pluto' write(*,*) 'enter new value:' 118 read(*,*,iostat=ierr) aphelie if(ierr.ne.0) goto 118 write(*,*) write(*,*) ' aphelie (new value):',aphelie else if (modif(1:len_trim(modif)) .eq. 'volcapa') then write(*,*) 'current value:',volcapa write(*,*) 'enter new value:' 119 read(*,*,iostat=ierr) volcapa if(ierr.ne.0) goto 119 write(*,*) write(*,*) ' volcapa (new value):',volcapa else if (modif(1:len_trim(modif)) .eq. 'year_day') then write(*,*) 'current value:',year_day write(*,*) 'enter new value:' 120 read(*,*,iostat=ierr) year_day if(ierr.ne.0) goto 120 write(*,*) write(*,*) ' year_day (new value):',year_day else if (modif(1:len_trim(modif)).eq.'rad') then write(*,*) 'current value:',rad write(*,*) 'enter new value:' 121 read(*,*,iostat=ierr) rad if(ierr.ne.0) goto 121 write(*,*) write(*,*) ' rad (new value):',rad else if (modif(1:len_trim(modif)).eq.'omeg') then write(*,*) 'current value:',omeg write(*,*) 'enter new value:' 122 read(*,*,iostat=ierr) omeg if(ierr.ne.0) goto 122 write(*,*) write(*,*) ' omeg (new value):',omeg else if (modif(1:len_trim(modif)).eq.'g') then write(*,*) 'current value:',g write(*,*) 'enter new value:' 123 read(*,*,iostat=ierr) g if(ierr.ne.0) goto 123 write(*,*) write(*,*) ' g (new value):',g else if (modif(1:len_trim(modif)).eq.'mugaz') then write(*,*) 'current value:',mugaz write(*,*) 'enter new value:' 124 read(*,*,iostat=ierr) mugaz if(ierr.ne.0) goto 124 write(*,*) write(*,*) ' mugaz (new value):',mugaz r=8.314/(mugaz/1000.0) rm=r write(*,*) ' r (new value):',r else if (modif(1:len_trim(modif)).eq.'rcp') then write(*,*) 'current value:',rcp write(*,*) 'enter new value:' 125 read(*,*,iostat=ierr) rcp if(ierr.ne.0) goto 125 write(*,*) write(*,*) ' rcp (new value):',rcp cpp=rm/rcp write(*,*) ' cpp (new value):',cpp else if (modif(1:len_trim(modif)).eq.'daysec') then write(*,*) 'current value:',daysec write(*,*) 'enter new value:' 126 read(*,*,iostat=ierr) daysec if(ierr.ne.0) goto 126 write(*,*) write(*,*) ' daysec (new value):',daysec write(*,*) ' want to update omeg (y/n)?' 128 read(*,*,iostat=ierr) choice if(ierr.ne.0) goto 128 if (choice.eq.'y') then omeg=2.*3.14159265/daysec write(*,*) ' omeg (new value):',omeg endif else if (modif(1:len_trim(modif)).eq.'lmax') then write(*,*) 'current value:',lmax write(*,*) 'enter new value:' 127 read(*,*,iostat=ierr) lmax if(ierr.ne.0) goto 127 write(*,*) write(*,*) ' lmax (new value):',lmax else if (modif(1:len_trim(modif)).eq.'tpal') then write(*,*) 'current value:',tpal write(*,*) 'enter new value:' 129 read(*,*,iostat=ierr) tpal if(ierr.ne.0) goto 129 write(*,*) write(*,*) ' tpal (new value):',tpal write(*,*) 'Modify obli / peri_day / ecc accordingly (y/n)?' 130 read(*,*,iostat=ierr) choice if(ierr.ne.0) goto 130 if (choice.eq.'y') then c Initialising orbital parameters: call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) c Updating obliquity and orbital parameters: call calcmilank(tpal,obliquit,peri_day,eccpal) sma=(periheli+aphelie)/2. periheli=sma*(1-eccpal) aphelie=sma*(1+eccpal) endif else if (modif(1:len_trim(modif)).eq.'adjust') then write(*,*) 'current value:',adjust write(*,*) 'enter new value:' 131 read(*,*,iostat=ierr) adjust if(ierr.ne.0) goto 131 write(*,*) write(*,*) ' adjust (new value):',adjust endif enddo ! of do while(modif(1:1).ne.'hello') 999 continue c----------------------------------------------------------------------- c Write values of physical constants after modifications c----------------------------------------------------------------------- write(*,*) '*****************************************************' write(*,*) 'Reading tab_cntrl when calling tabfi AFTER changes' write(*,*) '*****************************************************' write(*,*) ' OLD VALUE IN NC NEW VALUE IN NC' write(*,5) '(1) = ngridmx?',tab_cntrl(tab0+1),float(ngridmx) write(*,5) '(2) lmax',tab_cntrl(tab0+2),float(lmax) write(*,5) '(3) day_ini',tab_cntrl(tab0+3),float(day_ini) write(*,5) '(5) rad',tab_cntrl(tab0+5),rad write(*,5) '(10) daysec',tab_cntrl(tab0+10),daysec write(*,6) '(6) omeg',tab_cntrl(tab0+6),omeg write(*,5) '(7) g',tab_cntrl(tab0+7),g write(*,5) '(8) mugaz',tab_cntrl(tab0+8),mugaz write(*,5) '(9) rcp',tab_cntrl(tab0+9),rcp write(*,6) '(11) dtphys?',tab_cntrl(tab0+11),dtphys write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day write(*,5) '(15) periheli',tab_cntrl(tab0+15),periheli write(*,5) '(16) aphelie',tab_cntrl(tab0+16),aphelie write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit write(*,6) '(19) z0',tab_cntrl(tab0+19),z0 write(*,5) '(20) tpal',tab_cntrl(tab0+20),tpal write(*,5) '(21) adjust',tab_cntrl(tab0+21),adjust write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa write(*,*) write(*,*) ENDIF ! of if (Lmodif == 1) c----------------------------------------------------------------------- c Save some constants for later use (as routine arguments) c----------------------------------------------------------------------- p_omeg = omeg p_g = g rm= 8.314511E+0 *1000.E+0/mugaz p_cpp = rm/rcp p_mugaz = mugaz p_daysec = daysec p_rad=rad write(*,*),'NOW in the nc file: rad, g, omeg = ',p_rad,p_g,p_omeg write(*,*),'r is = ',rm,' mugaz is ',mugaz write(*,*),'cpp calculated from r and rcp is : ',p_cpp write(*,*),'daysec is : ',daysec end