! $Id: inigeom.f90 5159 2024-08-02 19:58:25Z fairhead $ SUBROUTINE inigeom ! Auteur : P. Le Van ! ............ Version du 01/04/2001 ........................ ! Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en- ! endroits que les aires aireij1,..aireij4 . ! Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol. USE fxhyp_m, ONLY: fxhyp USE fyhyp_m, ONLY: fyhyp USE comconst_mod, ONLY: pi, g, omeg, rad USE logic_mod, ONLY: fxyhypb, ysinus USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, & alphax,alphay,taux,tauy,transx,transy,pxo,pyo USE lmdz_ssum_scopy, ONLY: ssum USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, & tetagrot, tetatemp, coefdis, vert_prof_dissip USE lmdz_comgeom2 USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! !----------------------------------------------------------------------- ! .... Variables locales .... INTEGER :: i,j,itmax,itmay,iter REAL :: cvu(iip1,jjp1),cuv(iip1,jjm) REAL :: ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm 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 :: 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 ! ------------------------------------------------------------------ ! - - ! - calcul des coeff. ( cu, cv , 1./cu**2, 1./cv**2 ) - ! - - ! ------------------------------------------------------------------ ! les coef. ( cu, cv ) permettent de passer des vitesses naturelles ! aux vitesses covariantes et contravariantes , ou vice-versa ... ! on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu ! v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv ! on en tire : u(covariant) = cu * cu * u(contravariant) ! v(covariant) = cv * cv * 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(i,j) = rad * COS(y) * dx/dX . ! . cv( j ) = rad * dy/dY . ! . aire(i,j) = cu(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 . . . . . . . . . . . cv ! ************** aux points u, v, scalaires, et z **************** ! cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu ! 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) WRITE(6,3) 3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ', & //5x,' Calcul des elongations cu et cv comme sommes des 4 ' / & 5 x,' elong. cuij1, .. 4 , cvij1,.. 4 qui les entourent , aux& & '/ 5x,' memes endroits que les aires aireij1,...j4 . ' / ) IF( nitergdiv/=2 ) THEN gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. ) ELSE gamdi_gdiv = 0. ENDIF IF( nitergrot/=2 ) THEN gamdi_grot = coefdis/ ( REAL(nitergrot) -2. ) ELSE gamdi_grot = 0. ENDIF IF( niterh/=2 ) THEN gamdi_h = coefdis/ ( REAL(niterh) -2. ) ELSE gamdi_h = 0. ENDIF WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis, & nitergdiv,nitergrot,niterh pi = 2.* ASIN(1.) WRITE(6,990) ! ---------------------------------------------------------------- IF( .NOT.fxyhypb ) THEN IF( ysinus ) THEN WRITE(6,*) ' *** 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 WRITE(6,*) '*** 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. ! .... determination de transx ( pour le zoom ) par Newton-Raphson ... itmax = 10 eps = .1e-7 xo1 = 0. DO iter = 1, itmax x1 = xo1 f = x1+ alphax *SIN(x1-pxo) df = 1.+ alphax *COS(x1-pxo) x1 = x1 - f/df xdm = ABS( x1- xo1 ) IF( xdm<=eps )GO TO 11 xo1 = x1 END DO 11 CONTINUE transx = xo1 itmay = 10 eps = .1e-7 yo1 = 0. DO iter = 1,itmay y1 = yo1 f = y1 + alphay* SIN(y1-pyo) df = 1. + alphay* COS(y1-pyo) y1 = y1 -f/df ydm = ABS(y1-yo1) IF(ydm<=eps) GO TO 17 yo1 = y1 END DO 17 CONTINUE transy = yo1 CALL fxy ( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1, & rlatu2,yprimu2, & rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025) ENDIF ELSE ! .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol. ! ..................................................................... WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) ENDIF ! ------------------------------------------------------------------- rlatu(1) = ASIN(1.) rlatu(jjp1) = - rlatu(1) ! .... calcul aux poles .... yprimu(1) = 0. yprimu(jjp1) = 0. un4rad2 = 0.25 * rad * rad ! -------------------------------------------------------------------- ! -------------------------------------------------------------------- ! - - ! - calcul des aires ( aire,aireu,airev, 1./aire, 1./airez ) - ! - et de fext , 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 (i,j) ! ( i + 1/4, j + 1/4 ) : aireij2 (i,j) ! ( i - 1/4, j + 1/4 ) : aireij3 (i,j) ! ( i - 1/4, j - 1/4 ) : aireij4 (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,aireij2,aireij3,aireij4 qui sont affectees au ! point (i,j) . ! On definit en outre les coefficients alpha comme etant egaux a ! (aireij / aire), c.a.d par exp. alpha1(i,j)=aireij1(i,j)/aire(i,j) ! De meme, toute aire centree en 1 point U est egale a la somme des ! 4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U . ! Idem pour airev, airez . ! On a ,pour chaque maille : dX = dY = 1 ! . V ! aireij4 . . aireij1 ! U . . P . U ! aireij3 . . aireij2 ! . V ! .................................................................... ! Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4 ! qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen ! taires 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( i,1 ) = un4rad2 * coslatm * xprp * yprm aireij3( 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( i,1 ) = 0. aireij4( i,1 ) = 0. cuij1 ( i,1 ) = 0. cuij4 ( i,1 ) = 0. cvij1 ( i,1 ) = 0. cvij4 ( i,1 ) = 0. ENDDO END IF IF ( j == jjp1 ) THEN yprp = yprimu2(j-1) rlatp = rlatu2 (j-1) !cc yprp = fyprim( REAL(j) - 0.25 ) !cc rlatp = fy ( REAL(j) - 0.25 ) coslatp = COS( rlatp ) radclatp = 0.5* rad * coslatp DO i = 1,iim xprp = xprimp025( i ) xprm = xprimm025( i ) aireij1( i,jjp1 ) = un4rad2 * coslatp * xprp * yprp aireij4( 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( i,jjp1 ) = 0. aireij3( i,jjp1 ) = 0. cvij2 ( i,jjp1 ) = 0. cvij3 ( i,jjp1 ) = 0. cuij2 ( i,jjp1 ) = 0. cuij3 ( i,jjp1 ) = 0. ENDDO END IF ! IF ( j > 1 .AND. j < jjp1 ) THEN rlatp = rlatu2 ( j-1 ) yprp = yprimu2( j-1 ) rlatm = rlatu1 ( j ) yprm = yprimu1( j ) !c rlatp = fy ( REAL(j) - 0.25 ) !c yprp = fyprim( REAL(j) - 0.25 ) !c rlatm = fy ( REAL(j) + 0.25 ) !c yprm = fyprim( REAL(j) + 0.25 ) coslatm = COS( rlatm ) coslatp = COS( rlatp ) radclatp = 0.5* rad * coslatp radclatm = 0.5* rad * coslatm ai14 = un4rad2 * coslatp * yprp ai23 = un4rad2 * coslatm * yprm DO i = 1,iim xprp = xprimp025( i ) xprm = xprimm025( i ) aireij1 ( i,j ) = ai14 * xprp aireij2 ( i,j ) = ai23 * xprp aireij3 ( i,j ) = ai23 * xprm aireij4 ( i,j ) = ai14 * xprm cuij1 ( i,j ) = radclatp * xprp cuij2 ( i,j ) = radclatm * xprp cuij3 ( i,j ) = radclatm * xprm cuij4 ( i,j ) = radclatp * xprm cvij1 ( i,j ) = 0.5* rad * yprp cvij2 ( i,j ) = 0.5* rad * yprm cvij3 ( i,j ) = cvij2(i,j) cvij4 ( i,j ) = cvij1(i,j) END DO END IF ! ........ periodicite ............ cvij1 (iip1,j) = cvij1 (1,j) cvij2 (iip1,j) = cvij2 (1,j) cvij3 (iip1,j) = cvij3 (1,j) cvij4 (iip1,j) = cvij4 (1,j) cuij1 (iip1,j) = cuij1 (1,j) cuij2 (iip1,j) = cuij2 (1,j) cuij3 (iip1,j) = cuij3 (1,j) cuij4 (iip1,j) = cuij4 (1,j) aireij1 (iip1,j) = aireij1 (1,j ) aireij2 (iip1,j) = aireij2 (1,j ) aireij3 (iip1,j) = aireij3 (1,j ) aireij4 (iip1,j) = aireij4 (1,j ) END DO ! .............................................................. DO j = 1, jjp1 DO i = 1, iim aire ( i,j ) = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) + & aireij4(i,j) alpha1 ( i,j ) = aireij1(i,j) / aire(i,j) alpha2 ( i,j ) = aireij2(i,j) / aire(i,j) alpha3 ( i,j ) = aireij3(i,j) / aire(i,j) alpha4 ( i,j ) = aireij4(i,j) / aire(i,j) alpha1p2( i,j ) = alpha1 (i,j) + alpha2 (i,j) alpha1p4( i,j ) = alpha1 (i,j) + alpha4 (i,j) alpha2p3( i,j ) = alpha2 (i,j) + alpha3 (i,j) alpha3p4( i,j ) = alpha3 (i,j) + alpha4 (i,j) END DO aire (iip1,j) = aire (1,j) alpha1 (iip1,j) = alpha1 (1,j) alpha2 (iip1,j) = alpha2 (1,j) alpha3 (iip1,j) = alpha3 (1,j) alpha4 (iip1,j) = alpha4 (1,j) alpha1p2(iip1,j) = alpha1p2(1,j) alpha1p4(iip1,j) = alpha1p4(1,j) alpha2p3(iip1,j) = alpha2p3(1,j) alpha3p4(iip1,j) = alpha3p4(1,j) END DO ! DO j = 1,jjp1 DO i = 1,iim aireu (i,j)= aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) + & aireij3(i+1,j) unsaire ( i,j)= 1./ aire(i,j) unsair_gam1( i,j)= unsaire(i,j)** ( - gamdi_gdiv ) unsair_gam2( i,j)= unsaire(i,j)** ( - gamdi_h ) airesurg ( i,j)= aire(i,j)/ g END DO aireu (iip1,j) = aireu (1,j) unsaire (iip1,j) = unsaire(1,j) unsair_gam1(iip1,j) = unsair_gam1(1,j) unsair_gam2(iip1,j) = unsair_gam2(1,j) airesurg (iip1,j) = airesurg(1,j) END DO DO j = 1,jjm DO i=1,iim airev (i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) + & aireij4(i,j+1) ENDDO DO i=1,iim airez = aireij2(i,j)+aireij1(i,j+1)+aireij3(i+1,j) + & aireij4(i+1,j+1) unsairez(i,j) = 1./ airez unsairz_gam(i,j)= unsairez(i,j)** ( - gamdi_grot ) fext (i,j) = airez * SIN(rlatv(j))* 2.* omeg ENDDO airev (iip1,j) = airev(1,j) unsairez (iip1,j) = unsairez(1,j) fext (iip1,j) = fext(1,j) unsairz_gam(iip1,j) = unsairz_gam(1,j) END DO ! ..... Calcul des elongations cu,cv, cvu ......... DO j = 1, jjm DO i = 1, iim cv(i,j) = 0.5 *( cvij2(i,j)+cvij3(i,j)+cvij1(i,j+1)+cvij4(i,j+1)) cvu(i,j)= 0.5 *( cvij1(i,j)+cvij4(i,j)+cvij2(i,j) +cvij3(i,j) ) cuv(i,j)= 0.5 *( cuij2(i,j)+cuij3(i,j)+cuij1(i,j+1)+cuij4(i,j+1)) unscv2(i,j) = 1./ ( cv(i,j)*cv(i,j) ) ENDDO DO i = 1, iim cuvsurcv (i,j) = airev(i,j) * unscv2(i,j) cvsurcuv (i,j) = 1./cuvsurcv(i,j) cuvscvgam1(i,j) = cuvsurcv (i,j) ** ( - gamdi_gdiv ) cuvscvgam2(i,j) = cuvsurcv (i,j) ** ( - gamdi_h ) cvscuvgam(i,j) = cvsurcuv (i,j) ** ( - gamdi_grot ) ENDDO cv (iip1,j) = cv (1,j) cvu (iip1,j) = cvu (1,j) unscv2 (iip1,j) = unscv2 (1,j) cuv (iip1,j) = cuv (1,j) cuvsurcv (iip1,j) = cuvsurcv (1,j) cvsurcuv (iip1,j) = cvsurcuv (1,j) cuvscvgam1(iip1,j) = cuvscvgam1(1,j) cuvscvgam2(iip1,j) = cuvscvgam2(1,j) cvscuvgam(iip1,j) = cvscuvgam(1,j) ENDDO DO j = 2, jjm DO i = 1, iim cu(i,j) = 0.5*(cuij1(i,j)+cuij4(i+1,j)+cuij2(i,j)+cuij3(i+1,j)) unscu2 (i,j) = 1./ ( cu(i,j) * cu(i,j) ) cvusurcu (i,j) = aireu(i,j) * unscu2(i,j) cusurcvu (i,j) = 1./ cvusurcu(i,j) cvuscugam1 (i,j) = cvusurcu(i,j) ** ( - gamdi_gdiv ) cvuscugam2 (i,j) = cvusurcu(i,j) ** ( - gamdi_h ) cuscvugam (i,j) = cusurcvu(i,j) ** ( - gamdi_grot ) ENDDO cu (iip1,j) = cu(1,j) unscu2 (iip1,j) = unscu2(1,j) cvusurcu (iip1,j) = cvusurcu(1,j) cusurcvu (iip1,j) = cusurcvu(1,j) cvuscugam1(iip1,j) = cvuscugam1(1,j) cvuscugam2(iip1,j) = cvuscugam2(1,j) cuscvugam (iip1,j) = cuscvugam(1,j) ENDDO ! .... calcul aux poles .... DO i = 1, iip1 cu ( i, 1 ) = 0. unscu2( i, 1 ) = 0. cvu ( i, 1 ) = 0. cu (i, jjp1) = 0. unscu2(i, jjp1) = 0. cvu (i, jjp1) = 0. ENDDO ! .............................................................. DO j = 1, jjm DO i= 1, iim airvscu2 (i,j) = airev(i,j)/ ( cuv(i,j) * cuv(i,j) ) aivscu2gam(i,j) = airvscu2(i,j)** ( - gamdi_grot ) ENDDO airvscu2 (iip1,j) = airvscu2(1,j) aivscu2gam(iip1,j) = aivscu2gam(1,j) ENDDO DO j=2,jjm DO i=1,iim airuscv2 (i,j) = aireu(i,j)/ ( cvu(i,j) * cvu(i,j) ) aiuscv2gam (i,j) = airuscv2(i,j)** ( - gamdi_grot ) ENDDO airuscv2 (iip1,j) = airuscv2 (1,j) aiuscv2gam(iip1,j) = aiuscv2gam(1,j) ENDDO ! calcul des aires aux poles : ! ----------------------------- apoln = SSUM(iim,aire(1,1),1) apols = SSUM(iim,aire(1,jjp1),1) unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) ) unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) ) unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) ) unsapolsga2 = 1./ ( apols ** ( - gamdi_h ) ) !----------------------------------------------------------------------- ! gtitre='Coriolis version ancienne' ! gfichier='fext1' ! CALL writestd(fext,iip1*jjm) ! changement F. Hourdin calcul conservatif pour fext ! constang contient le produit a * cos ( latitude ) * omega DO i=1,iim constang(i,1) = 0. ENDDO DO j=1,jjm-1 DO i=1,iim constang(i,j+1) = rad*omeg*cu(i,j+1)*COS(rlatu(j+1)) ENDDO ENDDO DO i=1,iim constang(i,jjp1) = 0. ENDDO ! periodicite en longitude DO j=1,jjm fext(iip1,j) = fext(1,j) ENDDO DO j=1,jjp1 constang(iip1,j) = constang(1,j) ENDDO ! fin du changement !----------------------------------------------------------------------- WRITE(6,*) ' *** Coordonnees de la grille *** ' WRITE(6,995) WRITE(6,*) ' LONGITUDES aux pts. V ( degres ) ' WRITE(6,995) DO i=1,iip1 rlonvv(i) = rlonv(i)*180./pi ENDDO WRITE(6,400) rlonvv WRITE(6,995) WRITE(6,*) ' LATITUDES aux pts. V ( degres ) ' WRITE(6,995) DO i=1,jjm rlatuu(i)=rlatv(i)*180./pi ENDDO WRITE(6,400) (rlatuu(i),i=1,jjm) DO i=1,iip1 rlonvv(i)=rlonu(i)*180./pi ENDDO WRITE(6,995) WRITE(6,*) ' LONGITUDES aux pts. U ( degres ) ' WRITE(6,995) WRITE(6,400) rlonvv WRITE(6,995) WRITE(6,*) ' LATITUDES aux pts. U ( degres ) ' WRITE(6,995) DO i=1,jjp1 rlatuu(i)=rlatu(i)*180./pi ENDDO WRITE(6,400) (rlatuu(i),i=1,jjp1) WRITE(6,995) 444 format(f10.3,f6.0) 400 FORMAT(1x,8f8.2) 990 FORMAT(//) 995 FORMAT(/) RETURN END SUBROUTINE inigeom