Index: trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90 (revision 3936) @@ -164,7 +164,7 @@ flym_ipm(1)=flym_ipm(1)*pfluxuv/5. ELSE IF (.not.diurnal) THEN - mu_sol(:) = mu0(:) + mu_sol(:) = max(mu0(:),0.) mu_ipm(:) = max(mu_sol(:), 0.5) - flym_ipm(:)= mu0(:)*72.5e10 + flym_ipm(:)= mu_sol(:)*72.5e10 ELSE ! case with full fit to Gladstone et al. results ! 1) get longitude/latitude (in radian) of anti-subsolar point (max de mu0 - 180) @@ -182,7 +182,7 @@ DO ig=1,ngrid ! calculation of cosinus of incident angle for IPM flux - mu_sol(ig) = mu0(ig) + mu_sol(ig) = max(mu0(ig),0.) mu_ipm(ig) = max(mu_sol(ig), 0.5) - IF (mu0(ig).LT.1.e-4) THEN ! Daytime + IF (mu_sol(ig).LE.0.) THEN ! Nighttime ! Distance to subsolar point dist=acos(sin(latitude(ig))*sin(latit)+cos(latitude(ig))* & @@ -191,6 +191,6 @@ flym_ipm(ig)=(valmin_dl-valmin)/(90.**puis)*(dist)**puis & +valmin - ELSE ! Nightime - flym_ipm(ig)= mu0(ig)*(valmax-valmin_dl)+valmin_dl + ELSE ! Daytime + flym_ipm(ig)= mu_sol(ig)*(valmax-valmin_dl)+valmin_dl ENDIF ! Proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) ? Index: trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90 (revision 3936) @@ -43,6 +43,4 @@ ! rho_q(nq) ! tracer densities (kg.m-3) ! qext(nq) ! Single Scat. Extinction coeff at 0.67 um -! alpha_lift(nq) ! saltation vertical flux/horiz flux ratio (m-1) -! alpha_devil(nq) ! lifting coeeficient by dust devil ! rho_dust ! Mars dust density ! rho_ice ! Water ice density @@ -110,9 +108,6 @@ IF (.NOT.ALLOCATED(rho_q)) ALLOCATE(rho_q(nq)) IF (.NOT.ALLOCATED(qext)) ALLOCATE(qext(nq)) - IF (.NOT.ALLOCATED(alpha_lift)) ALLOCATE(alpha_lift(nq)) - IF (.NOT.ALLOCATED(alpha_devil)) ALLOCATE(alpha_devil(nq)) IF (.NOT.ALLOCATED(qextrhor)) ALLOCATE(qextrhor(nq)) ! IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq)) - IF (.NOT.ALLOCATED(is_chim)) ALLOCATE(is_chim(nqtot)) IF (.NOT.ALLOCATED(is_rad)) ALLOCATE(is_rad(nqtot)) IF (.NOT.ALLOCATED(is_recomb)) ALLOCATE(is_recomb(nqtot)) @@ -125,10 +120,7 @@ !! initialization - alpha_lift(:) = 0. - alpha_devil(:) = 0. mmol(:) = 0. aki(:) = 0. cpi(:) = 0. - is_chim(:) = 0 is_rad(:) = 0 is_recomb(:) = 0 @@ -383,8 +375,4 @@ if (is_master) close(407) - - ! Calculate number of species in the chemistry - nesp = sum(is_chim) - write(*,*) 'Number of species in the chemistry nesp = ',nesp ! Calculate number of microphysical tracer @@ -451,6 +439,4 @@ write(*,*) Write(*,*) '******** initracer : dust transport parameters :' - write(*,*) 'alpha_lift = ', alpha_lift - write(*,*) 'alpha_devil = ', alpha_devil write(*,*) 'radius = ', radius write(*,*) 'Qext = ', qext @@ -504,14 +490,4 @@ cpi(iq) end if - ! option is_chim - if (index(tracline,'is_chim=' ) /= 0) then - read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) & - is_chim(iq) - write(*,*) ' Parameter value (traceur.def) : is_chim=', & - is_chim(iq) - else - write(*,*) ' Parameter value (default) : is_chim=', & - is_chim(iq) - end if ! option is_rad if (index(tracline,'is_rad=') /= 0) then Index: trunk/LMDZ.PLUTO/libf/phypluto/interpolateCO2CH4.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/interpolateCO2CH4.F90 (revision 3935) +++ (revision ) @@ -1,140 +1,0 @@ -subroutine interpolateN2CH4(wn,temp,presN2,presCH4,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the N2-CH4 CIA opacity, using a lookup table from - ! Turbet et al. 2020, Icarus, Volume 346, article id. 113762 - ! - ! Authors - ! ------- - ! M. Turbet (2023) - ! - !================================================================== - - use datafile_mod, only: datadir - implicit none - - ! input - double precision wn ! wavenumber (cm^-1) - double precision temp ! temperature (Kelvin) - double precision presN2 ! N2 partial pressure (Pascals) - double precision presCH4 ! CH4 partial pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=240) - parameter(nT=6) - double precision, parameter :: losch = 2.6867774e19 - ! Loschmit's number (molecule cm^-3 at STP) - ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 - ! see Richard et al. 2011, JQSRT for details - - double precision amagatN2 - double precision amagatCH4 - double precision wn_arr(nS) - double precision temp_arr(nT) - double precision abs_arr(nS,nT) - - integer k,iT - logical firstcall - - save wn_arr, temp_arr, abs_arr !read by master - - character*100 dt_file - integer strlen,ios - - character(LEN=*), parameter :: fmat1 = "(A21,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" - - character*21 bleh - double precision blah, Ttemp - integer nres - integer ind - - if(temp.gt.600)then - print*,'Your temperatures are too high for this N2-CH4 CIA dataset.' - print*,'Please run mixed N2-CH4 atmospheres below T = 600 K.' - stop - endif - - if(temp.lt.100)then - print*,'Your temperatures are too low for this N2-CH4 CIA dataset.' - print*,'Please run mixed N2-CH4 atmospheres above T = 100 K.' - stop - endif - - amagatN2 = (273.15/temp)*(presN2/101325.0) - amagatCH4 = (273.15/temp)*(presCH4/101325.0) - - if(firstcall)then ! called by sugas_corrk only - print*,'----------------------------------------------------' - print*,'Initialising N2-CH4 continuum from Turbet et al. 2020' - - ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/N2-CH4_TURBET2020.cia' - - -!$OMP MASTER - open(33,file=dt_file,form='formatted',status='old',iostat=ios) - if (ios.ne.0) then ! file not found - write(*,*) 'Error from interpolateN2CH4' - write(*,*) 'Data file ',trim(dt_file),' not found.' - write(*,*) 'Check that your path to datagcm:',trim(datadir) - write(*,*) 'is correct. You can change it in callphys.def with:' - write(*,*) 'datadir = /absolute/path/to/datagcm' - write(*,*) 'Also check that the continuum data continuum_data/N2-CH4_TURBET2020.cia is there.' - call abort - else - - do iT=1,nT - - read(33,fmat1) bleh,blah,blah,nres,Ttemp - if(nS.ne.nres)then - print*,'Resolution given in file: ',trim(dt_file) - print*,'is ',nres,', which does not match nS.' - print*,'Please adjust nS value in interpolateN2CH4.F90' - stop - endif - temp_arr(iT)=Ttemp - - do k=1,nS - read(33,*) wn_arr(k),abs_arr(k,it) - write(*,*) 'for k = ', k, ' we have ', wn_arr(k),abs_arr(k,it) - end do - - end do - - endif - close(33) -!$OMP END MASTER -!$OMP BARRIER - - print*,'interpolateN2CH4: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' N2 partial pressure ',presN2,' Pa' - print*,' and CH4 partial pressure ',presCH4,' Pa' - - endif - - call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) - -! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' -! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' - - abcoef=abcoef*100.0*amagatN2*amagatCH4 ! convert to m^-1 - -! print*,'We have ',amagatN2,' amagats of N2' -! print*,'and ',amagaCH4,' amagats of CH4' -! print*,'So the absorption is ',abcoef,' m^-1' - - -! unlike for Rayleigh scattering, we do not currently weight by the BB function -! however our bands are normally thin, so this is no big deal. - - - return -end subroutine interpolateN2CH4 - Index: trunk/LMDZ.PLUTO/libf/phypluto/interpolateN2CH4.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/interpolateN2CH4.F90 (revision 3936) +++ trunk/LMDZ.PLUTO/libf/phypluto/interpolateN2CH4.F90 (revision 3936) @@ -0,0 +1,140 @@ +subroutine interpolateN2CH4(wn,temp,presN2,presCH4,abcoef,firstcall,ind) + + !================================================================== + ! + ! Purpose + ! ------- + ! Calculates the N2-CH4 CIA opacity, using a lookup table from + ! Turbet et al. 2020, Icarus, Volume 346, article id. 113762 + ! + ! Authors + ! ------- + ! M. Turbet (2023) + ! + !================================================================== + + use datafile_mod, only: datadir + implicit none + + ! input + double precision wn ! wavenumber (cm^-1) + double precision temp ! temperature (Kelvin) + double precision presN2 ! N2 partial pressure (Pascals) + double precision presCH4 ! CH4 partial pressure (Pascals) + + ! output + double precision abcoef ! absorption coefficient (m^-1) + + integer nS,nT + parameter(nS=240) + parameter(nT=6) + double precision, parameter :: losch = 2.6867774e19 + ! Loschmit's number (molecule cm^-3 at STP) + ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 + ! see Richard et al. 2011, JQSRT for details + + double precision amagatN2 + double precision amagatCH4 + double precision wn_arr(nS) + double precision temp_arr(nT) + double precision abs_arr(nS,nT) + + integer k,iT + logical firstcall + + save wn_arr, temp_arr, abs_arr !read by master + + character*100 dt_file + integer strlen,ios + + character(LEN=*), parameter :: fmat1 = "(A21,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" + + character*21 bleh + double precision blah, Ttemp + integer nres + integer ind + + if(temp.gt.600)then + print*,'Your temperatures are too high for this N2-CH4 CIA dataset.' + print*,'Please run mixed N2-CH4 atmospheres below T = 600 K.' + stop + endif + + if(temp.lt.100)then + print*,'Your temperatures are too low for this N2-CH4 CIA dataset.' + print*,'Please run mixed N2-CH4 atmospheres above T = 100 K.' + stop + endif + + amagatN2 = (273.15/temp)*(presN2/101325.0) + amagatCH4 = (273.15/temp)*(presCH4/101325.0) + + if(firstcall)then ! called by sugas_corrk only + print*,'----------------------------------------------------' + print*,'Initialising N2-CH4 continuum from Turbet et al. 2020' + + ! 1.1 Open the ASCII files + dt_file=TRIM(datadir)//'/continuum_data/N2-CH4_TURBET2020.cia' + + +!$OMP MASTER + open(33,file=dt_file,form='formatted',status='old',iostat=ios) + if (ios.ne.0) then ! file not found + write(*,*) 'Error from interpolateN2CH4' + write(*,*) 'Data file ',trim(dt_file),' not found.' + write(*,*) 'Check that your path to datagcm:',trim(datadir) + write(*,*) 'is correct. You can change it in callphys.def with:' + write(*,*) 'datadir = /absolute/path/to/datagcm' + write(*,*) 'Also check that the continuum data continuum_data/N2-CH4_TURBET2020.cia is there.' + call abort + else + + do iT=1,nT + + read(33,fmat1) bleh,blah,blah,nres,Ttemp + if(nS.ne.nres)then + print*,'Resolution given in file: ',trim(dt_file) + print*,'is ',nres,', which does not match nS.' + print*,'Please adjust nS value in interpolateN2CH4.F90' + stop + endif + temp_arr(iT)=Ttemp + + do k=1,nS + read(33,*) wn_arr(k),abs_arr(k,it) + write(*,*) 'for k = ', k, ' we have ', wn_arr(k),abs_arr(k,it) + end do + + end do + + endif + close(33) +!$OMP END MASTER +!$OMP BARRIER + + print*,'interpolateN2CH4: At wavenumber ',wn,' cm^-1' + print*,' temperature ',temp,' K' + print*,' N2 partial pressure ',presN2,' Pa' + print*,' and CH4 partial pressure ',presCH4,' Pa' + + endif + + call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) + +! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' +! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' + + abcoef=abcoef*100.0*amagatN2*amagatCH4 ! convert to m^-1 + +! print*,'We have ',amagatN2,' amagats of N2' +! print*,'and ',amagaCH4,' amagats of CH4' +! print*,'So the absorption is ',abcoef,' m^-1' + + +! unlike for Rayleigh scattering, we do not currently weight by the BB function +! however our bands are normally thin, so this is no big deal. + + + return +end subroutine interpolateN2CH4 + Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3936) @@ -33,6 +33,5 @@ igcm_n2,igcm_ch4_gas,igcm_ch4_ice,igcm_haze,& igcm_co_gas,igcm_co_ice,igcm_prec_haze,lw_n2,lw_ch4,lw_co,& - alpha_lift, alpha_devil, qextrhor, & - nesp, is_chim + qextrhor use time_phylmdz_mod, only: diagfi_output_rate, startfi_output_rate, nday use phyetat0_mod, only: phyetat0 @@ -535,6 +534,4 @@ zdtsw(:,:) = 0.0 zdtlw(:,:) = 0.0 - zdqc(:,:,:)=0. - zdqsc(:,:)=0. ! Initialize tracer names, indexes and properties. @@ -1447,4 +1444,7 @@ endif + zdqc(:,:,:)=0.0 + zdqsc(:,:)=0.0 + call condense_n2(ngrid,nlayer,nq,ptimestep, & capcal,pplay,pplev,tsurf,pt, & @@ -2276,4 +2276,5 @@ call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw) call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn) + call write_output("mu0","cos zenith anlge","",mu0) call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw) call write_output("GND","heat flux from ground","W m-2",fluxgrd) Index: trunk/LMDZ.PLUTO/libf/phypluto/stelang.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/stelang.F (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/stelang.F (revision 3936) @@ -85,5 +85,5 @@ c----- AF24: replaced flat param by 0 - rap = 1./((1.-0)**2) +c rap = 1./((1.-0)**2) 100 CONTINUE @@ -98,9 +98,9 @@ C DO jl=1,kgrid - ztim1=psilat(jl)*ptim1*rap + ztim1=psilat(jl)*ptim1 !*rap ztim2=pcolat(jl)*ptim2 ztim3=pcolat(jl)*ptim3 pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl) - pmu0(jl)=pmu0(jl)/SQRT(pcolat(jl)**2+(rap**2)*(psilat(jl)**2)) +c pmu0(jl)=pmu0(jl)/SQRT(pcolat(jl)**2+(rap**2)*(psilat(jl)**2)) ENDDO @@ -112,8 +112,8 @@ IF (pmu0(jl).gt.0.) THEN pfract(jl)=1. -c pmu0(jl)=sqrt(1224.*pmu0(jl)*pmu0(jl)+1.)/35. - ELSE -c pmu0(jl)=0. - pfract(jl)=0. +c pmu0(jl)=sqrt(1224.*pmu0(jl)*pmu0(jl)+1.)/35. + ELSE +c pmu0(jl)=0. + pfract(jl)=0. ENDIF ENDDO Index: trunk/LMDZ.PLUTO/libf/phypluto/surface_nature.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/surface_nature.F (revision 3935) +++ (revision ) @@ -1,64 +1,0 @@ - SUBROUTINE surface_nature(ngrid,nq,obliquit,qsurf,qsurfliquid - & ,qsurfsnow,rnat,oceanarea) - - USE surfdat_h - USE comsoil_h - USE geometry_mod, ONLY: cell_area - USE tracer_h - - IMPLICIT none - -!================================================================== -! -! Purpose -! ------- -! Defines a few things -! -! Authors -! ------- -! B. Charnay (2010) -! -! Called by -! --------- -! physiq.F -! -! Calls -! ----- -! none -! -! Notes -! ----- -! rnat is terrain type: 0-ocean; 1-continent; 2-continental ice -! -!================================================================== - - integer ngrid,nq - - REAL qsurf(ngrid,nq),ps(ngrid) - REAL qsurfliquid(ngrid) - REAL qsurfsnow(ngrid) - INTEGER iq, ig - INTEGER rnat(ngrid) - REAL oceanarea - REAL obliquit - - do ig=1,ngrid - rnat(ig)=1 - dryness(ig)=1 !(coefficient for evaporation) - if (inertiedat(ig,1).gt.1E4) then - rnat(ig)=0 - end if - end do - -! surface of all the oceans - - oceanarea=0. - do ig=1,ngrid - if (rnat(ig).eq.0)then - oceanarea=oceanarea+cell_area(ig) - end if - enddo - - return - end - Index: trunk/LMDZ.PLUTO/libf/phypluto/tabfi_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/tabfi_mod.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/tabfi_mod.F90 (revision 3936) @@ -153,4 +153,6 @@ dtemisice(:)=0 !time scale for snow metamorphism volcapa=1000000 ! volumetric heat capacity of subsurface + tpal=0. + adjust=0. ELSE Index: trunk/LMDZ.PLUTO/libf/phypluto/testconserv.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/testconserv.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/testconserv.F90 (revision 3936) @@ -12,5 +12,5 @@ ! Purpose ! ------- -! Test conservation of tracers +! Test conservation of gas trace tracers (CO and CH4 for now) ! ! Inputs @@ -96,9 +96,4 @@ enddo - iq = igcm1 - dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig) - Wtots = Wtots + pqcs(ig,iq)*cell_area(ig) - ncons(ig)=ncons(ig)+zdqs(ig,iq) - IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN iq = igcm2 Index: trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90 (revision 3936) @@ -12,6 +12,5 @@ integer, save :: nqtot ! total number of tracers - integer, save :: nesp ! number of species in the chemistry -!$OMP THREADPRIVATE(nqtot,nesp) +!$OMP THREADPRIVATE(nqtot) logical :: moderntracdef=.false. ! Standard or modern traceur.def @@ -25,6 +24,4 @@ real, save, allocatable :: rho_q(:) ! tracer densities (kg.m-3) real, save, allocatable :: qext(:) ! Single Scat. Extinction coeff at 0.67 um - real, save, allocatable :: alpha_lift(:) ! saltation vertical flux/horiz flux ratio (m-1) - real, save, allocatable :: alpha_devil(:) ! lifting coeeficient by dust devil real, save, allocatable :: qextrhor(:) ! Intermediate for computing opt. depth from q @@ -41,10 +38,9 @@ integer,save :: nmono real,save :: ref_r0 ! for computing reff=ref_r0*r0 (in log.n. distribution) -!$OMP THREADPRIVATE(noms,mmol,aki,cpi,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor, & +!$OMP THREADPRIVATE(noms,mmol,aki,cpi,radius,rho_q,qext,qextrhor, & !$OMP varian,r3n_q,rho_dust,rho_ice,rho_n2,lw_n2,ref_r0) - integer, save, allocatable :: is_chim(:) ! 1 if tracer used in chemistry, else 0 integer, save, allocatable :: is_rad(:) ! 1 if "" "" in radiative transfer, else 0 -!$OMP THREADPRIVATE(is_chim,is_rad) +!$OMP THREADPRIVATE(is_rad) integer, save, allocatable :: is_recomb(:) ! 1 if tracer used in recombining scheme, else 0 (if 1, must have is_rad=1)