Index: trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F (revision 3232) +++ trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F (revision 3258) @@ -39,9 +39,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" @@ -82,5 +82,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,tsurf(1) REAL u(llm),v(llm) ! zonal, meridional wind REAL gru,grv ! prescribed "geostrophic" background wind @@ -90,7 +90,16 @@ REAL,ALLOCATABLE :: tsoil(:) ! subsurface soil temperature (K) ! REAL n2ice ! n2ice layer (kg.m-2) !not used anymore - integer :: i_n2_ice=0 ! tracer index of n2 ice - integer :: i_h2o_ice=0 ! tracer index of h2o ice - integer :: i_h2o_vap=0 ! tracer index of h2o vapor + integer :: i_n2=0 ! tracer index of n2 ice + integer :: i_ch4_ice=0 ! tracer index of ch4 ice + integer :: i_ch4_gas=0 ! tracer index of ch4 gas + integer :: i_co_ice=0 ! tracer index of co ice + integer :: i_co_gas=0 ! tracer index of co gas + integer :: i_prec_haze=0 ! tracer index of haze + integer :: i_haze=0 ! tracer index of haze + integer :: i_haze10=0 ! tracer index of haze + integer :: i_haze30=0 ! tracer index of haze + integer :: i_haze50=0 ! tracer index of haze + integer :: i_haze100=0 ! tracer index of haze + REAL emis(1) ! surface layer REAL q2(llm+1) ! Turbulent Kinetic Energy @@ -100,5 +109,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(:,:) @@ -112,5 +121,5 @@ REAL tmp1(0:llm),tmp2(0:llm) integer :: nq !=1 ! number of tracers - + character*2 str2 character (len=7) :: str7 @@ -123,4 +132,11 @@ real Hscale, Hmax, rho, dz +! pluto specific + real ch4ref ! % ch4 + real coref ! % co + real hazeref ! haze kg/kg + real prechazeref ! prec haze kg/kg + + ! added by RW for autozlevs computation logical autozlevs @@ -144,5 +160,5 @@ c INITIALISATION c======================================================================= -! 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) @@ -203,15 +219,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. @@ -223,15 +239,15 @@ 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 emin_turb = 1.e-6 ! energie minimale ~1.e-8 lmixmin = 30 ! longueur de melange ~100 - + 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 n2 snow (north) @@ -242,5 +258,5 @@ c ------------------------------------------------------ -c Load parameters from "run.def" and "gases.def" +c Load parameters from "run.def" and "gases.def" c ------------------------------------------------------ @@ -303,5 +319,5 @@ stop endif - + do iq=1,nq ! minimal version, just read in the tracer names, 1 per line @@ -313,15 +329,17 @@ endif enddo !of do iq=1,nq -! check for n2_ice / h2o_ice tracers: - i_n2_ice=0 - i_h2o_ice=0 - i_h2o_vap=0 +! check for n2 / h2o_ice tracers: + i_n2=0 do iq=1,nq if (tname(iq)=="n2") then - i_n2_ice=iq - elseif (tname(iq)=="h2o_ice") then - i_h2o_ice=iq - elseif (tname(iq)=="h2o_vap") then - i_h2o_vap=iq + i_n2=iq + elseif (tname(iq)=="ch4_ice") then + i_ch4_ice=iq + elseif (tname(iq)=="co_ice") then + i_co_ice=iq + elseif (tname(iq)=="ch4_gas") then + i_ch4_gas=iq + elseif (tname(iq)=="co_gas") then + i_co_gas=iq endif enddo @@ -339,9 +357,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) @@ -361,5 +379,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) ?' @@ -416,5 +434,5 @@ PRINT *,"STOP. I NEED year_day IN RCM1D.DEF." STOP - ELSE + ELSE PRINT *,"--> year_day = ",year_day ENDIF @@ -449,9 +467,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) @@ -461,5 +479,5 @@ ELSE PRINT *,"--> obliquit = ",obliquit - ENDIF + ENDIF psurf = -99999. @@ -519,6 +537,6 @@ nday=ndt - ndt=ndt*day_step - dtphys=daysec/day_step + ndt=ndt*day_step + dtphys=daysec/day_step write(*,*)"-------------------------------------" write(*,*)"-------------------------------------" @@ -545,5 +563,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 @@ -608,103 +626,59 @@ 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" do iq=1,nq - - txt="" - write(txt,"(a)") tname(iq) - write(*,*)" tracer:",trim(txt) - - ! n2 - if (txt.eq."n2_ice") then - q(:,iq)=0. ! kg/kg of atmosphere - qsurf(iq)=0. ! kg/m2 at the surface - ! Look for a "profile_n2_ice" input file - open(91,file='profile_n2_ice',status='old', - & form='formatted',iostat=ierr) - if (ierr.eq.0) then - read(91,*) qsurf(iq) - do ilayer=1,nlayer - read(91,*) q(ilayer,iq) - enddo - else - write(*,*) "No profile_n2_ice file!" - endif - close(91) - endif ! of if (txt.eq."n2") - - ! 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 - open(91,file='profile_h2o_vap',status='old', - & form='formatted',iostat=ierr) - if (ierr.eq.0) then - read(91,*) qsurf(iq) - do ilayer=1,nlayer - read(91,*) q(ilayer,iq) - enddo - else - write(*,*) "No profile_h2o_vap file!" - endif - close(91) - endif ! of if (txt.eq."h2o_vap") - - ! WATER ICE - if (txt.eq."h2o_ice") then - q(:,iq)=0. ! kg/kg of atmosphere - qsurf(iq)=0. ! kg/m2 at the surface - ! Look for a "profile_h2o_ice" input file - open(91,file='profile_h2o_ice',status='old', - & form='formatted',iostat=ierr) - if (ierr.eq.0) then - read(91,*) qsurf(iq) - do ilayer=1,nlayer - read(91,*) q(ilayer,iq) - enddo - else - write(*,*) "No profile_h2o_ice file!" - endif - close(91) - endif ! of if (txt.eq."h2o_ice") - - !_vap - if((txt .ne. 'h2o_vap') .and. - & (index(txt,'_vap' ) .ne. 0)) then - q(:,iq)=0. !kg/kg of atmosphere - qsurf(iq) = 0. !kg/kg of atmosphere - ! Look for a "profile_...._vap" input file - tracer_profile_file_name="" - tracer_profile_file_name='profile_'//txt - open(91,file=tracer_profile_file_name,status='old', - & form="formatted",iostat=ierr) - if (ierr .eq. 0) then - read(91,*)qsurf(iq) - do ilayer=1,nlayer - read(91,*)q(ilayer,iq) - enddo - else - write(*,*) "No initial profile " - write(*,*) " for this tracer :" - write(*,*) txt - endif - close(91) - endif ! (txt .eq. "_vap") - !_ice - if((txt.ne."h2o_ice") .and. - & (index(txt,'_ice' ) /= 0)) then - q(:,iq)=0. !kg/kg of atmosphere - qsurf(iq) = 0. !kg/kg of atmosphere - endif ! we only initialize the solid at 0 + + if (iq.eq.i_n2) then + DO ilayer=1,nlayer + q(ilayer,iq) = 1 + ENDDO + else if (iq.eq.i_ch4_gas) then + ch4ref=0.5 ! default value for vmr ch4 + PRINT *,'vmr CH4 (%) ?' + call getin("ch4ref",ch4ref) + write(*,*) " CH4 ref (%) = ",ch4ref + DO ilayer=1,nlayer + q(ilayer,iq) = 0.01*ch4ref*16./28. + ENDDO + else if (iq.eq.i_co_gas) then + coref=0.05 ! default value for vmr ch4 + PRINT *,'vmr CO (%) ?' + call getin("coref",coref) + write(*,*) " CO ref (%) = ",coref + DO ilayer=1,nlayer + q(ilayer,iq) = 0.01*coref*28./28. + ENDDO + ! else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30).or.(iq.eq.i_haze).or.(iq.eq.i_haze50).or.(iq.eq.i_haze100)) then + ! hazeref=0. ! default value for haze mmr + ! PRINT *,'qhaze (kg/kg) ?' + ! call getin("hazeref",hazeref) + ! write(*,*) " haze ref (kg/kg) = ",hazeref + ! DO ilayer=1,nlayer + ! q(ilayer,iq) = hazeref + ! ENDDO + ! else if (iq.eq.i_prec_haze) then + ! prechazeref=0. ! default value for vmr ch4 + ! PRINT *,'qprechaze (kg/kg) ?' + ! call getin("prechazeref",prechazeref) + ! write(*,*) " prec haze ref (kg/kg) = ",prechazeref + ! DO ilayer=1,nlayer + ! q(ilayer,iq) = prechazeref + ! ENDDO + + else + DO ilayer=1,nlayer + q(ilayer,iq) = 0. + ENDDO + endif enddo ! of do iq=1,nq - + endif ! of tracer @@ -712,5 +686,5 @@ c ------------------------------------------------------ ptif=2.E+0*omeg*sinlat(1) - + c vent geostrophique @@ -748,10 +722,10 @@ ! value if there is no snow - if(i_n2_ice.gt.0)then - qsurf(i_n2_ice)=0 ! default value for n2ice + if(i_n2.gt.0)then + qsurf(i_n2)=0 ! default value for n2ice print*,'Initial n2 ice on the surface (kg.m-2)' - call getin("n2ice",qsurf(i_n2_ice)) - write(*,*) " n2ice = ",qsurf(i_n2_ice) - IF (qsurf(i_n2_ice).ge.1.E+0) THEN + call getin("n2ice",qsurf(i_n2)) + write(*,*) " n2ice = ",qsurf(i_n2) + IF (qsurf(i_n2).ge.1.E+0) THEN ! if we have some n2 ice on the surface, change emissivity if (latitude(1).ge.0) then ! northern hemisphere @@ -793,5 +767,5 @@ if(autozlevs)then - + open(91,file="z2sig.def",form='formatted') read(91,*) Hscale @@ -800,6 +774,6 @@ enddo close(91) - - + + print*,'Hmax = ',Hmax,' km' print*,'Auto-shifting Hscale to:' @@ -807,7 +781,7 @@ Hscale = Hmax / log(psurf/pceil) print*,'Hscale = ',Hscale,' km' - + ! none of this matters if we dont care about zlay - + endif @@ -833,9 +807,9 @@ plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) ENDDO - + DO ilayer=1,nlayer play(ilayer)=aps(ilayer)+psurf*bps(ilayer) ENDDO - + @@ -914,5 +888,5 @@ c======================================================================= -c BOUCLE TEMPORELLE DU MODELE 1D +c BOUCLE TEMPORELLE DU MODELE 1D c======================================================================= @@ -933,5 +907,5 @@ ENDIF -c calcul du geopotentiel +c calcul du geopotentiel c ~~~~~~~~~~~~~~~~~~~~~ @@ -968,5 +942,5 @@ , day,time,dtphys, , plev,play,phi, - , u, v,temp, q, + , u, v,temp, q, , w, C - sorties @@ -976,5 +950,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. @@ -992,5 +966,5 @@ ENDDO c end if -c +c c c Calcul du temps au pas de temps suivant @@ -1056,4 +1030,4 @@ c ======================================================== end !rcm1d - - + +