c*********************************************************************** SUBROUTINE eq_diff_rota3(fait,nu,y,dt,ad,as,bd,bs,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 Talon & Zahn 1997, Krot=3 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 : d_conv, 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 REAL (kind=dp), INTENT(out), DIMENSION(:,:,0:) :: ad, as REAL (kind=dp), INTENT(out), DIMENSION(:) :: bd, bs REAL (kind=dp), DIMENSION(3,2,0:1) :: dfrl REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: yd 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 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 nulles dans les ZC c dfrl(1,1:2,0:1) : d frl(4) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(2,1:2,0:1) : d frl(11) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(3,1:2,0:1) : d frl(20) / (d Omega : d dOmega) : / (d U : d dU) c les coefficients CALL coeff_rota3(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ées dfrl IF(lmix(nu))THEN c dans les ZC rotation solide bs(1)=frl(1)*y(1,0)-frl(2) as(1,1,0)=frl(1) bd(1)=frl(4)*y(1,1) ad(1,1,1)=frl(4) c flux de Uf fictif c conditions limites de JPZh IF(lim_jpz)THEN bs(2)=frl(8)*y(1,0)-frl(2) !formulation complète as(2,1,0)=frl(8) c bs(2)=frl(1)*y(1,0)-frl(2) !formulation simplifiée c as(2,1,0)=frl(1) bd(2)=frl(3)*y(1,0)*y(2,0) ad(2,1,0)=frl(3)*y(2,0) ad(2,2,0)=frl(3)*y(1,0) c conditions limites de PM ELSE bs(2)=y(2,0)-frl(9)*y(1,1) as(2,1,1)=-frl(9) as(2,2,0)=1.d0 ENDIF c théta = 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 psi = 0 bs(5)=y(5,0) ; as(5,5,0)=1.d0 ELSE c dans les ZR advection + diffusion du moment cinétique bs(1)=frl(1)*y(1,0)-frl(2) as(1,1,0)=frl(1) bd(1)=frl(3)*y(1,0)*y(2,0)+frl(4)*y(1,1) c dfrl(1,1:2,0:1) : d frl(4) / (d Omega : d dOmega) : (d / d U) ad(1,1,0)=frl(3)*y(2,0)+y(1,1)*dfrl(1,1,0) !d / d Omega ad(1,1,1)=frl(4)+y(1,1)*dfrl(1,1,1) !d / d dOmega ad(1,2,0)=frl(3)*y(1,0)+y(1,1)*dfrl(1,2,0) !d / d U ad(1,2,1)=y(1,1)*dfrl(1,2,1) !d / d dU c définition de U bs(2)=(frl(5)*y(1,0)**2-frl(7))*y(1,0)**2+frl(10)+frl(6)*y(2,0) 1 +frl(11)*y(3,0)+frl(15)*y(4,0)+frl(12)*y(5,0)-frl(13)*y(5,1) c dfrl(2,1,0:1) : d frl(11) / (d Omega : d dOmega) : (d / d U) as(2,1,0)=2.d0*(2.d0*frl(5)*y(1,0)**2-frl(7))*y(1,0) 1 +y(3,0)*dfrl(2,1,0) as(2,1,1)=y(3,0)*dfrl(2,1,1) as(2,2,0)=frl(6)+y(3,0)*dfrl(2,2,0) as(2,2,1)=y(3,0)*dfrl(2,2,1) as(2,3,0)=frl(11) as(2,4,0)=frl(15) as(2,5,0)=frl(12) as(2,5,1)=-frl(13) bd(2)=frl(17)*y(3,0)+frl(18)*y(4,0)+frl(16)*y(5,0)+frl(14)*y(5,1) ad(2,3,0)=frl(17) ad(2,4,0)=frl(18) ad(2,5,0)=frl(16) ad(2,5,1)=frl(14) c définition de theta bs(3)=y(3,0)-frl(19)*y(1,0)*y(1,1) as(3,1,0)=-frl(19)*y(1,1) as(3,1,1)=-frl(19)*y(1,0) as(3,3,0)=1.d0 c évolution de Lambda bs(4)=frl(20)*y(4,0)-frl(22)-frl(21)*y(2,0) c dfrl(3,1,0:1) : d frl(20) / (d Omega : d dOmega) : (d / d U) as(4,1,:)=dfrl(3,1,:)*y(4,0) as(4,2,0)=-frl(21)+dfrl(3,2,0)*y(4,0) as(4,2,1)=dfrl(3,2,1)*y(4,0) as(4,4,0)=frl(20) c définition de psi bs(5)=y(5,0)-frl(24)*y(4,0)+frl(23)*y(3,0) as(5,3,0)=frl(23) as(5,4,0)=-frl(24) as(5,5,0)=1.d0 ENDIF c-----------------vérification des dérivées------------------- c deriv=.TRUE. 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) IF(deriv)THEN IF(lmix(nu))THEN PRINT*,'convectif' ELSE PRINT*,'radiatif' ENDIF PRINT*,'nu=',nu PRINT*,'y(:,0)' ; WRITE(*,2000)y(:,0) PRINT*,'y(:,1)' ; WRITE(*,2000)y(:,1) PRINT*,'frl' ; WRITE(*,2000)frl PRINT*,'bs' ; WRITE(*,2000)bs PRINT*,'bd' ; WRITE(*,2000)bd ALLOCATE(bss(nrot),bds(nrot),frls(nrl),yd(nrot,0:1)) yd=y ; frls=frl DO id=0,1 DO i=1,nrot stor0=yd(i,id) ; stor=stor0*unpdx ; IF(stor == 0.d0)stor=dx dstor=stor-stor0 ; yd(i,id)=stor bss=0.d0 ; bds=0.d0 CALL coeff_rota3(dt,nu,yd,frl) IF(lmix(nu))THEN bss(1)=frl(1)*yd(1,0)-frl(2) bds(1)=frl(4)*yd(1,1) bss(2)=frl(3)*yd(1,0)*yd(2,0) bss(3)=yd(3,0) bss(4)=yd(4,0) bss(5)=yd(5,0) ELSE bss(1)=frl(1)*yd(1,0)-frl(2) bds(1)=frl(3)*yd(1,0)*yd(2,0)+frl(4)*yd(1,1) bss(2)=(frl(5)*yd(1,0)**2-frl(7))*yd(1,0)**2 1 +frl(10)+frl(6)*yd(2,0) 2 +frl(11)*yd(3,0)+frl(15)*yd(4,0)+frl(12)*yd(5,0)-frl(13)*yd(5,1) bds(2)=frl(17)*yd(3,0)+frl(18)*yd(4,0)+frl(16)*yd(5,0) 1 +frl(14)*yd(5,1) bss(3)=yd(3,0)-frl(19)*yd(1,0)*yd(1,1) bss(4)=frl(20)*yd(4,0)-frl(22)-frl(21)*yd(2,0) bss(5)=yd(5,0)-frl(24)*yd(4,0)+frl(23)*yd(3,0) ENDIF 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) IF(i <=2 .AND. id <=1 .AND. .FALSE.)THEN c WRITE(*,*)frl(11),frls(11),yd(i,id),y(i,id) WRITE(*,2000)dfrl(1,i,id),(frl(4)-frls(4))/dstor, 1 dfrl(2,i,id),(frl(11)-frls(11))/dstor, 2 dfrl(3,i,id),(frl(20)-frls(20))/dstor ENDIF yd(i,id)=stor0 PAUSE'eq_diff_rot test dérivées' ENDDO ENDDO DEALLOCATE(bds,bss,frls,yd) ; deriv=.FALSE. ENDIF c------------------fin des tests -------------------------- c parties intégrées CASE(1) c flux de U fictif sur la couche externe, signe - (si lim_jpz=.TRUE.) bs(2)=frl(3)*y(1,0)*y(2,0) as(2,1,0)=frl(3)*y(2,0) ; as(2,2,0)=frl(3)*y(1,0) CASE(2) c parties intégrées, au signe près sur la face radiative d'une limite ZR/ZC c pour l'équation de diffusion bs(1)=frl(3)*y(1,0)*y(2,0)+frl(4)*y(1,1) c dfrl(1,1:2,0:1) : d frl(4) / (d Omega : d dOmega) : (d / d U) as(1,1,0)=frl(3)*y(2,0)+y(1,1)*dfrl(1,1,0) as(1,1,1)=frl(4)+y(1,1)*dfrl(1,1,1) as(1,2,0)=frl(3)*y(1,0)+y(1,1)*dfrl(1,2,0) as(1,2,1)=y(1,1)*dfrl(1,2,1) c pour U bs(2)=frl(17)*y(3,0)+frl(18)*y(4,0)+frl(16)*y(5,0)+frl(14)*y(5,1) as(2,3,0)=frl(17) ; as(2,4,0)=frl(18) ; as(2,5,0)=frl(16) as(2,5,1)=frl(14) c Omega (U+Ulim) + Dv dOmega/dnu sur la face radiative d'une limite ZR/ZC cas PM CASE(3) bs(2)=y(1,0)*(y(2,0)+Ulim)+frl(9)*y(1,1) as(2,1,0)=y(2,0)+Ulim as(2,1,1)=frl(9) as(2,2,0)=y(1,0) END SELECT RETURN END SUBROUTINE eq_diff_rota3