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	SUBROUTINE dissip(vcov,ucov,teta,p,dv,du,dh)
2
3	! 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
8	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
13	IMPLICIT NONE
14
15	! 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
20	! 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
26	INTEGER :: l, ij
27
28	!-----------------------------------------------------------------------
29
30	! initialisations:
31
32	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
41	! Calcul de la dissipation:
42
43	! Calcul de la partie grad ( div ) :
44
45	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
51	DO l = 1, llm
52
53	DO ij = 1, iip1
54	gdx(ij,l) = 0.
55	gdx(ij+ip1jm,l) = 0.
56	END DO
57
58	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
66	! calcul de la partie n X grad ( rot ):
67
68	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
74
75	DO l = 1, llm
76	DO ij = 1, iip1
77	grx(ij,l) = 0.
78	END DO
79
80	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
88	! calcul de la partie div ( grad ):
89
90	IF (lstardis) THEN
91
92	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
98	CALL divgrad2(llm,teta,deltapres,niterh,gdx)
99	ELSE
100	CALL divgrad(llm,teta,niterh,gdx)
101	END IF
102
103	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
109	END SUBROUTINE dissip