Changeset 1086 for LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg.F
- Timestamp:
- Feb 3, 2009, 10:50:31 AM (15 years ago)
- File:
-
- 1 edited
-
LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg.F (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg.F
r524 r1086 3 3 ! 4 4 SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire, 5 . griscal ,iter) 6 5 & griscal ,iter) 6 7 USE filtre 8 7 9 IMPLICIT NONE 8 10 c======================================================================= … … 46 48 #include "dimensions.h" 47 49 #include "paramet.h" 48 #include "parafilt.h"49 50 #include "coefils.h" 50 c 51 INTEGER nlat,nbniv,ifiltre,iter 52 INTEGER i,j,l,k 53 INTEGER iim2,immjm 54 INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil 55 56 REAL champ( iip1,nlat,nbniv) 57 REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs 58 COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus) 59 , , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs) 60 , , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus) 61 REAL eignq(iim), sdd1(iim),sdd2(iim) 51 52 INTEGER nlat,nbniv,ifiltre,iter 53 INTEGER i,j,l,k 54 INTEGER iim2,immjm 55 INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil 56 57 REAL champ( iip1,nlat,nbniv) 58 59 REAL eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim) 62 60 LOGICAL griscal 63 61 INTEGER hemisph, iaire 64 c 62 63 LOGICAL,SAVE :: first=.TRUE. 64 65 REAL, SAVE :: sdd12(iim,4) 66 67 INTEGER, PARAMETER :: type_sddu=1 68 INTEGER, PARAMETER :: type_sddv=2 69 INTEGER, PARAMETER :: type_unsddu=3 70 INTEGER, PARAMETER :: type_unsddv=4 71 72 INTEGER :: sdd1_type, sdd2_type 73 74 IF (first) THEN 75 sdd12(1:iim,type_sddu) = sddu(1:iim) 76 sdd12(1:iim,type_sddv) = sddv(1:iim) 77 sdd12(1:iim,type_unsddu) = unsddu(1:iim) 78 sdd12(1:iim,type_unsddv) = unsddv(1:iim) 79 80 first=.FALSE. 81 ENDIF 65 82 66 83 IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) 67 *STOP'Pas de transformee simple dans cette version'68 84 & STOP'Pas de transformee simple dans cette version' 85 69 86 IF( iter.EQ. 2 ) THEN 70 PRINT *,' Pas d iteration du filtre dans cette version !'71 *, ' Utiliser old_filtreg et repasser !'72 73 ENDIF 74 87 PRINT *,' Pas d iteration du filtre dans cette version !' 88 & , ' Utiliser old_filtreg et repasser !' 89 STOP 90 ENDIF 91 75 92 IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN 76 PRINT *,' Cette routine ne calcule le filtre inverse que ',77 *' sur la grille des scalaires !'78 79 ENDIF 80 93 PRINT *,' Cette routine ne calcule le filtre inverse que ' 94 & , ' sur la grille des scalaires !' 95 STOP 96 ENDIF 97 81 98 IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN 82 PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' 83 *,' corriger et repasser !' 84 STOP 85 ENDIF 86 c 87 99 PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' 100 & , ' corriger et repasser !' 101 STOP 102 ENDIF 103 88 104 iim2 = iim * iim 89 105 immjm = iim * jjm 90 c 91 c 106 92 107 IF( griscal ) THEN 93 108 IF( nlat. NE. jjp1 ) THEN 94 PRINT 1111 95 STOP 96 ELSE 97 c 98 IF( iaire.EQ.1 ) THEN 99 CALL SCOPY( iim, sddv, 1, sdd1, 1 ) 100 CALL SCOPY( iim, unsddv, 1, sdd2, 1 ) 101 ELSE 102 CALL SCOPY( iim, unsddv, 1, sdd1, 1 ) 103 CALL SCOPY( iim, sddv, 1, sdd2, 1 ) 104 END IF 105 c 106 jdfil1 = 2 107 jffil1 = jfiltnu 108 jdfil2 = jfiltsu 109 jffil2 = jjm 110 END IF 109 PRINT 1111 110 STOP 111 ELSE 112 113 IF( iaire.EQ.1 ) THEN 114 sdd1_type = type_sddu 115 sdd2_type = type_unsddu 116 ELSE 117 sdd1_type = type_unsddu 118 sdd2_type = type_sddu 119 ENDIF 120 121 c IF( iaire.EQ.1 ) THEN 122 c CALL SCOPY( iim, sddv, 1, sdd1, 1 ) 123 c CALL SCOPY( iim, unsddv, 1, sdd2, 1 ) 124 c ELSE 125 c CALL SCOPY( iim, unsddv, 1, sdd1, 1 ) 126 c CALL SCOPY( iim, sddv, 1, sdd2, 1 ) 127 c END IF 128 129 jdfil1 = 2 130 jffil1 = jfiltnu 131 jdfil2 = jfiltsu 132 jffil2 = jjm 133 END IF 111 134 ELSE 112 IF( nlat.NE.jjm ) THEN 113 PRINT 2222 114 STOP 115 ELSE 116 c 117 IF( iaire.EQ.1 ) THEN 118 CALL SCOPY( iim, sddu, 1, sdd1, 1 ) 119 CALL SCOPY( iim, unsddu, 1, sdd2, 1 ) 120 ELSE 121 CALL SCOPY( iim, unsddu, 1, sdd1, 1 ) 122 CALL SCOPY( iim, sddu, 1, sdd2, 1 ) 123 END IF 124 c 125 jdfil1 = 1 126 jffil1 = jfiltnv 127 jdfil2 = jfiltsv 128 jffil2 = jjm 129 END IF 135 IF( nlat.NE.jjm ) THEN 136 PRINT 2222 137 STOP 138 ELSE 139 140 IF( iaire.EQ.1 ) THEN 141 sdd1_type = type_sddu 142 sdd2_type = type_unsddu 143 ELSE 144 sdd1_type = type_unsddu 145 sdd2_type = type_sddu 146 ENDIF 147 148 c IF( iaire.EQ.1 ) THEN 149 c CALL SCOPY( iim, sddu, 1, sdd1, 1 ) 150 c CALL SCOPY( iim, unsddu, 1, sdd2, 1 ) 151 c ELSE 152 c CALL SCOPY( iim, unsddu, 1, sdd1, 1 ) 153 c CALL SCOPY( iim, sddu, 1, sdd2, 1 ) 154 c END IF 155 156 jdfil1 = 1 157 jffil1 = jfiltnv 158 jdfil2 = jfiltsv 159 jffil2 = jjm 160 END IF 130 161 END IF 131 c 132 c 133 DO 100 hemisph = 1, 2 134 c 135 IF ( hemisph.EQ.1 ) THEN 136 jdfil = jdfil1 137 jffil = jffil1 138 ELSE 139 jdfil = jdfil2 140 jffil = jffil2 141 END IF 142 143 144 DO 50 l = 1, nbniv 145 DO 30 j = jdfil,jffil 146 147 148 DO 5 i = 1, iim 149 champ(i,j,l) = champ(i,j,l) * sdd1(i) 150 5 CONTINUE 151 c 152 153 IF( hemisph. EQ. 1 ) THEN 154 155 IF( ifiltre. EQ. -2 ) THEN 156 #ifdef CRAY 157 CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , 158 * 1, iim, iim ) 159 #else 160 #ifdef BLAS 161 CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim, 162 . champ(1,j,l), 1, 0.0, eignq, 1) 163 #else 164 DO k = 1, iim 165 eignq(k) = 0.0 162 163 DO hemisph = 1, 2 164 165 IF ( hemisph.EQ.1 ) THEN 166 jdfil = jdfil1 167 jffil = jffil1 168 ELSE 169 jdfil = jdfil2 170 jffil = jffil2 171 END IF 172 173 DO l = 1, nbniv 174 DO j = jdfil,jffil 175 DO i = 1, iim 176 champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i) 177 END DO 178 END DO 179 END DO 180 181 IF( hemisph. EQ. 1 ) THEN 182 183 IF( ifiltre. EQ. -2 ) THEN 184 185 DO j = jdfil,jffil 186 CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 187 & matrinvn(1,1,j), 188 & iim, champ(1,j,1), iip1*nlat, 0.0, 189 & eignq(1,j-jdfil+1,1), iim*nlat) 190 END DO 191 192 ELSE IF ( griscal ) THEN 193 194 DO j = jdfil,jffil 195 CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 196 & matriceun(1,1,j), 197 & iim, champ(1,j,1), iip1*nlat, 0.0, 198 & eignq(1,j-jdfil+1,1), iim*nlat) 199 END DO 200 201 ELSE 202 203 DO j = jdfil,jffil 204 CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 205 & matricevn(1,1,j), 206 & iim, champ(1,j,1), iip1*nlat, 0.0, 207 & eignq(1,j-jdfil+1,1), iim*nlat) 208 END DO 209 210 ENDIF 211 212 ELSE 213 214 IF( ifiltre. EQ. -2 ) THEN 215 216 DO j = jdfil,jffil 217 CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 218 & matrinvs(1,1,j-jfiltsu+1), 219 & iim, champ(1,j,1), iip1*nlat, 0.0, 220 & eignq(1,j-jdfil+1,1), iim*nlat) 221 END DO 222 223 224 ELSE IF ( griscal ) THEN 225 226 DO j = jdfil,jffil 227 CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 228 & matriceus(1,1,j-jfiltsu+1), 229 & iim, champ(1,j,1), iip1*nlat, 0.0, 230 & eignq(1,j-jdfil+1,1), iim*nlat) 231 END DO 232 233 ELSE 234 235 DO j = jdfil,jffil 236 CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 237 & matricevs(1,1,j-jfiltsv+1), 238 & iim, champ(1,j,1), iip1*nlat, 0.0, 239 & eignq(1,j-jdfil+1,1), iim*nlat) 240 END DO 241 242 ENDIF 243 244 ENDIF 245 246 IF( ifiltre.EQ. 2 ) THEN 247 248 DO l = 1, nbniv 249 DO j = jdfil,jffil 250 DO i = 1, iim 251 champ( i,j,l ) = 252 & (champ(i,j,l) + eignq(i,j-jdfil+1,l)) 253 & * sdd12(i,sdd2_type) ! sdd2(i) 254 END DO 255 END DO 256 END DO 257 258 ELSE 259 260 DO l = 1, nbniv 261 DO j = jdfil,jffil 262 DO i = 1, iim 263 champ( i,j,l ) = 264 & (champ(i,j,l) - eignq(i,j-jdfil+1,l)) 265 & * sdd12(i,sdd2_type) ! sdd2(i) 266 END DO 267 END DO 268 END DO 269 270 ENDIF 271 272 DO l = 1, nbniv 273 DO j = jdfil,jffil 274 champ( iip1,j,l ) = champ( 1,j,l ) 275 END DO 276 END DO 277 278 166 279 ENDDO 167 DO k = 1, iim 168 DO i = 1, iim 169 eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l) 170 ENDDO 171 ENDDO 172 #endif 173 #endif 174 ELSE IF ( griscal ) THEN 175 #ifdef CRAY 176 CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq , 177 * 1, iim, iim ) 178 #else 179 #ifdef BLAS 180 CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim, 181 . champ(1,j,l), 1, 0.0, eignq, 1) 182 #else 183 DO k = 1, iim 184 eignq(k) = 0.0 185 ENDDO 186 DO i = 1, iim 187 DO k = 1, iim 188 eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l) 189 ENDDO 190 ENDDO 191 #endif 192 #endif 193 ELSE 194 #ifdef CRAY 195 CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , 196 * 1, iim, iim ) 197 #else 198 #ifdef BLAS 199 CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim, 200 . champ(1,j,l), 1, 0.0, eignq, 1) 201 #else 202 DO k = 1, iim 203 eignq(k) = 0.0 204 ENDDO 205 DO i = 1, iim 206 DO k = 1, iim 207 eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l) 208 ENDDO 209 ENDDO 210 #endif 211 #endif 212 ENDIF 213 214 ELSE 215 216 IF( ifiltre. EQ. -2 ) THEN 217 #ifdef CRAY 218 CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , 219 * eignq, 1, iim, iim ) 220 #else 221 #ifdef BLAS 222 CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim, 223 . champ(1,j,l), 1, 0.0, eignq, 1) 224 #else 225 DO k = 1, iim 226 eignq(k) = 0.0 227 ENDDO 228 DO i = 1, iim 229 DO k = 1, iim 230 eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) 231 ENDDO 232 ENDDO 233 #endif 234 #endif 235 ELSE IF ( griscal ) THEN 236 #ifdef CRAY 237 CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , 238 * eignq, 1, iim, iim ) 239 #else 240 #ifdef BLAS 241 CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim, 242 . champ(1,j,l), 1, 0.0, eignq, 1) 243 #else 244 DO k = 1, iim 245 eignq(k) = 0.0 246 ENDDO 247 DO i = 1, iim 248 DO k = 1, iim 249 eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) 250 ENDDO 251 ENDDO 252 #endif 253 #endif 254 ELSE 255 #ifdef CRAY 256 CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 , 257 * eignq, 1, iim, iim ) 258 #else 259 #ifdef BLAS 260 CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim, 261 . champ(1,j,l), 1, 0.0, eignq, 1) 262 #else 263 DO k = 1, iim 264 eignq(k) = 0.0 265 ENDDO 266 DO i = 1, iim 267 DO k = 1, iim 268 eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) 269 ENDDO 270 ENDDO 271 #endif 272 #endif 273 ENDIF 274 275 ENDIF 276 c 277 IF( ifiltre.EQ. 2 ) THEN 278 DO 15 i = 1, iim 279 champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i) 280 15 CONTINUE 281 ELSE 282 DO 16 i=1,iim 283 champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i) 284 16 CONTINUE 285 ENDIF 286 c 287 champ( iip1,j,l ) = champ( 1,j,l ) 288 c 289 30 CONTINUE 290 c 291 50 CONTINUE 292 c 293 100 CONTINUE 294 c 280 295 281 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a 296 *filtrer, sur la grille des scalaires'/)282 & filtrer, sur la grille des scalaires'/) 297 283 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi 298 *ltrer, sur la grille de V ou de Z'/)284 & ltrer, sur la grille de V ou de Z'/) 299 285 RETURN 300 286 END
Note: See TracChangeset
for help on using the changeset viewer.
