|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 subroutine gradient_bicubic(NX1, NY1, src_array, sou_mask, 00002 $ src_latitudes, src_longitudes, 00003 $ id_per, cd_per, 00004 $ gradient_i, gradient_j, gradient_ij) 00005 00006 C**** 00007 C ***************************** 00008 C * OASIS ROUTINE - LEVEL ? * 00009 C * ------------- ------- * 00010 C ***************************** 00011 C 00012 C**** *gradient_bicubic* - calculate gradients for bicubic remapping 00013 C 00014 C Purpose: 00015 C ------- 00016 C Calculation of gradients for bicubic interpolation. In contrast to 00017 C the gradients of conservative remapping, these gradients are 00018 C calculated with respect to grid rows and grid lines. 00019 C 00020 C** Interface: 00021 C --------- 00022 C *CALL* *gradient_bicubic*(NX1, NY1, src_array, sou_mask, 00023 C src_latitudes, src_longitudes, 00024 C gradient_i, gradient_j, gradient_ij) 00025 C 00026 C Input: 00027 C ----- 00028 C NX1 : grid dimension in x-direction (integer) 00029 C NY1 : grid dimension in y-direction (integer) 00030 C src_array : array on source grid (real 2D) 00031 C sou_mask : source grid mask (integer 2D) 00032 C src_latitudes : latitudes on source grid (real 2D) 00033 C src_longitudes : longitudes on source grid (real 2D) 00034 C id_per : number of overlapping points for source grid 00035 C cd_per : grip periodicity type 00036 C 00037 C Output: 00038 C ------ 00039 C gradient_i : gradient in i-direction (real 2D) 00040 C gradient_j : gradient in j-direction (real 2D) 00041 C gradient_ij : gradient in ij-direction (real 2D) 00042 C 00043 C History: 00044 C ------- 00045 C Version Programmer Date Description 00046 C ------- ---------- ---- ----------- 00047 C 2.5 V. Gayler 2002/04/05 created 00048 C 00049 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00050 USE constants 00051 USE kinds_mod 00052 USE mod_unit 00053 USE mod_printing 00054 00055 IMPLICIT NONE 00056 !----------------------------------------------------------------------- 00057 ! INTENT(IN) 00058 !----------------------------------------------------------------------- 00059 INTEGER (int_kind), INTENT(IN) :: 00060 $ NX1, NY1, ! source grid dimensiones 00061 $ id_per ! nbr of overlapping grid points 00062 00063 CHARACTER*8, INTENT(IN) :: 00064 $ cd_per ! grip periodicity type 00065 00066 REAL (kind = real_kind), DIMENSION(NX1,NY1), INTENT(IN) :: 00067 $ src_array ! array on source grid 00068 00069 INTEGER (int_kind), DIMENSION(NX1,NY1), INTENT(IN) :: 00070 $ sou_mask ! source grid mask 00071 00072 REAL (kind = real_kind), DIMENSION(NX1,NY1), INTENT(IN) :: 00073 $ src_latitudes, ! source grid latitudes 00074 $ src_longitudes ! source grid longitudes 00075 00076 !----------------------------------------------------------------------- 00077 ! INTENT(OUT) 00078 !----------------------------------------------------------------------- 00079 REAL (kind = real_kind), DIMENSION(NX1,NY1), INTENT(OUT) :: 00080 $ gradient_i, ! gradient in i-direction (real 2D) 00081 $ gradient_j, ! gradient in j-direction (real 2D) 00082 $ gradient_ij ! gradient in ij-direction (real 2D) 00083 00084 !----------------------------------------------------------------------- 00085 ! LOCAL VARIABLES 00086 !----------------------------------------------------------------------- 00087 INTEGER (int_kind) :: 00088 $ i, j, ! looping indicees 00089 $ ip1, jp1, im1, jm1 00090 00091 REAL (kind = real_kind) :: 00092 $ di, dj, ! factor depending on grid cell distance 00093 $ gradient_ij1, ! gradient needed to calculate gradient_ij 00094 $ gradient_ij2 ! gradient needed to calculate gradient_ij 00095 00096 REAL (kind = real_kind), DIMENSION(NX1,NY1) :: 00097 $ src_lon, ! source grid longitudes [radiants] 00098 $ src_lat, ! source grid latitudes [radiants] 00099 $ pi180 ! conversion factor: deg -> rad 00100 00101 INTEGER (int_kind) :: il_maskval= 0_int_kind 00102 00103 !----------------------------------------------------------------------- 00104 ! 00105 IF (nlogprt .GE. 2) THEN 00106 WRITE (UNIT = nulou,FMT = *)' ' 00107 WRITE (UNIT = nulou,FMT = *)'Entering routine gradient_bicubic' 00108 WRITE (UNIT = nulou,FMT = *)' ' 00109 CALL FLUSH(nulou) 00110 ENDIF 00111 ! 00112 ! Transformation from degree to radiant 00113 ! ------------------------------------- 00114 pi180 = 1.74532925199432957692e-2 ! =PI/180 00115 00116 src_lon = src_longitudes * pi180 00117 src_lat = src_latitudes * pi180 00118 00119 ! Initialization 00120 ! -------------- 00121 gradient_i = 0. 00122 gradient_j = 0. 00123 gradient_ij = 0. 00124 00125 ! calculate gradients 00126 ! ------------------- 00127 DO i = 1, NX1 00128 DO j = 1, NY1 00129 00130 IF (sou_mask (i,j) /= il_maskval) THEN 00131 00132 di = half 00133 dj = half 00134 00135 ip1 = i + 1 00136 im1 = i - 1 00137 IF (i == NX1) THEN 00138 IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian 00139 IF (cd_per == 'R') ip1 = NX1 00140 ENDIF 00141 IF (i == 1 ) THEN 00142 IF (cd_per == 'P') im1 = NX1 - id_per 00143 IF (cd_per == 'R') im1 = 1 00144 ENDIF 00145 jp1 = j + 1 00146 jm1 = j - 1 00147 IF (j == NY1) THEN ! treatment of the last.. 00148 jp1 = NY1 00149 dj = one 00150 ENDIF 00151 IF (j == 1 ) THEN ! .. and the first grid-row 00152 jm1 = 1 00153 dj = one 00154 ENDIF 00155 00156 00157 ! gradient i 00158 ! ---------- 00159 IF (sou_mask(ip1,j) /= il_maskval .OR. 00160 $ sou_mask(im1,j) /= il_maskval) THEN 00161 IF (sou_mask(ip1,j) == il_maskval) THEN 00162 ip1 = i 00163 di = one 00164 ELSE IF (sou_mask(im1,j) == il_maskval) THEN 00165 im1 = i 00166 di = one 00167 ENDIF 00168 gradient_i(i,j) = 00169 $ di * (src_array(ip1,j) - src_array(im1,j)) 00170 ENDIF 00171 00172 ! gradient j 00173 ! ---------- 00174 IF (sou_mask(i,jp1) /= il_maskval .OR. 00175 $ sou_mask(i,jm1) /= il_maskval) THEN 00176 IF (sou_mask(i,jp1) == il_maskval) THEN 00177 jp1 = j 00178 dj = one 00179 ELSE IF (sou_mask(i,jm1) == il_maskval) THEN 00180 jm1 = j 00181 dj = one 00182 ENDIF 00183 gradient_j(i,j) = 00184 $ dj * (src_array(i,jp1) - src_array(i,jm1)) 00185 ENDIF 00186 ! 00187 ! gradient ij 00188 ! ----------- 00189 di = half 00190 dj = half 00191 ip1 = i + 1 00192 im1 = i - 1 00193 IF (i == NX1) THEN 00194 IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian 00195 IF (cd_per == 'R') ip1 = NX1 00196 ENDIF 00197 IF (i == 1 ) THEN 00198 IF (cd_per == 'P') im1 = NX1 - id_per 00199 IF (cd_per == 'R') im1 = 1 00200 ENDIF 00201 jp1 = j + 1 00202 jm1 = j - 1 00203 IF (j == NY1) THEN ! treatment of the last.. 00204 jp1 = NY1 00205 dj = one 00206 ENDIF 00207 IF (j == 1 ) THEN ! .. and the first grid-row 00208 jm1 = 1 00209 dj = one 00210 ENDIF 00211 00212 gradient_ij1 = zero 00213 IF (sou_mask(ip1,jp1) /= il_maskval .OR. 00214 $ sou_mask(im1,jp1) /= il_maskval) THEN 00215 IF (sou_mask(ip1,jp1) == il_maskval .AND. 00216 $ sou_mask(i,jp1) /= il_maskval) THEN 00217 ip1 = i 00218 di = one 00219 ELSE IF (sou_mask(im1,jp1) == il_maskval .AND. 00220 $ sou_mask(i,jp1) /= il_maskval) THEN 00221 im1 = i 00222 di = one 00223 ELSE 00224 di = zero 00225 ENDIF 00226 gradient_ij1 = 00227 $ di * (src_array(ip1,jp1) - src_array(im1,jp1)) 00228 ENDIF 00229 00230 di = half 00231 ip1 = i + 1 00232 im1 = i - 1 00233 IF (i == NX1) THEN 00234 IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian 00235 IF (cd_per == 'R') ip1 = NX1 00236 ENDIF 00237 IF (i == 1) THEN 00238 IF (cd_per == 'P') im1 = NX1 - id_per 00239 IF (cd_per == 'R') im1 = 1 00240 ENDIF 00241 gradient_ij2 = zero 00242 IF (sou_mask(ip1,jm1) /= il_maskval .OR. 00243 $ sou_mask(im1,jm1) /= il_maskval) THEN 00244 IF (sou_mask(ip1,jm1) == il_maskval .AND. 00245 $ sou_mask(i,jm1) /= il_maskval) THEN 00246 ip1 = i 00247 di = one 00248 ELSE IF (sou_mask(im1,jm1) == il_maskval .AND. 00249 $ sou_mask(i,jm1) /= il_maskval) THEN 00250 im1 = i 00251 di = one 00252 ELSE 00253 di = zero 00254 ENDIF 00255 gradient_ij2 = 00256 $ di * (src_array(ip1,jm1) - src_array(im1,jm1)) 00257 ENDIF 00258 00259 IF (gradient_ij1 /= zero .AND. gradient_ij2 /= zero) THEN 00260 gradient_ij(i,j) = dj * (gradient_ij1 - gradient_ij2) 00261 ENDIF 00262 ENDIF 00263 00264 ENDDO 00265 ENDDO 00266 ! 00267 IF (nlogprt .GE. 2) THEN 00268 WRITE (UNIT = nulou,FMT = *)' ' 00269 WRITE (UNIT = nulou,FMT = *)'Leaving routine gradient_bicubic' 00270 WRITE (UNIT = nulou,FMT = *)' ' 00271 CALL FLUSH(nulou) 00272 ENDIF 00273 ! 00274 RETURN 00275 00276 END SUBROUTINE gradient_bicubic
1.7.4 
