SUBROUTINE aerave_new ( ndata, & longdata,epdata,omegdata,gdata, & longref,epref,temp,nir,longir & ,epir,omegir,gir,qref,omegaref ) IMPLICIT NONE c....................................................................... c c R.Fournier 02/1996 c (modif F.Forget 02/1996) c le spectre est decoupe en "nir" bandes et cette routine calcule c les donnees radiatives moyenne sur chaque bande : l'optimisation c est faite pour une temperature au sol "temp" et une epaisseur c optique de l'atmosphere "epref" a la longueur d'onde "longref" c c dans la version actuelle, les ponderations sont independantes de c l'epaisseur optique : c'est a dire que "omegir", "gir" c et "epir/epre" sont independants de "epref". c en effet les ponderations sont choisies pour une solution exacte c en couche mince et milieu isotherme. c c entree c c ndata : taille des champs data c longdata,epdata,omegdata,gdata : proprietes radiative de l'aerosol c (longdata longueur d'onde en METRES) c * longref : longueur d'onde a laquelle l'epaisseur optique c est connue c * epref : epaisseur optique a longref c * temp : temperature choisie pour la ponderation (Planck) c * nir : nombre d'intervals dans la discretisation spectrale c du GCM c * longir : longueurs d'onde definissant ces intervals c c sortie c c * epir : epaisseur optique moyenne pour chaque interval c * omegir : "scattering albedo" moyen pour chaque interval c * gir : "assymetry factor" moyen pour chaque interval c * qref : extinction coefficient at reference wavelength c * omegaref : single scat. albedo at reference wavelength c c....................................................................... c REAL longref REAL epref REAL temp INTEGER nir REAL*8 longir(nir+1) REAL epir(nir) REAL omegir(nir) REAL gir(nir) c c....................................................................... c INTEGER iir,nirmx PARAMETER (nirmx=100) INTEGER idata,ndata c c....................................................................... c REAL emit REAL totalemit(nirmx) REAL longdata(ndata),epdata(ndata) & ,omegdata(ndata),gdata(ndata) REAL qextcorrdata(ndata) INTEGER ibande,nbande PARAMETER (nbande=1000) REAL long,deltalong INTEGER ilong INTEGER i1,i2 REAL c1,c2 REAL factep,qextcorr,omeg,g REAL qref,omegaref c c....................................................................... c DOUBLE PRECISION tmp1 REAL tmp2,tmp3 c c long=longref !if(nir.eq.27)then !print*,'long',long !print*,'longdata',longdata !print*,'epdata',epdata !print*,'omegdata',omegdata !print*,'gdata',gdata !print*,'data looks aok!' !print*,'ndata=',ndata !print*,'longdata',shape(longdata) !print*,'epdata',shape(epdata) !print*,'omegdata',shape(omegdata) !print*,'gdata',shape(gdata) ! print*,'longref',longref !print*,'epref',epref !print*,'temp',temp !print*,'nir',nir !print*,'longir',longir !print*,'epir',epir !print*,'omegir',gir !print*,'qref',qref !print*,'longir=',longir !stop !endif c******************************************************** c interpolation ilong=1 DO idata=2,ndata IF (long.gt.longdata(idata)) ilong=idata ENDDO i1=ilong i2=ilong+1 IF (i2.gt.ndata) i2=ndata IF (long.lt.longdata(1)) i2=1 IF (i1.eq.i2) THEN c1=1.E+0 c2=0.E+0 ELSE c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) ENDIF c******************************************************** c qref=c1*epdata(i1)+c2*epdata(i2) omegaref=c1*omegdata(i1)+c2*omegdata(i2) factep=qref/epref DO idata=1,ndata qextcorrdata(idata)=epdata(idata)/factep ENDDO c c....................................................................... c c write(*,*) 'TB17 : nb wl nir,nbande=',nir,nbande DO iir=1,nir c c....................................................................... c deltalong=(longir(iir+1)-longir(iir)) / nbande totalemit(iir)=0.E+0 epir(iir)=0.E+0 omegir(iir)=0.E+0 gir(iir)=0.E+0 c c....................................................................... c DO ibande=1,nbande c c....................................................................... c long=longir(iir) + (ibande-0.5E+0) * deltalong CALL blackl(DBLE(long),DBLE(temp),tmp1) emit=REAL(tmp1) c c....................................................................... c c******************************************************** c interpolation ilong=1 DO idata=2,ndata IF (long.gt.longdata(idata)) ilong=idata ENDDO i1=ilong i2=ilong+1 IF (i2.gt.ndata) i2=ndata IF (long.lt.longdata(1)) i2=1 IF (i1.eq.i2) THEN c1=1.E+0 c2=0.E+0 ELSE c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) ENDIF c******************************************************** c qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2) omeg=c1*omegdata(i1)+c2*omegdata(i2) g=c1*gdata(i1)+c2*gdata(i2) c c....................................................................... c c write(*,*) 'TB17:dtlong,emit,qcorr=',deltalong, emit,qextcorr totalemit(iir)=totalemit(iir)+deltalong*emit epir(iir)=epir(iir)+deltalong*emit*qextcorr omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g c c....................................................................... c ENDDO c c....................................................................... c gir(iir)=gir(iir)/omegir(iir) omegir(iir)=omegir(iir)/epir(iir) epir(iir)=epir(iir)/totalemit(iir) c c....................................................................... c ENDDO c c...................................................................... c c Diagnostic de controle si on moyenne sur tout le spectre vis ou IR : c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c tmp2=0.E+0 c DO iir=1,nir c tmp2=tmp2+totalemit(iir) c ENDDO c tmp3=5.67E-8 * temp**4 c IF (abs((tmp2-tmp3)/tmp3).gt.0.05E+0) THEN c PRINT *,'!!!! <---> il manque du Planck (voir moyenne.F)' c PRINT *,'somme des bandes :',tmp2,'--- Planck:',tmp3 c ENDIF c c...................................................................... c RETURN END