Changeset 232 for LMDZ.3.3/branches/rel-LF/libf/dyn3d/fyhyp.F
- Timestamp:
- Jun 20, 2001, 3:53:15 PM (23 years ago)
- File:
-
- 1 edited
-
LMDZ.3.3/branches/rel-LF/libf/dyn3d/fyhyp.F (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/branches/rel-LF/libf/dyn3d/fyhyp.F
r2 r232 1 SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau,deltay , 2 , rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ) 3 1 c 2 c $Header$ 3 c 4 SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau , 5 , rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 , 6 , champmin,champmax ) 7 8 cc ... Version du 01/04/2001 .... 4 9 5 10 IMPLICIT NONE … … 13 18 c 14 19 c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc) 15 c dzoom etant la distance totale de la zone du zoom 16 c tau la transition , normalement = 1 . 17 18 c N.B : on doit avoir : grossism * dzoom < 1 . 19 c ************** 20 c 21 c Pour Indoex , on a pris : 22 c ******* 23 c grossism = 2.5 , dzoom = 7/24 en x et y , pour iim = 128 et jjm=64 24 c yzoomdeg = 0. , tau = 1. et delaty = 10. 20 c dzoom etant la distance totale de la zone du zoom ( en radians ) 21 c tau la raideur de la transition de l'interieur a l'exterieur du zoom 22 c 23 c 24 c N.B : Il vaut mieux avoir : grossism * dzoom < pi/2 (radians) ,en lati. 25 c ******************************************************************** 25 26 c 26 27 c … … 29 30 30 31 INTEGER nmax , nmax2 31 PARAMETER ( nmax = 50000, nmax2 = 2*nmax )32 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 32 33 c 33 34 c 34 35 c ....... arguments d'entree ....... 35 36 c 36 REAL yzoomdeg, grossism,dzoom,tau , deltay 37 REAL yzoomdeg, grossism,dzoom,tau 38 c ( rentres par run.def ) 37 39 38 40 c ....... arguments de sortie ....... … … 42 44 43 45 c 44 c ..... Champs locaux .....46 c ..... champs locaux ..... 45 47 c 46 48 47 REAl ylat(jjp1), yprim(jjp1) 48 REAL yuv 49 REAL ytild(0:nmax2) 50 REAL fhyp(0:nmax),ffdx(0:nmax),beta,Ytprim(0:nmax2) 51 SAVE Ytprim, ytild,Yf 52 REAL Yf(0:nmax2),yypr(0:nmax2) 53 REAL yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) 54 REAL pi,depi,pis2,epsilon,yzoom 55 REAL yo1,yi,ylon2,fxm,ymoy,yint,Yprimin 56 REAL ypn,deply,y00 49 REAL*8 ylat(jjp1), yprim(jjp1) 50 REAL*8 yuv 51 REAL*8 yt(0:nmax2) 52 REAL*8 fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2) 53 SAVE Ytprim, yt,Yf 54 REAL*8 Yf(0:nmax2),yypr(0:nmax2) 55 REAL*8 yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) 56 REAL*8 pi,depi,pis2,epsilon,y0,pisjm 57 REAL*8 yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax 58 REAL*8 yfi,Yf1,ffdy 59 REAL*8 ypn,deply,y00 57 60 SAVE y00, deply 58 61 … … 60 63 INTEGER jpn,jjpn 61 64 SAVE jpn 62 65 REAL*8 a0,a1,a2,a3,yi2,heavyy0,heavyy0m 66 REAL*8 fa(0:nmax2),fb(0:nmax2) 67 REAL y0min,y0max 68 69 REAL*8 heavyside 70 EXTERNAL heavyside 63 71 64 72 pi = 2. * ASIN(1.) 65 73 depi = 2. * pi 66 74 pis2 = pi/2. 67 epsilon = 1.e-6 68 yzoom = yzoomdeg * pi/180. 69 70 75 pisjm = pi/ FLOAT(jjm) 76 epsilon = 1.e-3 77 y0 = yzoomdeg * pi/180. 78 79 IF( dzoom.LT.1.) THEN 80 dzoom = dzoom * pi 81 ELSEIF( dzoom.LT. 12. ) THEN 82 WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug 83 ,menter et relancer ! ' 84 STOP 1 85 ELSE 86 dzoom = dzoom * pi/180. 87 ENDIF 88 89 WRITE(6,18) 90 WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)' 91 WRITE(6,24) y0,grossism,tau,dzoom 71 92 72 93 DO i = 0, nmax2 73 ytild(i) = FLOAT(i) /nmax2 74 IF( ytild(i).EQ.0.5 ) ytild(i) = ytild(i) + 1.e-6 75 ENDDO 76 77 DO i = 1, nmax 78 fhyp(i) = TANH ( ( ytild(i) - 0.5*(1.- dzoom) ) / 79 , ( tau * ytild(i) * ( 0.5 -ytild(i))) ) 80 ENDDO 81 82 fhyp( 0 ) = - 1. 83 fhyp( nmax ) = 1. 94 yt(i) = - pis2 + FLOAT(i)* pi /nmax2 95 ENDDO 96 97 heavyy0m = heavyside( -y0 ) 98 heavyy0 = heavyside( y0 ) 99 y0min = 2.*y0*heavyy0m - pis2 100 y0max = 2.*y0*heavyy0 + pis2 101 102 DO i = 0, nmax2 103 IF( yt(i).LT.y0 ) THEN 104 fa (i) = tau* (yt(i)-y0+dzoom/2. ) 105 fb(i) = (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) ) 106 ELSEIF ( yt(i).GT.y0 ) THEN 107 fa(i) = tau *(y0-yt(i)+dzoom/2. ) 108 fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) 109 ENDIF 110 111 IF( 200.* fb(i) .LT. - fa(i) ) THEN 112 fhyp ( i) = - 1. 113 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 114 fhyp ( i) = 1. 115 ELSE 116 fhyp(i) = TANH ( fa(i)/fb(i) ) 117 ENDIF 118 119 IF( yt(i).EQ.y0 ) fhyp(i) = 1. 120 IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1. 121 122 ENDDO 84 123 85 124 cc .... Calcul de beta .... 86 125 c 87 ffdx( 0 ) = 0. 88 89 DO i = 1, nmax 90 ymoy = 0.5 * ( ytild(i-1) + ytild( i ) ) 91 fxm = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) ) / 92 , ( tau * ymoy * ( 0.5 -ymoy)) ) 93 ffdx(i) = ffdx(i-1) + fxm * ( ytild(i) - ytild(i-1) ) 94 ENDDO 95 96 beta = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 ) 126 ffdy = 0. 127 128 DO i = 1, nmax2 129 ymoy = 0.5 * ( yt(i-1) + yt( i ) ) 130 IF( ymoy.LT.y0 ) THEN 131 fa(i)= tau * ( ymoy-y0+dzoom/2.) 132 fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy ) 133 ELSEIF ( ymoy.GT.y0 ) THEN 134 fa(i)= tau * ( y0-ymoy+dzoom/2. ) 135 fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 ) 136 ENDIF 137 138 IF( 200.* fb(i) .LT. - fa(i) ) THEN 139 fxm ( i) = - 1. 140 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 141 fxm ( i) = 1. 142 ELSE 143 fxm(i) = TANH ( fa(i)/fb(i) ) 144 ENDIF 145 IF( ymoy.EQ.y0 ) fxm(i) = 1. 146 IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1. 147 ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) ) 148 149 ENDDO 150 151 beta = ( grossism * ffdy - pi ) / ( ffdy - pi ) 152 153 IF( 2.*beta - grossism.LE. 0.) THEN 154 155 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 156 ,tine fyhyp est mauvaise ! ' 157 WRITE(6,*)'Modifier les valeurs de grossismy ,tauy ou dzoomy', 158 , ' et relancer ! *** ' 159 CALL ABORT 160 161 ENDIF 97 162 c 98 163 c ..... calcul de Ytprim ..... 99 164 c 100 165 101 DO i = 0, nmax 166 DO i = 0, nmax2 102 167 Ytprim(i) = beta + ( grossism - beta ) * fhyp(i) 103 168 ENDDO 104 c105 DO i = 0, nmax106 Ytprim( nmax2 - i ) = Ytprim( i )107 ENDDO108 c109 169 110 170 c ..... Calcul de Yf ........ 111 171 112 Yf(0) = 0. 113 DO i = 1, nmax 114 ymoy = 0.5 * ( ytild(i-1) + ytild( i ) ) 115 fxm = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) ) / 116 , ( tau * ymoy * ( 0.5 -ymoy)) ) 117 yypr(i) = beta + ( grossism - beta ) * fxm 118 ENDDO 119 120 DO i = 1,nmax 121 yypr(nmax2-i+1) = yypr(i) 122 ENDDO 123 124 DO i=1,nmax2 125 Yf(i) = Yf(i-1) + yypr(i) * ( ytild(i) - ytild(i-1) ) 126 ENDDO 127 c 128 c 172 Yf(0) = - pis2 173 DO i = 1, nmax2 174 yypr(i) = beta + ( grossism - beta ) * fxm(i) 175 ENDDO 176 177 DO i=1,nmax2 178 Yf(i) = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) ) 179 ENDDO 129 180 130 181 c **************************************************************** … … 133 184 c ..... yuv = 0.5 si calcul des latitudes aux pts. V ..... 134 185 c 186 WRITE(6,18) 135 187 c 136 188 DO 5000 ik = 1,4 … … 149 201 jlat = jjm 150 202 ENDIF 151 152 c 203 c 204 yo1 = 0. 153 205 DO 1500 j = 1,jlat 154 155 ylon2 = ( FLOAT(j) + yuv -1.) / FLOAT(jjm)156 157 206 yo1 = 0. 158 yi = ylon2 159 160 c 161 DO 500 iter = 1,300 162 207 ylon2 = - pis2 + pisjm * ( FLOAT(j) + yuv -1.) 208 yfi = ylon2 209 c 163 210 DO 250 it = nmax2,0,-1 164 IF( y i.GE.ytild(it)) GO TO 350211 IF( yfi.GE.Yf(it)) GO TO 350 165 212 250 CONTINUE 166 167 213 it = 0 168 yi = ytild(it)169 170 214 350 CONTINUE 171 215 216 yi = yt(it) 172 217 IF(it.EQ.nmax2) THEN 173 218 it = nmax2 -1 174 Yf(it+1) = 1.219 Yf(it+1) = pis2 175 220 ENDIF 176 221 c ................................................................. … … 179 224 c ................................................................. 180 225 181 yint = ( Yf(it+1)-Yf(it) ) / ( ytild(it+1)-ytild(it) ) * 182 + ( yi-ytild(it) ) + Yf(it) 183 Yprimin = ( Ytprim(it+1)-Ytprim(it) )/ ( ytild(it+1)-ytild(it) ) * 184 + ( yi-ytild(it) ) + Ytprim(it) 185 yi = yi - (yint-ylon2)/Yprimin 186 187 IF( ABS(yi-yo1).LE.epsilon) GO TO 550 188 yo1 = yi 189 c 226 CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1), 227 , yt(it),yt(it+1) , a0,a1,a2,a3 ) 228 229 Yf1 = Yf(it) 230 Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi 231 232 DO 500 iter = 1,300 233 yi = yi - ( Yf1 - yfi )/ Yprimin 234 235 IF( ABS(yi-yo1).LE.epsilon) GO TO 550 236 yo1 = yi 237 yi2 = yi * yi 238 Yf1 = a0 + a1 * yi + a2 * yi2 + a3 * yi2 * yi 239 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi2 190 240 500 CONTINUE 191 PRINT *,' *** PAS DE SOLUTION ****',j,ylon2,iter192 STOP 4241 WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter 242 STOP 2 193 243 550 CONTINUE 194 195 yprim(j) = pi /( FLOAT(jjm) * Yprimin) 196 yvrai(j) = pi * (yi - 0.5) + yzoom 244 c 245 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi* yi 246 yprim(j) = pi / ( jjm * Yprimin ) 247 yvrai(j) = yi 197 248 198 249 1500 CONTINUE 199 200 cc print *,' LAT avant reorgan '201 cc print 68,(yyvrai(j),j=1,jlat)202 250 203 251 DO j = 1, jlat -1 204 252 IF( yvrai(j+1). LT. yvrai(j) ) THEN 205 PRINT *,' PBS. avec rlat(',j+1,' plus petit que rlat(',j, 206 , ')' 207 STOP 208 ENDIF 209 ENDDO 210 211 PRINT 18 212 PRINT *,'Reorganisation des latitudes pour avoir entre - pi/2 ', 213 , ' et pi/2 ' 214 c 215 253 WRITE(6,*) ' PBS. avec rlat(',j+1,') plus petit que rlat(',j, 254 , ')' 255 STOP 3 256 ENDIF 257 ENDDO 258 259 WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2' 260 , ,' et pi/2 ' 261 c 216 262 IF( ik.EQ.1 ) THEN 217 ypn = pis2 - deltay * pi/180.263 ypn = pis2 218 264 DO j = jlat,1,-1 219 265 IF( yvrai(j).LE. ypn ) GO TO 1502 … … 239 285 ENDDO 240 286 241 242 243 287 c *********** Fin de la reorganisation ************* 244 288 c 245 289 1600 CONTINUE 246 247 248 290 249 291 DO j = 1, jlat … … 256 298 ENDDO 257 299 258 259 300 IF( ik.EQ.1 ) THEN 260 PRINT 18 261 PRINT *, ' YLAT en U apres ( en deg. ) ' 262 PRINT 68,(yvrai(j),j=1,jlat) 263 PRINT *,' YPRIM ' 264 PRINT 68,( yprim(j),j=1,jlat) 301 c WRITE(6,18) 302 c WRITE(6,*) ' YLAT en U apres ( en deg. ) ' 303 c WRITE(6,68) (yvrai(j),j=1,jlat) 304 cc WRITE(6,*) ' YPRIM ' 305 cc WRITE(6,445) ( yprim(j),j=1,jlat) 306 265 307 DO j = 1, jlat 266 308 rrlatu(j) = ylat( j ) 267 309 yyprimu(j) = yprim( j ) 268 310 ENDDO 269 c 311 270 312 ELSE IF ( ik.EQ. 2 ) THEN 271 PRINT 18 272 PRINT *, ' YLAT en V apres ( en deg. ) ' 273 PRINT 68,(yvrai(j),j=1,jlat) 274 PRINT *,' YPRIM ' 275 PRINT 68,( yprim(j),j=1,jlat) 313 c WRITE(6,18) 314 c WRITE(6,*) ' YLAT en V apres ( en deg. ) ' 315 c WRITE(6,68) (yvrai(j),j=1,jlat) 316 cc WRITE(6,*)' YPRIM ' 317 cc WRITE(6,445) ( yprim(j),j=1,jlat) 318 276 319 DO j = 1, jlat 277 320 rrlatv(j) = ylat( j ) 278 321 yyprimv(j) = yprim( j ) 279 322 ENDDO 280 c 323 281 324 ELSE IF ( ik.EQ. 3 ) THEN 282 PRINT 18 283 PRINT *, ' YLAT en U + 0.75 apres ( en deg. ) ' 284 PRINT 68,(yvrai(j),j=1,jlat) 285 PRINT *,' YPRIM ' 286 PRINT 68,( yprim(j),j=1,jlat) 325 c WRITE(6,18) 326 c WRITE(6,*) ' YLAT en U + 0.75 apres ( en deg. ) ' 327 c WRITE(6,68) (yvrai(j),j=1,jlat) 328 cc WRITE(6,*) ' YPRIM ' 329 cc WRITE(6,445) ( yprim(j),j=1,jlat) 330 287 331 DO j = 1, jlat 288 332 rlatu2(j) = ylat( j ) … … 291 335 292 336 ELSE IF ( ik.EQ. 4 ) THEN 293 PRINT 18 294 PRINT *, ' YLAT en U + 0.25 apres ( en deg. ) ' 295 PRINT 68,(yvrai(j),j=1,jlat) 296 PRINT *,' YPRIM ' 297 PRINT 68,( yprim(j),j=1,jlat) 337 c WRITE(6,18) 338 c WRITE(6,*) ' YLAT en U + 0.25 apres ( en deg. ) ' 339 c WRITE(6,68)(yvrai(j),j=1,jlat) 340 cc WRITE(6,*) ' YPRIM ' 341 cc WRITE(6,68) ( yprim(j),j=1,jlat) 342 298 343 DO j = 1, jlat 299 344 rlatu1(j) = ylat( j ) 300 345 yprimu1(j) = yprim( j ) 301 346 ENDDO 347 302 348 ENDIF 303 349 304 350 5000 CONTINUE 305 351 c 352 WRITE(6,18) 353 c 306 354 c ..... fin de la boucle do 5000 ..... 307 355 356 DO j = 1, jjm 357 ylat(j) = rrlatu(j) - rrlatu(j+1) 358 ENDDO 359 champmin = 1.e12 360 champmax = -1.e12 361 DO j = 1, jjm 362 champmin = MIN( champmin, ylat(j) ) 363 champmax = MAX( champmax, ylat(j) ) 364 ENDDO 365 champmin = champmin * 180./pi 366 champmax = champmax * 180./pi 367 368 24 FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3) 308 369 18 FORMAT(/) 309 370 68 FORMAT(1x,7f9.2)
Note: See TracChangeset
for help on using the changeset viewer.
