--- trunk/libf/dyn3d/inigeom.f90 2010/03/05 16:43:45 25 +++ trunk/libf/dyn3d/inigeom.f90 2013/06/24 15:39:52 70 @@ -7,102 +7,48 @@ SUBROUTINE inigeom ! Auteur : P. Le Van - ! Version du 01/04/2001 ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes ! endroits que les aires aireij1_2d, ..., aireij4_2d. - ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à dérivée - ! tangente hyperbolique - ! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2) + ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à + ! dérivée tangente hyperbolique. Calcul des coefficients cu_2d, + ! cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les coefficients cu_2d et cv_2d + ! permettent de passer des vitesses naturelles aux vitesses + ! covariantes et contravariantes, ou vice-versa. + + ! On a : + ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d + ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d + + ! On en tire : + ! u(covariant) = cu_2d * cu_2d * u(contravariant) + ! v(covariant) = cv_2d * cv_2d * v(contravariant) + + ! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2 + ! et - jm / 2 <= Y <= jm / 2 + + ! x est la longitude du point en radians. + ! y est la latitude du point en radians. + ! + ! On a : cu_2d(i, j) = rad * cos(y) * dx / dX + ! cv(j) = rad * dy / dY + ! aire_2d(i, j) = cu_2d(i, j) * cv(j) + ! + ! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que + ! dépendant de j uniquement, sera ici indicé aussi en i pour un + ! adressage plus facile en ij. + + ! xprimu et xprimv sont respectivement les valeurs de dx / dX aux + ! points u et v. yprimu et yprimv sont respectivement les valeurs + ! de dy / dY aux points u et v. rlatu et rlatv sont respectivement + ! les valeurs de la latitude aux points u et v. cvu et cv_2d sont + ! respectivement les valeurs de cv_2d aux points u et v. - ! les coef. ( cu_2d, cv_2d ) permettent de passer des vitesses naturelles - ! aux vitesses covariantes et contravariantes , ou vice-versa ... + ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d + ! aux points u, v, scalaires, et z. Cf. "inigeom.txt". - ! on a : u (covariant) = cu_2d * u (naturel) , u(contrav)= u(nat)/cu_2d - ! v (covariant) = cv_2d * v (naturel) , v(contrav)= v(nat)/cv_2d - - ! on en tire : u(covariant) = cu_2d * cu_2d * u(contravariant) - ! v(covariant) = cv_2d * cv_2d * v(contravariant) - - ! on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2 - ! = = - ! et - jm/2 < Y < jm/2 - ! = = - - ! . x est la longitude du point en radians . - ! . y est la latitude du point en radians . - ! . . - ! . on a : cu_2d(i, j) = rad * COS(y) * dx/dX . - ! . cv( j ) = rad * dy/dY . - ! . aire_2d(i, j) = cu_2d(i, j) * cv(j) . - ! . . - ! . y, dx/dX, dy/dY calcules aux points concernes . - ! , - ! cv , bien que dependant de j uniquement, sera ici indice aussi en i - ! pour un adressage plus facile en ij . - - ! ************** aux points u et v , ***************** - ! xprimu et xprimv sont respectivement les valeurs de dx/dX - ! yprimu et yprimv . . . . . . . . . . . dy/dY - ! rlatu et rlatv . . . . . . . . . . .la latitude - ! cvu et cv_2d . . . . . . . . . . . cv_2d - - ! ************** aux points u, v, scalaires, et z **************** - ! cu_2d, cuv, cuscal, cuz sont respectiv. les valeurs de cu_2d - - ! Exemple de distribution de variables sur la grille dans le - ! domaine de travail ( X, Y ) . - ! DX=DY= 1 - - ! + represente un point scalaire ( p.exp la pression ) - ! > represente la composante zonale du vent - ! V represente la composante meridienne du vent - ! o represente la vorticite - - ! ---- , car aux poles , les comp.zonales covariantes sont nulles - - ! i -> - - ! 1 2 3 4 5 6 7 8 - ! j - ! v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- - - ! V o V o V o V o V o V o V o V o - - ! 2 + > + > + > + > + > + > + > + > - - ! V o V o V o V o V o V o V o V o - - ! 3 + > + > + > + > + > + > + > + > - - ! V o V o V o V o V o V o V o V o - - ! 4 + > + > + > + > + > + > + > + > - - ! V o V o V o V o V o V o V o V o - - ! 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- - - ! Ci-dessus, on voit que le nombre de pts.en longitude est egal - ! a IM = 8 - ! De meme , le nombre d'intervalles entre les 2 poles est egal - ! a JM = 4 - - ! Les points scalaires ( + ) correspondent donc a des valeurs - ! entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) . - - ! Les vents U ( > ) correspondent a des valeurs semi- - ! entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1) - - ! Les vents V ( V ) correspondent a des valeurs entieres - ! de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5) - - USE dimens_m, ONLY : iim, jjm - USE paramet_m, ONLY : iip1, jjp1 - USE comconst, ONLY : g, omeg, pi, rad - USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh - USE logic, ONLY : fxyhypb, ysinus + USE comconst, ONLY : g, omeg, rad USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, & alpha1p2_2d, alpha1p4_2d, alpha1_2d, & alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, & @@ -112,34 +58,34 @@ rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, & unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, & unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv + USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh + use conf_gcm_m, ONLY : fxyhypb, ysinus + USE dimens_m, ONLY : iim, jjm + use fxy_m, only: fxy + use fxyhyper_m, only: fxyhyper + use jumble, only: new_unit + use nr_util, only: pi + USE paramet_m, ONLY : iip1, jjp1 USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, & grossismy, pxo, pyo, taux, tauy, transx, transy ! Variables locales - INTEGER i, j, itmax, itmay, iter + INTEGER i, j, itmax, itmay, iter, unit REAL cvu(iip1, jjp1), cuv(iip1, jjm) - REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm + REAL ai14, ai23, airez, un4rad2 REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm REAL coslatm, coslatp, radclatm, radclatp - REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), & - cuij4(iip1, jjp1) - REAL cvij1(iip1, jjp1), cvij2(iip1, jjp1), cvij3(iip1, jjp1), & - cvij4(iip1, jjp1) - REAL rlonvv(iip1), rlatuu(jjp1) - REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm), yprimv(jjm), & - yprimu(jjp1) + REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m + REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m + REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm) + real yprimv(jjm), yprimu(jjp1) REAL gamdi_gdiv, gamdi_grot, gamdi_h - REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1) - SAVE rlatu1, yprimu1, rlatu2, yprimu2, yprimv, yprimu - SAVE rlonm025, xprimm025, rlonp025, xprimp025 - - real aireij1_2d(iim + 1, jjm + 1) - real aireij2_2d(iim + 1, jjm + 1) - real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1) + real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, & + aireij4_2d ! in m2 real airuscv2_2d(iim + 1, jjm) - real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm) + real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm) real aivscu2gam_2d(iim + 1, jjm) !------------------------------------------------------------------ @@ -147,17 +93,17 @@ PRINT *, 'Call sequence information: inigeom' IF (nitergdiv/=2) THEN - gamdi_gdiv = coefdis/(real(nitergdiv)-2.) + gamdi_gdiv = coefdis / (real(nitergdiv)-2.) ELSE gamdi_gdiv = 0. END IF IF (nitergrot/=2) THEN - gamdi_grot = coefdis/(real(nitergrot)-2.) + gamdi_grot = coefdis / (real(nitergrot)-2.) ELSE gamdi_grot = 0. END IF IF (niterh/=2) THEN - gamdi_h = coefdis/(real(niterh)-2.) + gamdi_h = coefdis / (real(niterh)-2.) ELSE gamdi_h = 0. END IF @@ -166,25 +112,21 @@ print *, "gamdi_grot = ", gamdi_grot print *, "gamdi_h = ", gamdi_h - WRITE (6, 990) - - IF ( .NOT. fxyhypb) THEN + IF (.NOT. fxyhypb) THEN IF (ysinus) THEN - print *, ' *** Inigeom , Y = Sinus ( Latitude ) *** ' - - ! utilisation de f(x, y ) avec y = sinus de la latitude ... - + print *, ' Inigeom, Y = Sinus (Latitude) ' + ! utilisation de f(x, y) avec y = sinus de la latitude CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, & rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, & xprimm025, rlonp025, xprimp025) ELSE - print *, '*** Inigeom , Y = Latitude , der. sinusoid . ***' - ! utilisation de f(x, y) a tangente sinusoidale , y etant la latit + print *, 'Inigeom, Y = Latitude, der. sinusoid .' + ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit - pxo = clon*pi/180. - pyo = 2.*clat*pi/180. + pxo = clon * pi / 180. + pyo = 2. * clat * pi / 180. - ! determination de transx ( pour le zoom ) par Newton-Raphson . + ! determination de transx (pour le zoom) par Newton-Raphson itmax = 10 eps = .1E-7 @@ -192,9 +134,9 @@ xo1 = 0. DO iter = 1, itmax x1 = xo1 - f = x1 + alphax*sin(x1-pxo) - df = 1. + alphax*cos(x1-pxo) - x1 = x1 - f/df + f = x1 + alphax * sin(x1-pxo) + df = 1. + alphax * cos(x1-pxo) + x1 = x1 - f / df xdm = abs(x1-xo1) IF (xdm<=eps) EXIT xo1 = x1 @@ -208,9 +150,9 @@ yo1 = 0. DO iter = 1, itmay y1 = yo1 - f = y1 + alphay*sin(y1-pyo) - df = 1. + alphay*cos(y1-pyo) - y1 = y1 - f/df + f = y1 + alphay * sin(y1-pyo) + df = 1. + alphay * cos(y1-pyo) + y1 = y1 - f / df ydm = abs(y1-yo1) IF (ydm<=eps) EXIT yo1 = y1 @@ -223,190 +165,113 @@ rlonp025, xprimp025) END IF ELSE - ! .... Utilisation de fxyhyper , f(x, y) a derivee tangente hyperbol. - print *, '*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' + ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique + print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique' CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, & taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, & yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, & rlonp025, xprimp025) END IF - rlatu(1) = asin(1.) + rlatu(1) = pi / 2. rlatu(jjp1) = -rlatu(1) - ! .... calcul aux poles .... + ! Calcul aux pôles yprimu(1) = 0. yprimu(jjp1) = 0. - un4rad2 = 0.25*rad*rad - - ! calcul des aires ( aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez ) - ! - et de fext_2d , force de coriolis extensive . - - ! A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y) , sont - ! affectees 4 aires entourant P , calculees respectivement aux points - ! ( i + 1/4, j - 1/4 ) : aireij1_2d (i, j) - ! ( i + 1/4, j + 1/4 ) : aireij2_2d (i, j) - ! ( i - 1/4, j + 1/4 ) : aireij3_2d (i, j) - ! ( i - 1/4, j - 1/4 ) : aireij4_2d (i, j) - - ! , - ! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y). - ! Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme - ! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont - ! affectees au - ! point (i, j) . - ! On definit en outre les coefficients alpha comme etant egaux a - ! (aireij / aire_2d), c.a.d par exp. - ! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j) - - ! De meme, toute aire centree en 1 point U est egale a la somme des - ! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant - ! le point U. - ! Idem pour airev_2d, airez . + un4rad2 = 0.25 * rad * rad - ! On a , pour chaque maille : dX = dY = 1 - - ! . V - - ! aireij4_2d . . aireij1_2d - - ! U . . P . U - - ! aireij3_2d . . aireij2_2d - - ! . V - - ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d, - ! aireij3_2d, aireij4_2d - ! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations - ! elementaires - ! cuij et les 4 elongat. cvij qui sont calculees aux memes - ! endroits que les aireij . - - ! ....... do 35 : boucle sur les jjm + 1 latitudes ..... - - DO j = 1, jjp1 - - IF (j==1) THEN - - yprm = yprimu1(j) - rlatm = rlatu1(j) - - coslatm = cos(rlatm) - radclatm = 0.5*rad*coslatm - - DO i = 1, iim - xprp = xprimp025(i) - xprm = xprimm025(i) - aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm - aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm - cuij2(i, 1) = radclatm*xprp - cuij3(i, 1) = radclatm*xprm - cvij2(i, 1) = 0.5*rad*yprm - cvij3(i, 1) = cvij2(i, 1) - END DO - - DO i = 1, iim - aireij1_2d(i, 1) = 0. - aireij4_2d(i, 1) = 0. - cuij1(i, 1) = 0. - cuij4(i, 1) = 0. - cvij1(i, 1) = 0. - cvij4(i, 1) = 0. - END DO - - END IF - - IF (j==jjp1) THEN - yprp = yprimu2(j-1) - rlatp = rlatu2(j-1) - - coslatp = cos(rlatp) - radclatp = 0.5*rad*coslatp - - DO i = 1, iim - xprp = xprimp025(i) - xprm = xprimm025(i) - aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp - aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp - cuij1(i, jjp1) = radclatp*xprp - cuij4(i, jjp1) = radclatp*xprm - cvij1(i, jjp1) = 0.5*rad*yprp - cvij4(i, jjp1) = cvij1(i, jjp1) - END DO - - DO i = 1, iim - aireij2_2d(i, jjp1) = 0. - aireij3_2d(i, jjp1) = 0. - cvij2(i, jjp1) = 0. - cvij3(i, jjp1) = 0. - cuij2(i, jjp1) = 0. - cuij3(i, jjp1) = 0. - END DO - - END IF - - IF (j>1 .AND. j