Changeset 5246 for LMDZ6/trunk/libf/dyn3d_common/extrapol.f90
- Timestamp:
- Oct 21, 2024, 2:58:45 PM (23 hours ago)
- File:
-
- 1 moved
-
LMDZ6/trunk/libf/dyn3d_common/extrapol.f90 (moved) (moved from LMDZ6/trunk/libf/dyn3d_common/extrapol.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3d_common/extrapol.f90
r5245 r5246 2 2 ! $Id$ 3 3 ! 4 C 5 C 6 SUBROUTINE extrapol (pfild, kxlon, kylat, pmask, 7 .norsud, ldper, knbor, pwork)8 9 c 10 cOASIS routine (Adaptation: Laurent Li, le 14 mars 1997)11 cFill up missed values by using the neighbor points12 c 13 INTEGERkxlon, kylat ! longitude and latitude dimensions (Input)14 INTEGERknbor ! minimum neighbor number (Input)15 LOGICALnorsud ! True if field is from North to South (Input)16 LOGICALldper ! True if take into account the periodicity (Input)17 REALpmask ! mask value (Input)18 REALpfild(kxlon,kylat) ! field to be extrapolated (Input/Output)19 REALpwork(kxlon,kylat) ! working space20 c 21 REALzwmsk22 INTEGERincre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat23 INTEGERix(9), jy(9) ! index arrays for the neighbors coordinates24 REALzmask(9)25 C 26 CWe search over the eight closest neighbors27 C 28 Cj+1 7 8 929 Cj 4 5 6 Current point 5 --> (i,j)30 Cj-1 1 2 331 Ci-1 i i+132 c 33 c 34 35 36 37 38 39 40 41 42 43 44 45 46 c 47 48 c 49 50 51 52 53 54 c 55 C* To avoid problems in floating point tests56 57 c 4 ! 5 ! 6 SUBROUTINE extrapol (pfild, kxlon, kylat, pmask, & 7 norsud, ldper, knbor, pwork) 8 IMPLICIT none 9 ! 10 ! OASIS routine (Adaptation: Laurent Li, le 14 mars 1997) 11 ! Fill up missed values by using the neighbor points 12 ! 13 INTEGER :: kxlon, kylat ! longitude and latitude dimensions (Input) 14 INTEGER :: knbor ! minimum neighbor number (Input) 15 LOGICAL :: norsud ! True if field is from North to South (Input) 16 LOGICAL :: ldper ! True if take into account the periodicity (Input) 17 REAL :: pmask ! mask value (Input) 18 REAL :: pfild(kxlon,kylat) ! field to be extrapolated (Input/Output) 19 REAL :: pwork(kxlon,kylat) ! working space 20 ! 21 REAL :: zwmsk 22 INTEGER :: incre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat 23 INTEGER :: ix(9), jy(9) ! index arrays for the neighbors coordinates 24 REAL :: zmask(9) 25 ! 26 ! We search over the eight closest neighbors 27 ! 28 ! j+1 7 8 9 29 ! j 4 5 6 Current point 5 --> (i,j) 30 ! j-1 1 2 3 31 ! i-1 i i+1 32 ! 33 ! 34 IF (norsud) THEN 35 DO j = 1, kylat 36 DO i = 1, kxlon 37 pwork(i,j) = pfild(i,kylat-j+1) 38 ENDDO 39 ENDDO 40 DO j = 1, kylat 41 DO i = 1, kxlon 42 pfild(i,j) = pwork(i,j) 43 ENDDO 44 ENDDO 45 ENDIF 46 ! 47 incre = 0 48 ! 49 DO j = 1, kylat 50 DO i = 1, kxlon 51 pwork(i,j) = pfild(i,j) 52 ENDDO 53 ENDDO 54 ! 55 !* To avoid problems in floating point tests 56 zwmsk = pmask - 1.0 57 ! 58 58 200 CONTINUE 59 incre = incre + 1 60 DO 99999 j = 1, kylat 61 DO 99999 i = 1, kxlon 62 IF (pfild(i,j).GT. zwmsk) THEN 63 pwork(i,j) = pfild(i,j) 64 inbor = 0 65 ideb = 1 66 ifin = 9 67 C 68 C* Fill up ix array 69 ix(1) = MAX (1,i-1) 70 ix(2) = i 71 ix(3) = MIN (kxlon,i+1) 72 ix(4) = MAX (1,i-1) 73 ix(5) = i 74 ix(6) = MIN (kxlon,i+1) 75 ix(7) = MAX (1,i-1) 76 ix(8) = i 77 ix(9) = MIN (kxlon,i+1) 78 C 79 C* Fill up iy array 80 jy(1) = MAX (1,j-1) 81 jy(2) = MAX (1,j-1) 82 jy(3) = MAX (1,j-1) 83 jy(4) = j 84 jy(5) = j 85 jy(6) = j 86 jy(7) = MIN (kylat,j+1) 87 jy(8) = MIN (kylat,j+1) 88 jy(9) = MIN (kylat,j+1) 89 C 90 C* Correct latitude bounds if southernmost or northernmost points 91 IF (j .EQ. 1) ideb = 4 92 IF (j .EQ. kylat) ifin = 6 93 C 94 C* Account for periodicity in longitude 95 C 96 IF (ldper) THEN 97 IF (i .EQ. kxlon) THEN 98 ix(3) = 1 99 ix(6) = 1 100 ix(9) = 1 101 ELSE IF (i .EQ. 1) THEN 102 ix(1) = kxlon 103 ix(4) = kxlon 104 ix(7) = kxlon 105 ENDIF 106 ELSE 107 IF (i .EQ. 1) THEN 108 ix(1) = i 109 ix(2) = i + 1 110 ix(3) = i 111 ix(4) = i + 1 112 ix(5) = i 113 ix(6) = i + 1 114 ENDIF 115 IF (i .EQ. kxlon) THEN 116 ix(1) = i -1 117 ix(2) = i 118 ix(3) = i - 1 119 ix(4) = i 120 ix(5) = i - 1 121 ix(6) = i 122 ENDIF 123 C 124 IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN 125 jy(1) = MAX (1,j-1) 126 jy(2) = MAX (1,j-1) 127 jy(3) = j 128 jy(4) = j 129 jy(5) = MIN (kylat,j+1) 130 jy(6) = MIN (kylat,j+1) 131 C 132 ideb = 1 133 ifin = 6 134 IF (j .EQ. 1) ideb = 3 135 IF (j .EQ. kylat) ifin = 4 136 ENDIF 137 ENDIF ! end for ldper test 138 C 139 C* Find unmasked neighbors 140 C 141 DO 230 k = ideb, ifin 142 zmask(k) = 0. 143 ilon = ix(k) 144 jlat = jy(k) 145 IF (pfild(ilon,jlat) .LT. zwmsk) THEN 146 zmask(k) = 1. 147 inbor = inbor + 1 148 ENDIF 149 230 CONTINUE 150 C 151 C* Not enough points around point P are unmasked; interpolation on P 152 C will be done in a future call to extrap. 153 C 154 IF (inbor .GE. knbor) THEN 155 pwork(i,j) = 0. 156 DO k = ideb, ifin 157 ilon = ix(k) 158 jlat = jy(k) 159 pwork(i,j) = pwork(i,j) 160 $ + pfild(ilon,jlat) * zmask(k)/ REAL(inbor) 161 ENDDO 162 ENDIF 163 C 164 ENDIF 165 99999 CONTINUE 166 C 167 C* 3. Writing back unmasked field in pfild 168 C ------------------------------------ 169 C 170 C* pfild then contains: 171 C - Values which were not masked 172 C - Interpolated values from the inbor neighbors 173 C - Values which are not yet interpolated 174 C 175 idoit = 0 176 DO j = 1, kylat 177 DO i = 1, kxlon 178 IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1 179 pfild(i,j) = pwork(i,j) 180 ENDDO 181 ENDDO 182 c 183 IF (idoit .ne. 0) GOTO 200 184 ccc PRINT*, "Number of extrapolation steps incre =", incre 185 c 186 IF (norsud) THEN 187 DO j = 1, kylat 188 DO i = 1, kxlon 189 pwork(i,j) = pfild(i,kylat-j+1) 190 ENDDO 191 ENDDO 192 DO j = 1, kylat 193 DO i = 1, kxlon 194 pfild(i,j) = pwork(i,j) 195 ENDDO 196 ENDDO 197 ENDIF 198 c 199 RETURN 200 END 59 incre = incre + 1 60 DO j = 1, kylat 61 DO i = 1, kxlon 62 IF (pfild(i,j).GT. zwmsk) THEN 63 pwork(i,j) = pfild(i,j) 64 inbor = 0 65 ideb = 1 66 ifin = 9 67 ! 68 !* Fill up ix array 69 ix(1) = MAX (1,i-1) 70 ix(2) = i 71 ix(3) = MIN (kxlon,i+1) 72 ix(4) = MAX (1,i-1) 73 ix(5) = i 74 ix(6) = MIN (kxlon,i+1) 75 ix(7) = MAX (1,i-1) 76 ix(8) = i 77 ix(9) = MIN (kxlon,i+1) 78 ! 79 !* Fill up iy array 80 jy(1) = MAX (1,j-1) 81 jy(2) = MAX (1,j-1) 82 jy(3) = MAX (1,j-1) 83 jy(4) = j 84 jy(5) = j 85 jy(6) = j 86 jy(7) = MIN (kylat,j+1) 87 jy(8) = MIN (kylat,j+1) 88 jy(9) = MIN (kylat,j+1) 89 ! 90 !* Correct latitude bounds if southernmost or northernmost points 91 IF (j .EQ. 1) ideb = 4 92 IF (j .EQ. kylat) ifin = 6 93 ! 94 !* Account for periodicity in longitude 95 ! 96 IF (ldper) THEN 97 IF (i .EQ. kxlon) THEN 98 ix(3) = 1 99 ix(6) = 1 100 ix(9) = 1 101 ELSE IF (i .EQ. 1) THEN 102 ix(1) = kxlon 103 ix(4) = kxlon 104 ix(7) = kxlon 105 ENDIF 106 ELSE 107 IF (i .EQ. 1) THEN 108 ix(1) = i 109 ix(2) = i + 1 110 ix(3) = i 111 ix(4) = i + 1 112 ix(5) = i 113 ix(6) = i + 1 114 ENDIF 115 IF (i .EQ. kxlon) THEN 116 ix(1) = i -1 117 ix(2) = i 118 ix(3) = i - 1 119 ix(4) = i 120 ix(5) = i - 1 121 ix(6) = i 122 ENDIF 123 ! 124 IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN 125 jy(1) = MAX (1,j-1) 126 jy(2) = MAX (1,j-1) 127 jy(3) = j 128 jy(4) = j 129 jy(5) = MIN (kylat,j+1) 130 jy(6) = MIN (kylat,j+1) 131 ! 132 ideb = 1 133 ifin = 6 134 IF (j .EQ. 1) ideb = 3 135 IF (j .EQ. kylat) ifin = 4 136 ENDIF 137 ENDIF ! end for ldper test 138 ! 139 !* Find unmasked neighbors 140 ! 141 DO k = ideb, ifin 142 zmask(k) = 0. 143 ilon = ix(k) 144 jlat = jy(k) 145 IF (pfild(ilon,jlat) .LT. zwmsk) THEN 146 zmask(k) = 1. 147 inbor = inbor + 1 148 ENDIF 149 END DO 150 ! 151 !* Not enough points around point P are unmasked; interpolation on P 152 ! will be done in a future call to extrap. 153 ! 154 IF (inbor .GE. knbor) THEN 155 pwork(i,j) = 0. 156 DO k = ideb, ifin 157 ilon = ix(k) 158 jlat = jy(k) 159 pwork(i,j) = pwork(i,j) & 160 + pfild(ilon,jlat) * zmask(k)/ REAL(inbor) 161 ENDDO 162 ENDIF 163 ! 164 ENDIF 165 END DO 166 END DO 167 ! 168 !* 3. Writing back unmasked field in pfild 169 ! ------------------------------------ 170 ! 171 !* pfild then contains: 172 ! - Values which were not masked 173 ! - Interpolated values from the inbor neighbors 174 ! - Values which are not yet interpolated 175 ! 176 idoit = 0 177 DO j = 1, kylat 178 DO i = 1, kxlon 179 IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1 180 pfild(i,j) = pwork(i,j) 181 ENDDO 182 ENDDO 183 ! 184 IF (idoit .ne. 0) GOTO 200 185 !cc PRINT*, "Number of extrapolation steps incre =", incre 186 ! 187 IF (norsud) THEN 188 DO j = 1, kylat 189 DO i = 1, kxlon 190 pwork(i,j) = pfild(i,kylat-j+1) 191 ENDDO 192 ENDDO 193 DO j = 1, kylat 194 DO i = 1, kxlon 195 pfild(i,j) = pwork(i,j) 196 ENDDO 197 ENDDO 198 ENDIF 199 ! 200 RETURN 201 END SUBROUTINE extrapol
Note: See TracChangeset
for help on using the changeset viewer.
