!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
***s/p radir7
*
#include "phy_macros_f.h"
subroutine radir7 (taux, fd, fm, fnuage, t, q, sh, ts, ps, nn, 1,3
% nk, oz, nkmax, ni2, nl,ni1, nt, mx, mxx, no3,
$ ncx, nco2,g1,g2,g3,th2o,tro3,yg3,bcn,dbcn,bo3,dbo3,
% to3, uu, tt, f, ku, iopt, opnua,
% xin, v, su, itemp, xpro, s, suo3, del, au, az,
% nhaut, nmoy, nbas,
$ qco2,uco2j, tco2j, wkco2j,trmin,tmem,
% ozpak1,ozpak2,
x reduc,flsh,flss,fldel,flnk,flnn,RADFIX,RADFLTR)
*
#include "impnone.cdk"
*
c variables a une dimension de travail
integer nn,nkmax,ni2,nl,ni1,nt,mx,i1,i2,m,ny,j,jind,
x indx1,indx2,mxx,no3,ncx,iopt,nk,mc,ip,i,l,it,indx,ny2
integer iwing,icen,nco2
c variables dimensionnees
real to,t1,t2,ap,po,elsa,epsil,corfac,eo3,eh2o,
x borne,z,aa,bb,w,xnu,rovlap,xmaxlap,bz,a,tx,wc,zx,
$ r,y,wc2,qco2
c
real tt(nt),tro3(mxx),to3(no3),g1(mxx,nt),g2(mxx,nt),
x g3(mxx,nt),th2o(mxx,nco2),bcn(nt,nco2),dbcn(nt,nco2),
y uu(mxx),yg3(mxx,ncx,nco2),bo3(nt),dbo3(nt)
c
real trnuage,diftlim
real sh(ni2,nk)
real taumax,secday
real voigth2o,voigto3
integer flnk, flnn
c variables dependant de ni2 ou ni1
real taux(ni2,flnk),fnuage(ni2,nk),t(ni1,nk),q(ni1,nk),
1 fm(ni2,flnn),fd(ni2,flnn),oz(ni2,nk),
$ ts(nl),ps(nl), f(ni2,nkmax,nkmax)
integer ku(ni2,nkmax,nkmax) ,itemp(ni2,nn)
real xin(ni2,nn), v(ni2,nn), su(ni2,nn),
2 xpro(ni2,nn), s(ni2,nn), trmin(nl), tmem(nl),
3 suo3(ni2,nn), del(ni2,nk) ,au(nl), az(nl),
4 nhaut(nl), nmoy(nl), nbas(nl)
real uco2j(ni2,nkmax,2),tco2j(ni2,nkmax,nkmax),
1 wkco2j(ni2,nkmax,3)
*
real*8 ozpak1,ozpak2
c
*
logical opnua
logical reduc,radfix,radfltr
real flsh(ni2,flnk), flss(ni2,flnn), fldel(ni2,flnk)
*
*Author
*
* l.garand (june 1989)
*
*Revision
*
* 001 g.pellerin - (march 1990) - standard documentation
* 002 louis garand - add co2 wing bands
* 003 g.pellerin (july 1990)
* standardization of thermodynamic functions
* 004 n. brunet - (may 1991) -
* new version of thermodynamic functions
* and file of constants
* 005 b. bilodeau - (august 1991) - adaptation to unix
* 006 y. chartier - (march 1993) - optimization
* 007 b. bilodeau and c. girard - (march 1993) -
* cooling rates changed at the top of the model
* emissivities set equal to 1
* 008 r. benoit (aug 93) local sigma
* 009 B. Bilodeau (April 95) - w3264 initialized in phydebu
* 010 M. Gagnon (June 1995) - Reduction mode
* 011 L. Garand (April 1995) Temperature effect on CO2 absorption
* (see vco2info) . Effective cloud amount (CF * Emissivity) passed
* to routine (see cldoptx routine, undoes correction 007)
* correction to cloud overlap assumption
* 012 L. Garand (April 1996) Transition from Lorentz to Voigt line
* shape following Giorgetta and Morcrette, MWR 1995, p. 3381-3383.
* 013 B. Dugas (Sep 96) - Implanter RADFIX
* 014 L. Garand (June 1997): extended radiation tables;
* remove packing of tables (code can be exported)
* 015 P.-A. Michelangeli (March 1998): Add radiative code of
* Fomichev and Blanchet (new FOMIC option)
* 016 B. Bilodeau (Jan 2001) - Automatic arrays
* 017 B. Dugas (Jan 2002) - Move all of the Fomichev code to fomichev.ftn
* 018 B. Dugas (Mar 2002) - In reduc mode, no smoothing
* 019 B. Dugas (Sep 2002) - Add background CO2 concentration argument QCO2
* 020 B. Bilodeau, J. Mailhot and L. Garand (Feb 2003) -
* NHAUT, NMOY and NBAS diagnostics
* 021 M. Lepine (March 2003) - CVMG... Replacements
* 022 B. Bilodeau (May 2003) - IBM conversion
* - calls to vsexp routine (from massvp4 library)
* - calls to exponen4 routine to calculate exponentiations
* - removal of useless exponentiations
* - iopt option removed; allowed merge of loops 111 and 112
* - invert dimension of some radiation tables in order to
* reduce the cache flooding (see litblrad2)
* 023 A-M. LEDUC (JUN 2003) - Add RADFLTR switch (radir6 to radir7)
* 024 F. Lemay (Sept 2003) - if radfix, vertical coordinate independent
* of surface pressure.
* 025 B. Bilodeau (Aug 2003) - exponen4 replaced by vspown1
*
*Object
* to calculate infra-red radiation rate of cooling and flux
* according to garand jas 1983
*
*
*Arguments
*
* - output -
* taux rate of cooling/warming in k/s
* fd downward flux at the surface in w/m**2
* fm upward flux at the surface in w/m**2
*
* - input -
* fnuage cloud amount in each layer (0.0 - 1.0) times emissivity
* t temperature in kelvin
* q specific humidity (kg/kg)
* sh sigma levels for t, q, oz
* ts surface temperature in kelvins
* ps surface pressure (n/m**2)
* nn number of flux levels (nk+1)
* nk number of layers
* oz o3 mixing ratio in kg/kg
* nkmax maximum number of levels (set in main program)
* ni2 maximum number of profiles to process
* nl actual number of profiles to process
* ni1 1st dimension of t and q
* nt dimension of table for temperature
* mx number of values of u by decade in the table
* mxx dimension of u in the table
* no3 dimension of ozone table: g1, g2, g3, th2o, bcn, dbcn, uu,
* yg3, bo3, dbo3, yg3o3, tro3, to2, uu: table precalculated
* in fg123
* ncx dimension of co2 for yg3
* nco2 (parameter nco2); (pointer from subroutine pntg123)
* g1 5+1; (pointer from subroutine pntg123)
* g2 g1 + mxx*nt (pointer)
* g3 g2 + mxx*nt (pointer)
* th2o g3 + mxx*nt (pointer)
* tro3 th2o + mxx*nco2 (pointer)
* yg3 tro3 + mxx (pointer)
* bcn yg3 + nco2*mxx*ncx (pointer)
* dbcn bcn + nco2*nt (pointer)
* bo3 dbcn+ nt*nco2; (pointer)
* dbo3 bo3 + nt (pointer)
* to3 dbo3 + nt (pointer)
* uu to3 + no3 (pointer)
* tt uu + mxx (pointer)
* f work field: (2 parts)
* upper triangle contains cloud transmission
* lower triangle contains h2o-co2 transmission overlap
* ku work field: (2 parts)
* upper triangle contains ozone table indices
* lower triangle contains h2o table indices
* iopt not used
* opnua .true. for cloud transmissivity to be .5 to .95
* (normal mode of operation, but unavailable in actual code)
* .false. for cloud transmissivity to be 100% (clear sky)
* (for research purposes)
* xin work field
* v work field
* su work field
* itemp work field
* xpro work field
* suo3 work field 2d for local sigma
* au work field
* reduc .true. to use interpolation
* .false. means we are working on full levels
* radfix .true. to use simple curve fits instead of full
* physics, from the lower stratosphere upwards
* radfltr .true. to apply smoothing of net fluxes
* s sigma flux levels (may be reduced)
* flsh full sigma levels
* flss full sigma flux levels
* del thickness between flux levels (may be reduced)
* fldel full ...
* flnk full number of levels excluding ground
* flnn full number of levels including ground
*
* - input
* qco2 background CO2 concentration in ppmv
*
* - output -
* uco2j amount of co2 in each layer in kg/m**2
* tco2j transmission of co2
*
* - input -
* wkco2j work field
* trmin work field
* tmem work field
*
* -input -
* ozpak1 not used
* ozpak2 not used
*
* - output -
* az work field; also used as output for zonal diagnostics
* of total clouds (variable nt)
* nhaut integrated high-level clouds
* nmoy " mid- " "
* nbas " low- " "
*
*notes
* su and fm can share the same space
* xpro and fd " " " " "
*
* uco2 amount of co2 in each layer in kg/m**2
* tco2 precalculated transmission of co2 according to
* exp(-tco2*ps) which gives the matrix of co2 transmission.
* uco2 and tco2 are calculated in co2info. upper triangle
* for co2 central band. lower triangle for wings
* (average of the two wings)
*
*
**
*
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
*
c
c
c limites de temperature
parameter (t1=180., t2=320., to=260.)
c
c dependence en temperature de l'absorbant h2o
parameter (ap=0.021)
c
c autres constantes physiques
parameter ( po=101325.)
parameter (elsa=1.66)
c
c
parameter (epsil=1.e-12)
parameter(corfac=600.)
c corfac modifie la temperature au premier niveau de flux
c pour le flux descendant pour tenir compte du toit trop bas
c en principe a modifier pour une discredisation choisie
c
! Beware of eo3 and eh2o! If values are changed, code
! must be modified as well because exponentiations have
! been eliminated
parameter (eo3=1.00)
parameter (eh2o=1.00)
parameter (voigth2o=30., voigto3=400.)
c exposants pour calcul des masses d'ozone et d'h2o
parameter (diftlim=8.)
c diftlim limite ts a + ou - diftlim de t(nn) dans radir
c ajoutant une robustesse dans le cas de fortes erreurs de ts
c
#include "consphy.cdk"
*
*
real stefinv
real rwngdbcn,rcendbcn,rwngbcn,rcenbcn,rg1,rg2,rg3
real rtro3,rdbo3,rbo3,rth2o
integer idco2,ido3,idwat,ncen,nc,nwng
integer i_tempo
real xlco2,xlo3,xlwat,aco2wng,aco2cen
real rido3mx,ridwatmx
real ln10inv, pogravinv
real auzt(nl)
pointer (auzt_, auzt)
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
AUTOMATIC ( flfm , REAL , (ni2,flnn))
AUTOMATIC ( flfd , REAL , (ni2,flnn))
AUTOMATIC ( ih , INTEGER , (ni2 ))
AUTOMATIC ( ib , INTEGER , (ni2 ))
AUTOMATIC ( azzt , real , (nl*2 ))
AUTOMATIC ( ps2pograv , real , (nl ))
AUTOMATIC ( sudelq , real , (nl ,nn ))
AUTOMATIC ( sudeloz , real , (nl ,nn ))
AUTOMATIC ( z1d , real , (nl ))
AUTOMATIC ( z1d4 , real , (nl ))
AUTOMATIC ( ts4 , real , (nl ))
*
************************************************************************
external intchamps, vco2inf2
*
*
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
* azzt et auzt doivent etre contigus en memoire
* (voir appel a vslog apres la boucle 111)
auzt_ = loc(azzt(nl+1))
* normalisation par stefan
stefinv = 1./stefan
*
ln10inv=1.0/log(10.)
pogravinv=1.0/po/grav
icen=2
iwing=1
* bande centrale =2 et aile=1 pour l'indice des tableaux co2
trnuage=0.
if(.not. opnua) trnuage=1.
c transmissivite des nuages sera 100% (=clair) si opnua est false
c typical clouds have 97% emissivity
c clouds colder than 220k or above 100mb have 0.50 emissivity
c
c matrix of co2 absorber amount starts at 10**-idco2
c similarly for ido3,idwat
c log10 of absorber amount starts at xlco2,xl03,xlwat
c caution idco2,ido3,idwat therefore needs to be changed if tables
c are changed in there beginning decade
c below are values for original tables (IRTAB4)
idco2=5
ido3=7
idwat=6
if(mxx.eq.501)then
c below are values for extended tables
idco2=9
ido3= 10
idwat= 7
endif
xlco2=-float(idco2) +1.e-6
xlo3= -float(ido3) + 1.e-6
xlwat=-float(idwat)+ 1.e-6
borne=uu(mxx)-1.e-6
rido3mx =(ido3*mx +1 +0.5)
ridwatmx=(idwat*mx +1 +0.5)
c
call vco2inf2(uco2j,tco2j,nl,nn,nk,nl,ni1,nkmax,sh,t,ps,
$ wkco2j(1,1,1),wkco2j(1,1,2),wkco2j(1,1,3),qco2)
do 100 l=1,nl
! Beware of eo3 and eh2o! If values are changed,
! the following code must be modified.
! suo3(l,1)=(s(l,1)+voigto3/ps(l))**eo3
suo3(l,1)=s(l,1)+voigto3/ps(l)
suo3(l,nn)=1.+voigto3/ps(l)
! z=(s(l,1)+voigto3/ps(l))**eh2o
z=s(l,1)+voigto3/ps(l)
v(l,1)=t(l,1) - (sh(l,1)-s(l,1))/(sh(l,2)-sh(l,1))*
$ (t(l,2)-t(l,1))
v(l,nn)=ts(l)
c discontinuite ts et tair non tenue en compte
c mais ts limite par diftlim
c v(l,nn)=amax1(ts(l),t(l,nk)-diftlim)
c v(l,nn)=amin1(ts(l),t(l,nk)+diftlim)
su(l,1)= exp(ap*(v(l,1)-to)) *z
su(l,nn)= exp(ap*(ts(l)-to))
100 continue
c
do 20 i=2,nk
do 120 l=1,nl
! Beware of eo3 and eh2o! If values are changed,
! the following code must be modified.
! z=(s(l,i)+voigth2o/ps(l))**eh2o
! suo3(l,i)=(s(l,i)+voigto3/ps(l))**eo3
z=s(l,i)+voigth2o/ps(l)
suo3(l,i)=s(l,i)+voigto3/ps(l)
c
v(l,i)=(t(l,i)+t(l,i-1))*0.5
su(l,i)= exp(ap*(v(l,i)-to)) *z
sudelq (l,i)=(su(l,i)+su(l,i-1))*del(l,i-1)*0.5 *q(l,i-1)
sudeloz(l,i)=(suo3(l,i)+suo3(l,i-1))*del(l,i-1)*0.5 *oz(l,i-1)
120 continue
c
20 continue
do l=1,nl
sudelq (l,nn)=(su(l,nn)+su(l,nn-1))*del(l,nn-1)*0.5 *q(l,nn-1)
sudeloz(l,nn)=(suo3(l,nn)+suo3(l,nn-1))*del(l,nn-1)*0.5 *oz(l,nn-1)
ps2pograv(l)=ps(l)*ps(l)*pogravinv*elsa
enddo
aa=(nt-1)/(t2-t1)
bb=1.-aa*t1
c
do 1 i=1,nn
c
do 101 l=1,nl
ku(l,i,i)=1
! it=max0(1,nint(aa*v(l,i)+bb))
i_tempo=aa*v(l,i)+bb+0.5
itemp(l,i)=min(max(1,i_tempo),nt)
f(l,i,i)=1.
tmem(l)=1.
trmin(l)=1.
au(l)=0.
az(l)=0.
c maintenant v(l,i) remplace par v(l,i)**4
! v(l,i)=v(l,i)**4
101 continue
call vspown1 (v(1,i),v(1,i),4.,nl)
ip=i+1
if(i.eq.nn)go to 40
c
do 11 j=ip,nn
* indx1=int(0.1 / (s(j-1)+epsil) ) + 1
jind=j-2
jind=max0(jind,1)
c
do 111 l=1,nl
*
c transmission through cloud layer xnu
xnu=1.-fnuage(l,j-1)
if(fnuage(l,jind).lt.0.01)then
c entering a new cloud isolated from upper one
c keep in memory transmission down to top of new cloud tmem
c trmin is minimum transmission in cloud layer processed
c basic idea is random overlap (hence tmem X xnu = cloud transmittance)
c for cloud layers separated by clear ones; but maximum overlap
c (hence cloud transmittance is min (tmem,xnu))for adjacent cloud layers.
tmem(l)= f(l,i,j-1)
trmin(l)= xnu
else
c inside a cloud use maximum overlap
c compute minimum transmission between adjacent layers
trmin(l)=min(trmin(l),xnu)
endif
*
f(l,i,j)= tmem(l) * trmin(l)
c elimine nuages si opnua false
f(l,i,j)=max(f(l,i,j),trnuage)
c
au(l)=au(l)+sudelq (l,j)
az(l)=az(l)+sudeloz(l,j)
azzt(l)=(az(l)*ps2pograv(l)+epsil)
auzt(l)=(au(l)*ps2pograv(l)+epsil)
111 continue
call vslog(azzt,azzt,nl*2)
do l=1,nl
bz=azzt(l)*ln10inv
bz=amax1(bz,xlo3)
ku(l,i,j)=(mx*bz +rido3mx)
! ku(l,i,j)=nint(mx*bz +itmp)
c pour ozone indice 3 superieur a indice 2
a=auzt(l)*ln10inv
a=amax1(a,uu(1))
a=amin1(a,borne)
m=(mx*a +ridwatmx)
ku(l,j,i)=m
c pour transmission h2o-co2 indice 2 superieur a indice 3
f(l,j,i)=tco2j(l,i,j) * th2o(m,icen)
enddo
11 continue
c
1 continue
c
40 continue
c FLUX MONTANT *** FLUX UP
c
call vspown1 (ts4,ts,4.,nl)
DO 130 L=1,NL
! FM(L,NN)=STEFAN*TS(L)**4
FM(L,NN)=STEFAN*TS4(L)
FD(L,1)=0.
! IT=MAX0(1,NINT(AA*TS(L)+BB))
I_TEMPO=AA*TS(L)+BB+0.5
IT=MIN(MAX(1,I_TEMPO),NT)
AZ(L)=FLOAT(IT)
130 CONTINUE
C
I1=1
DO 2 I=I1,NK
C CALCUL DE L'INTEGRALE PROCHE.
C
DO 102 L=1,NL
aco2wng=uco2j(l,i,1)
aco2cen=uco2j(l,i,2)
TX=T(L,I)
! IT=MAX0(1,NINT(AA*TX+BB))
I_TEMPO=AA*TX+BB+0.5
IT=MIN(MAX(1,I_TEMPO),NT)
WC=ALOG10(ACO2WNG+EPSIL)
WC=AMAX1(WC,xlco2)
WC=AMIN1(WC,.9999)
WC2=ALOG10(ACO2CEN+EPSIL)
WC2=AMAX1(WC2,xlco2)
WC2=AMIN1(WC2,.9999)
M=KU(L,I+1,I)
! NWNG=NINT(MX*WC+ idco2*MX+1)
NWNG=MX*WC+ idco2*MX+1.5
! NCEN=NINT(MX*WC2+idco2*MX+1)
NCEN=MX*WC2+idco2*MX+1.5
A=DBCN(IT,ICEN) * YG3(M,NCEN,ICEN)
Z=DBCN(IT,IWING)*YG3(M,NWNG,IWING)
ZX=G3(M,IT)+A +DBO3(IT)*TO3(KU(L,I,I+1))*TRO3(M)+Z
XIN(L,I)=(V(L,I)-V(L,I+1))*ZX*F(L,I,I+1)
XPRO(L,I)=-XIN(L,I)
102 CONTINUE
C
C
IF(I.EQ.NK)GO TO 50
C
IP=I+1
DO 104 L=1,NL
M=KU(L,IP,I)
IT=ITEMP(L,IP)
A=DBCN(IT,ICEN)
Z=DBCN(IT,IWING)*TH2O(M,IWING) *TCO2J(L,IP,I)
AU(L)=G2(M,IT)+A*F(L,IP,I) +Z
X +DBO3(IT)*TRO3(M)*TO3(KU(L,I,IP))
104 CONTINUE
C
C CALCUL DE L'INTEGRALE LOINTAINE.
C
DO 3 J=IP,NK
C
DO 103 L=1,NL
R=AU(L)
IT=ITEMP(L,J+1)
M=KU(L,J+1,I)
A=DBCN(IT,ICEN)
Z=DBCN(IT,IWING)*TH2O(M,IWING)*TCO2J(L,J+1,I)
AU(L)=G2(M,IT) +A *F(L,J+1,I) +Z
X +DBO3(IT)*TRO3(M)*TO3(KU(L,I,J+1))
XIN(L,I)=XIN(L,I)+.5*(R+AU(L))*(V(L,J)-V(L,J+1))
x *F(L,I,J+1)
103 CONTINUE
C
3 CONTINUE
C
50 CONTINUE
C
DO 140 L=1,NL
FM(L,I)=STEFAN* V(L,I) - XIN(L,I)
140 CONTINUE
C
2 CONTINUE
C
C
C FLUX DESCENDANT
C
if (radfix) then
C CORRECTION TEMPERATURE AU TOIT
do L=1,NL
Z1D(L)= T(L,1)+CORFAC*S(L,1)
end do
call vspown1(z1d4,z1d,4.,nl)
DO 125 L=1,NL
I_TEMPO=AA*Z1D(L)+BB+0.5
IT=MIN(MAX(1,I_TEMPO),NT)
ITEMP(L,1)=IT
xpro(l,1)=xpro(l,1)*(z1d4(L)-v(l,2))/
$ sign(max(abs(v(l,1)-v(l,2)),epsil),v(l,1)-v(l,2))
* XPRO(L,1)=XPRO(L,1)*(Z**4-V(L,2))/(V(L,1)-V(L,2)+EPSIL)
V(L,1)=Z1D4(L)
125 CONTINUE
endif
C
C
I2=2
DO 55 I=I2,NN
C
DO 105 L=1,NL
XIN(L,I)=-XPRO(L,I-1)
105 CONTINUE
55 CONTINUE
C
DO 5 I=I2,NN
IF(I.EQ.I2)GO TO 60
C
DO 106 L=1,NL
IT=ITEMP(L,1)
M=KU(L,I,1)
A=DBCN(IT,ICEN)
Z=DBCN(IT,IWING)*TH2O(M,IWING) *TCO2J(L,I,1)
AU(L)=G2(M,IT)+A*F(L,I,1) +Z
X +DBO3(IT)* TRO3(M) *TO3(KU(L,1,I))
106 CONTINUE
C
DO 9 J=1,I-2
C
DO 109 L=1,NL
R=AU(L)
IT=ITEMP(L,J+1)
M=KU(L,I,J+1)
A=DBCN(IT,ICEN)
z=0.
if(tco2j(l,i,j+1).gt.1.e-4)then
Z=DBCN(IT,IWING)*TH2O(M,IWING)*TCO2J(L,I,J+1)
endif
AU(L)=G2(M,IT) +A*F(L,I,J+1) +Z
X +DBO3(IT) *TRO3(M) * TO3(KU(L,J+1,I))
XIN(L,I)=XIN(L,I)+.5*(R+AU(L))*(V(L,J)-V(L,J+1))
x *F(L,J+1,I)
109 CONTINUE
C
9 CONTINUE
C
60 CONTINUE
C
C CALCUL DU TERME DE REFROIDISSEMENT VERS L'ESPACE.
C
DO 160 L=1,NL
IT=ITEMP(L,1)
M=KU(L,I,1)
A=BCN(IT,ICEN)
Z=BCN(IT,IWING)*TH2O(M,IWING) *TCO2J(L,I,1)
Y=G1(M,IT) +A*F(L,I,1) + Z
1 +BO3(IT)*TRO3(M)*TO3(KU(L,1,I))
FD(L,I)=XIN(L,I)+ STEFAN*V(L,I) -V(L,1) *Y *F(L,1,I)
160 CONTINUE
C
5 CONTINUE
c
c
ccc if(radfix) then
c ajustement au toit fd(1) est extrapole et non zero
do 175 l=1,nl
fd(l,1)=amax1(0.,fd(l,2)-(fd(l,3)-fd(l,2))*del(l,1)/del(l,2))
fd(l,1)=amin1(fd(l,1),fd(l,2))
175 continue
ccc else
ccc do l=1,nl
* Utiliser la moyenne des deux approches - BD, 26 novembre 1993.
ccc fd(l,1)=amax1( 0.,
ccc + fd(l,2)-( fd(l,3)-fd(l,2) )*del(l,1)/del(l,2)
ccc + )
ccc + / 2.0
ccc fd(l,1)=amin1(fd(l,1),fd(l,2))
* Ou bien, mettre le flux descendant a zero - BD, 4 juin 1997.
ccc fd(l,1)=0.
ccc enddo
ccc endif
c
c recherche du jump dans transmission h2o (indice de tableau =1)
c elimination par interpolation du flux vers le bas a ce niveau
c ou les indices de tableau passent de 1 a > 1.
c ceci se produit entre 5 et 10 mb.
do 671 i=3,nk
do 672 l=1,nl
if(ku(l,i,1) .gt.1. .and. ku(l,i-1,1).eq.1)
+ fd(l,i)=fd(l,i-1)+(fd(l,i+1)-fd(l,i-1))/(s(l,i+1)-s(l,i-1))*
+ (s(l,i)-s(l,i-1))
672 continue
671 continue
c
c
* calcul des indices IH et IB pour nuages 2-D
* IH = niveau le plus pres de sigma=0.4
* IB = niveau le plus pres de sigma=0.7
do j=1,nk
do l=1,nl
if (sh(l,j).le.0.4) ih(l) = j
if (sh(l,j).le.0.7) ib(l) = j
end do
end do
*
do 208 l=1,nl
c nuages totaux
az(l)=1.-f(l,1,nn)
*
* diagnostics de nuages 2-D
nhaut(l) = 1. - f(l, 1 ,IH(l))
nmoy (l) = 1. - f(l,IH(l),IB(l))
nbas (l) = 1. - f(l,IB(l),nn )
*
208 continue
*
c En mode reduction, interpoler fd et fm aux niveaux de flux complets
if( reduc ) then
call intchamps(flfm,fm,flss,s,nl,flnn,nn)
call intchamps(flfd,fd,flss,s,nl,flnn,nn)
do i=1,flnn
do l=1,nl
fm(l,i) = flfm(l,i)
fd(l,i) = flfd(l,i)
enddo
enddo
endif
do 2005 i=1,flnn
do 2006 l=1,nl
xin(l,i)=fm(l,i)
v(l,i)=fd(l,i)
2006 continue
2005 continue
c smoothing 1/4-1/2-1/4 des flux
*************************************
*
if ( radfltr ) then
do i=2,flnk
do l=1,nl
fm(l,i)=0.25*(xin(l,i-1)+xin(l,i+1)) + 0.5*xin(l,i)
fd(l,i)=0.25*(v(l,i-1)+v(l,i+1)) + 0.5*v(l,i)
enddo
enddo
else
do i=2,flnk
do l=1,nl
fm(l,i)=xin(l,i)
fd(l,i)=v(l,i)
enddo
enddo
endif
***************************************
secday=1./86400.
c taux de refroidissement
do 8 i=1,flnk
c
* taumax = (( -877.193 * (flsh(i)-0.05)**2 ) - 0.5)*secday
*
if(radfix) then
* >>>>>> Seulement si RADFIX est vrai <<<<<<<
do 108 l=1,nl
! taumax = (( -877.193 * (flsh(l,i)*ps(l)*1.E-5-0.05)**2) - 0.5)*secday
taumax= (( -877.193 * (flsh(l,i)*ps(l)*1.E-5-0.05)*
x (flsh(l,i)*ps(l)*1.E-5-0.05)) - 0.5)*secday
taux(l,i)=-(fm(l,i)-fd(l,i)-fm(l,i+1)+fd(l,i+1)) /
x (ps(l)*cpd*fldel(l,i)) *grav
c fit du taux au dessus de sigma=0.011 en fonction de temperature
*
* legere correction des taux au-dessus de 100 mb de facon a
* obtenir un equilibre quasi-parfait avec les taux visibles
* de sun1. correction inferieure a 0.15 k/j. le taux au toit
* est fixe a -1.9 k/jour.
r = flsh(l,i)*ps(l)*0.01
if (r.lt.100.) taux(l,i)= taux(l,i)+0.15*secday
if (i.eq.1) taux(l,i)= -1.9*secday
* taux(l,i)=cvmgt((-1.9*secday),taux(l,i),(i.eq.1))
* taux(l,i)=cvmgt((-1.9*secday),taux(l,i),flsh(l,i).lt.0.011)
*
* limite sur les taux (en general trop froids) lorsque sigma <= 0.05
if ((flsh(l,i)*ps(l)*1.E-5.le.0.05).and.(i.ne.1)) then
taux(l,i) = max(taumax,taux(l,i))
* >>>>>> Seulement si RADFIX est vrai <<<<<<<
endif
108 continue
*
else
*
* calcul des taux radiatifs sans corrections
do l=1,nl
taux(l,i)=-(fm(l,i)-fd(l,i)-fm(l,i+1)+fd(l,i+1)) /
x (ps(l)*cpd*fldel(l,i)) *grav
enddo
*
endif
8 continue
c
c aux deux bouts taux calcules avec flux sans smoothing
do 1088 l=1,nl
y=-(xin(l,1)-v(l,1)-xin(l,2)+v(l,2)) /
x (ps(l)*cpd*fldel(l,1)) *grav
c si radfix est vrai, le taux est fixe par fit en temperature
if ((flsh(l,1).ge.0.011) .and. .not.radfix ) taux(l,1) = y
taux(l,flnk)=-(xin(l,flnk)-v(l,flnk)-xin(l,flnn)+v(l,flnn)) /
x (ps(l)*cpd*fldel(l,flnk)) *grav
1088 continue
c
return
end