CCC********************************************************************* CCC IN THIS VERSION A* := 3/2 [ Gamma/f_L ] * (R_EM / R_E )**2 CCC IS CALCULATED CCC calculates albedo taking one map at the time CCC********************************************************************* real mlat,mlon, dh character*42 fname integer year,day, daym, mon character cldname*23, icename*18, resname*10, mt*2, dt*2, tm*4 character yeardir*2, yeartime*4 character*3 mn(12) double precision dsind, dcosd, dasind, dacosd, datan2d dimension icv(8192), ibuf(128,64), buf(128,64), ratio(128,64) dimension tmp(128,64),temp2(128,64),temp3(128,64),temp4(128,64) dimension temp5(128,64) equivalence (ic(1,1),icv(1)) common/clouds/ic(128,64) common/ice/ice(128,64) common/contrib/contrib(128,64) common/albmap/albmap(128,64) common/cos_sun/cos_sun(128,64) common/cos_moon/cos_moon(128,64) common/area/area(128,64) external sind, cosd, dsind, dcosd, dasind, dacosd external datan2d, acosd, asind, atan2d data mn /'jan','feb','mar','apr','may','jun', *'jul','aug','sep','oct','nov','dec'/ C yeardir='99' open(1, file='albedo2.ini') read(1,*) daym, tm, mon, year read(1,*) icename read(1,*) hst, hen, dh close(1) C if (year.eq.1994) yeardir='94' C if (year.eq.1995) yeardir='95' C if (year.eq.1999) yeardir='99' C if (year.eq.2001) yeardir='01' C if (year.eq.2000) yeardir='00' C if (year.eq.2002) yeardir='02' C if (year.eq.2003) yeardir='03' C if (year.eq.2004) yeardir='04' C if (year.eq.2005) yeardir='05' if (year.eq.1983) yeartime='1983' if (year.eq.1984) yeartime='1984' if (year.eq.1985) yeartime='1985' if (year.eq.1986) yeartime='1986' if (year.eq.1987) yeartime='1987' if (year.eq.1988) yeartime='1988' if (year.eq.1989) yeartime='1989' if (year.eq.1990) yeartime='1990' if (year.eq.1991) yeartime='1991' if (year.eq.1992) yeartime='1992' if (year.eq.1993) yeartime='1993' if (year.eq.1994) yeartime='1994' if (year.eq.1995) yeartime='1995' if (year.eq.1996) yeartime='1996' if (year.eq.1997) yeartime='1997' if (year.eq.1998) yeartime='1998' if (year.eq.1999) yeartime='1999' if (year.eq.2001) yeartime='2001' if (year.eq.2000) yeartime='2000' if (year.eq.2002) yeartime='2002' if (year.eq.2003) yeartime='2003' if (year.eq.2004) yeartime='2004' if (year.eq.2005) yeartime='2005' if (year.eq.2006) yeartime='2006' if (year.eq.2007) yeartime='2007' if (year.eq.2008) yeartime='2008' if (year.eq.2009) yeartime='2009' if (year.eq.2010) yeartime='2010' call istr2(daym,dt) call istr2(mon,mt) C cldname=dt//mn(mon)//tm//'.map' C resname=dt//mn(mon)//tm C cldname='13jun1200.map' C resname=dt//mn(mon)//tm if (year.eq.1983) icename='sicanada1984'//mt//'.map' if (year.ge.2005) icename='sicanada2004'//mt//'.map' C If one wants to fix the ice maps C icename='sicanada200106.map' cldname= 'T42_'//yeartime//'.'//mt//'.'//dt//'_davg.dat' resname= yeartime//'.'//mt//'.'//dt if (year.ge.2005) cldname='T42_2004.'//mt//'.'//dt//'_davg.dat' write(*,*) cldname write(*,*) icename write(*,*) resname write(*,*) year write(*,*) mon, daym write(*,*) hst, hen, dh call getmod call getgrid rem0=60.336 do i=1,128 do j=1,64 contrib(i,j)=0. albmap(i,j)=0. cos_sun(i,j)=0. cos_moon(i,j)=0. area(i,j)=0. enddo enddo c -------Reading Snow/Ice Cover Data ------------- open(unit=89,file='/data/model/icedata/'//icename) c open(unit=89,file='/data2/epb/Albedo_model/icedata/clear.map') READ(89,*) a,b,c,d do i=1,64 read(89,*) (ibuf(j,i), j=1,128) enddo close(89) c write(*,*) 'Snow/Ice Data is OK' do i=1,64 do j=1,128 if (j.lt.65) jj=j+64 if (j.ge.65) jj=j-64 ice(j,i)=ibuf(jj,i) enddo enddo c -------Reading Fractional Cloud Cover Data ------------- C open(87,file='/data2/epb/Albedo_model/cloudata/WSI/'// C *mn(mon)//yeardir//'/'//cldname) open(87,file='/data/model/cloudata/'// * 'ISCCP/'//cldname) c READ THE FILES INVERTED IN LATITUDE. WE PUT IT BACK LATER do i=64,1,-1 read(87,*) (ibuf(j,i), j=1,128) enddo close(87) c--------- To remove Blank Places C do i=1,5 C do j=1,128 C ibuf(j,i)=ibuf(j,6) c enddo c enddo c do i=60,64 c do j=1,128 c ibuf(j,i)=ibuf(j,59) c enddo c enddo C do i=1,64 C ibuf(128,i)=ibuf(127,i) C enddo c--------------------------------- C WE PUT LONGITUDES FROM -180 TO 180 do i=1,64 do j=1,128 if (j.lt.65) jj=j+64 if (j.ge.65) jj=j-64 ic(j,i)=ibuf(jj,i) enddo enddo c write(*,*) 'Cloud CoverData is OK' open(unit=3,file= '/data/model/results/'//resname//'.dat') hour=hst do while (hour.le.hen) day=getday(mon,daym,year) call sunmoon(year,day,hour,slon,slat,mlon,mlat,phase,rem) c write(*,*) year,day,hour,slon,slat,mlon,mlat,phase,rem c write(*,*) 'Subroutine _sunmoon_ is OK' call sums(slon,slat,mlon,mlat,ae,etos) c write(*,*) 'Subroutine _sums_ is OK' fl = sind(abs(phase))/3.141592654-abs(phase)*cosd(phase)/180. if (fl .le. 1.e-10) fl=1.e-10 Astar = 1.5*etos/fl Astar_=-2.5*log10(Astar) write(*,200) year,day,hour,slon,slat,mlon,mlat, * rem,phase,ae,Astar,Astar_ 200 format(2(i6),7(1x,f6.1),1x,3f8.4) c --------- Save the Results ------------------ write(3,200) year,day,hour,slon,slat,mlon,mlat, * rem,phase,ae,Astar,Astar_ c ---------------------------------------------- hour=hour+dh enddo 100 format(24i3) 300 format(128i1) close(3) c ----------------------------------------------- C PUT BACK THE LONGITUDE FROM 0 TO 360 do i=1,64 do j=1,128 if (j.lt.65) jj=j+64 if (j.ge.65) jj=j-64 tmp(j,i)=contrib(jj,i) temp2(j,i)=cos_sun(jj,i) temp3(j,i)=cos_moon(jj,i) temp4(j,i)=area(jj,i) temp5(j,i)= 1.5 * albmap(jj,i) /fl C if (temp5(j,i).le.0) temp5(j,i)= 0. C if (temp5(j,i).gt.0) temp5(j,i)= -2.5*log10(temp5(j,i)) enddo enddo C WHEN WE WRITE THE FILES IF GOES AGAIN 64->1 START AT SOUTH POLE open(1,file='/data/model/results/'//resname//'.lit') write(*,*) resname//'.lit' do i=64,1,-1 write(1,'(128f5.1)') (tmp(j,i), j=1,128) enddo close(1) end subroutine istr2(il,str) integer ival,il,index(10) character*(10) str1 character*(2) str character*(1) num(0:9) data num/'0','1','2','3','4','5','6','7','8','9'/ ival=il if (ival.lt.0) ival=ival*(-1) do 10 j=10,1,-1 i=ival/10**(j-1) ival=ival-i*10**(j-1) 10 index(j)=i str1=num(index(10))//num(index(9))//num(index(8))// *num(index(7))//num(index(6))//num(index(5))//num(index(4))// *num(index(3))//num(index(2))//num(index(1)) str=str1(9:10) return end subroutine getfile(nfile,fname) character*42 fname character*7 strmd(37) data strmd /'cloud1','cloud2','87feb10','87mar02', + '87mar22','87apr11','87may01','87may21', + '87jun10','87jun30','87jul20','87aug09', + '87aug29','87sep18','87oct08','87oct28', + '87nov17','87dec07','87dec27','88jan16', + '88feb05','88feb25','88mar16','88apr05', + '88apr25','88may15','88jun04','88jun24', + '88jul14','88aug03','88aug23','88sep12', + '88oct02','88oct22','88nov11','88dec01', + '88dec21' / c fname='data/'//strmd(nfile) fname=strmd(nfile) return end function getday(monthn, daym, year) integer monthn, daym, year integer month(12) data month /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ if (year.eq.1984.or.year.eq.1988.or.year.eq.1992.or.year.eq.1996 *.or.year.eq.2000) month(2)=29 sum=0 do i=1, monthn-1 sum=sum+month(i) enddo getday=sum+daym return end subroutine sunmoon(year,day,hour,slon,slat,mlon,mlat,phase,rem) c*******Returns, for a given year and day, the sub-solar and sub-lunar c*******points on the earth, as well as the earth-moon distance c*******INPUT: c******* year (integer): the year (integer) c******* day (integer): day of the year c******* hour (single ): UT in hours and fractions c*******OUTPUT: c******* slon (single ) : sub-solar latitude (degrees) c******* slat (single ) : sub-solar longitude (degrees east) c******* mlon (single ) : sub-lunar latitude (degrees) c******* mlat (single ) : sub-lunar longitude (degrees east) c******* phase (single ) : lunar phase (degrees) c******* rem (single ) : earth-moon distance (mean earth radii) integer year,day real hour,slon,slat,mlon,mlat,phase,rem double precision thetas, phis, thetam, phim double precision xn,tu,gmst double precision epsil,xl,g,xlam,dec,ra double precision beta,c1,c2,c3,pi double precision dsind, dcosd, dasind, dacosd, datan2d external sind, cosd, dsind, dcosd, dasind, dacosd external datan2d, acosd, asind, atan2d xn=-3288.5d0-(1991.d0-dble(year))*365.25d0 !number of days from J2000.0 xn=xn+dble(day)+dble(hour)/24.d0 tu=xn/36525.d0 !number of Julian centuries gmst=24110.54841d0+8640184.812866d0*tu gmst=gmst+0.093104d0*tu**2 gmst=gmst-6.2d-6*tu**3 gmst=gmst/3600.d0+dble(hour) gmst=dmod(gmst,24.d0) if (gmst.lt.0) gmst=gmst+24.d0 c*******do the sun first epsil=23.439d0-4.0d-7*xn !obliquity of ecliptic xl=280.460d0+0.9856474d0*xn !mean long. of sun g= 357.528d0+0.9856003d0*xn !mean anomaly xl=dmod(xl,360.d0) !reduce to interval -360 to 360 g =dmod(g ,360.d0) if (xl.lt.0.d0)xl=xl+360.d0 !reduce to interval 0 to 360 if (g .lt.0.d0)g =g +360.d0 xlam=xl+1.915d0*dsind(g) !ecliptic longitude xlam=xlam+2.d0-2.d0*dsind(2.d0*g) dec=dasind(dsind(xlam)*dsind(epsil)) !solar declination, ascension ra=datan2d(dcosd(epsil)*dsind(xlam),dcosd(xlam)) if (ra.lt.0.d0) ra=ra+360.d0 phis=ra-gmst*15.d0 if (phis.lt.0.d0) phis=phis+360.d0 thetas=dec slon=phis slat=thetas c*******now do the moon xlam=218.32d0+481267.883d0*tu xlam=xlam+6.29d0*dsind(134.9d0+477198.85d0*tu) xlam=xlam-1.27d0*dsind(259.2d0-413335.38d0*tu) xlam=xlam+0.66d0*dsind(235.7d0+890534.23d0*tu) xlam=xlam+0.21d0*dsind(269.9d0+954397.70d0*tu) xlam=xlam-0.19d0*dsind(357.5d0+ 35999.05d0*tu) xlam=xlam-0.11d0*dsind(186.6d0+966404.05d0*tu) beta=5.13d0*dsind( 93.3d0+483202.03*tu) beta=beta+0.28d0*dsind(228.2d0+960400.87d0*tu) beta=beta-0.28d0*dsind(318.3d0+ 6003.18d0*tu) beta=beta-0.17d0*dsind(217.6d0-407332.20d0*tu) c1=dcosd(beta)*dcosd(xlam) c2=0.9175d0*dcosd(beta)*dsind(xlam)-0.3978*dsind(beta) c3=0.3978d0*dcosd(beta)*dsind(xlam)+0.9175*dsind(beta) dec=dasind(c3) ra=datan2d(c2,c1) if (ra.lt.0.d0)ra=ra+360.d0 phim=ra-gmst*15.d0 if (phim.lt.0.d0) phim=phim+360.d0 thetam=dec mlon=phim mlat=thetam pi=0.9508d0 pi=pi+0.0518d0*dcosd(134.9d0+477198.85d0*tu) pi=pi+0.0095d0*dcosd(259.2d0-413335.38d0*tu) pi=pi+0.0078d0*dcosd(235.7d0+890534.23d0*tu) pi=pi+0.0028d0*dcosd(269.9d0+954397.70d0*tu) C distance in earth radii rem=1./dsind(pi) dotem= cosd(slat)*cosd(slon)*cosd(mlat)*cosd(mlon) dotem=dotem+cosd(slat)*sind(slon)*cosd(mlat)*sind(mlon) dotem=dotem+sind(slat)*sind(mlat) phase=acosd(-dotem) cross= cosd(slat)*cosd(slon)*cosd(mlat)*sind(mlon) cross=cross-cosd(slat)*sind(slon)*cosd(mlat)*cosd(mlon) if (cross.ge.0.) phase=-phase return end subroutine sums(slon,slat,mlon,mlat,ae,etos) real mlat,mlon common/grid/gridlon(128),gridlat(64),clat(64),mask(128,64) 1 ,re(3,128,64),dlat(64) common/clouds/ic(128,64) common/contrib/contrib(128,64) common/albmap/albmap(128,64) common/cos_sun/cos_sun(128,64) common/cos_moon/cos_moon(128,64) common/area/area(128,64) common/ice/ice(128,64) dimension rs(3),rm(3) double precision dsind, dcosd, dasind, dacosd, datan2d external sind, cosd, dsind, dcosd, dasind, dacosd external datan2d, acosd, asind, atan2d rs(3)=sind(slat) rm(3)=sind(mlat) rs(1)=cosd(slat)*cosd(slon) rm(1)=cosd(mlat)*cosd(mlon) rs(2)=cosd(slat)*sind(slon) rm(2)=cosd(mlat)*sind(mlon) dotsm=dot(rs,rm) ae=0. etos=0. do j=1,64 wt=clat(j)*dlat(j) C write(*,*) clat(j), dlat(j) do i=1,128 area(i,j)=wt dotes=dot(re(1,i,j),rs) if (dotes.ge.0) then if (dotes.ge.0.999999) dotes=0.999999 zenith=acosd(dotes) itype1=mask(i,j) icl=ic(i,j) if (icl.ge.95) then itype=12 else if (icl.ge.50) then if (itype1.ge.2.and.itype1.le.4) then itype=10 else if(itype1.eq.1) then itype=9 else itype=11 endif else if (icl.ge.5) then if (itype1.ge.2.and.itype1.le.4) then itype=7 else if(itype1.eq.1) then itype=6 else itype=8 endif else itype=itype1 if (ice(i,j).ne.0) itype=3 endif c******************************* c write(*,*) 'Vixa 2' c******************************* itype=2 asav=albedo(zenith,itype) ae=ae+asav*dotes*wt dotem=dot(re(1,i,j),rm) if (dotem.ge.0.) then if (dotem.ge.0.999999) dotem=0.999999 contrib(i,j)=1 cos_sun(i,j)= dotes cos_moon(i,j)= dotem C area(i,j)=wt view=acosd(dotem) azimuth=dotem*dotes-dotsm azimuth=azimuth/sqrt(1.-dotem**2) azimuth=azimuth/sqrt(1.-dotes**2) if (abs(azimuth).le.0.999999) then azimuth=acosd(azimuth) else if (azimuth.ge.0.999999) then azimuth=0. else azimuth=180. endif asav=asav*anifac(azimuth,view,zenith,itype) asym=albedo(view,itype)*anifac(azimuth,zenith,view,itype) asym=(asym+asav)/2. etos=etos+asym*wt*dotem*dotes albmap(i,j)= (asym*wt*dotem*dotes)/64.0 C write(*,*) 'Hello',j , i, wt, dotem, dotes endif endif enddo C write(*,*) 'Hello',j , i, wt, dotem, dotes enddo ae=ae/64. etos=etos/64. return end subroutine getmod common/model/alb(10,12),ani(8,7,10,12),zen(10),view(7),azi(8) data (zen(i) ,i=1,10)/25.84,36.87,45.57,53.13,60.00, 1 66.42,72.54,78.46,84.26,90.00/ data (view(i) ,i=1, 7)/15.,27.,39.,51.,63.,75.,90./ data (azi (i) ,i=1, 8)/9.,30.,60.,90.,120.,150.,171.,180./ open(unit=10,file='model.dat',status='old') do itype=1,12 do izen=1,10 read(10,*) alb(izen, itype) c write(*,*) alb(izen, itype) do iview=1,7 read(10,*) (ani(iazi,iview,izen,itype),iazi=1,8) enddo enddo enddo close(unit=10) return end subroutine getgrid external sind, cosd common/grid/gridlon(128),gridlat(64),clat(64),mask(128,64) 1 ,re(3,128,64),dlat(64) data (gridlat(i),i= 1, 4)/87.86379,85.09652,82.31291,79.52560/ data (gridlat(i),i= 5, 8)/76.73689,73.94751,71.15775,68.36775/ data (gridlat(i),i= 9,12)/65.57761,62.78735,59.99702,57.20663/ data (gridlat(i),i=13,16)/54.41620,51.62573,48.83524,46.04472/ data (gridlat(i),i=17,20)/43.25419,40.46365,37.67309,34.88252/ data (gridlat(i),i=21,24)/32.09194,29.30136,26.51077,23.72017/ data (gridlat(i),i=25,28)/20.92957,18.13897,15.34836,12.55776/ data (gridlat(i),i=29,32)/9.767145,6.976533, 4.18592,1.395307/ do i=1,128 gridlon(i)=(i-1)*2.8125 enddo do i=33,64 gridlat(i)=-gridlat(65-i) enddo do i=1,64 clat(i)=cosd(gridlat(i)) enddo dlat(1)=abs(90.-gridlat(2))/2.*3.141592654/180. do i=2,63 dlat(i)=abs(gridlat(i+1)-gridlat(i-1))/2.*3.141592654/180. enddo dlat(64)=dlat(1) write(*,*), dlat open (unit=10,file='t42grid.dat',status='old') do i=1,64 read(10,*) (mask(j,i),j=1,128) c write(7,100)(mask(j,i),j=1,128) enddo c 100 format(128I1) close(unit=10) do j=1,64 temps=sind(gridlat(j)) tempc=cosd(gridlat(j)) do i=1,128 re(1,i,j)=tempc*cosd(gridlon(i)) re(2,i,j)=tempc*sind(gridlon(i)) re(3,i,j)=temps enddo enddo return end function albedo(angle,itype) common/model/alb(10,12),ani(8,7,10,12),zen(10),view(7),azi(8) call sort(angle,zen,10,index) albedo=alb(index,itype) return end function anifac(azimuth,aview,zenith,itype) common/model/alb(10,12),ani(8,7,10,12),zen(10),view(7),azi(8) call sort(azimuth,azi,8,iazi) call sort(aview,view,7,iview) call sort(zenith,zen,10,izen) anifac=ani(iazi,iview,izen,itype) return end subroutine sort(val,array,idim,index) dimension array(1) do j=1,idim if (array(j).ge.val) go to 10 enddo write (6,*)'value exceeded array' stop 10 index=j return end function dot(v1,v2) dimension v1(3),v2(3) dot=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) return end function sind(degree) real degree sind=sin(degree*3.141592654/180.) return end function cosd(degree) real degree cosd=cos(degree*3.141592654/180.) return end double precision function dsind(degree) double precision degree dsind=dsin(degree*3.141592654d0/180.d0) return end double precision function dcosd(degree) double precision degree dcosd=dcos(degree*3.141592654/180.d0) return end double precision function dasind(argum) double precision argum dasind=(180.d0/3.141592654)*asin(argum) return end double precision function dacosd(argum) double precision argum dacosd=(180.d0/3.141592654)*acos(argum) return end double precision function datan2d(x,y) double precision x, y datan2d=(180.d0/3.141592654)*datan2(x,y) return end function asind(argum) real argum asind=(180./3.141592654)*asin(argum) return end function acosd(argum) real argum acosd=(180./3.141592654)*acos(argum) return end function atan2d(x,y) real argum atan2d=(180./3.141592654)*atan2(x,y) return end