|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 subroutine gradient(NX1, NY1, src_array, sou_mask, 00002 $ src_latitudes, src_longitudes, 00003 $ id_per, cd_per, 00004 $ grad_lat, grad_lon) 00005 00006 C**** 00007 C ***************************** 00008 C * OASIS ROUTINE - LEVEL ? * 00009 C * ------------- ------- * 00010 C ***************************** 00011 C 00012 C**** *gradient* - calculate gradients for conservative remapping 00013 C 00014 C Purpose: 00015 C ------- 00016 C Calculation of gradients in latitudinal and longitudinal direction. 00017 C In a first step the gradients in direction of source-grid rows 00018 C and lines are calculated. Then they are rotated to longitudinal 00019 C and latitudinal direction, using the scalar product. 00020 C This routine works for logically rectangular grids, only. 00021 C 00022 C** Interface: 00023 C --------- 00024 C *CALL* *gradient*(NX1, NY1, src_array, sou_mask, src_latitudes, 00025 C $ src_longitudes, grad_lat, grad_lon) 00026 C 00027 C Input: 00028 C ----- 00029 C NX1 : grid dimension in x-direction (integer) 00030 C NY1 : grid dimension in y-direction (integer) 00031 C src_array : array on source grid (real 2D) 00032 C sou_mask : source grid mask (integer 2D) 00033 C src_latitudes : latitudes on source grid (real 2D) 00034 C src_longitudes : longitudes on source grid (real 2D) 00035 C id_per : number of overlapping points for source grid 00036 C cd_per : grip periodicity type 00037 C 00038 C Output: 00039 C ------ 00040 C grad_lat : gradient in latitudinal direction (real 2D) 00041 C grad_lon : gradient in longitudinal direction (real 2D) 00042 C 00043 C History: 00044 C ------- 00045 C Version Programmer Date Description 00046 C ------- ---------- ---- ----------- 00047 C 2.5 V. Gayler 2001/09/20 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 dimensions 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 $ grad_lat, ! gradient in latitudinal direction 00081 $ grad_lon ! gradient in longitudinal direction 00082 00083 !----------------------------------------------------------------------- 00084 ! LOCAL VARIABLES 00085 !----------------------------------------------------------------------- 00086 INTEGER (int_kind) :: 00087 $ i, j, ! looping indicees 00088 $ ip1, jp1, im1, jm1 00089 00090 REAL (kind = real_kind) :: 00091 $ distance ! distance in rad 00092 00093 REAL (kind = real_kind) :: 00094 $ dVar_i, dVar_j, ! difference of Var in i / j direction 00095 $ dlat_i, dlat_j, ! difference in lat in i / j direction 00096 $ dlon_i, dlon_j, ! difference in lon in i / j direction 00097 $ dist_i, dist_j, ! distance in i / j direction 00098 $ grad_i, grad_j, ! gradient in i / j direction 00099 $ ABSold, ABSnew, lat_factor 00100 00101 REAL (kind = real_kind), DIMENSION(NX1,NY1) :: 00102 $ src_lon, ! source grid longitudes [radiants] 00103 $ src_lat, ! source grid latitudes [radiants] 00104 $ pi180 ! conversion factor: deg -> rad 00105 00106 INTEGER (int_kind) :: il_maskval= 0_int_kind 00107 00108 !----------------------------------------------------------------------- 00109 ! 00110 IF (nlogprt .GE. 2) THEN 00111 WRITE (UNIT = nulou,FMT = *)' ' 00112 WRITE (UNIT = nulou,FMT = *)' Entering routine gradient ' 00113 WRITE (UNIT = nulou,FMT = *)' ' 00114 CALL FLUSH(nulou) 00115 ENDIF 00116 ! 00117 ! Transformation from degree to radiant 00118 ! ------------------------------------- 00119 pi180 = 1.74532925199432957692e-2 ! =PI/180 00120 00121 src_lon = src_longitudes * pi180 00122 src_lat = src_latitudes * pi180 00123 00124 !----------------------------------------------------------------------- 00125 00126 DO i = 1, NX1 00127 00128 DO j = 1, NY1 00129 00130 IF (sou_mask(i,j) /= il_maskval) THEN 00131 00132 ip1 = i + 1 00133 im1 = i - 1 00134 IF (i == NX1) THEN 00135 IF (cd_per == 'P') ip1 = 1 + id_per ! the 0-meridian 00136 IF (cd_per == 'R') ip1 = NX1 00137 ENDIF 00138 IF (i == 1 ) THEN 00139 IF (cd_per == 'P') im1 = NX1 - id_per 00140 IF (cd_per == 'R') im1 = 1 00141 ENDIF 00142 jp1 = j + 1 00143 jm1 = j - 1 00144 IF (j == NY1) jp1 = NY1 ! treatment of the last.. 00145 IF (j == 1 ) jm1 = 1 ! .. and the first grid-row 00146 00147 IF (sou_mask(ip1,j) == il_maskval) ip1 = i 00148 IF (sou_mask(im1,j) == il_maskval) im1 = i 00149 IF (sou_mask(i,jp1) == il_maskval) jp1 = j 00150 IF (sou_mask(i,jm1) == il_maskval) jm1 = j 00151 00152 ! difference between neighbouring datapoints 00153 dVar_i = src_array(ip1,j) - src_array(im1,j) 00154 dVar_j = src_array(i,jp1) - src_array(i,jm1) 00155 00156 ! difference in latitudes 00157 dlat_i = src_lat(ip1,j) - src_lat(im1,j) 00158 dlat_j = src_lat(i,jp1) - src_lat(i,jm1) 00159 00160 ! difference in longitudes 00161 dlon_i = src_lon(ip1,j) - src_lon(im1,j) 00162 IF (dlon_i > PI) dlon_i = dlon_i - PI2 00163 IF (dlon_i < (-PI)) dlon_i = dlon_i + PI2 00164 dlon_j = src_lon(i,jp1) - src_lon(i,jm1) 00165 IF (dlon_j > PI) dlon_j = dlon_j - PI2 00166 IF (dlon_j < (-PI)) dlon_j = dlon_j + PI2 00167 lat_factor = cos(src_lat(i,j)) 00168 dlon_i = dlon_i * lat_factor 00169 dlon_j = dlon_j * lat_factor 00170 00171 ! distance 00172 dist_i = distance(src_lon(ip1,j), src_lat(ip1,j), 00173 $ src_lon(im1,j), src_lat(im1,j)) 00174 dist_j = distance(src_lon(i,jp1), src_lat(i,jp1), 00175 $ src_lon(i,jm1), src_lat(i,jm1)) 00176 00177 ! gradients: dVar / distance (= vector lenght) 00178 IF (dist_i /= 0.) THEN 00179 grad_i = dVar_i / dist_i 00180 ELSE 00181 grad_i = 0 00182 ENDIF 00183 IF (dist_j /= 0.) THEN 00184 grad_j = dVar_j / dist_j 00185 ELSE 00186 grad_j = 0 00187 ENDIF 00188 00189 ! projection by scalar product 00190 ! ---------------------------- 00191 grad_lon(i,j) = grad_i * dlon_i + grad_j * dlat_i 00192 grad_lat(i,j) = grad_i * dlon_j + grad_j * dlat_j 00193 00194 if (dist_i /= 0) then 00195 grad_lon(i,j) = grad_lon(i,j) / dist_i 00196 else 00197 grad_lon(i,j) = 0 00198 endif 00199 if (dist_j /= 0) then 00200 grad_lat(i,j) = grad_lat(i,j) / dist_j 00201 else 00202 grad_lat(i,j) = 0. 00203 endif 00204 00205 ! correct skale 00206 ! ------------- 00207 ABSold = SQRT(grad_i**2 + grad_j**2) 00208 ABSnew = SQRT(grad_lon(i,j)**2 + grad_lat(i,j)**2) 00209 IF (ABSnew > 1.E-10) THEN 00210 C grad_lon(i,j) = grad_lon(i,j)*ABSold/ABSnew 00211 grad_lon(i,j) = grad_lon(i,j) 00212 ELSE 00213 grad_lon(i,j) = 0.0 00214 ENDIF 00215 00216 ! test orthogonality 00217 ! ------------------ 00218 IF ((dlon_i*dlon_j+dlat_j*dlat_i) > 0.1) THEN 00219 print*, 'ORTHOGONAL? ', i, j, 00220 $ (dlon_i*dlon_j+dlat_j*dlat_i) 00221 ENDIF 00222 00223 ELSE 00224 00225 grad_lat(i,j) = 0. 00226 grad_lon(i,j) = 0. 00227 00228 ENDIF 00229 00230 ENDDO 00231 ENDDO 00232 IF (nlogprt .GE. 2) THEN 00233 WRITE (UNIT = nulou,FMT = *)' ' 00234 WRITE (UNIT = nulou,FMT = *)' Leaving routine gradient ' 00235 WRITE (UNIT = nulou,FMT = *)' ' 00236 CALL FLUSH(nulou) 00237 ENDIF 00238 RETURN 00239 00240 END SUBROUTINE gradient 00241 00242 00243 00244 00245 00246 00247
1.7.4 
