libf/dyn3d/dissip.f90

SUBROUTINE dissip(vcov,ucov,teta,p,dv,du,dh)

  ! From dyn3d/dissip.F,v 1.1.1.1 2004/05/19 12:53:05
  ! Avec nouveaux operateurs star :  gradiv2 , divgrad2, nxgraro2  ...
  ! Auteur:  P. Le Van                                                  
  ! Objet: dissipation horizontale                                             

  USE dimens_m, ONLY : llm
  USE paramet_m, ONLY : iip1, iip2, ip1jm, ip1jmp1, llmp1
  USE comdissnew, ONLY : lstardis, nitergdiv, nitergrot, niterh
  USE inidissip_m, ONLY : dtdiss, tetah, tetaudiv, tetaurot

  IMPLICIT NONE

  !   Arguments:                                                          
  REAL :: vcov(ip1jm,llm), ucov(ip1jmp1,llm), teta(ip1jmp1,llm)
  REAL, INTENT (IN) :: p(ip1jmp1,llmp1)
  REAL :: dv(ip1jm,llm), du(ip1jmp1,llm), dh(ip1jmp1,llm)

  !   Local:                                                              
  REAL :: gdx(ip1jmp1,llm), gdy(ip1jm,llm)
  REAL :: grx(ip1jmp1,llm), gry(ip1jm,llm)
  REAL :: te1dt(llm), te2dt(llm), te3dt(llm)
  REAL :: deltapres(ip1jmp1,llm)

  INTEGER :: l, ij

  !-----------------------------------------------------------------------

  !   initialisations:                                                    

  DO l = 1, llm
     te1dt(l) = tetaudiv(l)*dtdiss
     te2dt(l) = tetaurot(l)*dtdiss
     te3dt(l) = tetah(l)*dtdiss
  END DO
  du = 0.
  dv = 0.
  dh = 0.

  !   Calcul de la dissipation:                                           

  !   Calcul de la partie   grad  ( div ) :                               

  IF (lstardis) THEN
     CALL gradiv2(llm,ucov,vcov,nitergdiv,gdx,gdy)
  ELSE
     CALL gradiv(llm,ucov,vcov,nitergdiv,gdx,gdy)
  END IF

  DO l = 1, llm

     DO ij = 1, iip1
        gdx(ij,l) = 0.
        gdx(ij+ip1jm,l) = 0.
     END DO

     DO ij = iip2, ip1jm
        du(ij,l) = du(ij,l) - te1dt(l)*gdx(ij,l)
     END DO
     DO ij = 1, ip1jm
        dv(ij,l) = dv(ij,l) - te1dt(l)*gdy(ij,l)
     END DO
  END DO

  !   calcul de la partie   n X grad ( rot ):                             

  IF (lstardis) THEN
     CALL nxgraro2(llm,ucov,vcov,nitergrot,grx,gry)
  ELSE
     CALL nxgrarot(llm,ucov,vcov,nitergrot,grx,gry)
  END IF


  DO l = 1, llm
     DO ij = 1, iip1
        grx(ij,l) = 0.
     END DO

     DO ij = iip2, ip1jm
        du(ij,l) = du(ij,l) - te2dt(l)*grx(ij,l)
     END DO
     DO ij = 1, ip1jm
        dv(ij,l) = dv(ij,l) - te2dt(l)*gry(ij,l)
     END DO
  END DO

  !   calcul de la partie   div ( grad ):                                 

  IF (lstardis) THEN

     DO l = 1, llm
        DO ij = 1, ip1jmp1
           deltapres(ij,l) = amax1(0.,p(ij,l)-p(ij,l+1))
        END DO
     END DO

     CALL divgrad2(llm,teta,deltapres,niterh,gdx)
  ELSE
     CALL divgrad(llm,teta,niterh,gdx)
  END IF

  DO l = 1, llm
     DO ij = 1, ip1jmp1
        dh(ij,l) = dh(ij,l) - te3dt(l)*gdx(ij,l)
     END DO
  END DO

END SUBROUTINE dissip
1	guez	26	SUBROUTINE dissip(vcov,ucov,teta,p,dv,du,dh)
2	guez	3
3	guez	26	! From dyn3d/dissip.F,v 1.1.1.1 2004/05/19 12:53:05
4			! Avec nouveaux operateurs star : gradiv2 , divgrad2, nxgraro2 ...
5			! Auteur: P. Le Van
6			! Objet: dissipation horizontale
7	guez	3
8	guez	26	USE dimens_m, ONLY : llm
9			USE paramet_m, ONLY : iip1, iip2, ip1jm, ip1jmp1, llmp1
10			USE comdissnew, ONLY : lstardis, nitergdiv, nitergrot, niterh
11			USE inidissip_m, ONLY : dtdiss, tetah, tetaudiv, tetaurot
12	guez	3
13	guez	26	IMPLICIT NONE
14	guez	3
15	guez	26	! Arguments:
16			REAL :: vcov(ip1jm,llm), ucov(ip1jmp1,llm), teta(ip1jmp1,llm)
17			REAL, INTENT (IN) :: p(ip1jmp1,llmp1)
18			REAL :: dv(ip1jm,llm), du(ip1jmp1,llm), dh(ip1jmp1,llm)
19	guez	3
20	guez	26	! Local:
21			REAL :: gdx(ip1jmp1,llm), gdy(ip1jm,llm)
22			REAL :: grx(ip1jmp1,llm), gry(ip1jm,llm)
23			REAL :: te1dt(llm), te2dt(llm), te3dt(llm)
24			REAL :: deltapres(ip1jmp1,llm)
25	guez	3
26	guez	26	INTEGER :: l, ij
27	guez	3
28	guez	26	!-----------------------------------------------------------------------
29	guez	3
30	guez	26	! initialisations:
31	guez	3
32	guez	26	DO l = 1, llm
33			te1dt(l) = tetaudiv(l)*dtdiss
34			te2dt(l) = tetaurot(l)*dtdiss
35			te3dt(l) = tetah(l)*dtdiss
36			END DO
37			du = 0.
38			dv = 0.
39			dh = 0.
40	guez	3
41	guez	26	! Calcul de la dissipation:
42	guez	3
43	guez	26	! Calcul de la partie grad ( div ) :
44	guez	3
45	guez	26	IF (lstardis) THEN
46			CALL gradiv2(llm,ucov,vcov,nitergdiv,gdx,gdy)
47			ELSE
48			CALL gradiv(llm,ucov,vcov,nitergdiv,gdx,gdy)
49			END IF
50	guez	3
51	guez	26	DO l = 1, llm
52	guez	3
53	guez	26	DO ij = 1, iip1
54			gdx(ij,l) = 0.
55			gdx(ij+ip1jm,l) = 0.
56			END DO
57	guez	3
58	guez	26	DO ij = iip2, ip1jm
59			du(ij,l) = du(ij,l) - te1dt(l)*gdx(ij,l)
60			END DO
61			DO ij = 1, ip1jm
62			dv(ij,l) = dv(ij,l) - te1dt(l)*gdy(ij,l)
63			END DO
64			END DO
65	guez	3
66	guez	26	! calcul de la partie n X grad ( rot ):
67	guez	3
68	guez	26	IF (lstardis) THEN
69			CALL nxgraro2(llm,ucov,vcov,nitergrot,grx,gry)
70			ELSE
71			CALL nxgrarot(llm,ucov,vcov,nitergrot,grx,gry)
72			END IF
73	guez	3
74
75	guez	26	DO l = 1, llm
76			DO ij = 1, iip1
77			grx(ij,l) = 0.
78			END DO
79	guez	3
80	guez	26	DO ij = iip2, ip1jm
81			du(ij,l) = du(ij,l) - te2dt(l)*grx(ij,l)
82			END DO
83			DO ij = 1, ip1jm
84			dv(ij,l) = dv(ij,l) - te2dt(l)*gry(ij,l)
85			END DO
86			END DO
87	guez	3
88	guez	26	! calcul de la partie div ( grad ):
89	guez	3
90	guez	26	IF (lstardis) THEN
91	guez	3
92	guez	26	DO l = 1, llm
93			DO ij = 1, ip1jmp1
94			deltapres(ij,l) = amax1(0.,p(ij,l)-p(ij,l+1))
95			END DO
96			END DO
97	guez	3
98	guez	26	CALL divgrad2(llm,teta,deltapres,niterh,gdx)
99			ELSE
100			CALL divgrad(llm,teta,niterh,gdx)
101			END IF
102	guez	3
103	guez	26	DO l = 1, llm
104			DO ij = 1, ip1jmp1
105			dh(ij,l) = dh(ij,l) - te3dt(l)*gdx(ij,l)
106			END DO
107			END DO
108	guez	3
109	guez	26	END SUBROUTINE dissip