Changeset 5246 for LMDZ6/trunk/libf/filtrez
- Timestamp:
- Oct 21, 2024, 2:58:45 PM (10 months ago)
- Location:
- LMDZ6/trunk/libf/filtrez
- Files:
-
- 5 moved
-
acc.f90 (moved) (moved from LMDZ6/trunk/libf/filtrez/acc.F) (1 diff)
-
eigen.f90 (moved) (moved from LMDZ6/trunk/libf/filtrez/eigen.F) (1 diff)
-
eigen_sort.f90 (moved) (moved from LMDZ6/trunk/libf/filtrez/eigen_sort.F) (1 diff)
-
filtreg.F90 (moved) (moved from LMDZ6/trunk/libf/filtrez/filtreg.F) (1 diff)
-
inifgn.F90 (moved) (moved from LMDZ6/trunk/libf/filtrez/inifgn.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/filtrez/acc.f90
r5245 r5246 2 2 ! $Header$ 3 3 ! 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 end 4 subroutine acc(vec,d,im) 5 implicit none 6 integer :: im 7 real :: vec(im,im),d(im) 8 integer :: i,j 9 real ::sum 10 real,external :: ssum 11 do j=1,im 12 do i=1,im 13 d(i)=vec(i,j)*vec(i,j) 14 enddo 15 sum=ssum(im,d,1) 16 sum=sqrt(sum) 17 do i=1,im 18 vec(i,j)=vec(i,j)/sum 19 enddo 20 enddo 21 return 22 end subroutine acc -
LMDZ6/trunk/libf/filtrez/eigen.f90
r5245 r5246 2 2 ! $Header$ 3 3 ! 4 5 6 7 8 9 10 11 c 12 DO 48i = 1,im13 14 48 CONTINUE15 DO 49i = 1,iim16 17 49 CONTINUE18 c 19 cPRINT 70,d4 SUBROUTINE eigen( e,d) 5 IMPLICIT NONE 6 INCLUDE "dimensions.h" 7 real :: e( iim,iim ), d( iim ) 8 real :: asm( iim ) 9 integer :: im,i,j 10 im=iim 11 ! 12 DO i = 1,im 13 asm( i ) = d( im-i+1 ) 14 END DO 15 DO i = 1,iim 16 d( i ) = asm( i ) 17 END DO 18 ! 19 ! PRINT 70,d 20 20 70 FORMAT(5x,'Valeurs propres',/,8(1x,8f10.4,/),/) 21 22 c 23 DO 51i = 1,im24 DO 52j = 1,im25 26 52 CONTINUE27 DO 50j = 1,im28 29 50 CONTINUE30 51 CONTINUE21 print * 22 ! 23 DO i = 1,im 24 DO j = 1,im 25 asm( j ) = e( i , im-j+1 ) 26 END DO 27 DO j = 1,im 28 e( i,j ) = asm( j ) 29 END DO 30 END DO 31 31 32 33 END 32 RETURN 33 END SUBROUTINE eigen -
LMDZ6/trunk/libf/filtrez/eigen_sort.f90
r5245 r5246 2 2 ! $Header$ 3 3 ! 4 5 INTEGERn,np6 REALd(np),v(np,np)7 INTEGERi,j,k8 REALp4 SUBROUTINE eigen_sort(d,v,n,np) 5 INTEGER :: n,np 6 REAL :: d(np),v(np,np) 7 INTEGER :: i,j,k 8 REAL :: p 9 9 10 DO i=1,n-1 11 k=i 12 p=d(i) 13 DO j=i+1,n 14 IF(d(j).ge.p) THEN 15 k=j 16 p=d(j) 17 ENDIF 18 ENDDO 19 20 IF(k.ne.i) THEN 21 d(k)=d(i) 22 d(i)=p 23 DO j=1,n 24 p=v(j,i) 25 v(j,i)=v(j,k) 26 v(j,k)=p 27 ENDDO 28 ENDIF 29 ENDDO 10 DO i=1,n-1 11 k=i 12 p=d(i) 13 DO j=i+1,n 14 IF(d(j).ge.p) THEN 15 k=j 16 p=d(j) 17 ENDIF 18 ENDDO 30 19 31 RETURN 32 END 20 IF(k.ne.i) THEN 21 d(k)=d(i) 22 d(i)=p 23 DO j=1,n 24 p=v(j,i) 25 v(j,i)=v(j,k) 26 v(j,k)=p 27 ENDDO 28 ENDIF 29 ENDDO 30 31 RETURN 32 END SUBROUTINE eigen_sort -
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 -
LMDZ6/trunk/libf/filtrez/inifgn.F90
r5245 r5246 2 2 ! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifgn.F,v 1.1.1.1 2004-05-19 12:53:09 lmdzadmin Exp $ 3 3 ! 4 5 c 6 c ... H.Upadyaya , O.Sharma ... 7 c 8 9 c 10 11 12 4 SUBROUTINE inifgn(dv) 5 ! 6 ! ... H.Upadyaya , O.Sharma ... 7 ! 8 IMPLICIT NONE 9 ! 10 include "dimensions.h" 11 include "paramet.h" 12 include "comgeom.h" 13 13 14 c 15 REALvec(iim,iim),vec1(iim,iim)16 REALdlonu(iim),dlonv(iim)17 REALdu(iim),dv(iim),d(iim)18 REALpi19 INTEGERi,j,k,imm1,nrot20 C 21 22 c 23 24 REALSSUM25 c 14 ! 15 REAL :: vec(iim,iim),vec1(iim,iim) 16 REAL :: dlonu(iim),dlonv(iim) 17 REAL :: du(iim),dv(iim),d(iim) 18 REAL :: pi 19 INTEGER :: i,j,k,imm1,nrot 20 ! 21 include "coefils.h" 22 ! 23 EXTERNAL SSUM, acc,eigen,jacobi 24 REAL :: SSUM 25 ! 26 26 27 28 29 C 30 DO 5i=1,iim31 32 33 5 CONTINUE27 imm1 = iim -1 28 pi = 2.* ASIN(1.) 29 ! 30 DO i=1,iim 31 dlonu(i)= xprimu( i ) 32 dlonv(i)= xprimv( i ) 33 END DO 34 34 35 DO 12 i=1,iim 36 sddv(i) = SQRT(dlonv(i)) 37 sddu(i) = SQRT(dlonu(i)) 38 unsddu(i) = 1./sddu(i) 39 unsddv(i) = 1./sddv(i) 40 12 CONTINUE 41 C 42 DO 17 j=1,iim 43 DO 17 i=1,iim 44 vec(i,j) = 0. 45 vec1(i,j) = 0. 46 eignfnv(i,j) = 0. 47 eignfnu(i,j) = 0. 48 17 CONTINUE 49 c 50 c 51 eignfnv(1,1) = -1. 52 eignfnv(iim,1) = 1. 53 DO 20 i=1,imm1 54 eignfnv(i+1,i+1)= -1. 55 eignfnv(i,i+1) = 1. 56 20 CONTINUE 57 DO 25 j=1,iim 58 DO 25 i=1,iim 59 eignfnv(i,j) = eignfnv(i,j)/(sddu(i)*sddv(j)) 60 25 CONTINUE 61 DO 30 j=1,iim 62 DO 30 i=1,iim 63 eignfnu(i,j) = -eignfnv(j,i) 64 30 CONTINUE 65 c 35 DO i=1,iim 36 sddv(i) = SQRT(dlonv(i)) 37 sddu(i) = SQRT(dlonu(i)) 38 unsddu(i) = 1./sddu(i) 39 unsddv(i) = 1./sddv(i) 40 END DO 41 ! 42 DO j=1,iim 43 DO i=1,iim 44 vec(i,j) = 0. 45 vec1(i,j) = 0. 46 eignfnv(i,j) = 0. 47 eignfnu(i,j) = 0. 48 END DO 49 END DO 50 ! 51 ! 52 eignfnv(1,1) = -1. 53 eignfnv(iim,1) = 1. 54 DO i=1,imm1 55 eignfnv(i+1,i+1)= -1. 56 eignfnv(i,i+1) = 1. 57 END DO 58 DO j=1,iim 59 DO i=1,iim 60 eignfnv(i,j) = eignfnv(i,j)/(sddu(i)*sddv(j)) 61 END DO 62 END DO 63 DO j=1,iim 64 DO i=1,iim 65 eignfnu(i,j) = -eignfnv(j,i) 66 END DO 67 END DO 68 ! 66 69 #ifdef CRAY 67 68 70 CALL MXM(eignfnu,iim,eignfnv,iim,vec ,iim) 71 CALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim) 69 72 #else 70 71 72 73 74 75 76 77 78 79 73 DO j = 1, iim 74 DO i = 1, iim 75 vec (i,j) = 0.0 76 vec1(i,j) = 0.0 77 DO k = 1, iim 78 vec (i,j) = vec(i,j) + eignfnu(i,k) * eignfnv(k,j) 79 vec1(i,j) = vec1(i,j) + eignfnv(i,k) * eignfnu(k,j) 80 ENDDO 81 ENDDO 82 ENDDO 80 83 #endif 81 84 82 c 83 84 85 86 c 87 88 89 85 ! 86 CALL jacobi(vec,iim,iim,dv,eignfnv,nrot) 87 CALL acc(eignfnv,d,iim) 88 CALL eigen_sort(dv,eignfnv,iim,iim) 89 ! 90 CALL jacobi(vec1,iim,iim,du,eignfnu,nrot) 91 CALL acc(eignfnu,d,iim) 92 CALL eigen_sort(du,eignfnu,iim,iim) 90 93 91 cc ancienne version avec appels IMSL92 c 93 cCALL MXM(eignfnu,iim,eignfnv,iim,vec,iim)94 cCALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim)95 cCALL EVCSF(iim,vec,iim,dv,eignfnv,iim)96 cCALL acc(eignfnv,d,iim)97 cCALL eigen(eignfnv,dv)98 c 99 cCALL EVCSF(iim,vec1,iim,du,eignfnu,iim)100 cCALL acc(eignfnu,d,iim)101 cCALL eigen(eignfnu,du)94 !c ancienne version avec appels IMSL 95 ! 96 ! CALL MXM(eignfnu,iim,eignfnv,iim,vec,iim) 97 ! CALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim) 98 ! CALL EVCSF(iim,vec,iim,dv,eignfnv,iim) 99 ! CALL acc(eignfnv,d,iim) 100 ! CALL eigen(eignfnv,dv) 101 ! 102 ! CALL EVCSF(iim,vec1,iim,du,eignfnu,iim) 103 ! CALL acc(eignfnu,d,iim) 104 ! CALL eigen(eignfnu,du) 102 105 103 104 END 106 RETURN 107 END SUBROUTINE inifgn 105 108
Note: See TracChangeset
for help on using the changeset viewer.
