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

Extra checks around detection of emission grid point

parent 648f764a
camelot : F90=mpif90
camelot : F90FLAGS=-I ../dynamico/inc -i4 -r8 -auto -align all -assume realloc_lhs
camelot : F90FLAGS=-I ../dynamico/inc -i4 -r8 -auto -align all -assume realloc_lhs -g -traceback -check bounds -fp-model strict
camelot : NCFLAGS=-L/opt/netcdf4/4.4.1.1-parallel/ifort/lib -lnetcdff -Wl,-rpath -Wl,/opt/netcdf4/4.4.1.1-parallel/ifort/lib -lnetcdf
camelot : ICOFLAGS=-L ../dynamico/lib -licosa -L${MKLROOT}/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lm
camelot : XIOSFLAGS=-L ../XIOS/lib -licosa -lxios -lstdc++
camelot : LDFLAGS=$(ICOFLAGS) $(XIOSFLAGS) $(NCFLAGS)
camelot : LDFLAGS=$(ICOFLAGS) $(XIOSFLAGS) $(NCFLAGS)
camelot : all
all :
......
......@@ -9,7 +9,9 @@ MODULE emission
PRIVATE
SAVE
REAL, PARAMETER :: oneday = 86400. ! hard-coded
REAL, PARAMETER :: deg_to_rad=ASIN(1.)/90., &
oneday = 86400. ! hard-coded
INTEGER, PARAMETER :: log_unit = 15
INTEGER :: emission_ij ! index of cell where volcanic emissions occur
......@@ -33,47 +35,67 @@ CONTAINS
unjours = 86400.
CALL getin('unjours', unjours)
CALL init_emission(inout%lon, inout%lat)
CALL init_emission(inout%ngrid, inout%lon, inout%lat, emission_ij)
END SUBROUTINE init_physics
SUBROUTINE init_emission(lon,lat)
SUBROUTINE init_emission(ngrid, lon,lat, ij_min)
USE mpipara, ONLY : mpi_rank
USE icosa, ONLY : rstd
USE getin_mod, ONLY : getin
USE spherical_geom_mod, ONLY : dist_lonlat ! distance on the sphere
REAL(rstd), INTENT(IN) :: lon(:), lat(:)
REAL(rstd) :: emission_lon, emission_lat, dist, dist_min, dist_max
INTEGER :: ij, ij_min
USE spherical_geom_mod, ONLY : dist_lonlat, lonlat2xyz ! distance on the sphere
INTEGER, INTENT(IN) :: ngrid
REAL(rstd), INTENT(IN) :: lon(:), lat(:)
INTEGER, INTENT(OUT) :: ij_min
REAL(rstd) :: emission_lon, emission_lat, dist, dist_min, dist_max, &
xyz(3), emission_xyz(3)
INTEGER :: ij
WRITE(*,*) 'init_emission called'
emission_lon = 0.
CALL getin('emission_lon', emission_lon)
emission_lon = emission_lon * deg_to_rad
emission_lat = 0.
CALL getin('emission_lat', emission_lat)
dist_max = 1e-4
emission_lat = emission_lat * deg_to_rad
dist_max = 2e-2
CALL getin('emission_dist_max', dist_max)
dist_min = 4. ! distances on the unit sphere are always less than pi
ij_min = 1
DO ij=1, SIZE(lon)
CALL dist_lonlat(emission_lon, emission_lat, lon(ij), lat(ij), dist)
CALL lonlat2xyz(emission_lon, emission_lat, emission_xyz)
DO ij=1, ngrid
CALL lonlat2xyz(lon(ij), lat(ij), xyz)
xyz = xyz - emission_xyz
dist = SQRT( SUM(xyz*xyz) )
IF(dist<dist_min) THEN
dist_min = dist
ij_min = ij
END IF
END DO
PRINT *, 'emission lon, lat :', emission_lon, emission_lat
PRINT *, 'closest lon, lat :', lon(ij_min), lat(ij_min)
PRINT *, 'dist to closest cell :', ij, dist_min
IF(mpi_rank==0) THEN
PRINT *, 'emission ALL lon', lon(:)
PRINT *, 'emission ALL lat', lat(:)
END IF
! if the closest point is too far from the target, this means that
! this point is on another MPI process
IF(dist_min < dist_max) THEN
emission_ij = ij_min
IF(dist_min > dist_max) THEN
ij_min = -1
ELSE
! if the closest point is too far from the target, this means that
! this points is on another MPI process
emission_ij = -1
PRINT *, 'emission MIN(lon) MAX(lon)', SIZE(lon), ngrid, MINVAL(lon), MAXVAL(lon)
PRINT *, 'emission MIN(lat) MAX(lat)', SIZE(lat), ngrid, MINVAL(lat), MAXVAL(lat)
PRINT *, 'emission lon, lat :', mpi_rank, emission_lon, emission_lat
PRINT *, 'closest lon, lat :', mpi_rank, lon(ij_min), lat(ij_min)
PRINT *, 'closest dist, ij :', mpi_rank, dist_min, ij_min, ngrid
PRINT *, 'dist to closest cell :', mpi_rank, ij_min, dist_min
END IF
END SUBROUTINE init_emission
SUBROUTINE physics
......@@ -108,8 +130,11 @@ CONTAINS
! & 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
inout%dq(:,:,:)=0.
CALL emit(llm, time, inout%p(emission_ij, :), inout%geopot(emission_ij,:), &
& inout%dq(emission_ij,:,1) )
END IF
......
......@@ -25,6 +25,8 @@ disvert=read_apbp
dt=360
nqtot=5
itau_adv=1
time_scheme=euler
#----------- Dissipation ------------
nitergdiv=1
tau_graddiv=18000
......@@ -46,9 +48,14 @@ nudging_lat_end=70
#-------------- Physics -------------
physics=plugin
itau_physics=1
# in DEGREES
emission_lon = -72.2
emission_lat = -40.5
#---------------- Run ---------------
run_length=86400
# run_length=777600
# run_length=86400
run_length=777600
write_period=3600
etat0=dcmip2016_baroclinic_wave
#------------ Diagnostics -----------
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