!This oblate spheroid code is used to calculate a global mean from zonal !means given that the Earth is not a true sphere. ! !Based on the shape of the Earth, each 1-deg latitude mean taken at ![89.5, 88.5, 87.5, ..., 1.5, 0.5, -0.5, -1.5, ..., -88.5, -89.5] !is given a weighting value more closely representing the area of the !earth at that latitude. ! !It is meant to be used in place of cos(latitude) weighting (for sphere). ! Calculate the geodetic latitude weighting factors. program oblate real wate_g(180) real*8 re,rp,r,p,a,b,beta,f re = 6378.137 !radius around the equator rp = 6356.752 !radius around the poles r = (re*rp)**0.5 dtor = 3.14159/180.0 B = rp/re ! b=minor/major radii ratio boa2 = b**2 ek2 = 1. - b**2 ek = sqrt(ek2) sum_g = 0.0 do ilat=1,180 clat0 = 91.0 - float(ilat) if(clat0 .gt. 89.99) then clat0 = 89.99 endif wate0 = boa2*ar(ek,clat0) clat1 = 90.0 - float(ilat) if(clat1 .lt. -89.99) then clat1 = -89.99 endif wate1 = boa2*ar(ek,clat1) wate_g(ilat) = wate0 - wate1 sum_g = sum_g + wate_g(ilat) end do print *, sum_g open(1,file='CERES_geodetic1') do ilat = 1,180 xlat = 90.5 - ilat write(1,'(f7.2,f10.6)') xlat,wate_g(ilat)/sum_g end do close(1) stop end function ar( ek, gmu) dtor = 3.14159/180.0 Rtod = 180.0/3.14159 u = sin(gmu*dtor) If ((ek*u) .lt. 1.0) then if (ek .gt. 0.01) then ar =0.25*(u/(1.+ek*u)+u/(1.-ek*u)+alog((1.+ek*u)/(1.-ek*u))/ek) else ar =0.25*(u/(1.+ek*u)+u/(1.-ek*u)+2.*u) endif else ar = -1000. endif return end function geodeticlat( clat) re = 6378.137 rp = 6356.752 dtor = 3.14159/180.0 rtod = 180.0/3.14159 cm = tan(dtor*clat) ! sqrt(a)/sqrt(b)=sqrt(a/b) rx = re*rp*sqrt((1.0+cm**2)/(rp**2+cm**2*re**2)) ! earthrad(rx) at clat geodeticlat = rtod*atan(cm*re**2/rp**2) return end real function geod_c(rlat,iway) ! IWAY = 1 geoDETIC --> geoCENTRIC ! IWAY =-1 geoCENTRIC -->geoDETIC parameter (re=6367.0, requa= 6378.137 , rpole = 6356.752) parameter (dtr=3.1415927/180.0) parameter (factor= rpole**2/requa**2) tand(x) = tan(x*dtr) atand(x) = atan(x)/dtr x = tand(rlat) if(iway > 0) then geod_c = atand(x * factor ) elseif(iway < 0)then geod_c = atand(x / factor) else Stop ' geod_c:IWay ==0' endif if ( abs( rlat) > 89.9 ) geod_c = rlat return end