Commit 1adf8bde authored by Antoine Berchet's avatar Antoine Berchet

Merge branch 'LSCE' into 'master'

Lsce

See merge request !17
parents 34bf41ea 2efa32f3
Pipeline #62 failed with stages
......@@ -88,7 +88,7 @@ model :
version : std
# H matrix
file_pg : ~/cif/model_sources/dummy_gauss/Pasquill-Gifford.txt
file_pg : ~/CIF/model_sources/dummy_gauss/Pasquill-Gifford.txt
# Chemical scheme
chemistry:
......
......@@ -73,7 +73,7 @@ model :
version : std
# H matrix
file_pg : ~/cif/model_sources/dummy_gauss/Pasquill-Gifford.txt
file_pg : ~/CIF/model_sources/dummy_gauss/Pasquill-Gifford.txt
# Chemical scheme
chemistry:
......
......@@ -8,9 +8,9 @@ LDFLAGS = -L${NETCDFLIB} -lnetcdf -lnetcdff
RM_LIST = *.a *.o *.e *.mod core
OBJS = checkcfl.o \
endchimere.o findspec.o gaus2.o gausbt.o \
endchimere.o findspec.o \
iniboun.o inicdf.o inichem.o inichimere.o \
iniconc.o inidepo.o inidepout.o inidoms.o iniemis.o iniend.o iniequi.o \
iniconc.o inidepo.o inidepout.o inidoms.o iniemis.o iniend.o \
inigeom.o iniout.o iniparam.o iniphys.o \
iniread.o outprint.o readhour.o renewhour.o \
write_depo.o
......
SUBROUTINE GAUS2(F,XL,XU,RELER,ABSER,UROUND,ANSWR,IER,&
EXTRA1,EXTRA2,EXTRA3,NEXTRA,INDIC)
IMPLICIT NONE
!!!
! this old fashioned fortran full of GOTOs does not work with
! EMACS auto indenter !!
! Do not try to auto indent it
!!!
INTEGER, PARAMETER :: N1MAX=21
real(kind=8) :: F,XL,XU,RELER,ABSER,UROUND,ANSWR,EXTRA1,&
EXTRA2,EXTRA3,RATIO,XSUM
INTEGER :: NEXTRA,N1,I,N,NMAX=N1MAX,INDIC,IER
real(kind=8), DIMENSION(N1MAX) :: X,H
real(kind=8), DIMENSION(2,N1MAX) :: A
INTEGER, DIMENSION(N1MAX):: ISIDE
real(kind=8) :: FUN,XD,HD
EXTERNAL F
!
FUN(XD,HD)=0.5*HD*(F(XD+.2113248654052*HD,EXTRA1,&
EXTRA2,EXTRA3,NEXTRA,INDIC)+F(XD+.788675134598*HD,EXTRA1,&
EXTRA2,EXTRA3,NEXTRA,INDIC))
!
H(1)=XU-XL
A(2,1)=FUN(XL,H(1))
IF(IER.NE.1) GOTO 2
IF(10.*ABS(H(1))/RELER.LT.DMAX1(ABS(XU),ABS(XL))) GOTO 7
2 IF(ABS(XU-XL).GT.4.*UROUND*DMAX1(ABS(XL),ABS(XU))) GOTO 8
ANSWR=A(2,1)
IER=-2
RETURN
8 RATIO=DMAX1(ABS(XU/H(1)),ABS(XL/H(1)))
N1=-INT(1.4427*LOG(RATIO*UROUND))
! N1=2-INT(1.4427*LOG(RATIO*UROUND))
NMAX=MIN0(NMAX,N1)
IF(NMAX.GT.1) GOTO 10
IER=-1
RETURN
10 ISIDE(1)=2
DO 1 I=2,NMAX
ISIDE(I)=2
1 H(I)=.5*H(I-1)
X(2)=XL
N=2
4 XSUM=0.
A(1,N)=FUN(X(N),H(N))
A(2,N)=FUN(X(N)+H(N),H(N))
XSUM=A(1,N)+A(2,N)
IF(ABS(XSUM-A(ISIDE(N),N-1))/RELER.LT.ABS(XSUM)+ABSER) GOTO 3
IF(N.EQ.NMAX) GOTO 9
N=N+1
ISIDE(N)=1
X(N)=X(N-1)
GOTO 4
3 A(ISIDE(N),N-1)=XSUM
IF(ISIDE(N).EQ.1) GOTO 5
6 IF(N.EQ.2) GOTO 7
N=N-1
A(ISIDE(N),N-1)=A(1,N)+A(2,N)
IF(ISIDE(N).EQ.2) GOTO 6
5 ISIDE(N)=2
X(N)=X(N-1)+H(N-1)
GOTO 4
9 IER=N-1
XL=X(N)
XU=X(N)+2.*H(N)
RELER=XSUM
ABSER=A(ISIDE(N),N-1)
RETURN
7 IER=0
ANSWR=A(2,1)
RETURN
END
SUBROUTINE GAUSBT(F,XL,XU,RELER,ABSER,UROUND,ANSWR,&
IER,FIXSZ,BASESZ,INNER,TGAS,PGAS)
IMPLICIT NONE
EXTERNAL F
INTEGER, PARAMETER :: N1MAX=21
real(kind=8), DIMENSION(N1MAX) :: X,H
real(kind=8), DIMENSION(2,N1MAX) :: A
INTEGER, DIMENSION(N1MAX) :: ISIDE
real(kind=8) :: F,XL,XU,RELER,ABSER,UROUND,ANSWR,FIXSZ,BASESZ,TGAS,PGAS
real(kind=8) :: FUN,XD,HD,RATIO,XSUM
INTEGER :: IER,INNER,NMAX=N1MAX,N1,I,N
!
FUN(XD,HD)=0.5*HD*(F(XD+.2113248654052*HD,RELER,ABSER,UROUND,&
FIXSZ,BASESZ,INNER,TGAS,PGAS)+&
F(XD+.788675134598*HD,RELER,ABSER,UROUND,&
FIXSZ,BASESZ,INNER,TGAS,PGAS))
!
H(1)=XU-XL
A(2,1)=FUN(XL,H(1))
IF(IER.NE.1) GOTO 2
IF(10.*ABS(H(1))/RELER.LT.DMAX1(ABS(XU),ABS(XL))) GOTO 7
2 IF(ABS(XU-XL).GT.4.*UROUND*DMAX1(ABS(XL),ABS(XU))) GOTO 8
ANSWR=A(2,1)
IER=-2
RETURN
8 RATIO=DMAX1(ABS(XU/H(1)),ABS(XL/H(1)))
! N1=2-INT(1.4427*LOG(RATIO*UROUND))
N1=-INT(1.4427*LOG(RATIO*UROUND))
NMAX=MIN0(NMAX,N1)
IF(NMAX.GT.1) GOTO 10
IER=-1
RETURN
10 ISIDE(1)=2
DO 1 I=2,NMAX
ISIDE(I)=2
1 H(I)=.5*H(I-1)
X(2)=XL
N=2
4 XSUM=0.
A(1,N)=FUN(X(N),H(N))
A(2,N)=FUN(X(N)+H(N),H(N))
XSUM=A(1,N)+A(2,N)
IF(ABS(XSUM-A(ISIDE(N),N-1))/RELER.LT.ABS(XSUM)+ABSER) GOTO 3
IF(N.EQ.NMAX) GOTO 9
N=N+1
ISIDE(N)=1
X(N)=X(N-1)
GOTO 4
3 A(ISIDE(N),N-1)=XSUM
IF(ISIDE(N).EQ.1) GOTO 5
6 IF(N.EQ.2) GOTO 7
N=N-1
A(ISIDE(N),N-1)=A(1,N)+A(2,N)
IF(ISIDE(N).EQ.2) GOTO 6
5 ISIDE(N)=2
X(N)=X(N-1)+H(N-1)
GOTO 4
9 IER=N-1
XL=X(N)
XU=X(N)+2.*H(N)
RELER=XSUM
ABSER=A(ISIDE(N),N-1)
RETURN
7 IER=0
ANSWR=A(2,1)
RETURN
END
......@@ -46,7 +46,7 @@ subroutine inidepo
if(ns.gt.0) then
ndepo = ndepo + 1
indepo(ns) = ndepo
dMx(ndepo) = rMx
dMx(ndepo) = real(rMx)
dHx(ndepo) = rHx
df0(ndepo) = rf0
else
......
subroutine iniequi
!
implicit none
! local variables
integer :: nsweep,ndiv,iacalci,maxit
real(kind=8) :: eps,epsact
eps=1.d-5
maxit=60
nsweep=10
epsact=1.d-2
ndiv=5
iacalci=1
call SETPARM (-1, IACALCI, EPS, MAXIT, NSWEEP, EPSACT, NDIV)
!
return
END subroutine iniequi
......@@ -143,8 +143,6 @@ subroutine outprint(iprint, print_title, nh, writeout, parout)
toprint2d(izo, ime) = topcloc(izo, ime)
case('clwc')
toprint2d(izo, ime) = clwcloc(izo, ime, 1) * an / (airmloc(izo, ime, 1) * 29.)
case('humi')
toprint2d(izo, ime) = 1.e2 * rhloc(izo, ime, 1)
case default
call exit1('outprint.f90 : error in output parameter list')
end select
......
......@@ -245,14 +245,8 @@ subroutine locvalues
do ime = 1, nmerid
do izo = 1, nzonal
! Calculate relative humidity : Teten's formula
! print*,rank,'WW',izo,ime,ivert,airmloc(izo,ime,ivert)
rml = sphuloc(izo, ime, ivert) / airmloc(izo, ime, ivert) / 1.60
psat = 0.611 * 10**(7.5 * (temploc(izo, ime, ivert) - t0k) / (temploc(izo, ime, ivert) - 35.85))
psat = psat * 1000.
rhloc(izo, ime, ivert) = min(rml * presloc(izo, ime, ivert) / (0.62197 + rml) / psat, 1.d0)
clwckgkg = an * clwcloc(izo, ime, ivert) / (airmloc(izo, ime, ivert) * 29d0)
if(clwckgkg.le.1d-6) then
if(clwckgkg.le.1d-6) then ! 1d-6 in unit of clwckgkg [kg/kg]
incloud(izo, ime, ivert) = 0
else
incloud(izo, ime, ivert) = 1
......
......@@ -39,14 +39,14 @@ subroutine rates
real(kind = 8) :: wl, kc, ve
real(kind = 8) :: cw, gama
real(kind = 8) :: rhoa, diff, He, ft
real(kind = 8) :: tho, pe, at, za, diam
real(kind = 8) :: tho, za, diam
real(kind = 8) :: ph, ph1
real(kind = 8), dimension(ntemps) :: dtemp
!*****************************************************************************************
tho = 1. / dtr2
tho = dun / dtr2
diff = 0.25d-4
! First finds NO2 deposition index for HONO formation reaction
......@@ -72,10 +72,8 @@ subroutine rates
! Physical parameters
te = temploc(izo, ime, ivert)
pe = presloc(izo, ime, ivert)
ai = airmloc(izo, ime, ivert)
hu = sphuloc(izo, ime, ivert)
at = atteloc(izo, ime)
cw = clwcloc(izo, ime, ivert)
ph = phloc(izo, ime, ivert)
......@@ -139,7 +137,6 @@ subroutine rates
! Special case of ozone photolysis
else if(ity.eq.13) then
rate(nr, izo, ime, ivert) = rate(nr, izo, ime, ivert) * at
factor = hu / (hu + ai * (0.02909d0 * exp(70d0 / te) + 0.06545d0 * exp(110d0 / te)))
rate(nr, izo, ime, ivert) = phrate(iphoto(nr), izo, ime, ivert) * factor
......@@ -158,7 +155,7 @@ subroutine rates
f2 = tabrate(3, nr) * exp(-tabrate(4, nr) / te)
f3 = tabrate(5, nr) * exp(-tabrate(6, nr) / te)
f4 = tabrate(7, nr) * exp(-tabrate(8, nr) / te)
rate(nr, izo, ime, ivert) = 2d0 * sqrt(f1 * f2 * f3 * f4 / ((1 + f3) * (1 + f4)))
rate(nr, izo, ime, ivert) = 2d0 * sqrt(f1 * f2 * f3 * f4 / ((dun + f3) * (dun + f4)))
else if(ity.eq.9) then
f1 = tabrate(1, nr) * exp(-tabrate(2, nr) / te)
f2 = tabrate(3, nr) * exp(-tabrate(4, nr) / te)
......@@ -166,14 +163,14 @@ subroutine rates
f4 = tabrate(7, nr) * exp(-tabrate(8, nr) / te)
f3 = f3 / (dun + f3)
f4 = f4 / (dun + f4)
rate(nr, izo, ime, ivert) = 2d0 * sqrt(f1 * f2) * (dun - sqrt(f3 * f4)) * (dun - f4) / (2 - f3 - f4)
rate(nr, izo, ime, ivert) = 2d0 * sqrt(f1 * f2) * (dun - sqrt(f3 * f4)) * (dun - f4) / (2d0 - f3 - f4)
else if(ity.eq.16) then
ft = exp(tabrate(3, nr) * (1. / te - 1. / 298.))
if(tabrate(4, nr).eq.0) then
c1 = 10.**(-ph)
c1 = 10d0**(-ph)
c1 = c1**(tabrate(1, nr))
else
c1 = 1. + 13. * (10.**(-ph))
c1 = 1d0 + 13d0 * (10d0**(-ph))
endif
rate(nr, izo, ime, ivert) = 1.d-20
......@@ -182,59 +179,59 @@ subroutine rates
* tabrate(5, nr) * tabrate(6, nr) / c1
endif
else if(ity.eq.23) then
ft = exp(tabrate(4, nr) * (1. / te - 1. / 298.))
ft = exp(tabrate(4, nr) * (1.d0 / te - 1.d0 / 298.d0))
ph1 = min(ph, 5.0)
c1 = 10.**(-ph1)
c1 = 10.d0**(-ph1)
rate(nr, izo, ime, ivert) = ft * cw * tabrate(1, nr) * tabrate(2, nr) * &
tabrate(3, nr) * 8.314 * te * 1.e6 / (1.013e5 * conv) &
tabrate(3, nr) * 8.314d0 * te * 1.d6 / (1.013d5 * conv) &
/ (c1**(tabrate(5, nr)))
else if(ity.eq.24) then
ft = exp(tabrate(4, nr) * (1. / te - 1. / 298.))
ft = exp(tabrate(4, nr) * (1.d0 / te - 1.d0 / 298.d0))
ph1 = min(ph, 6.0)
c1 = 10.**(-ph1)
c1 = 10.d0**(-ph1)
rate(nr, izo, ime, ivert) = ft * cw * tabrate(1, nr) * tabrate(2, nr) * &
tabrate(3, nr) * 8.314 * te * 1.e6 / (1.013e5 * conv) &
tabrate(3, nr) * 8.314d0 * te * 1.d6 / (1.013d5 * conv) &
/ (c1**(tabrate(5, nr)))
else if(ity.eq.20) then
rate(nr, izo, ime, ivert) = 1.d-10
ve = sqrt(8d+3 * 8.314 * te / (pi * tabrate(2, nr)))
ve = sqrt(8d+3 * 8.314d0 * te / (pi * tabrate(2, nr)))
if(tabrate(4, nr).eq.0.) then
gama = tabrate(1, nr)
else
za = (log(1.) - log(tabrate(1, nr))) / (log(273.) - log(298.))
gama = za * (log(te) - log(298.)) + log(tabrate(1, nr))
za = (log(1.d0) - log(tabrate(1, nr))) / (log(273.d0) - log(298.d0))
gama = za * (log(te) - log(298.d0)) + log(tabrate(1, nr))
gama = exp(gama)
if(gama.gt.1.) gama = 1.
if(gama.lt.0.001) gama = 0.001
endif
if(tabrate(5, nr).ne.0.) then ! particles
diam = 10.d-6 ! cloud droplet diameter
c1 = diam / diff / 2.
c2 = 4. / (gama * ve)
rate(nr, izo, ime, ivert) = 1.d-6 * (1.d+8 * cw / (diam * 1.e2)) / (c1 + c2)
c1 = diam / diff / 2.d0
c2 = 4.d0 / (gama * ve)
rate(nr, izo, ime, ivert) = 1.d-6 * (1.d+8 * cw / (diam * 1.d2)) / (c1 + c2)
endif
if(rate(nr, izo, ime, ivert).gt.tho) rate(nr, izo, ime, ivert) = tho
else if(ity.eq.21) then
rhoa = ai * 29. * 1000. / an
rhoa = ai * 29.d0 * 1000.d0 / an
if(tabrate(1, nr).eq.1) then
if(cw.gt.1.d-11) then
wl = an * cw / (ai * 29.)
wl = an * cw / (ai * 29.d0)
else
wl = 0.0
wl = dzero
endif
ve = sqrt(8d+3 * 8.314 * te / (pi * tabrate(3, nr)))
kc = tabrate(2, nr) / 2. / diff + 4. / ve / tabrate(4, nr)
rate(nr, izo, ime, ivert) = 6. * wl * rhoa / (1.d+3 * tabrate(2, nr) * kc)
ve = sqrt(8d+3 * 8.314d0 * te / (pi * tabrate(3, nr)))
kc = tabrate(2, nr) / 2.d0 / diff + 4.d0 / ve / tabrate(4, nr)
rate(nr, izo, ime, ivert) = 6.d0 * wl * rhoa / (1.d+3 * tabrate(2, nr) * kc)
else
print *, '*** RATES: ERROR, ITY=21 and C1 != 1'
stop
endif
if(rate(nr, izo, ime, ivert).gt.tho) rate(nr, izo, ime, ivert) = tho
else if(ity.eq.22) then
He = tabrate(2, nr) * exp(-tabrate(3, nr) * (1. / te - 1. / 298.))
ve = sqrt(8d+3 * 8.314 * te / (pi * tabrate(4, nr)))
kc = tabrate(1, nr) / 2. / diff + 4. / ve / tabrate(5, nr) !
rate(nr, izo, ime, ivert) = 6.d2 / (8.314 * He * te * tabrate(1, nr) * kc)
He = tabrate(2, nr) * exp(-tabrate(3, nr) * (1.d0 / te - 1.d0 / 298.d0))
ve = sqrt(8d+3 * 8.314d0 * te / (pi * tabrate(4, nr)))
kc = tabrate(1, nr) / 2.d0 / diff + 4.d0 / ve / tabrate(5, nr) !
rate(nr, izo, ime, ivert) = 6.d2 / (8.314d0 * He * te * tabrate(1, nr) * kc)
if(rate(nr, izo, ime, ivert).gt.tho) rate(nr, izo, ime, ivert) = tho
else
print *, '*** ERROR: Reaction ', nr, ': Undefined rate type:', ity
......
......@@ -65,14 +65,12 @@ subroutine transmix(ns, izo, ime, ivert, trpr, trlo)
trpr = trpr + vfluxi(izo, ime, ivert, nv) * conc(ns, izo, ime, nv)
endif
enddo
!if(ime==66.and.izo==57.and.ns==32) print 102,'tttttttttttttt ',trpr
102 format(a15,e65.56)
! 3: Transport from resolved vertical velocity
if(inpv.eq.1) then
! UPWIND
!down side
if(ivert.gt.1)then
! if(ime==66.and.izo==57.and.ns==32) print 102,'wwwwwwwwwwwwww ',winvloc(izo,ime,ivert-1)
if(winvloc(izo, ime, ivert - 1).gt.dzero)then
trpr = trpr + winvloc(izo, ime, ivert - 1) / thlayloc(izo, ime, ivert) * conc(ns, izo, ime, ivert - 1)
endif
......@@ -85,7 +83,6 @@ subroutine transmix(ns, izo, ime, ivert, trpr, trlo)
trlo = trlo + winvloc(izo, ime, ivert) / thlayloc(izo, ime, ivert) * conc(ns, izo, ime, ivert)
endif
if(winvloc(izo, ime, ivert).lt.dzero)then
! if(ime==66.and.izo==57.and.ns==32) print 102,'wwwwwwwwww2222 ',winvloc(izo,ime,ivert)
trpr = trpr - winvloc(izo, ime, ivert) / thlayloc(izo, ime, ivert) * conc(ns, izo, ime, ivert + 1)
endif
else
......@@ -180,7 +177,6 @@ subroutine transmix(ns, izo, ime, ivert, trpr, trlo)
print*, 'NOT AN EXISTING SCHEME'
endif
endif !inpv
! if(ime==66.and.izo==57.and.ns==32) print 102,'uuuuuuuuuuuuuu ',trpr
! DEEP CONVECTION IMPLEMENTATION
if(ideepconv.ne.0.and.ideep(izo, ime).eq.1)then
......@@ -264,12 +260,9 @@ subroutine transmix(ns, izo, ime, ivert, trpr, trlo)
endif
endif ! of ideepconv=1
! DEEP CONVECTION IMPLEMENTATION
! if(ime==66.and.izo==57.and.ns==32) print 102,'dddddddddddddd ',trpr
! Horizontal transport
! 1: Western side
! if(ime==66.and.izo==57.and.ns==32) print 101,'hhhhhhhhhhhhhh ',inp,fluxw(izo,ime,ivert),fluxe(izo,ime,ivert),fluxs(izo,ime,ivert),fluxn(izo,ime,ivert)
101 format(a15,i2,4 x(e65.56,' '))
if(uwest(izo, ime, ivert).gt.dzero) then
if(inp.eq.1) then
! PPM
......
......@@ -90,9 +90,7 @@ contains
nitgsef = nitgs
endif
do ni = 1, nitgsef
! print 102,'BBBBBBBBBBBBBB ',conc(32, 57,66, :)
102 format(a15,17 x(e65.56,' '))
do ni = 1, nitgsef
do ivert = 1, nverti
do ime = 1, nmerid
if(dom_i==1) then
......@@ -104,7 +102,6 @@ contains
ratloss = sink / conc(ns, izo, ime, ivert)
conc(ns, izo, ime, ivert) = (cconc(ns, izo, ime, ivert) + dtr2 * source) / &
(dun + dtr2 * ratloss)
! if(ime==66.and.izo==57.and.ns==32) print 102,'bbbbbbbbbbbbbb ',conc(32, 57,66, ivert),source,sink,ratloss,cconc(ns,izo,ime,ivert)
end do
end do
call send_first_end_of_line
......
......@@ -116,8 +116,6 @@ subroutine vtransport
add = hlow / airmloc(izo, ime, ivert + 1)
winvloc(izo, ime, ivert) = fluxbal * add
endif
! if(ime==66.and.izo==57) print 102,'zzzzzzzzzzzzzz ',fluxbal,add,hlow,winvloc(izo,ime,ivert)
!102 format(a15,4 x(e65.56,' '))
! Calculation of mixing fluxes between layer IVERT and IVERT+1
......
......@@ -11,18 +11,23 @@ subroutine wdeposition(ns, izo, ime, ivert, wdlo)
! subroutine arguments
integer, intent(in) :: ns
integer, intent(in) :: izo, ime, ivert
real(kind = 8) :: wdlo
real(kind = 8) :: wdlo,tho
tho=dun/dtr2
wdlo = dzero
! For gases
if(inwetd(ns, 1).ne.0) then
wdlo = conc(ns, izo, ime, ivert) * wetdr1(inwetd(ns, 1), izo, ime, ivert)
! wdlo = conc(ns, izo, ime, ivert) * wetdr1(inwetd(ns, 1), izo, ime, ivert)
! DEBUG version 2017?
wdlo = conc(ns, izo, ime, ivert) * min(wetdr1(inwetd(ns, 1), izo, ime, ivert), tho)
endif
if(inwetd(ns, 2).ne.0) then
wdlo = wdlo + conc(ns, izo, ime, ivert) * wetdr2(inwetd(ns, 2), izo, ime, ivert)
! wdlo = wdlo + conc(ns, izo, ime, ivert) * wetdr2(inwetd(ns, 2), izo, ime, ivert)
! DEBUG version 2017?
wdlo = wdlo + conc(ns, izo, ime, ivert) * min(wetdr2(inwetd(ns, 2), izo, ime, ivert), tho)
endif
wetdepi(ns, izo, ime, ivert) = wdlo * thlayloc(izo, ime, ivert)
......
......@@ -43,7 +43,6 @@ subroutine worker
call prep_outprint(1)
call mpi_barrier(mpi_comm_world, ierr)
if(ad.eq.1) call awrite_concs_hour(1, rank)
!print*,'MMMMMMMMMMMMMM',ntabuzen,nphot,nlevphot
!print*,'Loop on physical time steps'
do nh = 1, nhourrun
......@@ -67,10 +66,9 @@ subroutine worker
call mpi_barrier(mpi_comm_world, ierr)
do ir = 1, ichemstep
call twostep
! if(nh==1.and.np==1)print 102,ir,conc(32, 57,66, :)
call worker_update_halo ! between workers. Master is not involved.
enddo
102 format('AAAAAAAAAAAAA ',i2,' ',17 x(e65.56,' '))
! IP
! RQ: pour le moment, on ne descend pas au niveau
! du pas de temps chimique pour la comparaison des obs
......@@ -119,8 +117,6 @@ subroutine worker
tabobs(i, 10) = ddp !in Pa
tabobs(i, 11) = airmloc(izo, ime, ivert) !in molec.cm-3
tabobs(i, 12) = thlayloc(izo, ime, ivert) !in cm
! if(tabobs(i, 3)==57 .and. tabobs(i,4)==66)print 101,'CCCCCCCCCCCCCC ',ihourrun,np,ivert,conc(ns, izo, ime, ivert),presloc(izo, ime, ivert - 1), presloc(izo, ime, ivert),tabobs(i,11),tabobs(i,12)
!101 format(a15,i2,' ',i2,' ',i2,' ',e65.56,' ',e65.56,' ',e65.56,' ',e65.56,' ',e65.56)
endif ! ihourrun,np,species
enddo ! nobs
......
......@@ -77,6 +77,7 @@ module worker_common
integer :: imstart, imend, izstart, izend ! limites en im et iz par rapport a la grille complete
! fin IP
real(kind = 8), dimension(:, :), allocatable :: xlong ! Longitudes
real(kind = 8), dimension(:, :), allocatable :: xlati ! Latitudes
real(kind = 8), dimension(:, :), allocatable :: clati ! Cosine of latitudes
real(kind = 8), dimension(:, :), allocatable :: slati ! Sine of latitudes
real(kind = 8), dimension(:, :), allocatable :: xsize ! Zonal cell length
......@@ -154,7 +155,6 @@ module worker_common
real(kind = 8), dimension(:, :, :), allocatable :: dtenloc
real(kind = 8), dimension(:, :, :), allocatable :: clwcloc
real(kind = 8), dimension(:, :, :), allocatable :: presloc
real(kind = 8), dimension(:, :, :), allocatable :: rhloc
real(kind = 8), dimension(:, :, :), allocatable :: scinloc
real(kind = 8), dimension(:, :, :), allocatable :: phloc
real(kind = 8), dimension(:, :), allocatable :: hghtloc ! Current mixing height
......@@ -322,6 +322,7 @@ contains
allocate(slati(nzonal, nmerid))
allocate(clati(nzonal, nmerid))
allocate(xlong(nzonal, nmerid))
allocate(xlati(nzonal, nmerid))
allocate(ideep(nzonal, nmerid))
allocate(fveg(nzonal, nmerid, nvegtype, nlduse))
allocate(emisaloc(nemisa, nzonal, nmerid, nlevemis))
......@@ -356,7 +357,6 @@ contains
allocate(winzloc(0:nzonalmax + 1, 0:nmeridmax + 1, nverti))
allocate(winmloc(0:nzonalmax + 1, 0:nmeridmax + 1, nverti))
allocate(presloc(nzonal, nmerid, nverti))
allocate(rhloc(nzonal, nmerid, nverti))
allocate(incloud(nzonal, nmerid, nverti))
allocate(dtenloc(nzonal, nmerid, nverti))
allocate(temploc(nzonal, nmerid, nverti))
......
......@@ -994,6 +994,8 @@ subroutine recv_dbl_arrays
clati(1:nzonal,1:nmerid) = dbuf2(:,:)
call mpi_recv(dbuf2, (nzonal)*(nmerid), mpi_double_precision,0,ias_xlong, mpi_comm_world,mpi_status_ignore,ierr)
xlong(1:nzonal,1:nmerid) = dbuf2(:,:)
call mpi_recv(dbuf2, (nzonal)*(nmerid), mpi_double_precision,0,ias_xlati, mpi_comm_world,mpi_status_ignore,ierr)
xlati(1:nzonal,1:nmerid) = dbuf2(:,:)
deallocate(dbuf2)
allocate(dbuf4(nzonal,nmerid,nvegtype,nlduse))
......@@ -1169,9 +1171,6 @@ subroutine recv_dbl_arrays
dbuf3(:,:,:) = hlayloc(1:nzonal,1:nmerid,1:nverti)
call mpi_send(dbuf3, nzonal*nmerid*nverti, mpi_double_precision,0, iar_hlayloc, mpi_comm_world,ierr)
dbuf3(:,:,:) = rhloc(1:nzonal,1:nmerid,1:nverti)
call mpi_send(dbuf3, nzonal*nmerid*nverti, mpi_double_precision,0, iar_rhloc, mpi_comm_world,ierr)
dbuf3(:,:,:) = presloc(1:nzonal,1:nmerid,1:nverti)
call mpi_send(dbuf3, nzonal*nmerid*nverti, mpi_double_precision,0, iar_presloc, mpi_comm_world,ierr)
......
subroutine zenith
! Calculates current cosine of zenithal angles
! INPUT : SOLTIM Hour lag between years' win. solstice and run sta
! CLATI Cosines of latitudes
! SLATI Sines of latitudes
! IZO Zonal addresses of cells
! IME Meridional addresses of cells
! OUTPUT: ZENILOC C