! $Id: ropp_pp_bangle_BARO.f90 2264 2009-10-09 18:05:20Z frhl !****s* ModelRefraction/ropp_pp_bangle_BARO * ! ! NAME ! ropp_pp_bangle_BARO - Compute bending angle from BAROCLIM spherical harmonics ! ! SYNOPSIS ! call ropp_pp_bangle_BARO(file, month, lat, alt, bangle) ! ----> NO LON ! ! DESCRIPTION ! This subroutine calculates a climatological bending angle profile for ! a given month, latitude and [longitude] from BAROCLIM data. Profiles are ! computed using Chebyshev polynomials and zonal harmonics with ! coefficients read from file. ! ! INPUTS ! character(len=*) :: mfile Model coefficients filename ! integer, :: month Month of year ! real(wp), :: lat Latitude ! real(wp), dim(:) :: alt Altitude levels on which to find N ! ! OUTPUT ! character(len=*) :: file Model coefficients filename ! real(wp), dim(:) :: bangle Bending angle field ! ! AUTHOR ! Met Office, Exeter, UK. ! Any comments on this software should be given via the ROM SAF ! Helpdesk at http://www.romsaf.org ! ! COPYRIGHT ! Copyright (c) 1998 Stig Syndergaard . ! For further details please refer to the file COPYRIGHT ! which you should have received as part of this distribution. ! !**** SUBROUTINE ropp_pp_bangle_BARO(mfile, month, lat, alt, bangle) !------------------------------------------------------------------------------- ! 1. Declarations !------------------------------------------------------------------------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_pp_BARO, not_this => ropp_pp_bangle_BARO IMPLICIT NONE CHARACTER(len=*), INTENT(inout) :: mfile ! Coefficients file INTEGER, INTENT(in) :: month ! Month of year REAL(wp), INTENT(in) :: lat ! Latitude REAL(wp), DIMENSION(:), INTENT(in) :: alt ! Altitude (m) REAL(wp), DIMENSION(:), INTENT(out) :: bangle ! Bending angle (rad.) TYPE(BAROcoeff_ba), SAVE :: coeff ! Spherical harmonic coeff REAL(wp), DIMENSION(:), ALLOCATABLE :: n1 ! Constants REAL(wp), DIMENSION(:), ALLOCATABLE :: n2 REAL(wp) :: theta, costheta REAL(wp) :: alpha, beta, a0 , alpha0 REAL(wp), PARAMETER :: deg2rad = 0.0174532925_wp ! degrees to radians REAL(wp), PARAMETER :: alt0 = 100.0_wp ! INTEGER, PARAMETER :: klimCoef = 127 ! WHY should this be hard coded, if then ! INTEGER, PARAMETER :: nlimCoef = 17 ! why no test when using baro-coeff? REAL(wp),DIMENSION(:), ALLOCATABLE :: chjs INTEGER :: i, k, ncoeff, kcoeff LOGICAL :: found CHARACTER(len=*), PARAMETER :: & filepath(9) = (/ 'data/BAROCLIM_coeff.nc ', & '../data/BAROCLIM_coeff.nc ', & '../../data/BAROCLIM_coeff.nc ', & '../../../data/BAROCLIM_coeff.nc ', & '../../../../data/BAROCLIM_coeff.nc ', & '*/data/BAROCLIM_coeff.nc ', & '*/*/data/BAROCLIM_coeff.nc ', & '*/*/*/data/BAROCLIM_coeff.nc ', & '*/*/*/*/data/BAROCLIM_coeff.nc ' /) !------------------------------------------------------------------------------- ! 2. Read coefficients file !------------------------------------------------------------------------------- IF (.not. BARO_read) THEN INQUIRE(File=mfile, Exist=Found) IF (.NOT. Found) THEN DO i=1,9 INQUIRE(File=filepath(i), Exist=Found) IF (Found) THEN mfile = filepath(i) EXIT ENDIF ENDDO ENDIF CALL ropp_pp_read_BARO(mfile, month, coeff) BARO_read = .true. ENDIF !------------------------------------------------------------------------------- ! 3. Constants for Clenshaw's recurrence formula !------------------------------------------------------------------------------- ! FIXME: Test that ncoeff = nlimCoeff (=18) and kcoeff = klimCoef (=128) ! Or should it just be that ncoeff .le. nlimCoeff? HJS ncoeff = SIZE(coeff%Aa) ! = nlimCoeff = 18? kcoeff = SIZE(coeff%Ac,2) ! = klimCoeff = 128? ! ncoeff = nlimCoef ! nlimCoeff ! kcoeff = klimCoef ! klimCoeff ALLOCATE(n1(ncoeff)) ALLOCATE(n2(ncoeff)) ALLOCATE(chjs(kcoeff)) DO i=1,ncoeff-1 n1(i) = (2.0_wp*i + 1.0_wp)/(i + 1.0_wp) n2(i) = (i + 1.0_wp)/(i + 2.0_wp) ENDDO !SJH !------------------------------------------------------------------------------- ! 4. Calculate lat/[lon] dependent parameters !------------------------------------------------------------------------------- theta = (90.0_wp - lat) * deg2rad costheta = COS(theta) DO k=1,kcoeff chjs(k) = ZonalHarmonics( coeff%Ac(1:ncoeff ,k) , & ncoeff , & n1, & n2, & costheta ) END DO alpha = ZonalHarmonics( coeff%Aa, & ncoeff, & n1, & n2, & costheta ) beta = ZonalHarmonics( coeff%Ab, & ncoeff, & n1, & n2, & costheta ) alpha0 = ZonalHarmonics( coeff%Ad, & ncoeff, & n1, & n2, & costheta ) a0 = ZonalHarmonics( coeff%Ae, & ncoeff, & n1, & n2, & costheta ) !------------------------------------------------------------------------------- ! 5. Compute bending angle !------------------------------------------------------------------------------- CALL calc_babangle(alt,bangle,alpha,beta,a0,alpha0) !------------------------------------------------------------------------------- ! 6. Clean up !------------------------------------------------------------------------------- DEALLOCATE(n1) DEALLOCATE(n2) DEALLOCATE(chjs) CONTAINS !------------------------------------------------------------------------------- ! 7. Calculation of bending angle - Chebyshev polynomial expansion !------------------------------------------------------------------------------- SUBROUTINE calc_babangle(alt, bangle, alpha, beta, a0 , alpha0) USE typesizes, ONLY: wp => EightByteReal IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(in) :: alt REAL(wp), DIMENSION(:), INTENT(out) :: bangle REAL(wp), INTENT(in) :: alpha REAL(wp), INTENT(in) :: beta REAL(wp), INTENT(in) :: a0 REAL(wp), INTENT(in) :: alpha0 REAL(wp) :: x, z , alt_baro REAL(wp) :: dk, dk1, dk2 REAL(wp) :: fx, hs REAL(wp), PARAMETER :: zref = 100.0_wp INTEGER :: i DO i=1,SIZE(alt) ! 6.0 Change ... alt_baro = alt(i)/1000._wp - a0 ! impact height in km ! 6.1 Change of variables IF(alt_baro < 0.0_wp) THEN z = TANH(alt_baro)/ alt0 ELSE z = alt_baro / alt0 ENDIF x = 1.0_wp - 2.0_wp*EXP(-z) ! 5.2 Clenshaw's recurrence formula (Chebychev polynomials) dk1 = 0.0_wp dk2 = 0.0_wp DO k=kcoeff,2,-1 dk = 2.0_wp*x*dk1 - dk2 + chjs(k) dk2 = dk1 dk1 = dk ENDDO fx = x*dk1 - dk2 + 0.5_wp*chjs(1) ! 5.3 Compute bending angle scale height hs = fx * EXP(-(z/zref)**2) + (alpha*z + beta) ! 5.4 Compute bending angle bangle(i) = alpha0 * exp(-alt_baro/hs) ENDDO END SUBROUTINE calc_babangle !------------------------------------------------------------------------------- ! 8. Calculation of spherical harmonics - Clenshaw's recurrence formula !------------------------------------------------------------------------------- FUNCTION ZonalHarmonics(A,nmaxCoef,n1,n2,costheta) RESULT(sharm) USE typesizes, ONLY: wp => EightByteReal IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(in) :: A INTEGER, INTENT(in) :: nmaxCoef REAL(wp), DIMENSION(:), INTENT(in) :: n1 REAL(wp), DIMENSION(:), INTENT(in) :: n2 REAL(wp), INTENT(in) :: costheta INTEGER :: n REAL(wp) :: sharm REAL(wp) :: dn , dn1 , dn2 dn1 = 0.0_wp dn2 = 0.0_wp DO n=nmaxCoef-1,1,-1 dn = costheta * n1( n ) * dn1 - n2( n ) * dn2 + A( n+1 ) dn2 = dn1 dn1 = dn ENDDO sharm = costheta * dn1 - 0.5_wp * dn2 + A( 1 ) END FUNCTION ZonalHarmonics END SUBROUTINE ropp_pp_bangle_BARO