Index: /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F (revision 3671) +++ /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F (revision 3672) @@ -36,5 +36,5 @@ use mod_interface_dyn_phys, only: init_interface_dyn_phys use inifis_mod, only: inifis - use phys_state_var_mod, only: phys_state_var_init + use phys_state_var_mod, only: phys_state_var_init, tsurf, tsoil use physiq_mod, only: physiq use restart1D_mod, only: writerestart1D @@ -48,9 +48,9 @@ !================================================================== -! +! ! Purpose ! ------- ! Run the physics package of the universal model in a 1D column. -! +! ! It can be compiled with a command like (e.g. for 25 layers): ! "makegcm -p std -d 25 rcm1d" @@ -95,5 +95,5 @@ REAL play(llm) ! Pressure at the middle of the layers (Pa) REAL plev(llm+1) ! intermediate pressure levels (pa) - REAL psurf,tsurf(1) + REAL psurf REAL u(llm),v(llm) ! zonal, meridional wind REAL gru,grv ! prescribed "geostrophic" background wind @@ -101,5 +101,4 @@ REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg) REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2) - REAL,ALLOCATABLE :: tsoil(:) ! subsurface soil temperature (K) ! REAL co2ice ! co2ice layer (kg.m-2) !not used anymore integer :: i_co2_ice=0 ! tracer index of co2 ice @@ -114,5 +113,5 @@ REAL du(llm),dv(llm),dtemp(llm) REAL dudyn(llm),dvdyn(llm),dtempdyn(llm) - REAL dpsurf(1) + REAL dpsurf(1) REAL,ALLOCATABLE :: dq(:,:) REAL,ALLOCATABLE :: dqdyn(:,:) @@ -126,5 +125,5 @@ REAL tmp1(0:llm),tmp2(0:llm) integer :: nq !=1 ! number of tracers - + character*2 str2 character (len=7) :: str7 @@ -181,5 +180,5 @@ endif -! check if 'rcm1d.def' file is around +! check if 'rcm1d.def' file is around open(90,file='rcm1d.def',status='old',form='formatted', & iostat=ierr) @@ -257,15 +256,15 @@ endif close(90) - + ! Initialize dimphy module - call init_dimphy(1,llm) - + call init_dimphy(1,llm) + ! now initialize arrays using phys_state_var_init ! but first initialise naerkind (from callphys.def) naerkind=0 !default call getin("naerkind",naerkind) - + call phys_state_var_init(nq) - + saveprofile=.false. saveprofile=.true. @@ -277,13 +276,13 @@ pi=2.E+0*asin(1.E+0) -c Parametres Couche limite et Turbulence +c Parametres Couche limite et Turbulence c -------------------------------------- - z0 = 1.e-2 ! surface roughness (m) ~0.01 - + z0 = 1.e-2 ! surface roughness (m) ~0.01 + c propriete optiques des calottes et emissivite du sol c ---------------------------------------------------- emissiv= 0.95 ! Emissivite du sol martien ~.95 emisice(1)=0.95 ! Emissivite calotte nord - emisice(2)=0.95 ! Emissivite calotte sud + emisice(2)=0.95 ! Emissivite calotte sud iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north) @@ -294,5 +293,5 @@ c ------------------------------------------------------ -c Load parameters from "run.def" and "gases.def" +c Load parameters from "run.def" and "gases.def" c ------------------------------------------------------ @@ -355,5 +354,5 @@ stop endif - + do iq=1,nq ! minimal version, just read in the tracer names, 1 per line @@ -391,9 +390,9 @@ nq=0 ! still, make allocations for 1 dummy tracer - allocate(tname(1)) + allocate(tname(1)) allocate(qsurf(1)) allocate(q(llm,1)) allocate(dq(llm,1)) - + ! Check that tracer boolean is compliant with number of tracers ! -- otherwise there is an error (and more generally we have to be consistent) @@ -413,5 +412,5 @@ !!! GEOPOTENTIAL. useless in 1D because control by surface pressure phisfi(1)=0.E+0 - !!! LATITUDE. can be set. + !!! LATITUDE. can be set. latitude=0 ! default value for latitude PRINT *,'latitude (in degrees) ?' @@ -468,5 +467,5 @@ PRINT *,"STOP. I NEED year_day IN RCM1D.DEF." STOP - ELSE + ELSE PRINT *,"--> year_day = ",year_day ENDIF @@ -501,9 +500,9 @@ PRINT *,"STOP. peri_day.gt.year_day" STOP - ELSE + ELSE PRINT *,"--> peri_day = ", peri_day - ENDIF - - obliquit = -99999. + ENDIF + + obliquit = -99999. PRINT *,'OBLIQUITY in deg ?' call getin("obliquit",obliquit) @@ -513,5 +512,5 @@ ELSE PRINT *,"--> obliquit = ",obliquit - ENDIF + ENDIF if (restart) then @@ -552,5 +551,5 @@ c Date (en sols depuis le solstice de printemps) du debut du run !if (restart) then - ! ! day + ! ! day ! ierr=NF90_INQ_VARID(nid_restart1D,'day',did) ! if (ierr==NF90_NOERR) then @@ -609,6 +608,6 @@ nday=ndt - ndt=ndt*day_step - dtphys=daysec/day_step + ndt=ndt*day_step + dtphys=daysec/day_step write(*,*)"-------------------------------------" write(*,*)"-------------------------------------" @@ -635,5 +634,5 @@ !!! - read callphys.def !!! - calculate sine and cosine of longitude and latitude -!!! - calculate mugaz and cp +!!! - calculate mugaz and cp !!! NB: some operations are being done dummily in inifis in 1D !!! - physical constants: nevermind, things are done allright below @@ -645,9 +644,4 @@ nsoil=nsoilmx - allocate(tsoil(nsoilmx)) - !! those are defined in comsoil_h.F90 - IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx)) - IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1)) - IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx)) ! At this point, both getin() and getin_p() functions have been used, @@ -698,11 +692,11 @@ ENDDO ENDDO - + DO iq=1,nq qsurf(iq) = 0. ENDDO - - if (tracer) then ! Initialize tracers here. - + + if (tracer) then ! Initialize tracers here. + write(*,*) "rcm1d : initializing tracers profiles" @@ -738,13 +732,13 @@ endif else ! restart - + txt="" write(txt,"(a)") tname(iq) write(*,*)" tracer:",trim(txt) - + ! CO2 if (txt.eq."co2_ice") then q(:,iq)=0. ! kg/kg of atmosphere - qsurf(iq)=0. ! kg/m2 at the surface + qsurf(iq)=0. ! kg/m2 at the surface ! Look for a "profile_co2_ice" input file open(91,file='profile_co2_ice',status='old', @@ -760,10 +754,10 @@ close(91) endif ! of if (txt.eq."co2") - + ! WATER VAPOUR if (txt.eq."h2o_vap") then q(:,iq)=0. ! kg/kg of atmosphere qsurf(iq)=0. ! kg/m2 at the surface - ! Look for a "profile_h2o_vap" input file + ! Look for a "profile_h2o_vap" input file open(91,file='profile_h2o_vap',status='old', & form='formatted',iostat=ierr) @@ -778,5 +772,5 @@ close(91) endif ! of if (txt.eq."h2o_vap") - + ! WATER ICE if (txt.eq."h2o_ice") then @@ -807,10 +801,10 @@ open(91,file=tracer_profile_file_name,status='old', & form="formatted",iostat=ierr) - if (ierr .eq. 0) then + if (ierr .eq. 0) then read(91,*)qsurf(iq) - do ilayer=1,nlayer + do ilayer=1,nlayer read(91,*)q(ilayer,iq) - enddo - else + enddo + else write(*,*) "No initial profile " write(*,*) " for this tracer :" @@ -818,7 +812,7 @@ endif close(91) - endif ! (txt .eq. "_vap") - !_ice - if((txt.ne."h2o_ice") .and. + endif ! (txt .eq. "_vap") + !_ice + if((txt.ne."h2o_ice") .and. & (index(txt,'_ice' ) /= 0)) then q(:,iq)=0. !kg/kg of atmosphere @@ -827,5 +821,5 @@ endif ! restart enddo ! of do iq=1,nq - + endif ! of tracer @@ -833,5 +827,5 @@ c ------------------------------------------------------ ptif=2.E+0*omeg*sinlat(1) - + c vent geostrophique @@ -941,5 +935,5 @@ if(autozlevs)then - + open(91,file="z2sig.def",form='formatted') read(91,*) Hscale @@ -948,6 +942,6 @@ enddo close(91) - - + + print*,'Hmax = ',Hmax,' km' print*,'Auto-shifting Hscale to:' @@ -955,7 +949,7 @@ Hscale = Hmax / log(psurf/pceil) print*,'Hscale = ',Hscale,' km' - + ! none of this matters if we dont care about zlay - + endif @@ -981,9 +975,9 @@ plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) ENDDO - + DO ilayer=1,nlayer play(ilayer)=aps(ilayer)+psurf*bps(ilayer) ENDDO - + @@ -1061,5 +1055,5 @@ DO isoil=1,nsoil inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia - tsoil(isoil)=tsurf(1) ! soil temperature + tsoil(1,isoil)=tsurf(1) ! soil temperature ENDDO @@ -1139,9 +1133,9 @@ & dtphys,time, & tsurf,tsoil,emis,albedo,q2,qsurf, - & cloudfrac,totcloudfrac,hice, + & cloudfrac,totcloudfrac,hice, & rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice) endif c======================================================================= -c BOUCLE TEMPORELLE DU MODELE 1D +c BOUCLE TEMPORELLE DU MODELE 1D c======================================================================= @@ -1162,5 +1156,5 @@ ENDIF -c calcul du geopotentiel +c calcul du geopotentiel c ~~~~~~~~~~~~~~~~~~~~~ @@ -1197,5 +1191,5 @@ , day,time,dtphys, , plev,play,phi, - , u, v,temp, q, + , u, v,temp, q, , w, C - sorties @@ -1205,5 +1199,5 @@ c evolution du vent : modele 1D c ----------------------------- - + c la physique calcule les derivees temporelles de u et v. c on y rajoute betement un effet Coriolis. @@ -1221,5 +1215,5 @@ ENDDO c end if -c +c c c Calcul du temps au pas de temps suivant @@ -1281,5 +1275,5 @@ if (lastcall) then call writerestart1D('restart1D.nc',nlayer,nsoil,day,time,psurf, - & temp,tsoil,u,v,nq,q) + & temp,tsoil,u,v,nq,q) endif @@ -1290,4 +1284,4 @@ c ======================================================== end !rcm1d - - + + Index: /trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F (revision 3671) +++ /trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F (revision 3672) @@ -34,5 +34,5 @@ use mod_interface_dyn_phys, only: init_interface_dyn_phys use inifis_mod, only: inifis - use phys_state_var_mod, only: phys_state_var_init + use phys_state_var_mod, only: phys_state_var_init, tsurf, tsoil use physiq_mod, only: physiq implicit none @@ -84,5 +84,5 @@ REAL play(llm) ! Pressure at the middle of the layers (Pa) REAL plev(llm+1) ! intermediate pressure levels (pa) - REAL psurf,tsurf(1) + REAL psurf REAL u(llm),v(llm) ! zonal, meridional wind REAL gru,grv ! prescribed "geostrophic" background wind @@ -90,5 +90,4 @@ REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg) REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2) - REAL,ALLOCATABLE :: tsoil(:) ! subsurface soil temperature (K) ! REAL n2ice ! n2ice layer (kg.m-2) !not used anymore integer :: i_n2=0 ! tracer index of n2 ice @@ -577,9 +576,4 @@ nsoil=nsoilmx - allocate(tsoil(nsoilmx)) - !! those are defined in comsoil_h.F90 - IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx)) - IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1)) - IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx)) ! At this point, both getin() and getin_p() functions have been used, @@ -866,5 +860,5 @@ OPEN(14,file='proftsoil',status='old',form='formatted',err=101) DO isoil=1,nsoil - READ (14,*) tsoil(isoil) + READ (14,*) tsoil(1,isoil) inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia ENDDO @@ -874,8 +868,8 @@ CLOSE(14) - else + else DO isoil=1,nsoil inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia - tsoil(isoil)=tsurf(1) ! soil temperature + tsoil(1,isoil)=tsurf(1) ! soil temperature ENDDO endif @@ -1065,5 +1059,5 @@ endif ! save haze profile - if (haze.and.lecthaze.eq.1) then + if (haze.and.lecthaze.eq.1) then OPEN(16,file='profhaze.out',form='formatted') DO iq=1,nq