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