c*********************************************************************** SUBROUTINE eq_diff_rota4(fait,nu,y,dt,ad,as,bd,bs,ntour,Ulim) c routine private du module mod_evol c formation des équations à résoudre dans resout_rota c pour la rotation par la méthode des éléments finis Galerkin c formalisme de Mathis & Zahn 2004, Krot=4 c avec u fictif /= 0 et flux de moment cinétique c Auteur: P.Morel, Département Cassiopée, O.C.A. c CESAM2k c entrées c fait=0 : point courant, fait=1 : extérieur c nu: point de calcul en m**2/3 c y: solution c dt : pas temporel c sorties c as, ad: coefficients de la spline/dérivée c bs, bd: seconds membres c--------------------------------------------------------------------- USE mod_donnees, ONLY : lim_jpz, nom_rot, nrl, nrot, secon6 USE mod_kind USE mod_numerique, ONLY : bsp1dn, no_croiss USE mod_variables, ONLY : mrot, tot_conv IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(nrot,0:1) :: y REAL (kind=dp), INTENT(in), OPTIONAL :: Ulim REAL (kind=dp), INTENT(in) :: dt, nu INTEGER, INTENT(in) :: fait, ntour REAL (kind=dp), INTENT(out), DIMENSION(:,:,0:) :: ad, as REAL (kind=dp), INTENT(out), DIMENSION(:) :: bd, bs REAL (kind=dp), DIMENSION(4,2,0:1) :: dfrl REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: ys REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: bds, bss, frls REAL (kind=dp), DIMENSION(nrl) :: frl REAL (kind=dp), PARAMETER :: dx=1.d-5, unpdx=1.d0+dx REAL (kind=dp) :: dstor, stor, stor0 INTEGER :: i, id, j LOGICAL :: deriv=.FALSE. c---------------------------------------------------------------------- 2000 FORMAT(8es10.3) c le nombre d'inconnues est nrot=5 c y(1:nrot,der): dérivée d'ordre der=0, der=1 c (der=0 pour la fonction, der=1 pour dérivée première) c bs(i) : coefficient de S dans < bs . S > = 0 c as(i,j,k) : dérivée bs(i) / y(j,k) c i=1,...,nv, j=1,...,nv, k=0,1 c bd(i) : coefficient de dS dans < bd . dS > = 0 c ad(i,j,k) : dérivée bd(i) / y(j,k) c i=1,...,nv, j=1,...,nv, k=0,1 c les dérivées c dfrl(1,1:2,0:1) : d frl(13) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(2,1:2,0:1) : d frl(23) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(3,1:2,0:1) : d frl(32) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(4,1:2,0:1) : d frl(33) / (d Omega : d dOmega) : / (d U : d dU) c les coefficients CALL coeff_rota4(dt,nu,y,frl,dfrl) c mises à 0 ad=0.d0 ; as=0.d0 ; bd=0.d0 ; bs=0.d0 SELECT CASE(fait) CASE(0) c dans une zone de mélange pas de dérivée dfrl IF(lmix(nu))THEN c dans les ZC équation de diffusion, rotation solide bs(1)=frl(26)*y(1,0)-frl(30) as(1,1,0)=frl(26) bd(1)=frl(32)*y(1,1) ad(1,1,1)=frl(32) c flux de U fictif c conditions limites de JPZh IF(lim_jpz)THEN c bs(2)=frl(34)*y(1,0)-frl(30) !formulation complète c as(2,1,0)=frl(34) bs(2)=frl(26)*y(1,0)-frl(30) !formulation simplifiée as(2,1,0)=frl(26) bd(2)=frl(31)*y(1,0)*y(2,0) ad(2,1,0)=frl(31)*y(2,0) ad(2,2,0)=frl(31)*y(1,0) c conditions limites de PM ELSE bs(2)=y(2,0)-frl(33)*y(1,1) as(2,1,1)=-frl(33) as(2,2,0)=1.d0 ENDIF c Psi = 0 bs(3)=y(3,0) ; as(3,3,0)=1.d0 c Lambda = 0 bs(4)=y(4,0) ; as(4,4,0)=1.d0 c Phi bs(5)=frl(29)*y(1,0)**2+frl(27)*y(5,0)+frl(28)*y(1,0)*y(1,1) as(5,1,0)=2.d0*frl(29)*y(1,0)+frl(28)*y(1,1) as(5,1,1)=frl(28)*y(1,0) as(5,5,0)=frl(27) bd(5)=y(5,0)+frl(25)*y(5,1) ad(5,5,0)=1.0d0 ad(5,5,1)=frl(25) ELSE c dans les ZR advection + diffusion du moment cinétique bs(1)=frl(26)*y(1,0)-frl(30) as(1,1,0)=frl(26) bd(1)=frl(31)*y(1,0)*y(2,0)+frl(32)*y(1,1) c dfrl(3,1,0:1) : d frl(32) / (d Omega : d dOmega) : / (d U : d dU) ad(1,1,0)=frl(31)*y(2,0)+dfrl(3,1,0)*y(1,1) ad(1,1,1)=frl(32)+dfrl(3,1,1)*y(1,1) ad(1,2,0)=frl(31)*y(1,0)+dfrl(3,2,0)*y(1,1) ad(1,2,1)=dfrl(3,2,1)*y(1,1) c définition de U bs(2)=(frl(1)*y(1,0)**2+frl(3))*y(1,0)**2 1 +frl(16)*y(2,0)+frl(13)*y(3,0)+frl(14)*y(4,0) 2 +(frl(6)*y(1,0)**2+frl(9))*y(5,0)+frl(15) 2 +(frl(2)*y(1,0)**2+frl(8)*y(5,0)+frl(4))*y(1,0)*y(1,1) 3 +frl(7)*y(1,0)*y(1,1)*y(5,1)+(frl(5)*y(1,0)**2+frl(10))*y(5,1) 4 +frl(12)*y(3,1) c dfrl(1,1:2,0:1) : d frl(13) / (d Omega : d dOmega) : / (d U : d dU) as(2,1,0)=2.d0*(2.d0*frl(1)*y(1,0)**2+frl(3))*y(1,0) 1 +2.d0*frl(6)*y(1,0)*y(5,0)+(3.d0*frl(2)*y(1,0)**2 2 +frl(8)*y(5,0)+frl(4))*y(1,1)+frl(7)*y(1,1)*y(5,1) 3 +2.d0*frl(5)*y(1,0)*y(5,1)+y(3,0)*dfrl(1,1,0) as(2,1,1)=(frl(2)*y(1,0)**2+frl(8)*y(5,0)+frl(4))*y(1,0) 1 +frl(7)*y(1,0)*y(5,1)+y(3,0)*dfrl(1,1,1) as(2,2,0)=frl(16)+y(3,0)*dfrl(1,2,0) as(2,2,1)=y(3,0)*dfrl(1,2,1) as(2,3,0)=frl(13) as(2,3,1)=frl(12) as(2,4,0)=frl(14) as(2,5,0)=frl(6)*y(1,0)**2+frl(9)+frl(8)*y(1,0)*y(1,1) as(2,5,1)=frl(7)*y(1,0)*y(1,1)+frl(5)*y(1,0)**2+frl(10) bd(2)=frl(11)*y(4,0)+frl(17)*y(3,1)+frl(18)*y(3,0) ad(2,3,0)=frl(18) ad(2,3,1)=frl(17) ad(2,4,0)=frl(11) c évolution de Psi bs(3)=frl(19)*y(1,0)*y(1,1)+frl(20)*y(4,0)+frl(21)*y(3,0) as(3,1,0)=frl(19)*y(1,1) as(3,1,1)=frl(19)*y(1,0) as(3,3,0)=frl(21) as(3,4,0)=frl(20) c évolution de Lambda bs(4)=frl(22)*y(2,0)+frl(23)*y(4,0)+frl(24) c dfrl(2,1:2,0:1) : d frl(23) / (d Omega : d dOmega) : / (d U : d dU) as(4,1:2,:)=dfrl(2,1:2,:)*y(4,0) as(4,2,0)=frl(22)+dfrl(2,2,0)*y(4,0) as(4,2,1)=dfrl(2,2,1)*y(4,0) as(4,4,0)=frl(23) c Phi bs(5)=frl(29)*y(1,0)**2+frl(27)*y(5,0)+frl(28)*y(1,0)*y(1,1) as(5,1,0)=2.d0*frl(29)*y(1,0)+frl(28)*y(1,1) as(5,1,1)=frl(28)*y(1,0) as(5,5,0)=frl(27) bd(5)=y(5,0)+frl(25)*y(5,1) ad(5,5,0)=1.0 ad(5,5,1)=frl(25) ENDIF c parties intégrées limites externes, nulles au centre CASE(1) c flux de U fictif sur la couche externe, signe - mis dans resout_rota4 c sur Uf bs(2) ne sert qu'avec lim_jpz=.TRUE. bs(2)=frl(31)*y(1,0)*y(2,0) as(2,1,0)=frl(31)*y(2,0) ; as(2,2,0)=frl(31)*y(1,0) c condition sur Phi sur la couche externe, signe - mis dans resout_rota4 bs(5)=2.d0/3.d0*frl(25)*y(5,1) as(5,5,1)=2.d0/3.d0*frl(25) CASE(2) c parties intégrées, au signe près sur la face radiative des limites ZR/ZC c pour l'équation de diffusion bs(1)=frl(31)*y(1,0)*y(2,0)+frl(32)*y(1,1) c dfrl(3,1:2,0:1) : d frl(32) / (d Omega : d dOmega) : / (d U : d dU) as(1,1,0)=frl(31)*y(2,0)+dfrl(3,1,0)*y(1,1) as(1,1,1)=frl(32)+dfrl(3,1,1)*y(1,1) as(1,2,0)=frl(31)*y(1,0)+dfrl(3,2,0)*y(1,1) as(1,2,1)=dfrl(3,2,1)*y(1,1) c pour U bs(2)=frl(17)*y(3,1)+frl(11)*y(4,0)+frl(18)*y(3,0) as(2,3,0)=frl(18) as(2,3,1)=frl(17) as(2,4,0)=frl(11) c Pour Phi bs(5)=y(5,0)+frl(25)*y(5,1) as(5,5,0)=1.d0 as(5,5,1)=frl(25) c Omega (U+Ulim) + Dv dOmega/dnu sur la face radiative d'une limite ZR/ZC cas PM CASE(3) IF(PRESENT(Ulim))THEN bs(2)=y(1,0)*(y(2,0)+Ulim)+frl(33)*y(1,1) as(2,1,0)=y(2,0)+Ulim ELSE bs(2)=y(1,0)*y(2,0)+frl(33)*y(1,1) as(2,1,0)=y(2,0) ENDIF c dfrl(4,1:2,0:1) : d frl(33) / (d Omega : d dOmega) : / (d U : d dU) as(2,1,0)=as(2,1,0)+dfrl(4,1,0)*y(1,1) as(2,1,1)=frl(33)+dfrl(4,1,1)*y(1,1) as(2,2,0)=y(1,0)+dfrl(4,2,0)*y(1,1) as(2,2,1)=dfrl(4,2,1)*y(1,1) END SELECT c-----------------vérification des dérivées------------------- c deriv=.TRUE. c deriv=(nu >= 6.E-01 .AND. nu <= 6.03E-01) !ZR c 1 .OR. (nu >= 1.E-01 .AND. nu <= 1.03E-01) !ZC c 2 .OR. (fait > 0) !limite deriv=deriv .AND. (ntour > 2) c deriv= fait > 0 c deriv= nu > 0.8d0 .AND. nu < 0.9d0 c deriv= nu > 0.141d0 .AND. nu < 0.142d0 c deriv= deriv .OR. nu > 1.155d0 .AND. nu < 1.157d0 c deriv= nu > mrot(convf(1)-1) .AND. nu < mrot(convf(1)+1) c deriv= deriv .OR. nu > mrot(convd(2)-1) .AND. nu < mrot(convd(2)+1) c deriv= nu > mrot(convf(1)-1) .AND. nu < mrot(convf(1)+5) IF(deriv)THEN IF(lmix(nu))THEN PRINT*,'convectif' WRITE(*,2000)frl(27:29) ELSE PRINT*,'radiatif' ENDIF PRINT*,'nu=',nu,', fait=',fait PRINT*,'y(:,0)' ; WRITE(*,2000)y(:,0) PRINT*,'y(:,1)' ; WRITE(*,2000)y(:,1) ALLOCATE(bss(nrot),bds(nrot),frls(nrl),ys(nrot,0:1)) frls=frl DO id=0,1 DO i=1,nrot ys=y stor0=ys(i,id) ; stor=stor0*unpdx ; IF(stor == 0.d0)stor=dx dstor=stor-stor0 ; ys(i,id)=stor bss=0.d0 ; bds=0.d0 CALL coeff_rota4(dt,nu,ys,frls) SELECT CASE(fait) CASE(0) c convectif IF(lmix(nu))THEN bss(1)=frls(26)*ys(1,0)-frls(30) bds(1)=frls(32)*ys(1,1) IF(lim_jpz)THEN c bss(2)=frls(34)*ys(1,0)-frls(30) !formulation complète bss(2)=frls(26)*ys(1,0)-frls(30) !formulation simplifiée bds(2)=frls(31)*ys(1,0)*ys(2,0) ELSE bss(2)=ys(2,0)-frls(33)*ys(1,1) ENDIF bss(3)=ys(3,0) bss(4)=ys(4,0) bss(5)=frls(29)*ys(1,0)**2+frls(27)*ys(5,0) 1 +frls(28)*ys(1,0)*ys(1,1) bds(5)=ys(5,0)+frls(25)*ys(5,1) c radiatif ELSE bss(1)=frls(26)*ys(1,0)-frls(30) bds(1)=frls(31)*ys(1,0)*ys(2,0)+frls(32)*ys(1,1) bss(2)=(frls(1)*ys(1,0)**2+frls(3))*ys(1,0)**2 1 +frls(16)*ys(2,0)+frls(13)*ys(3,0)+frls(14)*ys(4,0) 2 +(frls(6)*ys(1,0)**2+frls(9))*ys(5,0)+frls(15) 3 +(frls(2)*ys(1,0)**2+frls(8)*ys(5,0)+frls(4))*ys(1,0)*ys(1,1) 4 +frls(7)*ys(1,0)*ys(1,1)*ys(5,1)+(frls(5)*ys(1,0)**2 5 +frls(10))*ys(5,1)+frls(12)*ys(3,1) bds(2)=frls(11)*ys(4,0)+frls(17)*ys(3,1)+frls(18)*ys(3,0) bss(3)=frls(19)*ys(1,0)*ys(1,1)+frls(20)*ys(4,0)+frls(21)*ys(3,0) bss(4)=frls(22)*ys(2,0)+frls(23)*ys(4,0)+frls(24) bss(5)=frls(29)*ys(1,0)**2+frls(27)*ys(5,0) 1 +frls(28)*ys(1,0)*ys(1,1) bds(5)=ys(5,0)+frls(25)*ys(5,1) ENDIF c parties intégrées limites externes, nulles au centre CASE(1) bss(2)=frls(31)*ys(1,0)*ys(2,0) bss(5)=2.d0/3.d0*frls(25)*ys(5,1) c parties intégrées, face radiative CASE(2) bss(1)=frls(31)*ys(1,0)*ys(2,0)+frls(32)*ys(1,1) bss(2)=frls(17)*ys(3,1)+frls(11)*ys(4,0)+frls(18)*ys(3,0) bss(5)=ys(5,0)+frls(25)*ys(5,1) c Omega (U+Ulim) + Dv dOmega/dnu face radiative cas PM CASE(3) IF(PRESENT(Ulim))THEN bss(2)=ys(1,0)*(ys(2,0)+Ulim)+frls(33)*ys(1,1) ELSE bss(2)=ys(1,0)*ys(2,0)+frls(33)*ys(1,1) ENDIF END SELECT c écritures IF(id == 0)THEN PRINT*,'eq_diff_rota dérivées ',nom_rot(i) ELSE PRINT*,'eq_diff_rota dérivées d',nom_rot(i) ENDIF c WRITE(*,2000)bss(:) ; WRITE(*,2000)bs(:) c WRITE(*,2000)bds(:) ; WRITE(*,2000)bd(:) WRITE(*,2000)(as(j,i,id),(bss(j)-bs(j))/dstor,j=1,nrot) WRITE(*,2000)(ad(j,i,id),(bds(j)-bd(j))/dstor,j=1,nrot) c dérivées des frl IF( i <= 2)THEN PRINT*,'vérif des dfrl' WRITE(*,2000)dfrl(1,i,id),(frls(13)-frl(13))/dstor, 1 dfrl(2,i,id),(frls(23)-frl(23))/dstor, 2 dfrl(3,i,id),(frls(32)-frl(32))/dstor, 3 dfrl(4,i,id),(frls(33)-frl(33))/dstor ENDIF PAUSE'eq_diff_rot test dérivées' ENDDO ENDDO DEALLOCATE(bds,bss,frls,ys) ; deriv=.FALSE. ENDIF c------------------fin des tests -------------------------- RETURN END SUBROUTINE eq_diff_rota4