Changeset 2218 for LMDZ5/trunk/libf/dyn3d_common/fyhyp_m.F90
- Timestamp:
- Mar 2, 2015, 5:11:07 PM (10 years ago)
- File:
-
- 1 moved
-
LMDZ5/trunk/libf/dyn3d_common/fyhyp_m.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/fyhyp.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/dyn3d_common/fyhyp_m.F90
r2217 r2218 1 ! 2 ! $Id$ 3 ! 4 c 5 c 6 SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau , 7 , rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 , 8 , champmin,champmax ) 9 10 cc ... Version du 01/04/2001 .... 11 12 IMPLICIT NONE 13 c 14 c ... Auteur : P. Le Van ... 15 c 16 c ....... d'apres formulations de R. Sadourny ....... 17 c 18 c Calcule les latitudes et derivees dans la grille du GCM pour une 19 c fonction f(y) a tangente hyperbolique . 20 c 21 c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc) 22 c dzoom etant la distance totale de la zone du zoom ( en radians ) 23 c tau la raideur de la transition de l'interieur a l'exterieur du zoom 24 c 25 c 26 c N.B : Il vaut mieux avoir : grossism * dzoom < pi/2 (radians) ,en lati. 27 c ******************************************************************** 28 c 29 c 30 #include "dimensions.h" 31 #include "paramet.h" 32 33 INTEGER nmax , nmax2 34 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 35 c 36 c 37 c ....... arguments d'entree ....... 38 c 39 REAL yzoomdeg, grossism,dzooma,tau 40 c ( rentres par run.def ) 41 42 c ....... arguments de sortie ....... 43 c 44 REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm), 45 , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm) 46 47 c 48 c ..... champs locaux ..... 49 c 50 51 REAL dzoom 52 REAL(KIND=8) ylat(jjp1), yprim(jjp1) 53 REAL(KIND=8) yuv 54 REAL(KIND=8) yt(0:nmax2) 55 REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2) 56 SAVE Ytprim, yt,Yf 57 REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2) 58 REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) 59 REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm 60 REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax 61 REAL(KIND=8) yfi,Yf1,ffdy 62 REAL(KIND=8) ypn,deply,y00 63 SAVE y00, deply 64 65 INTEGER i,j,it,ik,iter,jlat 66 INTEGER jpn,jjpn 67 SAVE jpn 68 REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m 69 REAL(KIND=8) fa(0:nmax2),fb(0:nmax2) 70 REAL y0min,y0max 71 72 REAL(KIND=8) heavyside 73 74 pi = 2. * ASIN(1.) 75 depi = 2. * pi 76 pis2 = pi/2. 77 pisjm = pi/ REAL(jjm) 78 epsilon = 1.e-3 79 y0 = yzoomdeg * pi/180. 80 81 IF( dzooma.LT.1.) THEN 82 dzoom = dzooma * pi 83 ELSEIF( dzooma.LT. 12. ) THEN 84 WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug 85 ,menter et relancer ! ' 86 STOP 1 1 module fyhyp_m 2 3 IMPLICIT NONE 4 5 contains 6 7 SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) 8 9 ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32 10 11 ! Author: P. Le Van, from analysis by R. Sadourny 12 13 ! Calcule les latitudes et dérivées dans la grille du GCM pour une 14 ! fonction f(y) à dérivée tangente hyperbolique. 15 16 ! Il vaut mieux avoir : grossismy * dzoom < pi / 2 17 18 use coefpoly_m, only: coefpoly 19 20 include "dimensions.h" 21 ! for jjm 22 23 include "serre.h" 24 ! for clat, grossismy, dzoomy, tauy 25 26 REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1) 27 REAL, intent(out):: rlatv(jjm) 28 real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm) 29 30 ! Local: 31 32 DOUBLE PRECISION champmin, champmax 33 INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax 34 REAL dzoom ! distance totale de la zone du zoom (en radians) 35 DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1) 36 DOUBLE PRECISION yuv 37 DOUBLE PRECISION, save:: yt(0:nmax2) 38 DOUBLE PRECISION fhyp(0:nmax2), beta 39 DOUBLE PRECISION, save:: ytprim(0:nmax2) 40 DOUBLE PRECISION fxm(0:nmax2) 41 DOUBLE PRECISION, save:: yf(0:nmax2) 42 DOUBLE PRECISION yypr(0:nmax2) 43 DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1) 44 DOUBLE PRECISION pi, pis2, epsilon, y0, pisjm 45 DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin 46 DOUBLE PRECISION yfi, yf1, ffdy 47 DOUBLE PRECISION ypn, deply, y00 48 SAVE y00, deply 49 50 INTEGER i, j, it, ik, iter, jlat 51 INTEGER jpn, jjpn 52 SAVE jpn 53 DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m 54 DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2) 55 REAL y0min, y0max 56 57 DOUBLE PRECISION heavyside 58 59 !------------------------------------------------------------------- 60 61 print *, "Call sequence information: fyhyp" 62 63 pi = 2.*asin(1.) 64 pis2 = pi/2. 65 pisjm = pi/real(jjm) 66 epsilon = 1e-3 67 y0 = clat*pi/180. 68 dzoom = dzoomy*pi 69 print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):' 70 print *, y0, grossismy, tauy, dzoom 71 72 DO i = 0, nmax2 73 yt(i) = -pis2 + real(i)*pi/nmax2 74 END DO 75 76 heavyy0m = heavyside(-y0) 77 heavyy0 = heavyside(y0) 78 y0min = 2.*y0*heavyy0m - pis2 79 y0max = 2.*y0*heavyy0 + pis2 80 81 fa = 999.999 82 fb = 999.999 83 84 DO i = 0, nmax2 85 IF (yt(i)<y0) THEN 86 fa(i) = tauy*(yt(i)-y0 + dzoom/2.) 87 fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i)) 88 ELSE IF (yt(i)>y0) THEN 89 fa(i) = tauy*(y0-yt(i) + dzoom/2.) 90 fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0) 91 END IF 92 93 IF (200.*fb(i)<-fa(i)) THEN 94 fhyp(i) = -1. 95 ELSE IF (200.*fb(i)<fa(i)) THEN 96 fhyp(i) = 1. 87 97 ELSE 88 dzoom = dzooma * pi/180. 89 ENDIF 90 91 WRITE(6,18) 92 WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)' 93 WRITE(6,24) y0,grossism,tau,dzoom 94 95 DO i = 0, nmax2 96 yt(i) = - pis2 + REAL(i)* pi /nmax2 97 ENDDO 98 99 heavyy0m = heavyside( -y0 ) 100 heavyy0 = heavyside( y0 ) 101 y0min = 2.*y0*heavyy0m - pis2 102 y0max = 2.*y0*heavyy0 + pis2 103 104 fa = 999.999 105 fb = 999.999 106 107 DO i = 0, nmax2 108 IF( yt(i).LT.y0 ) THEN 109 fa (i) = tau* (yt(i)-y0+dzoom/2. ) 110 fb(i) = (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) ) 111 ELSEIF ( yt(i).GT.y0 ) THEN 112 fa(i) = tau *(y0-yt(i)+dzoom/2. ) 113 fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) 114 ENDIF 115 116 IF( 200.* fb(i) .LT. - fa(i) ) THEN 117 fhyp ( i) = - 1. 118 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 119 fhyp ( i) = 1. 120 ELSE 121 fhyp(i) = TANH ( fa(i)/fb(i) ) 122 ENDIF 123 124 IF( yt(i).EQ.y0 ) fhyp(i) = 1. 125 IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1. 126 127 ENDDO 128 129 cc .... Calcul de beta .... 130 c 131 ffdy = 0. 132 133 DO i = 1, nmax2 134 ymoy = 0.5 * ( yt(i-1) + yt( i ) ) 135 IF( ymoy.LT.y0 ) THEN 136 fa(i)= tau * ( ymoy-y0+dzoom/2.) 137 fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy ) 138 ELSEIF ( ymoy.GT.y0 ) THEN 139 fa(i)= tau * ( y0-ymoy+dzoom/2. ) 140 fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 ) 141 ENDIF 142 143 IF( 200.* fb(i) .LT. - fa(i) ) THEN 144 fxm ( i) = - 1. 145 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 146 fxm ( i) = 1. 147 ELSE 148 fxm(i) = TANH ( fa(i)/fb(i) ) 149 ENDIF 150 IF( ymoy.EQ.y0 ) fxm(i) = 1. 151 IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1. 152 ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) ) 153 154 ENDDO 155 156 beta = ( grossism * ffdy - pi ) / ( ffdy - pi ) 157 158 IF( 2.*beta - grossism.LE. 0.) THEN 159 160 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 161 ,tine fyhyp est mauvaise ! ' 162 WRITE(6,*)'Modifier les valeurs de grossismy ,tauy ou dzoomy', 163 , ' et relancer ! *** ' 164 CALL ABORT_GCM("FYHYP", "", 1) 165 166 ENDIF 167 c 168 c ..... calcul de Ytprim ..... 169 c 170 171 DO i = 0, nmax2 172 Ytprim(i) = beta + ( grossism - beta ) * fhyp(i) 173 ENDDO 174 175 c ..... Calcul de Yf ........ 176 177 Yf(0) = - pis2 178 DO i = 1, nmax2 179 yypr(i) = beta + ( grossism - beta ) * fxm(i) 180 ENDDO 181 182 DO i=1,nmax2 183 Yf(i) = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) ) 184 ENDDO 185 186 c **************************************************************** 187 c 188 c ..... yuv = 0. si calcul des latitudes aux pts. U ..... 189 c ..... yuv = 0.5 si calcul des latitudes aux pts. V ..... 190 c 191 WRITE(6,18) 192 c 193 DO 5000 ik = 1,4 194 195 IF( ik.EQ.1 ) THEN 196 yuv = 0. 197 jlat = jjm + 1 198 ELSE IF ( ik.EQ.2 ) THEN 199 yuv = 0.5 200 jlat = jjm 201 ELSE IF ( ik.EQ.3 ) THEN 202 yuv = 0.25 203 jlat = jjm 204 ELSE IF ( ik.EQ.4 ) THEN 205 yuv = 0.75 206 jlat = jjm 207 ENDIF 208 c 209 yo1 = 0. 210 DO 1500 j = 1,jlat 211 yo1 = 0. 212 ylon2 = - pis2 + pisjm * ( REAL(j) + yuv -1.) 213 yfi = ylon2 214 c 215 DO 250 it = nmax2,0,-1 216 IF( yfi.GE.Yf(it)) GO TO 350 217 250 CONTINUE 218 it = 0 219 350 CONTINUE 220 221 yi = yt(it) 222 IF(it.EQ.nmax2) THEN 223 it = nmax2 -1 224 Yf(it+1) = pis2 225 ENDIF 226 c ................................................................. 227 c .... Interpolation entre yi(it) et yi(it+1) pour avoir Y(yi) 228 c ..... et Y'(yi) ..... 229 c ................................................................. 230 231 CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1), 232 , yt(it),yt(it+1) , a0,a1,a2,a3 ) 233 234 Yf1 = Yf(it) 235 Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi 236 237 DO 500 iter = 1,300 238 yi = yi - ( Yf1 - yfi )/ Yprimin 239 240 IF( ABS(yi-yo1).LE.epsilon) GO TO 550 241 yo1 = yi 242 yi2 = yi * yi 243 Yf1 = a0 + a1 * yi + a2 * yi2 + a3 * yi2 * yi 244 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi2 245 500 CONTINUE 246 WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter 247 STOP 2 248 550 CONTINUE 249 c 250 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi* yi 251 yprim(j) = pi / ( jjm * Yprimin ) 252 yvrai(j) = yi 253 254 1500 CONTINUE 255 256 DO j = 1, jlat -1 257 IF( yvrai(j+1). LT. yvrai(j) ) THEN 258 WRITE(6,*) ' PBS. avec rlat(',j+1,') plus petit que rlat(',j, 259 , ')' 260 STOP 3 261 ENDIF 262 ENDDO 263 264 WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2' 265 , ,' et pi/2 ' 266 c 267 IF( ik.EQ.1 ) THEN 268 ypn = pis2 269 DO j = jlat,1,-1 270 IF( yvrai(j).LE. ypn ) GO TO 1502 271 ENDDO 272 1502 CONTINUE 273 274 jpn = j 275 y00 = yvrai(jpn) 276 deply = pis2 - y00 277 ENDIF 278 279 DO j = 1, jjm +1 - jpn 280 ylatt (j) = -pis2 - y00 + yvrai(jpn+j-1) 281 yprimm(j) = yprim(jpn+j-1) 282 ENDDO 283 284 jjpn = jpn 285 IF( jlat.EQ. jjm ) jjpn = jpn -1 286 287 DO j = 1,jjpn 288 ylatt (j + jjm+1 -jpn) = yvrai(j) + deply 289 yprimm(j + jjm+1 -jpn) = yprim(j) 290 ENDDO 291 292 c *********** Fin de la reorganisation ************* 293 c 294 1600 CONTINUE 295 98 fhyp(i) = tanh(fa(i)/fb(i)) 99 END IF 100 101 IF (yt(i)==y0) fhyp(i) = 1. 102 IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1. 103 END DO 104 105 ! Calcul de beta 106 107 ffdy = 0. 108 109 DO i = 1, nmax2 110 ymoy = 0.5*(yt(i-1) + yt(i)) 111 IF (ymoy<y0) THEN 112 fa(i) = tauy*(ymoy-y0 + dzoom/2.) 113 fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy) 114 ELSE IF (ymoy>y0) THEN 115 fa(i) = tauy*(y0-ymoy + dzoom/2.) 116 fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0) 117 END IF 118 119 IF (200.*fb(i)<-fa(i)) THEN 120 fxm(i) = -1. 121 ELSE IF (200.*fb(i)<fa(i)) THEN 122 fxm(i) = 1. 123 ELSE 124 fxm(i) = tanh(fa(i)/fb(i)) 125 END IF 126 IF (ymoy==y0) fxm(i) = 1. 127 IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1. 128 ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1)) 129 END DO 130 131 beta = (grossismy*ffdy-pi)/(ffdy-pi) 132 133 IF (2. * beta - grossismy <= 0.) THEN 134 print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' & 135 // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' & 136 // 'dzoomy et relancer.' 137 STOP 1 138 END IF 139 140 ! calcul de Ytprim 141 142 DO i = 0, nmax2 143 ytprim(i) = beta + (grossismy-beta)*fhyp(i) 144 END DO 145 146 ! Calcul de Yf 147 148 yf(0) = -pis2 149 DO i = 1, nmax2 150 yypr(i) = beta + (grossismy-beta)*fxm(i) 151 END DO 152 153 DO i = 1, nmax2 154 yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1)) 155 END DO 156 157 ! yuv = 0. si calcul des latitudes aux pts. U 158 ! yuv = 0.5 si calcul des latitudes aux pts. V 159 160 loop_ik: DO ik = 1, 4 161 IF (ik==1) THEN 162 yuv = 0. 163 jlat = jjm + 1 164 ELSE IF (ik==2) THEN 165 yuv = 0.5 166 jlat = jjm 167 ELSE IF (ik==3) THEN 168 yuv = 0.25 169 jlat = jjm 170 ELSE IF (ik==4) THEN 171 yuv = 0.75 172 jlat = jjm 173 END IF 174 175 yo1 = 0. 296 176 DO j = 1, jlat 297 ylat(j) = ylatt( jlat +1 -j ) 298 yprim(j) = yprimm( jlat +1 -j ) 299 ENDDO 300 301 DO j = 1, jlat 302 yvrai(j) = ylat(j)*180./pi 303 ENDDO 304 305 IF( ik.EQ.1 ) THEN 306 c WRITE(6,18) 307 c WRITE(6,*) ' YLAT en U apres ( en deg. ) ' 308 c WRITE(6,68) (yvrai(j),j=1,jlat) 309 cc WRITE(6,*) ' YPRIM ' 310 cc WRITE(6,445) ( yprim(j),j=1,jlat) 311 312 DO j = 1, jlat 313 rrlatu(j) = ylat( j ) 314 yyprimu(j) = yprim( j ) 315 ENDDO 316 317 ELSE IF ( ik.EQ. 2 ) THEN 318 c WRITE(6,18) 319 c WRITE(6,*) ' YLAT en V apres ( en deg. ) ' 320 c WRITE(6,68) (yvrai(j),j=1,jlat) 321 cc WRITE(6,*)' YPRIM ' 322 cc WRITE(6,445) ( yprim(j),j=1,jlat) 323 324 DO j = 1, jlat 325 rrlatv(j) = ylat( j ) 326 yyprimv(j) = yprim( j ) 327 ENDDO 328 329 ELSE IF ( ik.EQ. 3 ) THEN 330 c WRITE(6,18) 331 c WRITE(6,*) ' YLAT en U + 0.75 apres ( en deg. ) ' 332 c WRITE(6,68) (yvrai(j),j=1,jlat) 333 cc WRITE(6,*) ' YPRIM ' 334 cc WRITE(6,445) ( yprim(j),j=1,jlat) 335 336 DO j = 1, jlat 337 rlatu2(j) = ylat( j ) 338 yprimu2(j) = yprim( j ) 339 ENDDO 340 341 ELSE IF ( ik.EQ. 4 ) THEN 342 c WRITE(6,18) 343 c WRITE(6,*) ' YLAT en U + 0.25 apres ( en deg. ) ' 344 c WRITE(6,68)(yvrai(j),j=1,jlat) 345 cc WRITE(6,*) ' YPRIM ' 346 cc WRITE(6,68) ( yprim(j),j=1,jlat) 347 348 DO j = 1, jlat 349 rlatu1(j) = ylat( j ) 350 yprimu1(j) = yprim( j ) 351 ENDDO 352 353 ENDIF 354 355 5000 CONTINUE 356 c 357 WRITE(6,18) 358 c 359 c ..... fin de la boucle do 5000 ..... 360 361 DO j = 1, jjm 362 ylat(j) = rrlatu(j) - rrlatu(j+1) 363 ENDDO 364 champmin = 1.e12 365 champmax = -1.e12 366 DO j = 1, jjm 367 champmin = MIN( champmin, ylat(j) ) 368 champmax = MAX( champmax, ylat(j) ) 369 ENDDO 370 champmin = champmin * 180./pi 371 champmax = champmax * 180./pi 372 373 24 FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3) 374 18 FORMAT(/) 375 68 FORMAT(1x,7f9.2) 376 377 RETURN 378 END 177 yo1 = 0. 178 ylon2 = -pis2 + pisjm*(real(j) + yuv-1.) 179 yfi = ylon2 180 181 it = nmax2 182 DO while (it >= 1 .and. yfi < yf(it)) 183 it = it - 1 184 END DO 185 186 yi = yt(it) 187 IF (it==nmax2) THEN 188 it = nmax2 - 1 189 yf(it + 1) = pis2 190 END IF 191 192 ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi) 193 ! et Y'(yi) 194 195 CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), & 196 yt(it), yt(it + 1), a0, a1, a2, a3) 197 198 yf1 = yf(it) 199 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi 200 201 iter = 1 202 DO 203 yi = yi - (yf1-yfi)/yprimin 204 IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit 205 yo1 = yi 206 yi2 = yi*yi 207 yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi 208 yprimin = a1 + 2.*a2*yi + 3.*a3*yi2 209 END DO 210 if (abs(yi-yo1) > epsilon) then 211 print *, 'Pas de solution.', j, ylon2 212 STOP 1 213 end if 214 215 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi 216 yprim(j) = pi/(jjm*yprimin) 217 yvrai(j) = yi 218 END DO 219 220 DO j = 1, jlat - 1 221 IF (yvrai(j + 1)<yvrai(j)) THEN 222 print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', & 223 j, ')' 224 STOP 1 225 END IF 226 END DO 227 228 print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2' 229 230 IF (ik==1) THEN 231 ypn = pis2 232 DO j = jjm + 1, 1, -1 233 IF (yvrai(j)<=ypn) exit 234 END DO 235 236 jpn = j 237 y00 = yvrai(jpn) 238 deply = pis2 - y00 239 END IF 240 241 DO j = 1, jjm + 1 - jpn 242 ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1) 243 yprimm(j) = yprim(jpn + j-1) 244 END DO 245 246 jjpn = jpn 247 IF (jlat==jjm) jjpn = jpn - 1 248 249 DO j = 1, jjpn 250 ylatt(j + jjm + 1-jpn) = yvrai(j) + deply 251 yprimm(j + jjm + 1-jpn) = yprim(j) 252 END DO 253 254 ! Fin de la reorganisation 255 256 DO j = 1, jlat 257 ylat(j) = ylatt(jlat + 1-j) 258 yprim(j) = yprimm(jlat + 1-j) 259 END DO 260 261 DO j = 1, jlat 262 yvrai(j) = ylat(j)*180./pi 263 END DO 264 265 IF (ik==1) THEN 266 DO j = 1, jjm + 1 267 rlatu(j) = ylat(j) 268 yyprimu(j) = yprim(j) 269 END DO 270 ELSE IF (ik==2) THEN 271 DO j = 1, jjm 272 rlatv(j) = ylat(j) 273 END DO 274 ELSE IF (ik==3) THEN 275 DO j = 1, jjm 276 rlatu2(j) = ylat(j) 277 yprimu2(j) = yprim(j) 278 END DO 279 ELSE IF (ik==4) THEN 280 DO j = 1, jjm 281 rlatu1(j) = ylat(j) 282 yprimu1(j) = yprim(j) 283 END DO 284 END IF 285 END DO loop_ik 286 287 DO j = 1, jjm 288 ylat(j) = rlatu(j) - rlatu(j + 1) 289 END DO 290 champmin = 1e12 291 champmax = -1e12 292 DO j = 1, jjm 293 champmin = min(champmin, ylat(j)) 294 champmax = max(champmax, ylat(j)) 295 END DO 296 champmin = champmin*180./pi 297 champmax = champmax*180./pi 298 299 DO j = 1, jjm 300 IF (rlatu1(j) <= rlatu2(j)) THEN 301 print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j 302 STOP 13 303 ENDIF 304 305 IF (rlatu2(j) <= rlatu(j+1)) THEN 306 print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j 307 STOP 14 308 ENDIF 309 310 IF (rlatu(j) <= rlatu1(j)) THEN 311 print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j 312 STOP 15 313 ENDIF 314 315 IF (rlatv(j) <= rlatu2(j)) THEN 316 print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j 317 STOP 16 318 ENDIF 319 320 IF (rlatv(j) >= rlatu1(j)) THEN 321 print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j 322 STOP 17 323 ENDIF 324 325 IF (rlatv(j) >= rlatu(j)) THEN 326 print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j 327 STOP 18 328 ENDIF 329 ENDDO 330 331 print *, 'Latitudes' 332 print 3, champmin, champmax 333 334 3 Format(1x, ' Au centre du zoom, la longueur de la maille est', & 335 ' d environ ', f0.2, ' degres ', /, & 336 ' alors que la maille en dehors de la zone du zoom est ', & 337 "d'environ ", f0.2, ' degres ') 338 339 END SUBROUTINE fyhyp 340 341 end module fyhyp_m
Note: See TracChangeset
for help on using the changeset viewer.
