Index: /trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F =================================================================== --- /trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F (revision 1669) +++ /trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F (revision 1670) @@ -35,5 +35,7 @@ USE temps_mod, ONLY: day_ini USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 + use tabfi_mod, only: tabfi use iniphysiq_mod, only: iniphysiq + use phyetat0_mod, only: phyetat0 implicit none @@ -337,5 +339,5 @@ write(*,*) 'Reading file STARTFI' fichnom = 'startfi.nc' - CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx, + CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx, . nqtot,day_ini,time, . tsurf,tsoil,emis,q2,qsurf) Index: /trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F =================================================================== --- /trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F (revision 1669) +++ /trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F (revision 1670) @@ -35,4 +35,5 @@ USE temps_mod, ONLY: day_ini USE iniphysiq_mod, ONLY: iniphysiq + use phyetat0_mod, only: phyetat0 implicit none @@ -214,5 +215,5 @@ - CALL phyetat0 (ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot, + CALL phyetat0(.true.,ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot, . day_ini_fi,timefi, . tsurf,tsoil,emis,q2,qsurf) Index: /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 1669) +++ /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 1670) @@ -74,3 +74,7 @@ !$OMP THREADPRIVATE(iscallphys) + ! do we read a startphy.nc file (default=.true.) + logical,save :: startphy_file=.true. +!$OMP THREADPRIVATE(startphy_file) + END MODULE callkeys_mod Index: /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 1669) +++ /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 1670) @@ -95,4 +95,7 @@ call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step + ! do we read a startphy.nc file? (default: .true.) + call getin_p("startphy_file",startphy_file) + ! -------------------------------------------------------------- ! Reading the "callphys.def" file controlling some key options Index: unk/LMDZ.TITAN/libf/phytitan/phyetat0.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F90 (revision 1669) +++ (revision ) @@ -1,289 +1,0 @@ -subroutine phyetat0 (ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, & - day_ini,time,tsurf,tsoil, & - emis,q2,qsurf) - - - USE tracer_h, ONLY: noms - USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe - use iostart, only: nid_start, open_startphy, close_startphy, & - get_field, get_var, inquire_field, & - inquire_dimension, inquire_dimension_length - - implicit none - -!====================================================================== -! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818 -! Adaptation à Mars : Yann Wanherdrick -! Objet: Lecture de l etat initial pour la physique -!====================================================================== -#include "netcdf.inc" - -!====================================================================== -! INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4 -! PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille -!====================================================================== -! Arguments: -! --------- -! inputs: - integer,intent(in) :: ngrid - integer,intent(in) :: nlayer - character*(*),intent(in) :: fichnom ! "startfi.nc" file - integer,intent(in) :: tab0 - integer,intent(in) :: Lmodif - integer,intent(in) :: nsoil ! # of soil layers - integer,intent(in) :: nq - integer,intent(in) :: day_ini - real,intent(in) :: time - -! outputs: - real,intent(out) :: tsurf(ngrid) ! surface temperature - real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature - real,intent(out) :: emis(ngrid) ! surface emissivity - real,intent(out) :: q2(ngrid,nlayer+1) ! - real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface - -!====================================================================== -! Local variables: - -! INTEGER radpas -! REAL solaire - - real xmin,xmax ! to display min and max of a field -! - INTEGER ig,iq,lmax - INTEGER nid, nvarid - INTEGER ierr, i, nsrf -! integer isoil -! INTEGER length -! PARAMETER (length=100) - CHARACTER*7 str7 - CHARACTER*2 str2 - CHARACTER*1 yes -! - REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec - INTEGER nqold - -! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...) -! logical :: oldtracernames=.false. - integer :: count - character(len=30) :: txt ! to store some text - - INTEGER :: indextime=1 ! index of selected time, default value=1 - logical :: found - -! -! ALLOCATE ARRAYS IN surfdat_h -! -IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid)) -IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid)) -IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid)) -IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid)) -IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid)) -IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid)) -IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid)) - - -! open physics initial state file: -call open_startphy(fichnom) - - -! possibility to modify tab_cntrl in tabfi -write(*,*) -write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0 -call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, & - p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) - -!c -!c Lecture des latitudes (coordonnees): -!c -! ierr = NF_INQ_VARID (nid, "latitude", nvarid) -! IF (ierr.NE.NF_NOERR) THEN -! PRINT*, 'phyetat0: Le champ est absent' -! CALL abort -! ENDIF -!#ifdef NC_DOUBLE -! ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati) -!#else -! ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati) -!#endif -! IF (ierr.NE.NF_NOERR) THEN -! PRINT*, 'phyetat0: Lecture echouee pour ' -! CALL abort -! ENDIF -!c -!c Lecture des longitudes (coordonnees): -!c -! ierr = NF_INQ_VARID (nid, "longitude", nvarid) -! IF (ierr.NE.NF_NOERR) THEN -! PRINT*, 'phyetat0: Le champ est absent' -! CALL abort -! ENDIF -!#ifdef NC_DOUBLE -! ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long) -!#else -! ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long) -!#endif -! IF (ierr.NE.NF_NOERR) THEN -! PRINT*, 'phyetat0: Lecture echouee pour ' -! CALL abort -! ENDIF -!c -!c Lecture des aires des mailles: -!c -! ierr = NF_INQ_VARID (nid, "area", nvarid) -! IF (ierr.NE.NF_NOERR) THEN -! PRINT*, 'phyetat0: Le champ est absent' -! CALL abort -! ENDIF -!#ifdef NC_DOUBLE -! ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area) -!#else -! ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area) -!#endif -! IF (ierr.NE.NF_NOERR) THEN -! PRINT*, 'phyetat0: Lecture echouee pour ' -! CALL abort -! ENDIF -! xmin = 1.0E+20 -! xmax = -1.0E+20 -! xmin = MINVAL(area) -! xmax = MAXVAL(area) -! PRINT*,'Aires des mailles :', xmin, xmax - -! Load surface geopotential: -call get_field("phisfi",phisfi,found) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: surface geopotential range:", & - minval(phisfi), maxval(phisfi) -endif - -! Load bare ground albedo: -call get_field("albedodat",albedodat,found) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: Bare ground albedo range:", & - minval(albedodat), maxval(albedodat) -endif - -! ZMEA -call get_field("ZMEA",zmea,found) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: range:", & - minval(zmea), maxval(zmea) -endif - -! ZSTD -call get_field("ZSTD",zstd,found) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: range:", & - minval(zstd), maxval(zstd) -endif - -! ZSIG -call get_field("ZSIG",zsig,found) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: range:", & - minval(zsig), maxval(zsig) -endif - -! ZGAM -call get_field("ZGAM",zgam,found) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: range:", & - minval(zgam), maxval(zgam) -endif - -! ZTHE -call get_field("ZTHE",zthe,found) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: range:", & - minval(zthe), maxval(zthe) -endif - -! Surface temperature : -call get_field("tsurf",tsurf,found,indextime) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: Surface temperature range:", & - minval(tsurf), maxval(tsurf) -endif - -! Surface emissivity -call get_field("emis",emis,found,indextime) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: Surface emissivity range:", & - minval(emis), maxval(emis) -endif - -! pbl wind variance -call get_field("q2",q2,found,indextime) -if (.not.found) then - write(*,*) "phyetat0: Failed loading " - call abort -else - write(*,*) "phyetat0: PBL wind variance range:", & - minval(q2), maxval(q2) -endif - -! tracer on surface -if (nq.ge.1) then - do iq=1,nq - txt=noms(iq) - - !! There was a bug here. MT2015. - - !if (txt.eq."h2o_vap") then - ! There is no surface tracer for h2o_vap; - ! "h2o_ice" should be loaded instead - ! txt="h2o_ice" - ! write(*,*) 'phyetat0: loading surface tracer', & - ! ' h2o_ice instead of h2o_vap' - !endif - - call get_field(txt,qsurf(:,iq),found,indextime) - if (.not.found) then - write(*,*) "phyetat0: Failed loading <",trim(txt),">" - write(*,*) " ",trim(txt)," is set to zero" - qsurf(:,iq) = 0. - else - write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", & - minval(qsurf(:,iq)), maxval(qsurf(:,iq)) - endif - enddo -endif ! of if (nq.ge.1) - - -! Call to soil_settings, in order to read soil temperatures, -! as well as thermal inertia and volumetric heat capacity -call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime) -! -! close file: -! -call close_startphy - -END SUBROUTINE phyetat0 Index: /trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90 (revision 1670) +++ /trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90 (revision 1670) @@ -0,0 +1,253 @@ +module phyetat0_mod + +implicit none + +contains + +subroutine phyetat0 (startphy_file, & + ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, & + day_ini,time,tsurf,tsoil, & + emis,q2,qsurf) + + + use tabfi_mod, only: tabfi + USE tracer_h, ONLY: noms + USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe + use iostart, only: nid_start, open_startphy, close_startphy, & + get_field, get_var, inquire_field, & + inquire_dimension, inquire_dimension_length + + implicit none + +!====================================================================== +! Arguments: +! --------- +! inputs: + logical,intent(in) :: startphy_file ! .true. if reading start file + integer,intent(in) :: ngrid + integer,intent(in) :: nlayer + character*(*),intent(in) :: fichnom ! "startfi.nc" file + integer,intent(in) :: tab0 + integer,intent(in) :: Lmodif + integer,intent(in) :: nsoil ! # of soil layers + integer,intent(in) :: nq + integer,intent(out) :: day_ini + real,intent(out) :: time + +! outputs: + real,intent(out) :: tsurf(ngrid) ! surface temperature + real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature + real,intent(out) :: emis(ngrid) ! surface emissivity + real,intent(out) :: q2(ngrid,nlayer+1) ! + real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface + +!====================================================================== +! Local variables: + +! INTEGER radpas +! REAL solaire + + real xmin,xmax ! to display min and max of a field +! + INTEGER ig,iq,lmax + INTEGER nid, nvarid + INTEGER ierr, i, nsrf +! integer isoil +! INTEGER length +! PARAMETER (length=100) + CHARACTER*7 str7 + CHARACTER*2 str2 + CHARACTER*1 yes +! + REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec + INTEGER nqold + +! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...) +! logical :: oldtracernames=.false. + integer :: count + character(len=30) :: txt ! to store some text + + INTEGER :: indextime=1 ! index of selected time, default value=1 + logical :: found + + character(len=8) :: modname="phyetat0" + +! +! ALLOCATE ARRAYS IN surfdat_h +! +IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid)) +IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid)) +IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid)) +IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid)) +IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid)) +IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid)) +IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid)) + +if (startphy_file) then + ! open physics initial state file: + call open_startphy(fichnom) + + ! possibility to modify tab_cntrl in tabfi + write(*,*) + write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0 + call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, & + p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) + +else ! "academic" initialization of planetary parameters via tabfi + call tabfi (ngrid,0,0,0,day_ini,lmax,p_rad, & + p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) +endif ! of if (startphy_file) + +if (startphy_file) then + ! Load surface geopotential: + call get_field("phisfi",phisfi,found) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + phisfi(:)=0 +endif ! of if (startphy_file) +write(*,*) "phyetat0: surface geopotential range:", & + minval(phisfi), maxval(phisfi) + +if (startphy_file) then + ! Load bare ground albedo: + call get_field("albedodat",albedodat,found) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + albedodat(:)=0.5 ! would be better to read value from def file... +endif ! of if (startphy_file) +write(*,*) "phyetat0: Bare ground albedo range:", & + minval(albedodat), maxval(albedodat) + +! ZMEA +if (startphy_file) then + call get_field("ZMEA",zmea,found) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + zmea(:)=0 +endif ! of if (startphy_file) +write(*,*) "phyetat0: range:", & + minval(zmea), maxval(zmea) + +! ZSTD +if (startphy_file) then + call get_field("ZSTD",zstd,found) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + zstd(:)=0 +endif ! of if (startphy_file) +write(*,*) "phyetat0: range:", & + minval(zstd), maxval(zstd) + +! ZSIG +if (startphy_file) then + call get_field("ZSIG",zsig,found) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + zsig(:)=0 +endif ! of if (startphy_file) +write(*,*) "phyetat0: range:", & + minval(zsig), maxval(zsig) + +! ZGAM +if (startphy_file) then + call get_field("ZGAM",zgam,found) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + zgam(:)=0 +endif ! of if (startphy_file) +write(*,*) "phyetat0: range:", & + minval(zgam), maxval(zgam) + +! ZTHE +if (startphy_file) then + call get_field("ZTHE",zthe,found) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + zthe(:)=0 +endif ! of if (startphy_file) +write(*,*) "phyetat0: range:", & + minval(zthe), maxval(zthe) + +! Surface temperature : +if (startphy_file) then + call get_field("tsurf",tsurf,found,indextime) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + tsurf(:)=0 ! will be updated afterwards in physiq ! +endif ! of if (startphy_file) +write(*,*) "phyetat0: Surface temperature range:", & + minval(tsurf), maxval(tsurf) + +! Surface emissivity +if (startphy_file) then + call get_field("emis",emis,found,indextime) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + emis(:)=1 ! would be better to read value from def file... +endif ! of if (startphy_file) +write(*,*) "phyetat0: Surface emissivity range:", & + minval(emis), maxval(emis) + +! pbl wind variance +if (startphy_file) then + call get_field("q2",q2,found,indextime) + if (.not.found) then + call abort_physic(modname,"Failed loading ",1) + endif +else + q2(:,:)=0 +endif ! of if (startphy_file) +write(*,*) "phyetat0: PBL wind variance range:", & + minval(q2), maxval(q2) + +! tracer on surface +if (nq.ge.1) then + do iq=1,nq + txt=noms(iq) + if (startphy_file) then + call get_field(txt,qsurf(:,iq),found,indextime) + if (.not.found) then + write(*,*) "phyetat0: Failed loading <",trim(txt),">" + write(*,*) " ",trim(txt)," is set to zero" + qsurf(:,iq) = 0. + endif + else + qsurf(:,iq)=0 + endif ! of if (startphy_file) + write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", & + minval(qsurf(:,iq)), maxval(qsurf(:,iq)) + enddo! of do iq=1,nq +endif ! of if (nq.ge.1) + + +if (startphy_file) then + ! Call to soil_settings, in order to read soil temperatures, + ! as well as thermal inertia and volumetric heat capacity + call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime) +endif ! of if (startphy_file) +! +! close file: +! +if (startphy_file) call close_startphy + +end subroutine phyetat0 + +end module phyetat0_mod Index: /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 1669) +++ /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 1670) @@ -25,4 +25,5 @@ alpha_lift, alpha_devil, qextrhor use time_phylmdz_mod, only: ecritphy, iphysiq, nday + use phyetat0_mod, only: phyetat0 use phyredem, only: physdem0, physdem1 use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval @@ -146,5 +147,5 @@ ! ------------------ -#include "netcdf.inc" +include "netcdf.inc" ! Arguments : @@ -225,5 +226,5 @@ real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!) - integer l,ig,ierr,iq,nw + integer l,ig,ierr,iq,nw,isoil ! FOR DIAGNOSTIC : @@ -432,6 +433,18 @@ ! Read 'startfi.nc' file. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, & + call phyetat0(startphy_file, & + ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, & day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf) + if (.not.startphy_file) then + ! additionnal "academic" initialization of physics + if (is_master) write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!" + tsurf(:)=pt(:,1) + if (is_master) write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!" + do isoil=1,nsoilmx + tsoil(1:ngrid,isoil)=tsurf(1:ngrid) + enddo + if (is_master) write(*,*) "Physiq: initializing day_ini to pdat !" + day_ini=pday + endif if (pday.ne.day_ini) then Index: unk/LMDZ.TITAN/libf/phytitan/tabfi.F =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/tabfi.F (revision 1669) +++ (revision ) @@ -1,522 +1,0 @@ -c======================================================================= - SUBROUTINE tabfi(ngrid,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 nlayer) -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 , only: getin - - use surfdat_h, only: emisice, iceradius, dtemisice, - & emissiv - use comsoil_h, only: volcapa - use iostart, only: get_var - use mod_phys_lmdz_para, only: is_parallel - use planete_mod, only: year_day, periastr, apoastr, peri_day, - & obliquit, z0, lmixmin, emin_turb - use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r - use time_phylmdz_mod, only: dtphys, daysec - use callkeys_mod, only: check_cpp_match,force_cpp - implicit none - -#include "netcdf.inc" - -c----------------------------------------------------------------------- -c Declarations -c----------------------------------------------------------------------- - -c Arguments -c --------- - INTEGER,INTENT(IN) :: ngrid,nid,tab0 - INTEGER*4,INTENT(OUT) :: day_ini - INTEGER,INTENT(IN) :: Lmodif - INTEGER,INTENT(OUT) :: lmax - REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time - -c Variables -c --------- - INTEGER,PARAMETER :: length=100 - REAL tab_cntrl(length) ! array in which are stored the run's parameters - INTEGER ierr,nvarid - INTEGER size - CHARACTER modif*20 - LOGICAL :: found - - write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif - - IF (nid.eq.0) then -c----------------------------------------------------------------------- -c Initialization of various physical constants to defaut values (nid = 0 case) -c----------------------------------------------------------------------- - ELSE -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 - - call get_var("controle",tab_cntrl,found) - if (.not.found) then - write(*,*)"tabfi: Failed reading array" - call abort - else - write(*,*)'tabfi: tab_cntrl',tab_cntrl - endif -c -c Initialization of some physical constants -c informations on physics grid -! if(ngrid.ne.tab_cntrl(tab0+1)) then -! print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid' -! print*,tab_cntrl(tab0+1),ngrid -! endif - lmax = nint(tab_cntrl(tab0+2)) - day_ini = tab_cntrl(tab0+3) - time = tab_cntrl(tab0+4) - write (*,*) 'IN tabfi day_ini=',day_ini -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) - cpp=(8.314511/(mugaz/1000.0))/rcp - daysec = tab_cntrl(tab0+10) - dtphys = tab_cntrl(tab0+11) -c Informations about planet for the physics only - year_day = tab_cntrl(tab0+14) - periastr = tab_cntrl(tab0+15) - apoastr = 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) - lmixmin = tab_cntrl(tab0+20) - emin_turb = tab_cntrl(tab0+21) -c optical properties of polar caps and ground emissivity - emisice(1) = tab_cntrl(tab0+24) - emisice(2) = tab_cntrl(tab0+25) - emissiv = tab_cntrl(tab0+26) - iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north) - iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south) - dtemisice(1)= tab_cntrl(tab0+33) !time scale for snow metamorphism (north) - dtemisice(2)= tab_cntrl(tab0+34) !time scale for snow metamorphism (south) -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 - p_cpp = cpp - p_mugaz = mugaz - p_daysec = daysec - p_rad=rad - - ENDIF ! end of (nid = 0) - -c----------------------------------------------------------------------- -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(*,5) '(1) = ngrid?',tab_cntrl(tab0+1),float(ngrid) - 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) periastr',tab_cntrl(tab0+15),periastr - write(*,5) '(16) apoastr',tab_cntrl(tab0+16),apoastr - 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(*,6) '(21) emin_turb',tab_cntrl(tab0+21),emin_turb - write(*,5) '(20) lmixmin',tab_cntrl(tab0+20),lmixmin - - write(*,5) '(26) emissiv',tab_cntrl(tab0+26),emissiv - write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1) - write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2) - write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1) - write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2) - write(*,5) '(33) dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1) - write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2) - - write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa - - write(*,*) - write(*,*) 'Lmodif in tabfi!!!!!!!',Lmodif - -c----------------------------------------------------------------------- -c Modifications... -! NB: Modifying controls should only be done by newstart, and in seq mode - if ((Lmodif.eq.1).and.is_parallel) then - write(*,*) "tabfi: Error modifying tab_control should", - & " only happen in serial mode (eg: by newstart)" - stop - endif -c----------------------------------------------------------------------- - - IF(Lmodif.eq.1) then - - write(*,*) - write(*,*) 'Change values in tab_cntrl ? :' - write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' - write(*,*) '(Current values given above)' - write(*,*) - write(*,*) '(3) day_ini : Initial day (=0 at Ls=0)' - write(*,*) '(19) z0 : surface roughness (m)' - write(*,*) '(21) emin_turb : minimal energy (PBL)' - write(*,*) '(20) lmixmin : mixing length (PBL)' - write(*,*) '(26) emissiv : ground emissivity' - write(*,*) '(24 et 25) emisice : CO2 ice max emissivity ' - write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow' - write(*,*) '(33 et 34) dtemisice : time scale for snow', - & 'metamorphism' - write(*,*) '(35) volcapa : soil volumetric heat capacity' - write(*,*) '(18) obliquit : planet obliquity (deg)' - write(*,*) '(17) peri_day : periastron date (sols since Ls=0)' - write(*,*) '(15) periastr : min. star-planet dist (UA)' - write(*,*) '(16) apoastr : max. star-planet (UA)' - write(*,*) '(14) year_day : length of year (in sols)' - 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 ' - write(*,*) ' of the atmosphere (g/mol)' - write(*,*) '(9) rcp : r/Cp' - write(*,*) '(8)+(9) calc_cpp_mugaz : r/Cp and mugaz ' - write(*,*) ' computed from gases.def' - write(*,*) '(10) daysec : length of a sol (s)' - 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. 'emin_turb') then - write(*,*) 'current value:',emin_turb - write(*,*) 'enter new value:' - 103 read(*,*,iostat=ierr) emin_turb - if(ierr.ne.0) goto 103 - write(*,*) ' ' - write(*,*) ' emin_turb (new value):',emin_turb - - else if (modif(1:len_trim(modif)) .eq. 'lmixmin') then - write(*,*) 'current value:',lmixmin - write(*,*) 'enter new value:' - 104 read(*,*,iostat=ierr) lmixmin - if(ierr.ne.0) goto 104 - write(*,*) ' ' - write(*,*) ' lmixmin (new value):',lmixmin - - else if (modif(1:len_trim(modif)) .eq. 'emissiv') then - write(*,*) 'current value:',emissiv - write(*,*) 'enter new value:' - 105 read(*,*,iostat=ierr) emissiv - if(ierr.ne.0) goto 105 - write(*,*) ' ' - write(*,*) ' emissiv (new value):',emissiv - - else if (modif(1:len_trim(modif)) .eq. 'emisice') then - write(*,*) 'current value emisice(1) North:',emisice(1) - write(*,*) 'enter new value:' - 106 read(*,*,iostat=ierr) emisice(1) - if(ierr.ne.0) goto 106 - write(*,*) - write(*,*) ' emisice(1) (new value):',emisice(1) - write(*,*) - - write(*,*) 'current value emisice(2) South:',emisice(2) - write(*,*) 'enter new value:' - 107 read(*,*,iostat=ierr) emisice(2) - if(ierr.ne.0) goto 107 - write(*,*) - write(*,*) ' emisice(2) (new value):',emisice(2) - - else if (modif(1:len_trim(modif)) .eq. 'iceradius') then - write(*,*) 'current value iceradius(1) North:',iceradius(1) - write(*,*) 'enter new value:' - 110 read(*,*,iostat=ierr) iceradius(1) - if(ierr.ne.0) goto 110 - write(*,*) - write(*,*) ' iceradius(1) (new value):',iceradius(1) - write(*,*) - - write(*,*) 'current value iceradius(2) South:',iceradius(2) - write(*,*) 'enter new value:' - 111 read(*,*,iostat=ierr) iceradius(2) - if(ierr.ne.0) goto 111 - write(*,*) - write(*,*) ' iceradius(2) (new value):',iceradius(2) - - else if (modif(1:len_trim(modif)) .eq. 'dtemisice') then - write(*,*) 'current value dtemisice(1) North:',dtemisice(1) - write(*,*) 'enter new value:' - 112 read(*,*,iostat=ierr) dtemisice(1) - if(ierr.ne.0) goto 112 - write(*,*) - write(*,*) ' dtemisice(1) (new value):',dtemisice(1) - write(*,*) - - write(*,*) 'current value dtemisice(2) South:',dtemisice(2) - write(*,*) 'enter new value:' - 113 read(*,*,iostat=ierr) dtemisice(2) - if(ierr.ne.0) goto 113 - write(*,*) - write(*,*) ' dtemisice(2) (new value):',dtemisice(2) - - else if (modif(1:len_trim(modif)) .eq. 'obliquit') then - write(*,*) 'current value:',obliquit - write(*,*) 'obliquit should be 25.19 on current Mars' - 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 485 on current Mars' - 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. 'periastr') then - write(*,*) 'current value:',periastr - write(*,*) 'periastr should be 206.66 on present-day Mars' - write(*,*) 'enter new value:' - 117 read(*,*,iostat=ierr) periastr - if(ierr.ne.0) goto 117 - write(*,*) - write(*,*) ' periastr (new value):',periastr - - else if (modif(1:len_trim(modif)) .eq. 'apoastr') then - write(*,*) 'current value:',apoastr - write(*,*) 'apoastr should be 249.22 on present-day Mars' - write(*,*) 'enter new value:' - 118 read(*,*,iostat=ierr) apoastr - if(ierr.ne.0) goto 118 - write(*,*) - write(*,*) ' apoastr (new value):',apoastr - - 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.'rad') then - write(*,*) 'current value:',rad - write(*,*) 'enter new value:' - 120 read(*,*,iostat=ierr) rad - if(ierr.ne.0) goto 120 - write(*,*) - write(*,*) ' rad (new value):',rad - - else if (modif(1:len_trim(modif)).eq.'omeg') then - write(*,*) 'current value:',omeg - write(*,*) 'enter new value:' - 121 read(*,*,iostat=ierr) omeg - if(ierr.ne.0) goto 121 - write(*,*) - write(*,*) ' omeg (new value):',omeg - - else if (modif(1:len_trim(modif)).eq.'g') then - write(*,*) 'current value:',g - write(*,*) 'enter new value:' - 122 read(*,*,iostat=ierr) g - if(ierr.ne.0) goto 122 - write(*,*) - write(*,*) ' g (new value):',g - - else if (modif(1:len_trim(modif)).eq.'mugaz') then - write(*,*) 'current value:',mugaz - write(*,*) 'enter new value:' - 123 read(*,*,iostat=ierr) mugaz - if(ierr.ne.0) goto 123 - write(*,*) - write(*,*) ' mugaz (new value):',mugaz - r=8.314511/(mugaz/1000.0) - write(*,*) ' R (new value):',r - - else if (modif(1:len_trim(modif)).eq.'rcp') then - write(*,*) 'current value:',rcp - write(*,*) 'enter new value:' - 124 read(*,*,iostat=ierr) rcp - if(ierr.ne.0) goto 124 - write(*,*) - write(*,*) ' rcp (new value):',rcp - r=8.314511/(mugaz/1000.0) - cpp=r/rcp - write(*,*) ' cpp (new value):',cpp - - else if (modif(1:len_trim(modif)).eq.'calc_cpp_mugaz') then - write(*,*) 'current value rcp, mugaz:',rcp,mugaz - check_cpp_match=.false. - force_cpp=.false. - call su_gases - call calc_cpp_mugaz - write(*,*) - write(*,*) ' cpp (new value):',cpp - write(*,*) ' mugaz (new value):',mugaz - r=8.314511/(mugaz/1000.0) - rcp=r/cpp - write(*,*) ' rcp (new value):',rcp - - else if (modif(1:len_trim(modif)).eq.'daysec') then - write(*,*) 'current value:',daysec - write(*,*) 'enter new value:' - 125 read(*,*,iostat=ierr) daysec - if(ierr.ne.0) goto 125 - write(*,*) - write(*,*) ' daysec (new value):',daysec - -! added by RW! - else if (modif(1:len_trim(modif)).eq.'year_day') then - write(*,*) 'current value:',year_day - write(*,*) 'enter new value:' - 126 read(*,*,iostat=ierr) year_day - if(ierr.ne.0) goto 126 - write(*,*) - write(*,*) ' year_day (new value):',year_day - - 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(*,5) '(1) = ngrid?',tab_cntrl(tab0+1),float(ngrid) - 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) periastr',tab_cntrl(tab0+15),periastr - write(*,5) '(16) apoastr',tab_cntrl(tab0+16),apoastr - 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(*,6) '(21) emin_turb',tab_cntrl(tab0+21),emin_turb - write(*,5) '(20) lmixmin',tab_cntrl(tab0+20),lmixmin - - write(*,5) '(26) emissiv',tab_cntrl(tab0+26),emissiv - write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1) - write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2) - write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1) - write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2) - write(*,5) '(33) dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1) - write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2) - - 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 - p_cpp = cpp - p_mugaz = mugaz - p_daysec = daysec - p_rad=rad - - - end Index: /trunk/LMDZ.TITAN/libf/phytitan/tabfi_mod.F90 =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/tabfi_mod.F90 (revision 1670) +++ /trunk/LMDZ.TITAN/libf/phytitan/tabfi_mod.F90 (revision 1670) @@ -0,0 +1,585 @@ +MODULE tabfi_mod + +IMPLICIT NONE + +CONTAINS + +!======================================================================= + SUBROUTINE tabfi(ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad, & + p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) +!======================================================================= +! +! C. Hourdin 15/11/96 +! +! Object: Lecture du tab_cntrl physique dans un fichier +! ------ et initialisation des constantes physiques +! +! Arguments: +! ---------- +! +! Inputs: +! ------ +! +! - nid: unitne logique du fichier ou on va lire le tab_cntrl +! (ouvert dans le programme appellant) +! +! si nid=0: +! pas de lecture du tab_cntrl mais +! Valeurs par default des constantes physiques +! +! - tab0: Offset de tab_cntrl a partir duquel sont ranges +! les parametres physiques (50 pour start_archive) +! +! - Lmodif: si on souhaite modifier les constantes Lmodif = 1 = TRUE +! +! +! Outputs: +! -------- +! +! - day_ini: tab_cntrl(tab0+3) (Dans les cas ou l'on souhaite +! comparer avec le day_ini dynamique) +! +! - lmax: tab_cntrl(tab0+2) (pour test avec nlayer) +! +! - p_rad +! - p_omeg ! +! - p_g ! Constantes physiques ayant des +! - p_mugaz ! homonymes dynamiques +! - p_daysec ! +! +!======================================================================= +! to use 'getin_p' + use ioipsl_getin_p_mod, only: getin_p + + use surfdat_h, only: emisice, iceradius, dtemisice, & + emissiv + use comsoil_h, only: volcapa + use iostart, only: get_var + use mod_phys_lmdz_para, only: is_parallel + use planete_mod, only: year_day, periastr, apoastr, peri_day, & + obliquit, z0, lmixmin, emin_turb + use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r + use time_phylmdz_mod, only: dtphys, daysec + use callkeys_mod, only: check_cpp_match,force_cpp + implicit none + + include "netcdf.inc" + +!----------------------------------------------------------------------- +! Declarations +!----------------------------------------------------------------------- + +! Arguments +! --------- + INTEGER,INTENT(IN) :: ngrid,nid,tab0 + INTEGER*4,INTENT(OUT) :: day_ini + INTEGER,INTENT(IN) :: Lmodif + INTEGER,INTENT(OUT) :: lmax + REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time + +! Variables +! --------- + INTEGER,PARAMETER :: length=100 + REAL tab_cntrl(length) ! array in which are stored the run's parameters + INTEGER ierr,nvarid + INTEGER size + CHARACTER modif*20 + LOGICAL :: found + CHARACTER(len=5) :: modname="tabfi" + + write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif + + IF (nid.eq.0) then +!----------------------------------------------------------------------- +! Initialization of various physical constants to defaut values (nid = 0 case) +!----------------------------------------------------------------------- + tab_cntrl(:)=0 + lmax=0 ! not used anyways + !day_ini already set via inifis + time=0 +! Informations about planet for dynamics and physics + ! rad,cpp,g,r,rcp already initialized by inifis + omeg=-999. + call getin_p("omega",omeg) + if (omeg.eq.-999.) then + call abort_physic(modname,"Missing value for omega in def files!",1) + endif + mugaz=(8.3144621/r)*1.E3 + ! daysec already set by inifis + ! dtphys alread set by inifis +! Informations about planet for the physics only + year_day=-999. ! length of year, in standard days + call getin_p("year_day",year_day) + if (year_day.eq.-999.) then + call abort_physic(modname, & + "Missing value for year_day in def files!",1) + endif + periastr=-999. + call getin_p("periastron",periastr) + if (periastr.eq.-999.) then + call abort_physic(modname, & + "Missing value for periastron in def files!",1) + endif + apoastr=-999. + call getin_p("apoastron",apoastr) + if (apoastr.eq.-999.) then + call abort_physic(modname, & + "Missing value for apoastron in def files!",1) + endif + peri_day=-999. + call getin_p("periastron_day",peri_day) + if (peri_day.eq.-999.) then + call abort_physic(modname, & + "Missing value for periastron date in def files!",1) + endif + obliquit=-999. + call getin_p("obliquity",obliquit) + if (obliquit.eq.-999.) then + call abort_physic(modname, & + "Missing value for obliquity in def files!",1) + endif +! boundary layer and turbulence + z0=1.e-2 ! surface roughness length (m) + lmixmin=30 + emin_turb=1.e-6 +! optical properties of polar caps and ground emissivity + emisice(:)=0 + emissiv=0 + iceradius(:)=1.e-6 ! mean scat radius of CO2 snow + dtemisice(:)=0 !time scale for snow metamorphism + volcapa=1000000 ! volumetric heat capacity of subsurface + + ELSE +!----------------------------------------------------------------------- +! Initialization of physical constants by reading array tab_cntrl(:) +! which contains these parameters (nid != 0 case) +!----------------------------------------------------------------------- +! Read 'controle' array +! + + call get_var("controle",tab_cntrl,found) + if (.not.found) then + call abort_physic(modname,"Failed reading array",1) + else + write(*,*)'tabfi: tab_cntrl',tab_cntrl + endif +! +! Initialization of some physical constants +! informations on physics grid +! if(ngrid.ne.tab_cntrl(tab0+1)) then +! print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid' +! print*,tab_cntrl(tab0+1),ngrid +! endif + lmax = nint(tab_cntrl(tab0+2)) + day_ini = tab_cntrl(tab0+3) + time = tab_cntrl(tab0+4) + write (*,*) 'IN tabfi day_ini=',day_ini +! 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) + cpp=(8.314511/(mugaz/1000.0))/rcp + daysec = tab_cntrl(tab0+10) + dtphys = tab_cntrl(tab0+11) +! Informations about planet for the physics only + year_day = tab_cntrl(tab0+14) + periastr = tab_cntrl(tab0+15) + apoastr = tab_cntrl(tab0+16) + peri_day = tab_cntrl(tab0+17) + obliquit = tab_cntrl(tab0+18) +! boundary layer and turbulence + z0 = tab_cntrl(tab0+19) + lmixmin = tab_cntrl(tab0+20) + emin_turb = tab_cntrl(tab0+21) +! optical properties of polar caps and ground emissivity + emisice(1) = tab_cntrl(tab0+24) + emisice(2) = tab_cntrl(tab0+25) + emissiv = tab_cntrl(tab0+26) + iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north) + iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south) + dtemisice(1)= tab_cntrl(tab0+33) !time scale for snow metamorphism (north) + dtemisice(2)= tab_cntrl(tab0+34) !time scale for snow metamorphism (south) +! soil properties + volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity +!----------------------------------------------------------------------- +! Save some constants for later use (as routine arguments) +!----------------------------------------------------------------------- + p_omeg = omeg + p_g = g + p_cpp = cpp + p_mugaz = mugaz + p_daysec = daysec + p_rad=rad + + ENDIF ! end of (nid = 0) + +!----------------------------------------------------------------------- +! Write physical constants to output before modifying them +!----------------------------------------------------------------------- + + 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(*,5) '(1) = ngrid?',tab_cntrl(tab0+1),float(ngrid) + 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) periastr',tab_cntrl(tab0+15),periastr + write(*,5) '(16) apoastr',tab_cntrl(tab0+16),apoastr + 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(*,6) '(21) emin_turb',tab_cntrl(tab0+21),emin_turb + write(*,5) '(20) lmixmin',tab_cntrl(tab0+20),lmixmin + + write(*,5) '(26) emissiv',tab_cntrl(tab0+26),emissiv + write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1) + write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2) + write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1) + write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2) + write(*,5) '(33) dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1) + write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2) + + write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa + + write(*,*) + write(*,*) 'Lmodif in tabfi!!!!!!!',Lmodif + +!----------------------------------------------------------------------- +! Modifications... +! NB: Modifying controls should only be done by newstart, and in seq mode + if ((Lmodif.eq.1).and.is_parallel) then + write(*,*) "tabfi: Error modifying tab_control should", & + " only happen in serial mode (eg: by newstart)" + stop + endif +!----------------------------------------------------------------------- + + IF(Lmodif.eq.1) then + + write(*,*) + write(*,*) 'Change values in tab_cntrl ? :' + write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + write(*,*) '(Current values given above)' + write(*,*) + write(*,*) '(3) day_ini : Initial day (=0 at Ls=0)' + write(*,*) '(19) z0 : surface roughness (m)' + write(*,*) '(21) emin_turb : minimal energy (PBL)' + write(*,*) '(20) lmixmin : mixing length (PBL)' + write(*,*) '(26) emissiv : ground emissivity' + write(*,*) '(24 et 25) emisice : CO2 ice max emissivity ' + write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow' + write(*,*) '(33 et 34) dtemisice : time scale for snow metamorphism' + write(*,*) '(35) volcapa : soil volumetric heat capacity' + write(*,*) '(18) obliquit : planet obliquity (deg)' + write(*,*) '(17) peri_day : periastron date (sols since Ls=0)' + write(*,*) '(15) periastr : min. star-planet dist (UA)' + write(*,*) '(16) apoastr : max. star-planet (UA)' + write(*,*) '(14) year_day : length of year (in sols)' + 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 ' + write(*,*) ' of the atmosphere (g/mol)' + write(*,*) '(9) rcp : r/Cp' + write(*,*) '(8)+(9) calc_cpp_mugaz : r/Cp and mugaz ' + write(*,*) ' computed from gases.def' + write(*,*) '(10) daysec : length of a sol (s)' + 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. 'emin_turb') then + write(*,*) 'current value:',emin_turb + write(*,*) 'enter new value:' + 103 read(*,*,iostat=ierr) emin_turb + if(ierr.ne.0) goto 103 + write(*,*) ' ' + write(*,*) ' emin_turb (new value):',emin_turb + + else if (modif(1:len_trim(modif)) .eq. 'lmixmin') then + write(*,*) 'current value:',lmixmin + write(*,*) 'enter new value:' + 104 read(*,*,iostat=ierr) lmixmin + if(ierr.ne.0) goto 104 + write(*,*) ' ' + write(*,*) ' lmixmin (new value):',lmixmin + + else if (modif(1:len_trim(modif)) .eq. 'emissiv') then + write(*,*) 'current value:',emissiv + write(*,*) 'enter new value:' + 105 read(*,*,iostat=ierr) emissiv + if(ierr.ne.0) goto 105 + write(*,*) ' ' + write(*,*) ' emissiv (new value):',emissiv + + else if (modif(1:len_trim(modif)) .eq. 'emisice') then + write(*,*) 'current value emisice(1) North:',emisice(1) + write(*,*) 'enter new value:' + 106 read(*,*,iostat=ierr) emisice(1) + if(ierr.ne.0) goto 106 + write(*,*) + write(*,*) ' emisice(1) (new value):',emisice(1) + write(*,*) + + write(*,*) 'current value emisice(2) South:',emisice(2) + write(*,*) 'enter new value:' + 107 read(*,*,iostat=ierr) emisice(2) + if(ierr.ne.0) goto 107 + write(*,*) + write(*,*) ' emisice(2) (new value):',emisice(2) + + else if (modif(1:len_trim(modif)) .eq. 'iceradius') then + write(*,*) 'current value iceradius(1) North:',iceradius(1) + write(*,*) 'enter new value:' + 110 read(*,*,iostat=ierr) iceradius(1) + if(ierr.ne.0) goto 110 + write(*,*) + write(*,*) ' iceradius(1) (new value):',iceradius(1) + write(*,*) + + write(*,*) 'current value iceradius(2) South:',iceradius(2) + write(*,*) 'enter new value:' + 111 read(*,*,iostat=ierr) iceradius(2) + if(ierr.ne.0) goto 111 + write(*,*) + write(*,*) ' iceradius(2) (new value):',iceradius(2) + + else if (modif(1:len_trim(modif)) .eq. 'dtemisice') then + write(*,*) 'current value dtemisice(1) North:',dtemisice(1) + write(*,*) 'enter new value:' + 112 read(*,*,iostat=ierr) dtemisice(1) + if(ierr.ne.0) goto 112 + write(*,*) + write(*,*) ' dtemisice(1) (new value):',dtemisice(1) + write(*,*) + + write(*,*) 'current value dtemisice(2) South:',dtemisice(2) + write(*,*) 'enter new value:' + 113 read(*,*,iostat=ierr) dtemisice(2) + if(ierr.ne.0) goto 113 + write(*,*) + write(*,*) ' dtemisice(2) (new value):',dtemisice(2) + + else if (modif(1:len_trim(modif)) .eq. 'obliquit') then + write(*,*) 'current value:',obliquit + write(*,*) 'obliquit should be 25.19 on current Mars' + 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 485 on current Mars' + 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. 'periastr') then + write(*,*) 'current value:',periastr + write(*,*) 'periastr should be 206.66 on present-day Mars' + write(*,*) 'enter new value:' + 117 read(*,*,iostat=ierr) periastr + if(ierr.ne.0) goto 117 + write(*,*) + write(*,*) ' periastr (new value):',periastr + + else if (modif(1:len_trim(modif)) .eq. 'apoastr') then + write(*,*) 'current value:',apoastr + write(*,*) 'apoastr should be 249.22 on present-day Mars' + write(*,*) 'enter new value:' + 118 read(*,*,iostat=ierr) apoastr + if(ierr.ne.0) goto 118 + write(*,*) + write(*,*) ' apoastr (new value):',apoastr + + 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.'rad') then + write(*,*) 'current value:',rad + write(*,*) 'enter new value:' + 120 read(*,*,iostat=ierr) rad + if(ierr.ne.0) goto 120 + write(*,*) + write(*,*) ' rad (new value):',rad + + else if (modif(1:len_trim(modif)).eq.'omeg') then + write(*,*) 'current value:',omeg + write(*,*) 'enter new value:' + 121 read(*,*,iostat=ierr) omeg + if(ierr.ne.0) goto 121 + write(*,*) + write(*,*) ' omeg (new value):',omeg + + else if (modif(1:len_trim(modif)).eq.'g') then + write(*,*) 'current value:',g + write(*,*) 'enter new value:' + 122 read(*,*,iostat=ierr) g + if(ierr.ne.0) goto 122 + write(*,*) + write(*,*) ' g (new value):',g + + else if (modif(1:len_trim(modif)).eq.'mugaz') then + write(*,*) 'current value:',mugaz + write(*,*) 'enter new value:' + 123 read(*,*,iostat=ierr) mugaz + if(ierr.ne.0) goto 123 + write(*,*) + write(*,*) ' mugaz (new value):',mugaz + r=8.314511/(mugaz/1000.0) + write(*,*) ' R (new value):',r + + else if (modif(1:len_trim(modif)).eq.'rcp') then + write(*,*) 'current value:',rcp + write(*,*) 'enter new value:' + 124 read(*,*,iostat=ierr) rcp + if(ierr.ne.0) goto 124 + write(*,*) + write(*,*) ' rcp (new value):',rcp + r=8.314511/(mugaz/1000.0) + cpp=r/rcp + write(*,*) ' cpp (new value):',cpp + + else if (modif(1:len_trim(modif)).eq.'calc_cpp_mugaz') then + write(*,*) 'current value rcp, mugaz:',rcp,mugaz + check_cpp_match=.false. + force_cpp=.false. + call su_gases + call calc_cpp_mugaz + write(*,*) + write(*,*) ' cpp (new value):',cpp + write(*,*) ' mugaz (new value):',mugaz + r=8.314511/(mugaz/1000.0) + rcp=r/cpp + write(*,*) ' rcp (new value):',rcp + + else if (modif(1:len_trim(modif)).eq.'daysec') then + write(*,*) 'current value:',daysec + write(*,*) 'enter new value:' + 125 read(*,*,iostat=ierr) daysec + if(ierr.ne.0) goto 125 + write(*,*) + write(*,*) ' daysec (new value):',daysec + +! added by RW! + else if (modif(1:len_trim(modif)).eq.'year_day') then + write(*,*) 'current value:',year_day + write(*,*) 'enter new value:' + 126 read(*,*,iostat=ierr) year_day + if(ierr.ne.0) goto 126 + write(*,*) + write(*,*) ' year_day (new value):',year_day + + endif + enddo ! of do while(modif(1:1).ne.'hello') + + 999 continue + +!---------------------------------------------------------------------- +! Write values of physical constants after modifications +!----------------------------------------------------------------------- + + write(*,*) '*****************************************************' + write(*,*) 'Reading tab_cntrl when calling tabfi AFTER changes' + write(*,*) '*****************************************************' + write(*,5) '(1) = ngrid?',tab_cntrl(tab0+1),float(ngrid) + 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) periastr',tab_cntrl(tab0+15),periastr + write(*,5) '(16) apoastr',tab_cntrl(tab0+16),apoastr + 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(*,6) '(21) emin_turb',tab_cntrl(tab0+21),emin_turb + write(*,5) '(20) lmixmin',tab_cntrl(tab0+20),lmixmin + + write(*,5) '(26) emissiv',tab_cntrl(tab0+26),emissiv + write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1) + write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2) + write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1) + write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2) + write(*,5) '(33) dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1) + write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2) + + write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa + + write(*,*) + write(*,*) + + ENDIF ! of if (Lmodif == 1) + +!----------------------------------------------------------------------- +! Save some constants for later use (as routine arguments) +!----------------------------------------------------------------------- + p_omeg = omeg + p_g = g + p_cpp = cpp + p_mugaz = mugaz + p_daysec = daysec + p_rad=rad + + + END SUBROUTINE tabfi + +end module tabfi_mod