Commit 78d28175 authored by Thomas Dubos's avatar Thomas Dubos
Browse files

Polish RUN/ and create CHIMERE/

parent 0bd87daa
# echo 40 120 2 | gawk -f domains-geom.awk
# p = A.pref + B.ps
# pref=10^5 Pa
BEGIN {pref=1000}
{n=$1;ptop=$2;llthick=$3}
END {
dp1 = llthick
p1 = pref - llthick
a = 200 / llthick
#find lowest n200 such that layer thickness increases with altitude
n200min=4 #lowest reasonable number of layers in the lowest 200 hPa of the atmosphere
for(n200=n200min;n200<n;n200++) {
# Equation: a = (r^n - 1)/(r - 1) : find r to have n200 layers in lowest 200 hPa of the atmosphere
x = 1.01
for(i=1;i<=30;i++) {
f = 1 - a
for(k=1;k<=n200-1;k++) {
f = f + x^k
}
fp = 1
for(k=1;k<=n200-2;k++) {
fp = fp + (k+1)*x^k
}
x = x - f/fp
}
# compute thickness of highest layer in the lowest 200 hPa of the atmosphere
th = dp1 * x^(n200 - 1)
# thickness of layers abobe 200 hPa
dp2 = (800 - ptop) / (n - n200)
if(dp2 > th) {break}
}
#write n200 lowest layers in output file
p = pref
dp = dp1
for(k=1;k<=n200;k++) {
p = p - dp
dp = x*dp
ca = (ptop/pref)*(p1-p)/(p1-ptop)
cb = (p1/pref)*(p-ptop)/(p1-ptop)
if(ca<0) {ca = 0}
if(cb<0) {cb = 0}
printf "%9.5f%9.5f\n",ca,cb,p,dp/x
}
# create n - n200 layers between 800 hPa and top of model
dp = (800 - ptop) / (n - n200)
for(k=1;k<=(n-n200);k++) {
p = p - dp
ca = (ptop/pref)*(p1-p)/(p1-ptop)
cb = (p1/pref)*(p-ptop)/(p1-ptop)
if(ca<0) {ca = 0}
if(cb<0) {cb = 0}
printf "%9.5f%9.5f\n",ca,cb,p,dp
}
}
subroutine ini_temis
use chimere_common
use netcdf
implicit none
integer :: id, ierr, ifn_fntracer, k, new
integer, parameter :: manydata=1000
character(len=splen)tmpname
real(kind=iprec), allocatable,dimension(:) :: tmplon,tmplat
integer,external :: interdat
!------------------------------------------------------------------------------------------
print*,'prep_temis: Open and read tracer data file: ',trim(fntracer)
call opfi(ifn_fntracer,fntracer,'f','o')
nline_tracer=0
read(ifn_fntracer,*) ! first line is a comment line
do id=1,manydata
nline_tracer=nline_tracer+1
read(ifn_fntracer,*,end=1001)tmpname
enddo
1001 continue
close(ifn_fntracer)
nline_tracer=nline_tracer-1
if(chimverb.ge.5)print*,nline_tracer,' lines found in the raw tracer data file'
allocate(lonemis(nline_tracer))
allocate(latemis(nline_tracer))
allocate(izoemis(nline_tracer))
allocate(imeemis(nline_tracer))
allocate(trc_name(nline_tracer))
allocate(trc_index(nline_tracer))
allocate(tmptyp(nline_tracer))
allocate(tmplon(nline_tracer))
allocate(tmplat(nline_tracer))
allocate(tmpdate1(nline_tracer))
allocate(tmpdate2(nline_tracer))
allocate(tmpmass(nline_tracer))
allocate(tmpaltmin(nline_tracer))
allocate(tmpaltmax(nline_tracer))
allocate(tmpmolmass(nline_tracer))
!----------------------------------------------------------
! Second read of the data to have the values
if(chimverb.ge.5)print*,'prep_temis: re-open tracer file'
call opfi(ifn_fntracer,fntracer,'f','o')
read(ifn_fntracer,*) ! first line is a comment line
if(chimverb.ge.5)print*,'prep_temis reading the following lines from:'
if(chimverb.ge.5)print*,fntracer
if(chimverb.ge.5)print*,'==================================================='
do id=1,nline_tracer
read(ifn_fntracer,*,end=1002)trc_name(id),tmptyp(id),tmplon(id),tmplat(id), &
& tmpdate1(id),tmpdate2(id),tmpmass(id),tmpaltmin(id),tmpaltmax(id),tmpmolmass(id)
if(chimverb.ge.5)print*,'prep_temis :',trc_name(id),tmptyp(id),tmplon(id),tmplat(id), &
& tmpdate1(id),tmpdate2(id),tmpmass(id),tmpaltmin(id),tmpaltmax(id),tmpmolmass(id)
lonemis(id)=tmplon(id)
latemis(id)=tmplat(id)
! Necessary conversion in [cm] since routine AFTER met2chim
tmpaltmin(id)=tmpaltmin(id)*m2cm
tmpaltmax(id)=tmpaltmax(id)*m2cm
if((tmptyp(id).eq.2).and.(tmpmolmass(id).ne.100.))then
print*, 'WARNING: Molar mass must be specified to 100. for aerosol tracers'
print*, 'WARNING: Resetting molar mass to 100. for ',trc_name(id)
tmpmolmass(id)=100.
endif
enddo
1002 continue
close(ifn_fntracer)
if(chimverb.ge.5)print*,'==================================================='
print*,'prep_temis: successful reading of TRACERS.data'
!----------------------------------------------------------
! Search for the real number of different tracers
print*,'Search for the real number of different tracers'
ntrc=0
do id=1,nline_tracer
new=1
if(id.gt.1)then
do k=1,id-1
if(trc_name(id)(1:len_trim(trc_name(id))).eq.trc_name(k)(1:len_trim(trc_name(k))))new=0
enddo
endif
if(new.eq.1)then
ntrc=ntrc+1
endif
enddo
print*,ntrc,' different tracers found'
allocate(emittrc(ntrc))
allocate(emityp(ntrc))
allocate(molmass(ntrc))
! Fill the molar mass, type and name of the ntrc tracer species
ntrc=0
do id=1,nline_tracer
new=1
if(id.gt.1)then
do k=1,id-1
if(trc_name(id)(1:len_trim(trc_name(id))).eq.trc_name(k)(1:len_trim(trc_name(k))))new=0
enddo
endif
izoemis(id)=-999
imeemis(id)=-999
if(new.eq.1)then
ntrc=ntrc+1
emittrc(ntrc)=trc_name(id)
emityp(ntrc)=tmptyp(id)
molmass(ntrc)=tmpmolmass(id)
trc_index(id)=ntrc
else
do k=1,ntrc
! If the same tracer species appears on different lines, fill trc_index and check consistency
if(trc_name(id)(1:len_trim(trc_name(id))).eq.emittrc(k))then
trc_index(id)=k
if(emityp(k).ne.tmptyp(id))then
print*,'ERROR: Tracer '//emittrc(k)// ' appears on two different lines of TRACERS.data'
print*,'with two different types. Please check TRACERS.data'
stop 1
endif
if(molmass(k).ne.tmpmolmass(id))then
print*,'ERROR: Tracer '//emittrc(k)// ' appears on two different lines of TRACERS.data'
print*,'with two different molar masses. Please check TRACERS.data'
stop 1
endif
endif
enddo
endif
enddo
do k=1,ntrc
print*,' ',k,' ',emittrc(k)
enddo
! Calculate the emission duration for each line
allocate(emisduration(nline_tracer))
do id=1,nline_tracer
emisduration(id)=interdat(tmpdate1(id),tmpdate2(id))
if(emisduration(id).eq.0)emisduration(id)=1
print*,id,' : emissions during ',emisduration(id),' hours'
enddo
end subroutine ini_temis
!----------------------------------------------------------------------------------------
subroutine process_temis(ks,itimeslot)
! Altitude and layer depth must be here in [cm]
use chimere_common
use pncvar
use interp_tools
use gridvars, only : xlong, xlati, nzonal, nmerid, xlong_entire_domain, xlati_entire_domain
implicit none
integer,intent(in) :: ks,itimeslot
! Hypothesis:
! The tracer molar mass is arbitrarily set to 100 (parameter 'molmass')
integer, parameter :: manydata=1000
integer :: k,k2,id,is
character(len=4)extens
character(len=splen)tmpname
integer :: nv
real(kind=iprec) :: altd,altu
real(kind=iprec), allocatable,dimension(:) :: cumverti,totemis
logical, allocatable,dimension(:) :: isactive
logical :: species_update
real(kind=iprec), allocatable,dimension(:,:) :: vertemis,vertprox
real(kind=iprec), allocatable,dimension(:,:,:) :: bufw
real(kind=iprec), dimension(3) :: mdistrib
integer :: i,ime,izo,ip,ik,iv
logical :: update_emipa ! if particulate tracer emissions are present, emipa needs be updated
! Mostly for findcell
real(kind=iprec), allocatable,dimension(:,:,:,:):: corners_3d_temis
logical :: foundsm
integer :: lstix,lstiy
real(kind=iprec) :: xl,yl,surf,glat,ddx,ddy,x,y
real(kind=iprec), parameter :: degrad=0.0174533d0
real(kind=iprec), parameter :: radius=6371000.d0
if (.not. allocated(vertemis))allocate(vertemis(nline_tracer,nverti))
if (.not. allocated(cumverti))allocate(cumverti(nline_tracer))
if (.not. allocated(isactive))allocate(isactive(nline_tracer)) ! Is this tracer emission line active in the present worker and time
if (.not. allocated(bufw))allocate(bufw(nzonal,nmerid,nverti))
if (.not. allocated(totemis))allocate(totemis(nline_tracer)) ! Total mass emitted in one hour for one tracer line in molecules/cm2/hour for the enrire column
if (.not. allocated(vertprox))allocate(vertprox(nline_tracer,nverti))
if(chimverb.ge.7.and.rank.eq.0)print*,'entering temis',rank
! Initialisation of indices for each line at first pass only
do id=1,nline_tracer
if((izoemis(id).eq.-999).and.(id.eq.1))then
allocate(corners_3d_temis(nzonal,nmerid,4,3))
call calculate_corners(xlong(1:nzonal,1:nmerid),xlati(1:nzonal,1:nmerid),corners_3d_temis,nzonal,nmerid)
endif
if(izoemis(id).eq.-999)then
if(chimverb.ge.10)print*,'trac line',id,'rank',rank
izoemis(id) = -1
imeemis(id) = -1
lstix = 0
lstiy = 0
x=lonemis(id)
y=latemis(id)
call findcell(x,y,izo,ime,nzonal,nmerid,corners_3d_temis,lstix,lstiy,foundsm)
if(foundsm) then
izoemis(id)=izo !izoemis stays equal to -1 if not found
imeemis(id)=ime !imeemis stays equal to -1 if not found
glat=real(xlati(izo,ime))*degrad
ddx=real(xlong_entire_domain(domains%dom(rank+1)%izstart+izo, &
& domains%dom(rank+1)%imstart+ime-1)-xlong(izo,ime))
ddy=real(xlati_entire_domain(domains%dom(rank+1)%izstart+izo-1, &
& domains%dom(rank+1)%imstart+ime)-xlati(izo,ime))
xl=radius*cos(glat)*degrad*ddx
yl=radius*degrad*ddy
surf=xl*yl*1d4
emisfac(id)=1.d6*avogadro/(tmpmolmass(id)*surf*3600d0) ! conversion factor to molecule/cm2/s on 1 hour
endif
if(id.eq.nline_tracer)then
deallocate(corners_3d_temis)
endif
endif
if(chimverb.ge.10)print*,'end trac',id,'rank',rank
enddo
if(chimverb.ge.10)print*,'Exit first loop on id rank',rank
mdistrib(1)=0.3
mdistrib(2)=0.4
mdistrib(3)=0.3
! mdistrib(1)=0.0d0 !WARNING : only good for ash !!
! mdistrib(2)=0.2d0 !WARNING : only good for ash !!
! mdistrib(3)=0.8d0 !WARNING : only good for ash !!
! Calculation of mass for each line at current time step
totemis=0.
vertemis=0. ! emission values for each line and vertical level
vertprox=0. ! proxy for vertical distribution
cumverti=0.
isactive=.false.
if(chimverb.ge.10)print*,'begin second loop on id',rank
do id=1,nline_tracer
if(chimverb.ge.10)print*,'begin loop 2 temis',id,'rank',rank
if (izoemis(id) > 0 .and. imeemis(id) > 0) then
if((idate(itimeslot).ge.tmpdate1(id)).and.(idate(itimeslot).lt.tmpdate2(id)).or.((idate(itimeslot).eq.tmpdate1(id)).and.(tmpdate1(id).eq.tmpdate2(id))))then
if(tmpmass(id).gt.0.)then
isactive(id)=.true.
totemis(id)=(emisfac(id)*tmpmass(id)/emisduration(id)) !in mol/cm2/hour
do nv=1,nverti
if(nv.eq.1)then
altd=0.
altu=hlay(izoemis(id),imeemis(id),nv,ks)
else
altd=hlay(izoemis(id),imeemis(id),nv-1,ks)
altu=hlay(izoemis(id),imeemis(id),nv,ks)
endif
if (intersec(altd,altu,tmpaltmin(id),tmpaltmax(id)).gt.dzero)then
vertprox(id,nv)=intersec(altd,altu,tmpaltmin(id),tmpaltmax(id))
cumverti(id)=cumverti(id)+vertprox(id,nv)
endif
enddo
if(cumverti(id).gt.0.)then
do nv=1,nverti
vertemis(id,nv)=vertemis(id,nv)+totemis(id)*vertprox(id,nv)/cumverti(id) ! Affect the
enddo
else
print*,'WARNING: in prep_temis.F90'
print*,'A tracer emission mass is specified but the altitude range'
print*,tmpaltmin(id),tmpaltmax(id)
print*,'Does not intersect the CHIMERE altitude range',0.,hlay(izoemis(id),imeemis(id),nverti,ks)
print*,'This tracer emission will not be taken into account'
print*,'Check TRACERS.data'
endif
endif
endif
endif
enddo
if(chimverb.ge.10)print*,'temis 2 rank',rank
! Now update emisa (gaseous tracers)
do is=1,nemisa
! First deal with Gaseous species
species_update=.false.
bufw=0.0
do id=1,nline_tracer ! Loop on tracer lines to find all emissions in species is
if (tmptyp(id).eq.1)then
k=trc_index(id)
tmpname=emittrc(k)(1:len_trim(emittrc(k)))
if(emisa_species(is)%name(1:len_trim(emisa_species(is)%name)) &
& .eq.tmpname(1:len_trim(tmpname))) then
species_update=.true.
if(isactive(id))then
do nv=1,nverti
bufw(izoemis(id),imeemis(id),nv)=vertemis(id,nv)
enddo
endif
endif
endif
enddo
if(species_update)then
! Update emisa variable
emisa(is,1:nzonal,1:nmerid,1:nverti,ks)=emisa(is,1:nzonal,1:nmerid,1:nverti,ks) &
& +bufw(1:nzonal,1:nmerid,1:nverti)
if(chimverb.ge.10)print*,'rank ',rank,'begin update AEMISSIONS filefor',tmpname(1:len_trim(tmpname))
call pncvar_put_var(domains,var_emisa(is), &
& emisa(is,1:nzonal,1:nmerid,1:nverti,ks), rank, &
& (/nzonal,nmerid,nverti,1/), (/itimeslot+1/))
if(chimverb.ge.10)print*,'rank ',rank,'finish updating AEMISSIONS file for variable ',tmpname(1:len_trim(tmpname))
endif
if(chimverb.ge.10)print*,'temis 4 rank',rank
! Now deal with Particulate species
species_update=.false.
update_emipa=.false.
! aerosols tracers
bufw(:,:,:)=0.0
do id=1,nline_tracer
k=trc_index(id)
if(tmptyp(id).eq.2)then
do k2=1,3
if(k2.eq.1)extens='_fin'
if(k2.eq.2)extens='_coa'
if(k2.eq.3)extens='_big'
tmpname=emittrc(k)(1:len_trim(emittrc(k)))//extens
if(emisa_species(is)%name(1:len_trim(emisa_species(is)%name)) &
& .eq.tmpname(1:len_trim(tmpname)))then
update_emipa=.true. ! it is needed to update emipa at the end of the treatment of is
if (izoemis(id) > 0 .and. imeemis(id) > 0) then
do nv=1,nverti
bufw(izoemis(k),imeemis(k),nv)=vertemis(id,nv)*mdistrib(k2)
enddo
! Update emisa variable
emisa(is,1:nzonal,1:nmerid,1:nverti,ks)=emisa(is,1:nzonal,1:nmerid,1:nverti,ks) &
& +bufw(1:nzonal,1:nmerid,1:nverti)
end if
call pncvar_put_var(domains,var_emisa(is), &
& emisa(is,1:nzonal,1:nmerid,1:nverti,ks), rank, &
& (/nzonal,nmerid,nverti,1/), (/itimeslot+1/))
endif
enddo
endif
enddo
if(chimverb.ge.10)print*,'temis 5 rank',rank
! update emipa array for particulate emissions, if necessary
if(update_emipa)then
if(chimverb.ge.10)print*,'updating emipa...',rank
ip = iprofemipa(is)
if(ip.ne.0) then
ik = ispecemip(ip)
do i=1,ms
do ime=1,nmerid
do izo=1,nzonal
do iv=1,nverti
emipa(ik,i,izo,ime,iv,ks) = emisa(is,izo,ime,iv,ks)*emiprof(ip,i)
enddo
enddo
enddo
enddo
endif
if(chimverb.ge.10)print*,'end updating emipa...',rank
endif
if(chimverb.ge.10)print*,'temis 6 rank',rank
enddo ! of nemisa
! Synchronize disk to avoid data loss
if(rank.eq.0.and.chimverb.ge.7)print*,'temis begin sync rank',rank
print*,'Running pncvar_sync for tracers in ',emisa_nc%name
call pncvar_sync(emisa_nc%id)
print*,'Successful pncvar_sync for tracers in ',emisa_nc%name
if(rank.eq.0.and.chimverb.ge.7)print*,'temis end sync rank',rank
if((chimverb.ge.5).and.(rank==0))print*,'Leaving temis'
end subroutine process_temis
......@@ -3,7 +3,7 @@
#PBS -q std
#PBS -n
#PBS -l nodes=1:ppn=64
#PBS -l walltime=00:15:00
#PBS -l walltime=00:30:00
#PBS -l mem=240gb
#PBS -l vmem=240gb
......
......@@ -33,19 +33,17 @@ tau_divgrad=18000
guided_type=nudging
itau_nudging=30
# itau_nudging=60
# nudging_zone=circular
nudging_zone=lonlat
nudging_stiffness=8
nudging_lon_start=90
nudging_lon_end=180
nudging_lat_start=20
nudging_lat_end=70
# nudging_lon_center=40
# nudging_lat_center=40
# nudging_radius=3000
#-------------- Physics -------------
#---------------- Run ---------------
run_length=86400
# run_length=86400
run_length=777600
write_period=3600
etat0=dcmip2016_baroclinic_wave
#------------ Diagnostics -----------
......@@ -66,6 +66,7 @@ function cmd_merge()
{
cd $ROOT/NetCDF
date
rm -f merged.nc
cdo merge surf.nc uvtqz.nc merged.nc
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment