SURFEX v7.3: /home/dasprezs/EXPORT_v7_3/src/SURFEX/hor_interpol

Go to the documentation of this file.
00001 !     #########
00002 SUBROUTINE HOR_INTERPOL_ROTLATLON(KLUOUT,PFIELDIN,PFIELDOUT)
00003 !     #################################################################################
00004 !
00005 !!****  *HOR_INTERPOL_ROTLATLON * - Interpolation from a rotated lat/lon grid
00006 !!
00007 !!    PURPOSE
00008 !!    -------
00009 !!
00010 !!**  METHOD
00011 !!    ------
00012 !!
00013 !!    REFERENCE
00014 !!    ---------
00015 !!      
00016 !!
00017 !!    AUTHOR
00018 !!    ------
00019 !!     U. Andrae
00020 !!
00021 !!    MODIFICATIONS
00022 !!    -------------
00023 !!      Original    10/2007
00024 !!------------------------------------------------------------------
00025 !
00026 !
00027 !
00028 USE MODD_PREP,       ONLY : XLAT_OUT, XLON_OUT
00029 USE MODD_GRID_ROTLATLON
00030 USE MODD_SURF_PAR,   ONLY : XUNDEF
00031 USE MODD_GRID_GRIB,  ONLY : NNI
00032 
00033 !
00034 !
00035 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00036 USE PARKIND1  ,ONLY : JPRB
00037 !
00038 IMPLICIT NONE
00039 !
00040 !*      0.1    declarations of arguments
00041 !
00042 INTEGER,              INTENT(IN)  :: KLUOUT    ! logical unit of output listing
00043 REAL, DIMENSION(:,:), INTENT(IN)  :: PFIELDIN  ! field to interpolate horizontally
00044 REAL, DIMENSION(:,:), INTENT(OUT) :: PFIELDOUT ! interpolated field
00045 !
00046 !*      0.2    declarations of local variables
00047 !
00048 
00049  INTEGER, ALLOCATABLE :: ii(:),jj(:)
00050 
00051  REAL,    ALLOCATABLE :: XLAT_IND(:),XLON_IND(:),  
00052                            XRAT_OUT(:),XRON_OUT(:),  
00053                            w00(:),w01(:),            
00054                            w10(:),w11(:)  
00055 
00056  LOGICAL, ALLOCATABLE :: LMASK(:)
00057 
00058 INTEGER :: I,J,K,L,IJ,IJ00,IJ01,IJ10,IJ11,INO,JL
00059 REAL    :: WX,WY,WSUM
00060 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00061 !
00062 !-------------------------------------------------------------------------------------
00063 IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_ROTLATLON',0,ZHOOK_HANDLE)
00064 WRITE(KLUOUT,'(A)')' | Running rotated latlon interpolation'
00065 
00066 INO = SIZE(XLAT_OUT)
00067 
00068 !
00069 !*      1.    Allocations
00070 !
00071 ALLOCATE(XRAT_OUT(INO),       &
00072            XRON_OUT(INO),       &
00073            XLAT_IND(INO),       &
00074            XLON_IND(INO),       &
00075                  II(INO),       &
00076                  JJ(INO),       &
00077                 W00(INO),       &
00078                 W01(INO),       &
00079                 W10(INO),       &
00080                 W11(INO))  
00081 
00082 ALLOCATE(LMASK(NNI))
00083 !
00084 !*  Transformation of latitudes/longitudes into rotated coordinates
00085 
00086     WRITE(KLUOUT,*)'XLAT_OUT',XLAT_OUT(10:10)
00087     WRITE(KLUOUT,*)'XLON_OUT',XLON_OUT(10:10)
00088 
00089     CALL REGROT(XLON_OUT,XLAT_OUT,   &
00090                   XRON_OUT,XRAT_OUT,   &
00091                   INO,1,INO,1,         &
00092                   XRLOP,XRLAP,1)  
00093 
00094     WRITE(KLUOUT,*)'XRAT_OUT',XRAT_OUT(10:10)
00095     WRITE(KLUOUT,*)'XRON_OUT',XRON_OUT(10:10)
00096 
00097     DO IJ=1,INO
00098        XLAT_IND(IJ) = ( XRAT_OUT(IJ) - XRILA1) / XRDY + 1.
00099        XLON_IND(IJ) = ( XRON_OUT(IJ) - XRILO1) / XRDX + 1.
00100     ENDDO
00101 
00102     PFIELDOUT(:,:) = XUNDEF
00103 
00104     DO JL=1,SIZE(PFIELDIN,2)
00105 
00106     LMASK= .TRUE.
00107     WHERE ( ABS(PFIELDIN(:,JL)-XUNDEF) < 1.e-6 ) LMASK = .FALSE.
00108 
00109     DO IJ=1,INO
00110 
00111          II(IJ)  = INT(XLON_IND(IJ))
00112          JJ(IJ)  = INT(XLAT_IND(IJ))
00113 
00114          WX  = XLON_IND(IJ) - FLOAT(II(IJ))
00115          WY  = XLAT_IND(IJ) - FLOAT(JJ(IJ))
00116 
00117          W00(IJ) = (1.-WX)*(1.-WY)
00118          W01(IJ) = (1.-WX)*    WY
00119          W10(IJ) =     WX *(1.-WY)
00120          W11(IJ) =     WX *    WY
00121 
00122          K    = II(IJ)
00123          L    = JJ(IJ)
00124          IJ00 = k   + NRX*(l   -1)
00125          IJ01 = k   + NRX*(l+1 -1)
00126          IJ10 = k+1 + NRX*(l   -1)
00127          IJ11 = k+1 + NRX*(l+1 -1)
00128 
00129          IF (.NOT. LMASK(IJ00)) w00(IJ) = 0.
00130          IF (.NOT. LMASK(IJ01)) w01(IJ) = 0.
00131          IF (.NOT. LMASK(IJ10)) w10(IJ) = 0.
00132          IF (.NOT. LMASK(IJ11)) w11(IJ) = 0.
00133 
00134             wsum = w00(IJ) + w01(IJ) + &
00135                      w10(IJ) + w11(IJ)  
00136 
00137             IF ( ABS(wsum) < 1.e-6 ) CYCLE
00138 
00139             w00(IJ) = w00(IJ) / wsum
00140             w01(IJ) = w01(IJ) / wsum
00141             w10(IJ) = w10(IJ) / wsum
00142             w11(IJ) = w11(IJ) / wsum
00143 
00144     ENDDO
00145 
00146     !
00147     ! Bi linear
00148     !
00149 
00150     WRITE(KLUOUT,*)'NRX,NRY',NRX,NRY
00151 
00152        DO IJ=1,INO
00153 
00154           K    = II(IJ)
00155           L    = JJ(IJ)
00156           IJ00 = k   + NRX*(l   -1)
00157           IJ01 = k   + NRX*(l+1 -1)
00158           IJ10 = k+1 + NRX*(l   -1)
00159           IJ11 = k+1 + NRX*(l+1 -1)
00160 
00161           WRITE(KLUOUT,*)PFIELDIN(IJ00,JL)
00162 
00163           PFIELDOUT(IJ,JL) = w00(IJ)*PFIELDIN(IJ00,JL) +       &
00164                                w01(IJ)*PFIELDIN(IJ01,JL) +       &
00165                                w10(IJ)*PFIELDIN(IJ10,JL) +       &
00166                                w11(IJ)*PFIELDIN(IJ11,JL)  
00167  
00168        ENDDO
00169     ENDDO
00170 
00171 
00172 !
00173 !*      5.    Deallocations
00174 !
00175 DEALLOCATE(XRAT_OUT,XRON_OUT,    &
00176              XLAT_IND,XLON_IND,    &
00177              II,JJ,                &
00178              W00,W01,W10,W11,      &
00179              LMASK)  
00180 !
00181 IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_ROTLATLON',1,ZHOOK_HANDLE)
00182 CONTAINS
00183 !
00184       SUBROUTINE REGROT(PXREG,PYREG,PXROT,PYROT,KXDIM,KYDIM,KX,KY, &
00185                          PXCEN,PYCEN,KCALL)  
00186 !
00187 !
00188       IMPLICIT NONE
00189 !
00190 !-----------------------------------------------------------------------
00191 !
00192 !*    CONVERSION BETWEEN REGULAR AND ROTATED SPHERICAL COORDINATES.
00193 !*
00194 !*    PXREG     LONGITUDES OF THE REGULAR COORDINATES
00195 !*    PYREG     LATITUDES OF THE REGULAR COORDINATES
00196 !*    PXROT     LONGITUDES OF THE ROTATED COORDINATES
00197 !*    PYROT     LATITUDES OF THE ROTATED COORDINATES
00198 !*              ALL COORDINATES GIVEN IN DEGREES N (NEGATIVE FOR S)
00199 !*              AND DEGREES E (NEGATIVE VALUES FOR W)
00200 !*    KXDIM     DIMENSION OF THE GRIDPOINT FIELDS IN THE X-DIRECTION
00201 !*    KYDIM     DIMENSION OF THE GRIDPOINT FIELDS IN THE Y-DIRECTION
00202 !*    KX        NUMBER OF GRIDPOINT IN THE X-DIRECTION
00203 !*    KY        NUMBER OF GRIDPOINTS IN THE Y-DIRECTION
00204 !*    PXCEN     REGULAR LONGITUDE OF THE SOUTH POLE OF THE ROTATED GRID
00205 !*    PYCEN     REGULAR LATITUDE OF THE SOUTH POLE OF THE ROTATED GRID
00206 !*
00207 !*    KCALL=-1: FIND REGULAR AS FUNCTIONS OF ROTATED COORDINATES.
00208 !*    KCALL= 1: FIND ROTATED AS FUNCTIONS OF REGULAR COORDINATES.
00209 !*
00210 !*    J.E. HAUGEN   HIRLAM   JUNE -92
00211 !
00212 !-----------------------------------------------------------------------
00213 !
00214       INTEGER KXDIM,KYDIM,KX,KY,KCALL
00215       REAL PXREG(KXDIM,KYDIM),PYREG(KXDIM,KYDIM), 
00216             PXROT(KXDIM,KYDIM),PYROT(KXDIM,KYDIM), 
00217             PXCEN,PYCEN  
00218 !
00219 !-----------------------------------------------------------------------
00220 !
00221       REAL PI,ZRAD,ZSYCEN,ZCYCEN,ZXMXC,ZSXMXC,ZCXMXC,ZSYREG,ZCYREG, 
00222             ZSYROT,ZCYROT,ZCXROT,ZSXROT,ZRADI  
00223       INTEGER JY,JX
00224       REAL(KIND=JPRB) :: ZHOOK_HANDLE
00225 !
00226 !-----------------------------------------------------------------------
00227 !
00228       IF (LHOOK) CALL DR_HOOK('REGROT',0,ZHOOK_HANDLE)
00229       PI = 4.*ATAN(1.)
00230       ZRAD = PI/180.
00231       ZRADI = 1./ZRAD
00232       ZSYCEN = SIN(ZRAD*(PYCEN+90.))
00233       ZCYCEN = COS(ZRAD*(PYCEN+90.))
00234 !
00235       IF (KCALL.EQ.1) THEN
00236 !
00237       DO JY = 1,KY
00238       DO JX = 1,KX
00239 !
00240       ZXMXC  = ZRAD*(PXREG(JX,JY) - PXCEN)
00241       ZSXMXC = SIN(ZXMXC)
00242       ZCXMXC = COS(ZXMXC)
00243       ZSYREG = SIN(ZRAD*PYREG(JX,JY))
00244       ZCYREG = COS(ZRAD*PYREG(JX,JY))
00245       ZSYROT = ZCYCEN*ZSYREG - ZSYCEN*ZCYREG*ZCXMXC
00246       ZSYROT = MAX(ZSYROT,-1.0)
00247       ZSYROT = MIN(ZSYROT,+1.0)
00248 !
00249       PYROT(JX,JY) = ASIN(ZSYROT)*ZRADI
00250 !
00251       ZCYROT = COS(PYROT(JX,JY)*ZRAD)
00252       ZCXROT = (ZCYCEN*ZCYREG*ZCXMXC +     &
00253                  ZSYCEN*ZSYREG)/ZCYROT  
00254       ZCXROT = MAX(ZCXROT,-1.0)
00255       ZCXROT = MIN(ZCXROT,+1.0)
00256       ZSXROT = ZCYREG*ZSXMXC/ZCYROT
00257 !
00258       PXROT(JX,JY) = ACOS(ZCXROT)*ZRADI
00259 !
00260       IF (ZSXROT.LT.0.0) PXROT(JX,JY) = -PXROT(JX,JY)
00261 !
00262       ENDDO
00263       ENDDO
00264 !
00265       ELSEIF (KCALL.EQ.-1) THEN
00266 !
00267       DO JY = 1,KY
00268       DO JX = 1,KX
00269 !
00270       ZSXROT = SIN(ZRAD*PXROT(JX,JY))
00271       ZCXROT = COS(ZRAD*PXROT(JX,JY))
00272       ZSYROT = SIN(ZRAD*PYROT(JX,JY))
00273       ZCYROT = COS(ZRAD*PYROT(JX,JY))
00274       ZSYREG = ZCYCEN*ZSYROT + ZSYCEN*ZCYROT*ZCXROT
00275       ZSYREG = MAX(ZSYREG,-1.0)
00276       ZSYREG = MIN(ZSYREG,+1.0)
00277 !
00278       PYREG(JX,JY) = ASIN(ZSYREG)*ZRADI
00279 !
00280       ZCYREG = COS(PYREG(JX,JY)*ZRAD)
00281       ZCXMXC = (ZCYCEN*ZCYROT*ZCXROT -     &
00282                  ZSYCEN*ZSYROT)/ZCYREG  
00283       ZCXMXC = MAX(ZCXMXC,-1.0)
00284       ZCXMXC = MIN(ZCXMXC,+1.0)
00285       ZSXMXC = ZCYROT*ZSXROT/ZCYREG
00286       ZXMXC  = ACOS(ZCXMXC)*ZRADI
00287       IF (ZSXMXC.LT.0.0) ZXMXC = -ZXMXC
00288 !
00289       PXREG(JX,JY) = ZXMXC + PXCEN
00290 !
00291       ENDDO
00292       ENDDO
00293 !
00294       ELSE
00295       WRITE(6,'(1X,''INVALID KCALL IN REGROT'')')
00296       CALL ABORT
00297       ENDIF
00298       IF (LHOOK) CALL DR_HOOK('REGROT',1,ZHOOK_HANDLE)
00299 !
00300       END SUBROUTINE REGROT
00301 !
00302 !-------------------------------------------------------------------------------------
00303 END SUBROUTINE HOR_INTERPOL_ROTLATLON