Changeset 1441 for trunk/LMDZ.COMMON/libf/dyn3d_common
- Timestamp:
- Jun 4, 2015, 10:21:20 AM (10 years ago)
- Location:
- trunk/LMDZ.COMMON/libf/dyn3d_common
- Files:
-
- 2 added
- 2 edited
- 2 moved
-
fxhyp_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/fxhyp.F) (1 diff)
-
fyhyp_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/fyhyp.F) (1 diff)
-
grilles_gcm_netcdf_sub.F90 (modified) (2 diffs)
-
inigeom.F (modified) (2 diffs)
-
invert_zoom_x_m.F90 (added)
-
principal_cshift_m.F90 (added)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/dyn3d_common/fxhyp_m.F90
r1440 r1441 1 ! 2 ! $Id: fxhyp.F 1674 2012-10-29 16:27:03Z emillour $ 3 ! 4 c 5 c 6 SUBROUTINE fxhyp ( xzoomdeg,grossism,dzooma,tau , 7 , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025, 8 , champmin,champmax ) 9 10 c Auteur : P. Le Van 11 12 IMPLICIT NONE 13 14 c Calcule les longitudes et derivees dans la grille du GCM pour une 15 c fonction f(x) a tangente hyperbolique . 16 c 17 c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.) 18 c dzoom etant la distance totale de la zone du zoom 19 c tau la raideur de la transition de l'interieur a l'exterieur du zoom 20 c 21 c On doit avoir grossism x dzoom < pi ( radians ) , en longitude. 22 c ******************************************************************** 23 24 25 INTEGER nmax, nmax2 26 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 27 c 28 LOGICAL scal180 29 PARAMETER ( scal180 = .TRUE. ) 30 31 c scal180 = .TRUE. si on veut avoir le premier point scalaire pour 32 c une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres. 33 c sinon scal180 = .FALSE. 34 35 #include "dimensions.h" 36 #include "paramet.h" 37 38 c ...... arguments d'entree ....... 39 c 40 REAL xzoomdeg,dzooma,tau,grossism 41 42 c ...... arguments de sortie ...... 43 44 REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1), 45 , rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1) 46 47 c .... variables locales .... 48 c 49 REAL dzoom 50 REAL(KIND=8) xlon(iip1),xprimm(iip1),xuv 51 REAL(KIND=8) xtild(0:nmax2) 52 REAL(KIND=8) fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2) 53 REAL(KIND=8) Xf(0:nmax2),xxpr(0:nmax2) 54 REAL(KIND=8) xvrai(iip1),xxprim(iip1) 55 REAL(KIND=8) pi,depi,epsilon,xzoom,fa,fb 56 REAL(KIND=8) Xf1, Xfi , a0,a1,a2,a3,xi2 57 INTEGER i,it,ik,iter,ii,idif,ii1,ii2 58 REAL(KIND=8) xi,xo1,xmoy,xlon2,fxm,Xprimin 59 REAL(KIND=8) champmin,champmax,decalx 60 INTEGER is2 61 SAVE is2 62 63 REAL(KIND=8) heavyside 64 65 pi = 2. * ASIN(1.) 66 depi = 2. * pi 67 epsilon = 1.e-3 68 xzoom = xzoomdeg * pi/180. 69 c 70 if (iim==1) then 71 72 rlonm025(1)=-pi/2. 73 rlonv(1)=0. 74 rlonu(1)=pi 75 rlonp025(1)=pi/2. 76 rlonm025(2)=rlonm025(1)+depi 77 rlonv(2)=rlonv(1)+depi 78 rlonu(2)=rlonu(1)+depi 79 rlonp025(2)=rlonp025(1)+depi 80 81 xprimm025(:)=1. 82 xprimv(:)=1. 83 xprimu(:)=1. 84 xprimp025(:)=1. 85 champmin=depi 86 champmax=depi 87 return 88 89 endif 90 91 decalx = .75 92 IF( grossism.EQ.1..AND.scal180 ) THEN 93 decalx = 1. 94 ENDIF 95 96 WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx 97 c 98 IF( dzooma.LT.1.) THEN 99 dzoom = dzooma * depi 100 ELSEIF( dzooma.LT. 25. ) THEN 101 WRITE(6,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug 102 ,menter et relancer ! ' 103 STOP 1 104 ELSE 105 dzoom = dzooma * pi/180. 106 ENDIF 107 108 WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)' 109 WRITE(6,24) xzoom,grossism,tau,dzoom 110 111 DO i = 0, nmax2 112 xtild(i) = - pi + REAL(i) * depi /nmax2 113 ENDDO 114 115 DO i = nmax, nmax2 116 117 fa = tau* ( dzoom/2. - xtild(i) ) 118 fb = xtild(i) * ( pi - xtild(i) ) 119 120 IF( 200.* fb .LT. - fa ) THEN 121 fhyp ( i) = - 1. 122 ELSEIF( 200. * fb .LT. fa ) THEN 123 fhyp ( i) = 1. 124 ELSE 125 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 126 IF( 200.*fb + fa.LT.1.e-10 ) THEN 127 fhyp ( i ) = - 1. 128 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 129 fhyp ( i ) = 1. 130 ENDIF 131 ELSE 132 fhyp ( i ) = TANH ( fa/fb ) 133 ENDIF 134 ENDIF 135 IF ( xtild(i).EQ. 0. ) fhyp(i) = 1. 136 IF ( xtild(i).EQ. pi ) fhyp(i) = -1. 137 138 ENDDO 139 140 cc .... Calcul de beta .... 141 142 ffdx = 0. 143 144 DO i = nmax +1,nmax2 145 146 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 147 fa = tau* ( dzoom/2. - xmoy ) 148 fb = xmoy * ( pi - xmoy ) 149 150 IF( 200.* fb .LT. - fa ) THEN 151 fxm = - 1. 152 ELSEIF( 200. * fb .LT. fa ) THEN 153 fxm = 1. 154 ELSE 155 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 156 IF( 200.*fb + fa.LT.1.e-10 ) THEN 157 fxm = - 1. 158 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 159 fxm = 1. 160 ENDIF 161 ELSE 162 fxm = TANH ( fa/fb ) 163 ENDIF 164 ENDIF 165 166 IF ( xmoy.EQ. 0. ) fxm = 1. 167 IF ( xmoy.EQ. pi ) fxm = -1. 168 169 ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) ) 170 171 ENDDO 172 173 beta = ( grossism * ffdx - pi ) / ( ffdx - pi ) 174 175 IF( 2.*beta - grossism.LE. 0.) THEN 176 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 177 ,tine fxhyp est mauvaise ! ' 178 WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ', 179 , ' et relancer ! *** ' 180 CALL ABORT_GCM("FXHYP", "", 1) 181 ENDIF 182 c 183 c ..... calcul de Xprimt ..... 184 c 185 186 DO i = nmax, nmax2 187 Xprimt(i) = beta + ( grossism - beta ) * fhyp(i) 188 ENDDO 189 c 190 DO i = nmax+1, nmax2 191 Xprimt( nmax2 - i ) = Xprimt( i ) 192 ENDDO 193 c 194 195 c ..... Calcul de Xf ........ 196 197 Xf(0) = - pi 198 199 DO i = nmax +1, nmax2 200 201 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 202 fa = tau* ( dzoom/2. - xmoy ) 203 fb = xmoy * ( pi - xmoy ) 204 205 IF( 200.* fb .LT. - fa ) THEN 206 fxm = - 1. 207 ELSEIF( 200. * fb .LT. fa ) THEN 208 fxm = 1. 209 ELSE 210 fxm = TANH ( fa/fb ) 211 ENDIF 212 213 IF ( xmoy.EQ. 0. ) fxm = 1. 214 IF ( xmoy.EQ. pi ) fxm = -1. 215 xxpr(i) = beta + ( grossism - beta ) * fxm 216 217 ENDDO 218 219 DO i = nmax+1, nmax2 220 xxpr(nmax2-i+1) = xxpr(i) 221 ENDDO 222 223 DO i=1,nmax2 224 Xf(i) = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) ) 225 ENDDO 226 227 228 c ***************************************************************** 229 c 230 231 c ..... xuv = 0. si calcul aux pts scalaires ........ 232 c ..... xuv = 0.5 si calcul aux pts U ........ 233 c 234 WRITE(6,18) 235 c 236 DO 5000 ik = 1, 4 237 238 IF( ik.EQ.1 ) THEN 239 xuv = -0.25 240 ELSE IF ( ik.EQ.2 ) THEN 241 xuv = 0. 242 ELSE IF ( ik.EQ.3 ) THEN 243 xuv = 0.50 244 ELSE IF ( ik.EQ.4 ) THEN 245 xuv = 0.25 246 ENDIF 247 248 xo1 = 0. 249 250 ii1=1 251 ii2=iim 252 IF(ik.EQ.1.and.grossism.EQ.1.) THEN 253 ii1 = 2 254 ii2 = iim+1 255 ENDIF 256 DO 1500 i = ii1, ii2 257 258 xlon2 = - pi + (REAL(i) + xuv - decalx) * depi / REAL(iim) 259 260 Xfi = xlon2 261 c 262 DO 250 it = nmax2,0,-1 263 IF( Xfi.GE.Xf(it)) GO TO 350 264 250 CONTINUE 265 266 it = 0 267 268 350 CONTINUE 269 270 c ...... Calcul de Xf(xi) ...... 271 c 272 xi = xtild(it) 273 274 IF(it.EQ.nmax2) THEN 275 it = nmax2 -1 276 Xf(it+1) = pi 277 ENDIF 278 c ..................................................................... 279 c 280 c Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un 281 c polynome de degre 3 qui passe par les points (Xf(it),xtild(it) ) 282 c et (Xf(it+1),xtild(it+1) ) 283 284 CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1), 285 , xtild(it),xtild(it+1), a0, a1, a2, a3 ) 286 287 Xf1 = Xf(it) 288 Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi 289 290 DO 500 iter = 1,300 291 xi = xi - ( Xf1 - Xfi )/ Xprimin 292 293 IF( ABS(xi-xo1).LE.epsilon) GO TO 550 294 xo1 = xi 295 xi2 = xi * xi 296 Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi 297 Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2 298 500 CONTINUE 299 WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter 300 STOP 6 301 550 CONTINUE 302 303 xxprim(i) = depi/ ( REAL(iim) * Xprimin ) 304 xvrai(i) = xi + xzoom 305 306 1500 CONTINUE 307 308 309 310 IF(ik.EQ.1.and.grossism.EQ.1.) THEN 311 xvrai(1) = xvrai(iip1)-depi 312 xxprim(1) = xxprim(iip1) 313 ENDIF 314 315 DO i = 1 , iim 316 xlon(i) = xvrai(i) 317 xprimm(i) = xxprim(i) 318 ENDDO 319 DO i = 1, iim -1 320 IF( xvrai(i+1). LT. xvrai(i) ) THEN 321 WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i, 322 , ')' 323 STOP 7 324 ENDIF 325 ENDDO 326 c 327 c ... Reorganisation des longitudes pour les avoir entre - pi et pi .. 328 c ........................................................................ 329 330 champmin = 1.e12 331 champmax = -1.e12 332 DO i = 1, iim 333 champmin = MIN( champmin,xvrai(i) ) 334 champmax = MAX( champmax,xvrai(i) ) 335 ENDDO 336 337 IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 ) THEN 338 GO TO 1600 339 ELSE 340 WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi', 341 , ' et pi ' 342 c 343 IF( xzoom.LE.0.) THEN 344 IF( ik.EQ. 1 ) THEN 345 DO i = 1, iim 346 IF( xvrai(i).GE. - pi ) GO TO 80 347 ENDDO 348 WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! ' 349 STOP 8 350 80 CONTINUE 351 is2 = i 352 ENDIF 353 354 IF( is2.NE. 1 ) THEN 355 DO ii = is2 , iim 356 xlon (ii-is2+1) = xvrai(ii) 357 xprimm(ii-is2+1) = xxprim(ii) 358 ENDDO 359 DO ii = 1 , is2 -1 360 xlon (ii+iim-is2+1) = xvrai(ii) + depi 361 xprimm(ii+iim-is2+1) = xxprim(ii) 362 ENDDO 363 ENDIF 364 ELSE 365 IF( ik.EQ.1 ) THEN 366 DO i = iim,1,-1 367 IF( xvrai(i).LE. pi ) GO TO 90 368 ENDDO 369 WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! ' 370 STOP 9 371 90 CONTINUE 372 is2 = i 373 ENDIF 374 idif = iim -is2 375 DO ii = 1, is2 376 xlon (ii+idif) = xvrai(ii) 377 xprimm(ii+idif) = xxprim(ii) 378 ENDDO 379 DO ii = 1, idif 380 xlon (ii) = xvrai (ii+is2) - depi 381 xprimm(ii) = xxprim(ii+is2) 382 ENDDO 383 ENDIF 384 ENDIF 385 c 386 c ......... Fin de la reorganisation ............................ 387 388 1600 CONTINUE 389 390 391 xlon ( iip1) = xlon(1) + depi 392 xprimm( iip1 ) = xprimm (1 ) 393 394 DO i = 1, iim+1 395 xvrai(i) = xlon(i)*180./pi 396 ENDDO 397 398 IF( ik.EQ.1 ) THEN 399 c WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) ' 400 c WRITE(6,18) 401 c WRITE(6,68) xvrai 402 c WRITE(6,*) ' XPRIM k ',ik 403 c WRITE(6,566) xprimm 404 405 DO i = 1,iim +1 406 rlonm025(i) = xlon( i ) 407 xprimm025(i) = xprimm(i) 408 ENDDO 409 ELSE IF( ik.EQ.2 ) THEN 410 c WRITE(6,18) 411 c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) ' 412 c WRITE(6,68) xvrai 413 c WRITE(6,*) ' XPRIM k ',ik 414 c WRITE(6,566) xprimm 415 416 DO i = 1,iim + 1 417 rlonv(i) = xlon( i ) 418 xprimv(i) = xprimm(i) 419 ENDDO 420 421 ELSE IF( ik.EQ.3) THEN 422 c WRITE(6,18) 423 c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) ' 424 c WRITE(6,68) xvrai 425 c WRITE(6,*) ' XPRIM ik ',ik 426 c WRITE(6,566) xprimm 427 428 DO i = 1,iim + 1 429 rlonu(i) = xlon( i ) 430 xprimu(i) = xprimm(i) 431 ENDDO 432 433 ELSE IF( ik.EQ.4 ) THEN 434 c WRITE(6,18) 435 c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) ' 436 c WRITE(6,68) xvrai 437 c WRITE(6,*) ' XPRIM ik ',ik 438 c WRITE(6,566) xprimm 439 440 DO i = 1,iim + 1 441 rlonp025(i) = xlon( i ) 442 xprimp025(i) = xprimm(i) 443 ENDDO 444 445 ENDIF 446 447 5000 CONTINUE 448 c 449 WRITE(6,18) 450 c 451 c ........... fin de la boucle do 5000 ............ 452 453 DO i = 1, iim 454 xlon(i) = rlonv(i+1) - rlonv(i) 455 ENDDO 456 champmin = 1.e12 457 champmax = -1.e12 458 DO i = 1, iim 459 champmin = MIN( champmin, xlon(i) ) 460 champmax = MAX( champmax, xlon(i) ) 461 ENDDO 462 champmin = champmin * 180./pi 463 champmax = champmax * 180./pi 464 465 18 FORMAT(/) 466 24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3) 467 68 FORMAT(1x,7f9.2) 468 566 FORMAT(1x,7f9.4) 469 470 RETURN 471 END 1 module fxhyp_m 2 3 IMPLICIT NONE 4 5 contains 6 7 SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) 8 9 ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32 10 ! Author: P. Le Van, from formulas by R. Sadourny 11 12 ! Calcule les longitudes et dérivées dans la grille du GCM pour 13 ! une fonction f(x) à dérivée tangente hyperbolique. 14 15 ! Il vaut mieux avoir : grossismx \times dzoom < pi 16 17 ! Le premier point scalaire pour une grille regulière (grossismx = 18 ! 1., taux=0., clon=0.) est à - 180 degrés. 19 20 use arth_m, only: arth 21 use invert_zoom_x_m, only: invert_zoom_x, nmax 22 use nrtype, only: pi, pi_d, twopi, twopi_d, k8 23 use principal_cshift_m, only: principal_cshift 24 25 include "dimensions.h" 26 ! for iim 27 28 include "serre.h" 29 ! for clon, grossismx, dzoomx, taux 30 31 REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1) 32 real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1) 33 34 ! Local: 35 real rlonm025(iim + 1), rlonp025(iim + 1) 36 REAL dzoom, step 37 real d_rlonv(iim) 38 REAL(K8) xtild(0:2 * nmax) 39 REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax) 40 REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax) 41 REAL(K8) fa, fb 42 INTEGER i, is2 43 REAL(K8) xmoy, fxm 44 45 !---------------------------------------------------------------------- 46 47 print *, "Call sequence information: fxhyp" 48 49 test_iim: if (iim==1) then 50 rlonv(1)=0. 51 rlonu(1)=pi 52 rlonv(2)=rlonv(1)+twopi 53 rlonu(2)=rlonu(1)+twopi 54 55 xprimm025(:)=1. 56 xprimv(:)=1. 57 xprimu(:)=1. 58 xprimp025(:)=1. 59 else test_iim 60 test_grossismx: if (grossismx == 1.) then 61 step = twopi / iim 62 63 xprimm025(:iim) = step 64 xprimp025(:iim) = step 65 xprimv(:iim) = step 66 xprimu(:iim) = step 67 68 rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim) 69 rlonm025(:iim) = rlonv(:iim) - 0.25 * step 70 rlonp025(:iim) = rlonv(:iim) + 0.25 * step 71 rlonu(:iim) = rlonv(:iim) + 0.5 * step 72 else test_grossismx 73 dzoom = dzoomx * twopi_d 74 xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1) 75 76 ! Compute fhyp: 77 DO i = nmax, 2 * nmax 78 fa = taux * (dzoom / 2. - xtild(i)) 79 fb = xtild(i) * (pi_d - xtild(i)) 80 81 IF (200. * fb < - fa) THEN 82 fhyp(i) = - 1. 83 ELSE IF (200. * fb < fa) THEN 84 fhyp(i) = 1. 85 ELSE 86 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN 87 IF (200. * fb + fa < 1e-10) THEN 88 fhyp(i) = - 1. 89 ELSE IF (200. * fb - fa < 1e-10) THEN 90 fhyp(i) = 1. 91 END IF 92 ELSE 93 fhyp(i) = TANH(fa / fb) 94 END IF 95 END IF 96 97 IF (xtild(i) == 0.) fhyp(i) = 1. 98 IF (xtild(i) == pi_d) fhyp(i) = -1. 99 END DO 100 101 ! Calcul de beta 102 103 ffdx = 0. 104 105 DO i = nmax + 1, 2 * nmax 106 xmoy = 0.5 * (xtild(i-1) + xtild(i)) 107 fa = taux * (dzoom / 2. - xmoy) 108 fb = xmoy * (pi_d - xmoy) 109 110 IF (200. * fb < - fa) THEN 111 fxm = - 1. 112 ELSE IF (200. * fb < fa) THEN 113 fxm = 1. 114 ELSE 115 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN 116 IF (200. * fb + fa < 1e-10) THEN 117 fxm = - 1. 118 ELSE IF (200. * fb - fa < 1e-10) THEN 119 fxm = 1. 120 END IF 121 ELSE 122 fxm = TANH(fa / fb) 123 END IF 124 END IF 125 126 IF (xmoy == 0.) fxm = 1. 127 IF (xmoy == pi_d) fxm = -1. 128 129 ffdx = ffdx + fxm * (xtild(i) - xtild(i-1)) 130 END DO 131 132 print *, "ffdx = ", ffdx 133 beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d) 134 print *, "beta = ", beta 135 136 IF (2. * beta - grossismx <= 0.) THEN 137 print *, 'Bad choice of grossismx, taux, dzoomx.' 138 print *, 'Decrease dzoomx or grossismx.' 139 STOP 1 140 END IF 141 142 ! calcul de Xprimt 143 Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp 144 xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1) 145 146 ! Calcul de Xf 147 148 DO i = nmax + 1, 2 * nmax 149 xmoy = 0.5 * (xtild(i-1) + xtild(i)) 150 fa = taux * (dzoom / 2. - xmoy) 151 fb = xmoy * (pi_d - xmoy) 152 153 IF (200. * fb < - fa) THEN 154 fxm = - 1. 155 ELSE IF (200. * fb < fa) THEN 156 fxm = 1. 157 ELSE 158 fxm = TANH(fa / fb) 159 END IF 160 161 IF (xmoy == 0.) fxm = 1. 162 IF (xmoy == pi_d) fxm = -1. 163 xxpr(i) = beta + (grossismx - beta) * fxm 164 END DO 165 166 xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1) 167 168 Xf(0) = - pi_d 169 170 DO i=1, 2 * nmax - 1 171 Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1)) 172 END DO 173 174 Xf(2 * nmax) = pi_d 175 176 call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), & 177 xprimm025(:iim), xuv = - 0.25_k8) 178 call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), & 179 xuv = 0._k8) 180 call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), & 181 xuv = 0.5_k8) 182 call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), & 183 xprimp025(:iim), xuv = 0.25_k8) 184 end if test_grossismx 185 186 is2 = 0 187 188 IF (MINval(rlonm025(:iim)) < - pi - 0.1 & 189 .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN 190 IF (clon <= 0.) THEN 191 is2 = 1 192 193 do while (rlonm025(is2) < - pi .and. is2 < iim) 194 is2 = is2 + 1 195 end do 196 197 if (rlonm025(is2) < - pi) then 198 print *, 'Rlonm025 plus petit que - pi !' 199 STOP 1 200 end if 201 ELSE 202 is2 = iim 203 204 do while (rlonm025(is2) > pi .and. is2 > 1) 205 is2 = is2 - 1 206 end do 207 208 if (rlonm025(is2) > pi) then 209 print *, 'Rlonm025 plus grand que pi !' 210 STOP 1 211 end if 212 END IF 213 END IF 214 215 call principal_cshift(is2, rlonm025, xprimm025) 216 call principal_cshift(is2, rlonv, xprimv) 217 call principal_cshift(is2, rlonu, xprimu) 218 call principal_cshift(is2, rlonp025, xprimp025) 219 220 forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i) 221 print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, & 222 "degrees" 223 print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, & 224 "degrees" 225 226 ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu: 227 DO i = 1, iim + 1 228 IF (rlonp025(i) < rlonv(i)) THEN 229 print *, 'rlonp025(', i, ') = ', rlonp025(i) 230 print *, "< rlonv(", i, ") = ", rlonv(i) 231 STOP 1 232 END IF 233 234 IF (rlonv(i) < rlonm025(i)) THEN 235 print *, 'rlonv(', i, ') = ', rlonv(i) 236 print *, "< rlonm025(", i, ") = ", rlonm025(i) 237 STOP 1 238 END IF 239 240 IF (rlonp025(i) > rlonu(i)) THEN 241 print *, 'rlonp025(', i, ') = ', rlonp025(i) 242 print *, "> rlonu(", i, ") = ", rlonu(i) 243 STOP 1 244 END IF 245 END DO 246 end if test_iim 247 248 END SUBROUTINE fxhyp 249 250 end module fxhyp_m -
trunk/LMDZ.COMMON/libf/dyn3d_common/fyhyp_m.F90
r1440 r1441 1 ! 2 ! $Id: fyhyp.F 1403 2010-07-01 09:02:53Z fairhead $ 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 use nrtype, only: k8 20 21 include "dimensions.h" 22 ! for jjm 23 24 include "serre.h" 25 ! for clat, grossismy, dzoomy, tauy 26 27 REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1) 28 REAL, intent(out):: rlatv(jjm) 29 real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm) 30 31 ! Local: 32 33 REAL(K8) champmin, champmax 34 INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax 35 REAL dzoom ! distance totale de la zone du zoom (en radians) 36 REAL(K8) ylat(jjm + 1), yprim(jjm + 1) 37 REAL(K8) yuv 38 REAL(K8), save:: yt(0:nmax2) 39 REAL(K8) fhyp(0:nmax2), beta 40 REAL(K8), save:: ytprim(0:nmax2) 41 REAL(K8) fxm(0:nmax2) 42 REAL(K8), save:: yf(0:nmax2) 43 REAL(K8) yypr(0:nmax2) 44 REAL(K8) yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1) 45 REAL(K8) pi, pis2, epsilon, y0, pisjm 46 REAL(K8) yo1, yi, ylon2, ymoy, yprimin 47 REAL(K8) yfi, yf1, ffdy 48 REAL(K8) ypn, deply, y00 49 SAVE y00, deply 50 51 INTEGER i, j, it, ik, iter, jlat 52 INTEGER jpn, jjpn 53 SAVE jpn 54 REAL(K8) a0, a1, a2, a3, yi2, heavyy0, heavyy0m 55 REAL(K8) fa(0:nmax2), fb(0:nmax2) 56 REAL y0min, y0max 57 58 REAL(K8) heavyside 59 60 !------------------------------------------------------------------- 61 62 print *, "Call sequence information: fyhyp" 63 64 pi = 2.*asin(1.) 65 pis2 = pi/2. 66 pisjm = pi/real(jjm) 67 epsilon = 1e-3 68 y0 = clat*pi/180. 69 dzoom = dzoomy*pi 70 print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):' 71 print *, y0, grossismy, tauy, dzoom 72 73 DO i = 0, nmax2 74 yt(i) = -pis2 + real(i)*pi/nmax2 75 END DO 76 77 heavyy0m = heavyside(-y0) 78 heavyy0 = heavyside(y0) 79 y0min = 2.*y0*heavyy0m - pis2 80 y0max = 2.*y0*heavyy0 + pis2 81 82 fa = 999.999 83 fb = 999.999 84 85 DO i = 0, nmax2 86 IF (yt(i)<y0) THEN 87 fa(i) = tauy*(yt(i)-y0 + dzoom/2.) 88 fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i)) 89 ELSE IF (yt(i)>y0) THEN 90 fa(i) = tauy*(y0-yt(i) + dzoom/2.) 91 fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0) 92 END IF 93 94 IF (200.*fb(i)<-fa(i)) THEN 95 fhyp(i) = -1. 96 ELSE IF (200.*fb(i)<fa(i)) THEN 97 fhyp(i) = 1. 87 98 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 99 fhyp(i) = tanh(fa(i)/fb(i)) 100 END IF 101 102 IF (yt(i)==y0) fhyp(i) = 1. 103 IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1. 104 END DO 105 106 ! Calcul de beta 107 108 ffdy = 0. 109 110 DO i = 1, nmax2 111 ymoy = 0.5*(yt(i-1) + yt(i)) 112 IF (ymoy<y0) THEN 113 fa(i) = tauy*(ymoy-y0 + dzoom/2.) 114 fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy) 115 ELSE IF (ymoy>y0) THEN 116 fa(i) = tauy*(y0-ymoy + dzoom/2.) 117 fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0) 118 END IF 119 120 IF (200.*fb(i)<-fa(i)) THEN 121 fxm(i) = -1. 122 ELSE IF (200.*fb(i)<fa(i)) THEN 123 fxm(i) = 1. 124 ELSE 125 fxm(i) = tanh(fa(i)/fb(i)) 126 END IF 127 IF (ymoy==y0) fxm(i) = 1. 128 IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1. 129 ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1)) 130 END DO 131 132 beta = (grossismy*ffdy-pi)/(ffdy-pi) 133 134 IF (2. * beta - grossismy <= 0.) THEN 135 print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' & 136 // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' & 137 // 'dzoomy et relancer.' 138 STOP 1 139 END IF 140 141 ! calcul de Ytprim 142 143 DO i = 0, nmax2 144 ytprim(i) = beta + (grossismy-beta)*fhyp(i) 145 END DO 146 147 ! Calcul de Yf 148 149 yf(0) = -pis2 150 DO i = 1, nmax2 151 yypr(i) = beta + (grossismy-beta)*fxm(i) 152 END DO 153 154 DO i = 1, nmax2 155 yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1)) 156 END DO 157 158 ! yuv = 0. si calcul des latitudes aux pts. U 159 ! yuv = 0.5 si calcul des latitudes aux pts. V 160 161 loop_ik: DO ik = 1, 4 162 IF (ik==1) THEN 163 yuv = 0. 164 jlat = jjm + 1 165 ELSE IF (ik==2) THEN 166 yuv = 0.5 167 jlat = jjm 168 ELSE IF (ik==3) THEN 169 yuv = 0.25 170 jlat = jjm 171 ELSE IF (ik==4) THEN 172 yuv = 0.75 173 jlat = jjm 174 END IF 175 176 yo1 = 0. 296 177 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 178 yo1 = 0. 179 ylon2 = -pis2 + pisjm*(real(j) + yuv-1.) 180 yfi = ylon2 181 182 it = nmax2 183 DO while (it >= 1 .and. yfi < yf(it)) 184 it = it - 1 185 END DO 186 187 yi = yt(it) 188 IF (it==nmax2) THEN 189 it = nmax2 - 1 190 yf(it + 1) = pis2 191 END IF 192 193 ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi) 194 ! et Y'(yi) 195 196 CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), & 197 yt(it), yt(it + 1), a0, a1, a2, a3) 198 199 yf1 = yf(it) 200 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi 201 202 iter = 1 203 DO 204 yi = yi - (yf1-yfi)/yprimin 205 IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit 206 yo1 = yi 207 yi2 = yi*yi 208 yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi 209 yprimin = a1 + 2.*a2*yi + 3.*a3*yi2 210 END DO 211 if (abs(yi-yo1) > epsilon) then 212 print *, 'Pas de solution.', j, ylon2 213 STOP 1 214 end if 215 216 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi 217 yprim(j) = pi/(jjm*yprimin) 218 yvrai(j) = yi 219 END DO 220 221 DO j = 1, jlat - 1 222 IF (yvrai(j + 1)<yvrai(j)) THEN 223 print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', & 224 j, ')' 225 STOP 1 226 END IF 227 END DO 228 229 print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2' 230 231 IF (ik==1) THEN 232 ypn = pis2 233 DO j = jjm + 1, 1, -1 234 IF (yvrai(j)<=ypn) exit 235 END DO 236 237 jpn = j 238 y00 = yvrai(jpn) 239 deply = pis2 - y00 240 END IF 241 242 DO j = 1, jjm + 1 - jpn 243 ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1) 244 yprimm(j) = yprim(jpn + j-1) 245 END DO 246 247 jjpn = jpn 248 IF (jlat==jjm) jjpn = jpn - 1 249 250 DO j = 1, jjpn 251 ylatt(j + jjm + 1-jpn) = yvrai(j) + deply 252 yprimm(j + jjm + 1-jpn) = yprim(j) 253 END DO 254 255 ! Fin de la reorganisation 256 257 DO j = 1, jlat 258 ylat(j) = ylatt(jlat + 1-j) 259 yprim(j) = yprimm(jlat + 1-j) 260 END DO 261 262 DO j = 1, jlat 263 yvrai(j) = ylat(j)*180./pi 264 END DO 265 266 IF (ik==1) THEN 267 DO j = 1, jjm + 1 268 rlatu(j) = ylat(j) 269 yyprimu(j) = yprim(j) 270 END DO 271 ELSE IF (ik==2) THEN 272 DO j = 1, jjm 273 rlatv(j) = ylat(j) 274 END DO 275 ELSE IF (ik==3) THEN 276 DO j = 1, jjm 277 rlatu2(j) = ylat(j) 278 yprimu2(j) = yprim(j) 279 END DO 280 ELSE IF (ik==4) THEN 281 DO j = 1, jjm 282 rlatu1(j) = ylat(j) 283 yprimu1(j) = yprim(j) 284 END DO 285 END IF 286 END DO loop_ik 287 288 DO j = 1, jjm 289 ylat(j) = rlatu(j) - rlatu(j + 1) 290 END DO 291 champmin = 1e12 292 champmax = -1e12 293 DO j = 1, jjm 294 champmin = min(champmin, ylat(j)) 295 champmax = max(champmax, ylat(j)) 296 END DO 297 champmin = champmin*180./pi 298 champmax = champmax*180./pi 299 300 DO j = 1, jjm 301 IF (rlatu1(j) <= rlatu2(j)) THEN 302 print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j 303 STOP 13 304 ENDIF 305 306 IF (rlatu2(j) <= rlatu(j+1)) THEN 307 print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j 308 STOP 14 309 ENDIF 310 311 IF (rlatu(j) <= rlatu1(j)) THEN 312 print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j 313 STOP 15 314 ENDIF 315 316 IF (rlatv(j) <= rlatu2(j)) THEN 317 print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j 318 STOP 16 319 ENDIF 320 321 IF (rlatv(j) >= rlatu1(j)) THEN 322 print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j 323 STOP 17 324 ENDIF 325 326 IF (rlatv(j) >= rlatu(j)) THEN 327 print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j 328 STOP 18 329 ENDIF 330 ENDDO 331 332 print *, 'Latitudes' 333 print 3, champmin, champmax 334 335 3 Format(1x, ' Au centre du zoom, la longueur de la maille est', & 336 ' d environ ', f0.2, ' degres ', /, & 337 ' alors que la maille en dehors de la zone du zoom est ', & 338 "d'environ ", f0.2, ' degres ') 339 340 END SUBROUTINE fyhyp 341 342 end module fyhyp_m -
trunk/LMDZ.COMMON/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90
r1422 r1441 31 31 INTEGER out_latudim,out_latvdim,out_dim(3) 32 32 INTEGER out_levdim 33 34 INTEGER, PARAMETER :: longcles = 2035 REAL clesphy0(longcles)36 33 37 34 INTEGER start(4),COUNT(4) … … 59 56 pa= 50000. 60 57 61 CALL conf_gcm( 99, .TRUE. , clesphy0)58 CALL conf_gcm( 99, .TRUE. ) 62 59 CALL iniconst 63 60 CALL inigeom -
trunk/LMDZ.COMMON/libf/dyn3d_common/inigeom.F
r1422 r1441 20 20 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, 21 21 . alphax,alphay,taux,tauy,transx,transy,pxo,pyo 22 use fxhyp_m, only: fxhyp 23 use fyhyp_m, only: fyhyp 22 24 23 25 IMPLICIT NONE … … 266 268 WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' 267 269 268 CALL fxyhyper( clat, grossismy, dzoomy, tauy , 269 , clon, grossismx, dzoomx, taux , 270 , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 , 271 , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 ) 272 270 CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) 271 CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) 273 272 274 273 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
