Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit df428349 authored by Thomas Dubos's avatar Thomas Dubos Committed by Romain Pennel
Browse files

Time and height-dependent emissions

parent 332000f6
......@@ -8,6 +8,17 @@ MODULE emission
IMPLICIT NONE
PRIVATE
SAVE
! Read from run.def
REAL :: emission_time_unit, & ! time unit in seconds for following parameters
& zmin, & ! altitude in m
& zmax ! altitude in m
INTEGER, PARAMETER :: interval_nb=7
REAL, DIMENSION(interval_nb), PARAMETER :: &
& emission_start = (/ 19, 24, 31, 24+19, 48, 48+7, 48+19 /), & ! start of emission
& emission_end = (/ 24, 31, 24+19, 48, 48+7, 48+19, 72 /), & ! end of emission
& amount = (/ 5.21e7, 8.33e7, 5.50e7, 1.0e17, 2.0e7, 3.0e7, 1.0e7 /) ! in kg
REAL, PARAMETER :: deg_to_rad=ASIN(1.)/90., &
oneday = 86400. ! hard-coded
......@@ -64,6 +75,14 @@ CONTAINS
dist_max = 2e-2
CALL getin('emission_dist_max', dist_max)
emission_time_unit = 3600.
CALL getin('emission_time_unit', emission_time_unit)
zmin = 8000.
CALL getin('emission_zmin', zmin)
zmax = 11000.
CALL getin('emission_zmax', zmax)
dist_min = 4. ! distances on the unit sphere are always less than pi
ij_min = 1
CALL lonlat2xyz(emission_lon, emission_lat, emission_xyz)
......@@ -100,9 +119,10 @@ CONTAINS
SUBROUTINE physics
USE mpipara, ONLY : is_mpi_master
USE icosa, ONLY : llm
USE icosa, ONLY : llm, g
USE physics_interface_mod, ONLY : inout => physics_inout
REAL :: dps(inout%ngrid), play(inout%ngrid, llm), pphi(inout%ngrid, llm)
REAL :: pg(llm+1), z(llm+1), mass(llm), dQ(llm)
REAL :: timestep, time, jourvrai, gmtime
INTEGER :: l
IF(is_mpi_master) WRITE(log_unit,*) 'emission called', SHAPE(inout%p), SHAPE(inout%pk)
......@@ -122,34 +142,65 @@ CONTAINS
pphi(:,l) = pphi(:,l) - inout%geopot(:,1)
END DO
! ! go
! CALL emission(inout%ngrid,llm, &
! & firstcall,lastcall, &
! & jourvrai, gmtime, timestep, &
! & inout%p, play, pphi, &
! & inout%ulon, inout%ulat, inout%temp, &
! & inout%dulon, inout%dulat, inout%dtemp, dps)
inout%dtemp(:,:)=0.
inout%dulon(:,:)=0.
inout%dulat(:,:)=0.
inout%dq(:,:,:)=0.
IF( emission_ij>0 ) THEN
CALL emit(llm, time, inout%p(emission_ij, :), inout%geopot(emission_ij,:), &
& inout%dq(emission_ij,:,1) )
pg(:) = inout%p(emission_ij, :) / g
z(:) = inout%geopot(emission_ij,:) / g
! dQ in kg per level per time_unit
CALL emit(llm, time/emission_time_unit, pg, z, mass, dQ)
! mass is in kg of air per m^2
! we must divide dQ by cell area to obtain kg/m^2 per time_unit
! dq in kg/kg per level per second
dQ(:) = dQ(:)/mass(:)/emission_time_unit/inout%Ai(emission_ij)
inout%dq(emission_ij, :, 1) = dQ(:)
PRINT *, 'emission dQ in kg/kg/s', dQ(:)
END IF
END SUBROUTINE physics
SUBROUTINE emit(llm, time, p, geopot, dq)
FUNCTION intersec(a,b,c,d) ! intersection of [a b] and [c d]
REAL, INTENT(IN) :: a,b,c,d
REAL :: intersec
intersec = MIN(b,d)-MAX(a,c)
intersec = MAX(0., intersec)
END FUNCTION intersec
SUBROUTINE emit(llm, time, pg, z, mass, dQ)
USE icosa, ONLY : rstd
INTEGER, INTENT(IN) :: llm ! number of model layers
REAL(rstd), INTENT(IN) :: time ! current time in seconds since start of simulation
REAL(rstd), INTENT(IN) :: p(:), & ! pressure at interfaces
& geopot(:) ! geopotential at interfaces
REAL(rstd), INTENT(OUT) :: dq(:)
! uniform emission on the whole column
dq(:) = 1.
INTEGER, INTENT(IN) :: llm ! number of model layers
REAL(rstd), INTENT(IN) :: time ! current time in emission_time_unit since start of simulation
REAL(rstd), INTENT(IN) :: pg(llm+1), & ! pressure/g at interfaces
& z(llm+1) ! altitude above ground at interfaces
REAL(rstd), INTENT(OUT) :: mass(llm) ! mass in each full level in kg/m^2
REAL(rstd), INTENT(OUT) :: dQ(llm) ! added tracer mass in each full level in kg/emission_time_unit
INTEGER :: l, n
REAL :: rate
! compute mass in each full level
DO l=1, llm
mass(l)=pg(l)-pg(l+1)
END DO
IF(.FALSE.) THEN
! uniform emission on the whole column
dq(:) = mass(:)
ELSE
dQ(:)=0.
DO n=1, interval_nb
! emission rate in kg/m/time_unit
rate = amount(n) / (zmax-zmin) / (emission_end(n) - emission_start(n))
IF(time>=emission_start(n) .AND. time<emission_end(n)) THEN
DO l=1, llm
dQ(l) = dQ(l) + rate*intersec(z(l), z(l+1), zmin, zmax)
! PRINT *, 'emission : level, bottom, top, intersec ', l, z(l), z(l+1), dQ(l)
END DO
END IF
END DO
END IF
PRINT *, 'emission dQ(:)' , dQ
END SUBROUTINE emit
SUBROUTINE compute_play(ngrid, llm, plev, play)
......
......@@ -12,7 +12,12 @@ cd "$PBS_O_WORKDIR"
module load git python/3.6-anaconda50
./scripts/netcdf.sh unpack
# we now have symlinks Grib/*.zip instead of a big tarfile
# => unzip instead of unpack
# ./scripts/netcdf.sh unpack
./scripts/netcdf.sh unzip
./scripts/netcdf.sh extract
./scripts/netcdf.sh create
./scripts/netcdf.sh merge
......
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