Changeset 5246 for LMDZ6/trunk/libf/filtrez/filtreg.F90
- Timestamp:
- Oct 21, 2024, 2:58:45 PM (4 days ago)
- File:
-
- 1 moved
-
LMDZ6/trunk/libf/filtrez/filtreg.F90 (moved) (moved from LMDZ6/trunk/libf/filtrez/filtreg.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/filtrez/filtreg.F90
r5245 r5246 2 2 ! $Header$ 3 3 ! 4 SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire, 5 &griscal ,iter)6 7 8 9 10 c=======================================================================11 c 12 cAuteur: P. Le Van 07/10/9713 c------14 c 15 cObjet: filtre matriciel longitudinal ,avec les matrices precalculees16 cpour l'operateur Filtre .17 c------18 c 19 cArguments:20 c----------21 c 22 cnblat nombre de latitudes a filtrer23 cnbniv nombre de niveaux verticaux a filtrer24 cchamp(iip1,nblat,nbniv) en entree : champ a filtrer25 cen sortie : champ filtre26 cifiltre +1 Transformee directe27 c-1 Transformee inverse28 c+2 Filtre directe29 c-2 Filtre inverse30 c 31 ciaire 1 si champ intensif32 c2 si champ extensif (pondere par les aires)33 c 34 citer 1 filtre simple35 c 36 c=======================================================================37 c 38 c 39 cVariable Intensive40 cifiltre = 1 filtre directe41 cifiltre =-1 filtre inverse42 c 43 cVariable Extensive44 cifiltre = 2 filtre directe45 cifiltre =-2 filtre inverse46 c 47 c 48 49 50 51 52 INTEGERnlat,nbniv,ifiltre,iter53 INTEGERi,j,l,k54 INTEGERiim2,immjm55 INTEGERjdfil1,jdfil2,jffil1,jffil2,jdfil,jffil56 57 REALchamp( iip1,nlat,nbniv)58 59 REALeignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)60 LOGICALgriscal61 INTEGERhemisph, iaire62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)86 & STOP'Pas de transformee simple dans cette version'87 88 89 PRINT *,' Pas d iteration du filtre dans cette version !'90 & , ' Utiliser old_filtreg et repasser !'91 92 93 94 95 PRINT *,' Cette routine ne calcule le filtre inverse que '96 &, ' sur la grille des scalaires !'97 98 99 100 101 PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'102 &, ' corriger et repasser !'103 104 105 106 107 108 109 110 IF( nlat.NE. jjp1 ) THEN111 112 113 114 115 116 117 118 119 120 121 122 123 cIF( iaire.EQ.1 ) THEN124 c CALL SCOPY( iim, sddv, 1, sdd1, 1 ) 125 cCALL SCOPY( iim, unsddv, 1, sdd2, 1 )126 cELSE127 cCALL SCOPY( iim, unsddv, 1, sdd1, 1 )128 cCALL SCOPY( iim, sddv, 1, sdd2, 1 )129 cEND IF130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 cIF( iaire.EQ.1 ) THEN151 c CALL SCOPY( iim, sddu, 1, sdd1, 1 ) 152 cCALL SCOPY( iim, unsddu, 1, sdd2, 1 )153 cELSE154 cCALL SCOPY( iim, unsddu, 1, sdd1, 1 )155 cCALL SCOPY( iim, sddu, 1, sdd2, 1 )156 cEND IF157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 IF( hemisph.EQ. 1 ) THEN184 185 IF( ifiltre.EQ. -2 ) THEN186 187 188 #ifdef BLAS 189 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,190 & matrinvn(1,1,j),191 & iim, champ(1,j,1), iip1*nlat, 0.0,192 &eignq(1,j-jdfil+1,1), iim*nlat)193 #else 194 eignq(:,j-jdfil+1,:)195 $= matmul(matrinvn(:,:,j), champ(:iim,j,:))196 #endif 197 198 199 200 201 202 #ifdef BLAS 203 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,204 & matriceun(1,1,j),205 & iim, champ(1,j,1), iip1*nlat, 0.0,206 &eignq(1,j-jdfil+1,1), iim*nlat)207 #else 208 eignq(:,j-jdfil+1,:)209 $= matmul(matriceun(:,:,j), champ(:iim,j,:))210 #endif 211 212 213 ELSE214 215 216 #ifdef BLAS 217 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,218 & matricevn(1,1,j),219 & iim, champ(1,j,1), iip1*nlat, 0.0,220 &eignq(1,j-jdfil+1,1), iim*nlat)221 #else 222 eignq(:,j-jdfil+1,:)223 $= matmul(matricevn(:,:,j), champ(:iim,j,:))224 #endif 225 226 227 228 229 230 231 IF( ifiltre.EQ. -2 ) THEN232 233 234 #ifdef BLAS 235 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,236 & matrinvs(1,1,j-jfiltsu+1),237 & iim, champ(1,j,1), iip1*nlat, 0.0,238 &eignq(1,j-jdfil+1,1), iim*nlat)239 #else 240 eignq(:,j-jdfil+1,:)241 $ = matmul(matrinvs(:,:,j-jfiltsu+1),242 $champ(:iim,j,:))243 #endif 244 245 246 247 248 249 250 #ifdef BLAS 251 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,252 & matriceus(1,1,j-jfiltsu+1),253 & iim, champ(1,j,1), iip1*nlat, 0.0,254 &eignq(1,j-jdfil+1,1), iim*nlat)255 #else 256 eignq(:,j-jdfil+1,:)257 $ = matmul(matriceus(:,:,j-jfiltsu+1),258 $champ(:iim,j,:))259 #endif 260 261 262 ELSE263 264 265 #ifdef BLAS 266 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,267 & matricevs(1,1,j-jfiltsv+1),268 & iim, champ(1,j,1), iip1*nlat, 0.0,269 &eignq(1,j-jdfil+1,1), iim*nlat)270 #else 271 eignq(:,j-jdfil+1,:)272 $ = matmul(matricevs(:,:,j-jfiltsv+1),273 $champ(:iim,j,:))274 #endif 275 276 277 278 279 280 281 282 283 284 285 286 champ( i,j,l ) =287 & (champ(i,j,l) + eignq(i,j-jdfil+1,l))288 &* sdd12(i,sdd2_type) ! sdd2(i)289 290 291 292 293 294 295 296 297 298 champ( i,j,l ) =299 & (champ(i,j,l) - eignq(i,j-jdfil+1,l))300 &* sdd12(i,sdd2_type) ! sdd2(i)301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a317 & filtrer, sur la grille des scalaires'/)318 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi319 & ltrer, sur la grille de V ou de Z'/)320 321 END 4 SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire, & 5 griscal ,iter) 6 7 USE filtreg_mod 8 9 IMPLICIT NONE 10 !======================================================================= 11 ! 12 ! Auteur: P. Le Van 07/10/97 13 ! ------ 14 ! 15 ! Objet: filtre matriciel longitudinal ,avec les matrices precalculees 16 ! pour l'operateur Filtre . 17 ! ------ 18 ! 19 ! Arguments: 20 ! ---------- 21 ! 22 ! nblat nombre de latitudes a filtrer 23 ! nbniv nombre de niveaux verticaux a filtrer 24 ! champ(iip1,nblat,nbniv) en entree : champ a filtrer 25 ! en sortie : champ filtre 26 ! ifiltre +1 Transformee directe 27 ! -1 Transformee inverse 28 ! +2 Filtre directe 29 ! -2 Filtre inverse 30 ! 31 ! iaire 1 si champ intensif 32 ! 2 si champ extensif (pondere par les aires) 33 ! 34 ! iter 1 filtre simple 35 ! 36 !======================================================================= 37 ! 38 ! 39 ! Variable Intensive 40 ! ifiltre = 1 filtre directe 41 ! ifiltre =-1 filtre inverse 42 ! 43 ! Variable Extensive 44 ! ifiltre = 2 filtre directe 45 ! ifiltre =-2 filtre inverse 46 ! 47 ! 48 INCLUDE "dimensions.h" 49 INCLUDE "paramet.h" 50 INCLUDE "coefils.h" 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) 60 LOGICAL :: griscal 61 INTEGER :: hemisph, iaire 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 ( iim == 1 ) return ! no filtre in 2D y-z 75 76 IF (first) THEN 77 sdd12(1:iim,type_sddu) = sddu(1:iim) 78 sdd12(1:iim,type_sddv) = sddv(1:iim) 79 sdd12(1:iim,type_unsddu) = unsddu(1:iim) 80 sdd12(1:iim,type_unsddv) = unsddv(1:iim) 81 82 first=.FALSE. 83 ENDIF 84 85 IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) & 86 STOP 'Pas de transformee simple dans cette version' 87 88 IF( iter.EQ. 2 ) THEN 89 PRINT *,' Pas d iteration du filtre dans cette version !'& 90 & , ' Utiliser old_filtreg et repasser !' 91 STOP 92 ENDIF 93 94 IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN 95 PRINT *,' Cette routine ne calcule le filtre inverse que ' & 96 , ' sur la grille des scalaires !' 97 STOP 98 ENDIF 99 100 IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN 101 PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' & 102 , ' corriger et repasser !' 103 STOP 104 ENDIF 105 106 iim2 = iim * iim 107 immjm = iim * jjm 108 109 IF( griscal ) THEN 110 IF( nlat.NE. jjp1 ) THEN 111 PRINT 1111 112 STOP 113 ELSE 114 115 IF( iaire.EQ.1 ) THEN 116 sdd1_type = type_sddv 117 sdd2_type = type_unsddv 118 ELSE 119 sdd1_type = type_unsddv 120 sdd2_type = type_sddv 121 ENDIF 122 123 ! IF( iaire.EQ.1 ) THEN 124 ! CALL SCOPY( iim, sddv, 1, sdd1, 1 ) 125 ! CALL SCOPY( iim, unsddv, 1, sdd2, 1 ) 126 ! ELSE 127 ! CALL SCOPY( iim, unsddv, 1, sdd1, 1 ) 128 ! CALL SCOPY( iim, sddv, 1, sdd2, 1 ) 129 ! END IF 130 131 jdfil1 = 2 132 jffil1 = jfiltnu 133 jdfil2 = jfiltsu 134 jffil2 = jjm 135 END IF 136 ELSE 137 IF( nlat.NE.jjm ) THEN 138 PRINT 2222 139 STOP 140 ELSE 141 142 IF( iaire.EQ.1 ) THEN 143 sdd1_type = type_sddu 144 sdd2_type = type_unsddu 145 ELSE 146 sdd1_type = type_unsddu 147 sdd2_type = type_sddu 148 ENDIF 149 150 ! IF( iaire.EQ.1 ) THEN 151 ! CALL SCOPY( iim, sddu, 1, sdd1, 1 ) 152 ! CALL SCOPY( iim, unsddu, 1, sdd2, 1 ) 153 ! ELSE 154 ! CALL SCOPY( iim, unsddu, 1, sdd1, 1 ) 155 ! CALL SCOPY( iim, sddu, 1, sdd2, 1 ) 156 ! END IF 157 158 jdfil1 = 1 159 jffil1 = jfiltnv 160 jdfil2 = jfiltsv 161 jffil2 = jjm 162 END IF 163 END IF 164 165 DO hemisph = 1, 2 166 167 IF ( hemisph.EQ.1 ) THEN 168 jdfil = jdfil1 169 jffil = jffil1 170 ELSE 171 jdfil = jdfil2 172 jffil = jffil2 173 END IF 174 175 DO l = 1, nbniv 176 DO j = jdfil,jffil 177 DO i = 1, iim 178 champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i) 179 END DO 180 END DO 181 END DO 182 183 IF( hemisph.EQ. 1 ) THEN 184 185 IF( ifiltre.EQ. -2 ) THEN 186 187 DO j = jdfil,jffil 188 #ifdef BLAS 189 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & 190 matrinvn(1,1,j), & 191 iim, champ(1,j,1), iip1*nlat, 0.0, & 192 eignq(1,j-jdfil+1,1), iim*nlat) 193 #else 194 eignq(:,j-jdfil+1,:) & 195 = matmul(matrinvn(:,:,j), champ(:iim,j,:)) 196 #endif 197 END DO 198 199 ELSE IF ( griscal ) THEN 200 201 DO j = jdfil,jffil 202 #ifdef BLAS 203 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & 204 matriceun(1,1,j), & 205 iim, champ(1,j,1), iip1*nlat, 0.0, & 206 eignq(1,j-jdfil+1,1), iim*nlat) 207 #else 208 eignq(:,j-jdfil+1,:) & 209 = matmul(matriceun(:,:,j), champ(:iim,j,:)) 210 #endif 211 END DO 212 213 ELSE 214 215 DO j = jdfil,jffil 216 #ifdef BLAS 217 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & 218 matricevn(1,1,j), & 219 iim, champ(1,j,1), iip1*nlat, 0.0, & 220 eignq(1,j-jdfil+1,1), iim*nlat) 221 #else 222 eignq(:,j-jdfil+1,:) & 223 = matmul(matricevn(:,:,j), champ(:iim,j,:)) 224 #endif 225 END DO 226 227 ENDIF 228 229 ELSE 230 231 IF( ifiltre.EQ. -2 ) THEN 232 233 DO j = jdfil,jffil 234 #ifdef BLAS 235 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & 236 matrinvs(1,1,j-jfiltsu+1), & 237 iim, champ(1,j,1), iip1*nlat, 0.0, & 238 eignq(1,j-jdfil+1,1), iim*nlat) 239 #else 240 eignq(:,j-jdfil+1,:) & 241 = matmul(matrinvs(:,:,j-jfiltsu+1), & 242 champ(:iim,j,:)) 243 #endif 244 END DO 245 246 247 ELSE IF ( griscal ) THEN 248 249 DO j = jdfil,jffil 250 #ifdef BLAS 251 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & 252 matriceus(1,1,j-jfiltsu+1), & 253 iim, champ(1,j,1), iip1*nlat, 0.0, & 254 eignq(1,j-jdfil+1,1), iim*nlat) 255 #else 256 eignq(:,j-jdfil+1,:) & 257 = matmul(matriceus(:,:,j-jfiltsu+1), & 258 champ(:iim,j,:)) 259 #endif 260 END DO 261 262 ELSE 263 264 DO j = jdfil,jffil 265 #ifdef BLAS 266 CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, & 267 matricevs(1,1,j-jfiltsv+1), & 268 iim, champ(1,j,1), iip1*nlat, 0.0, & 269 eignq(1,j-jdfil+1,1), iim*nlat) 270 #else 271 eignq(:,j-jdfil+1,:) & 272 = matmul(matricevs(:,:,j-jfiltsv+1), & 273 champ(:iim,j,:)) 274 #endif 275 END DO 276 277 ENDIF 278 279 ENDIF 280 281 IF( ifiltre.EQ. 2 ) THEN 282 283 DO l = 1, nbniv 284 DO j = jdfil,jffil 285 DO i = 1, iim 286 champ( i,j,l ) = & 287 (champ(i,j,l) + eignq(i,j-jdfil+1,l)) & 288 * sdd12(i,sdd2_type) ! sdd2(i) 289 END DO 290 END DO 291 END DO 292 293 ELSE 294 295 DO l = 1, nbniv 296 DO j = jdfil,jffil 297 DO i = 1, iim 298 champ( i,j,l ) = & 299 (champ(i,j,l) - eignq(i,j-jdfil+1,l)) & 300 * sdd12(i,sdd2_type) ! sdd2(i) 301 END DO 302 END DO 303 END DO 304 305 ENDIF 306 307 DO l = 1, nbniv 308 DO j = jdfil,jffil 309 champ( iip1,j,l ) = champ( 1,j,l ) 310 END DO 311 END DO 312 313 314 ENDDO 315 316 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a& 317 & filtrer, sur la grille des scalaires'/) 318 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi& 319 & ltrer, sur la grille de V ou de Z'/) 320 RETURN 321 END SUBROUTINE filtreg
Note: See TracChangeset
for help on using the changeset viewer.
