dyn3d/Dissipation/inidissip.f

module inidissip_m

  use dimens_m, only: llm

  IMPLICIT NONE

  private llm

  REAL dtdiss
  integer idissip ! periode de la dissipation (en pas)
  real tetaudiv(llm),tetaurot(llm),tetah(llm)   
  real cdivu, crot, cdivh

contains

  SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
       tetagrot, tetatemp)

    ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
    ! Initialisation de la dissipation horizontale                        

    USE comconst, ONLY : dtvr
    USE comvert, ONLY : preff, presnivs
    USE conf_gcm_m, ONLY : iperiod
    USE dimens_m, ONLY : jjm, llm
    USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
    use new_unit_m, only: new_unit

    LOGICAL, intent(in):: lstardis
    INTEGER, intent(in):: nitergdiv, nitergrot, niterh
    REAL, intent(in):: tetagdiv, tetagrot, tetatemp

    ! Variables local to the procedure:
    REAL zvert(llm)
    REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm)
    REAL ullm, vllm, umin, vmin, zhmin, zhmax
    REAL zllm, z1llm
    INTEGER l, ij, idum, ii, unit
    REAL tetamin
    REAL ran1

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

    PRINT *, 'Call sequence information: inidissip'

    !   calcul des valeurs propres des operateurs par methode iterrative:   

    crot = -1.
    cdivu = -1.
    cdivh = -1.

    !   calcul de la valeur propre de divgrad:                              

    idum = 0
    DO l = 1, llm
       DO ij = 1, ip1jmp1
          deltap(ij, l) = 1.
       END DO
    END DO

    idum = -1
    zh(1) = ran1(idum) - .5
    idum = 0
    DO ij = 2, ip1jmp1
       zh(ij) = ran1(idum) - .5
    END DO

    CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)

    CALL minmax(iip1*jjp1, zh, zhmin, zhmax)

    IF (zhmin>=zhmax) THEN
       PRINT *, '  Inidissip  zh min max  ', zhmin, zhmax
       STOP 'probleme generateur alleatoire dans inidissip'
    END IF

    zllm = abs(zhmax)
    DO l = 1, 50
       IF (lstardis) THEN
          CALL divgrad2(1, zh, deltap, niterh, zh)
       ELSE
          CALL divgrad(1, zh, niterh, zh)
       END IF

       CALL minmax(iip1*jjp1, zh, zhmin, zhmax)

       zllm = abs(zhmax)
       z1llm = 1./zllm
       DO ij = 1, ip1jmp1
          zh(ij) = zh(ij)*z1llm
       END DO
    END DO

    IF (lstardis) THEN
       cdivh = 1./zllm
    ELSE
       cdivh = zllm**(-1./niterh)
    END IF

    !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)     

    PRINT *, 'calcul des valeurs propres'

    DO  ii = 1, 2

       DO ij = 1, ip1jmp1
          zu(ij) = ran1(idum) - .5
       END DO
       CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
       DO ij = 1, ip1jm
          zv(ij) = ran1(idum) - .5
       END DO
       CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)

       CALL minmax(iip1*jjp1, zu, umin, ullm)
       CALL minmax(iip1*jjm, zv, vmin, vllm)

       ullm = abs(ullm)
       vllm = abs(vllm)

       DO  l = 1, 50
          IF (ii==1) THEN
             IF (lstardis) THEN
                CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
             ELSE
                CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
             END IF
          ELSE
             IF (lstardis) THEN
                CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
             ELSE
                CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
             END IF
          END IF

          CALL minmax(iip1*jjp1, zu, umin, ullm)
          CALL minmax(iip1*jjm, zv, vmin, vllm)

          ullm = abs(ullm)
          vllm = abs(vllm)

          zllm = max(ullm, vllm)
          z1llm = 1./zllm
          DO ij = 1, ip1jmp1
             zu(ij) = zu(ij)*z1llm
          END DO
          DO ij = 1, ip1jm
             zv(ij) = zv(ij)*z1llm
          END DO
       end DO

       IF (ii==1) THEN
          IF (lstardis) THEN
             cdivu = 1./zllm
          ELSE
             cdivu = zllm**(-1./nitergdiv)
          END IF
       ELSE
          IF (lstardis) THEN
             crot = 1./zllm
          ELSE
             crot = zllm**(-1./nitergrot)
          END IF
       END IF

    END DO

    PRINT *, 'cdivu = ', cdivu
    PRINT *, 'crot = ', crot
    PRINT *, 'cdivh = ', cdivh

    ! Variation verticale du coefficient de dissipation :
    zvert = 2. - 1. / (1. + (1. - preff / presnivs)**2)

    tetaudiv = zvert / tetagdiv
    tetaurot = zvert / tetagrot
    tetah = zvert / tetatemp
    call new_unit(unit)
    open(unit, file="inidissip.csv", status="replace", action="write")
    write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
    do l = 1, llm
       write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
    end do
    close(unit)
    print *, 'Created file "inidissip.csv".'

    tetamin = min(1E6, minval(1. / tetaudiv), minval(1. / tetaurot), &
         minval(1. / tetah))
    PRINT *, 'tetamin = ', tetamin
    idissip = max(iperiod, int(tetamin / (2 * dtvr * iperiod)) * iperiod)
    PRINT *, 'idissip = ', idissip
    dtdiss = idissip * dtvr
    PRINT *, 'dtdiss = ', dtdiss

  END SUBROUTINE inidissip

end module inidissip_m
1	guez	26	module inidissip_m
2	guez	3
3	guez	26	use dimens_m, only: llm
4	guez	3
5	guez	26	IMPLICIT NONE
6	guez	3
7	guez	26	private llm
8	guez	3
9	guez	26	REAL dtdiss
10			integer idissip ! periode de la dissipation (en pas)
11			real tetaudiv(llm),tetaurot(llm),tetah(llm)
12			real cdivu, crot, cdivh
13	guez	3
14	guez	26	contains
15	guez	3
16	guez	26	SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
17			tetagrot, tetatemp)
18	guez	3
19	guez	26	! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
20			! Initialisation de la dissipation horizontale
21	guez	3
22	guez	26	USE comconst, ONLY : dtvr
23			USE comvert, ONLY : preff, presnivs
24			USE conf_gcm_m, ONLY : iperiod
25			USE dimens_m, ONLY : jjm, llm
26			USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
27			use new_unit_m, only: new_unit
28	guez	3
29	guez	26	LOGICAL, intent(in):: lstardis
30			INTEGER, intent(in):: nitergdiv, nitergrot, niterh
31			REAL, intent(in):: tetagdiv, tetagrot, tetatemp
32	guez	3
33	guez	26	! Variables local to the procedure:
34			REAL zvert(llm)
35			REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm)
36			REAL ullm, vllm, umin, vmin, zhmin, zhmax
37			REAL zllm, z1llm
38			INTEGER l, ij, idum, ii, unit
39			REAL tetamin
40			REAL ran1
41	guez	3
42	guez	26	!-----------------------------------------------------------------------
43	guez	3
44	guez	26	PRINT *, 'Call sequence information: inidissip'
45	guez	3
46	guez	26	! calcul des valeurs propres des operateurs par methode iterrative:
47	guez	3
48	guez	26	crot = -1.
49			cdivu = -1.
50			cdivh = -1.
51	guez	3
52	guez	26	! calcul de la valeur propre de divgrad:
53	guez	3
54	guez	26	idum = 0
55			DO l = 1, llm
56			DO ij = 1, ip1jmp1
57			deltap(ij, l) = 1.
58			END DO
59			END DO
60	guez	3
61	guez	26	idum = -1
62			zh(1) = ran1(idum) - .5
63			idum = 0
64			DO ij = 2, ip1jmp1
65			zh(ij) = ran1(idum) - .5
66			END DO
67	guez	3
68	guez	26	CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
69	guez	3
70	guez	26	CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
71	guez	3
72	guez	26	IF (zhmin>=zhmax) THEN
73			PRINT *, ' Inidissip zh min max ', zhmin, zhmax
74			STOP 'probleme generateur alleatoire dans inidissip'
75			END IF
76	guez	3
77	guez	26	zllm = abs(zhmax)
78			DO l = 1, 50
79			IF (lstardis) THEN
80			CALL divgrad2(1, zh, deltap, niterh, zh)
81			ELSE
82			CALL divgrad(1, zh, niterh, zh)
83			END IF
84	guez	3
85	guez	26	CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
86	guez	3
87	guez	26	zllm = abs(zhmax)
88			z1llm = 1./zllm
89			DO ij = 1, ip1jmp1
90			zh(ij) = zh(ij)*z1llm
91			END DO
92			END DO
93	guez	3
94	guez	26	IF (lstardis) THEN
95			cdivh = 1./zllm
96			ELSE
97			cdivh = zllm**(-1./niterh)
98			END IF
99	guez	3
100	guez	26	! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2)
101	guez	3
102	guez	26	PRINT *, 'calcul des valeurs propres'
103	guez	3
104	guez	26	DO ii = 1, 2
105	guez	3
106	guez	26	DO ij = 1, ip1jmp1
107			zu(ij) = ran1(idum) - .5
108			END DO
109			CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
110			DO ij = 1, ip1jm
111			zv(ij) = ran1(idum) - .5
112			END DO
113			CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
114	guez	3
115	guez	26	CALL minmax(iip1*jjp1, zu, umin, ullm)
116			CALL minmax(iip1*jjm, zv, vmin, vllm)
117	guez	3
118	guez	26	ullm = abs(ullm)
119			vllm = abs(vllm)
120	guez	3
121	guez	26	DO l = 1, 50
122			IF (ii==1) THEN
123			IF (lstardis) THEN
124			CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
125			ELSE
126			CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
127			END IF
128			ELSE
129			IF (lstardis) THEN
130			CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
131			ELSE
132			CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
133			END IF
134			END IF
135	guez	3
136	guez	26	CALL minmax(iip1*jjp1, zu, umin, ullm)
137			CALL minmax(iip1*jjm, zv, vmin, vllm)
138	guez	3
139	guez	26	ullm = abs(ullm)
140			vllm = abs(vllm)
141	guez	3
142	guez	26	zllm = max(ullm, vllm)
143			z1llm = 1./zllm
144			DO ij = 1, ip1jmp1
145			zu(ij) = zu(ij)*z1llm
146			END DO
147			DO ij = 1, ip1jm
148			zv(ij) = zv(ij)*z1llm
149			END DO
150			end DO
151	guez	3
152	guez	26	IF (ii==1) THEN
153			IF (lstardis) THEN
154			cdivu = 1./zllm
155			ELSE
156			cdivu = zllm**(-1./nitergdiv)
157			END IF
158			ELSE
159			IF (lstardis) THEN
160			crot = 1./zllm
161			ELSE
162			crot = zllm**(-1./nitergrot)
163			END IF
164			END IF
165	guez	3
166	guez	26	END DO
167	guez	3
168	guez	26	PRINT *, 'cdivu = ', cdivu
169			PRINT *, 'crot = ', crot
170			PRINT *, 'cdivh = ', cdivh
171	guez	3
172	guez	26	! Variation verticale du coefficient de dissipation :
173			zvert = 2. - 1. / (1. + (1. - preff / presnivs)**2)
174	guez	3
175	guez	26	tetaudiv = zvert / tetagdiv
176			tetaurot = zvert / tetagrot
177			tetah = zvert / tetatemp
178			call new_unit(unit)
179			open(unit, file="inidissip.csv", status="replace", action="write")
180			write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
181			do l = 1, llm
182			write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
183			end do
184			close(unit)
185			print *, 'Created file "inidissip.csv".'
186	guez	3
187	guez	26	tetamin = min(1E6, minval(1. / tetaudiv), minval(1. / tetaurot), &
188			minval(1. / tetah))
189			PRINT *, 'tetamin = ', tetamin
190			idissip = max(iperiod, int(tetamin / (2 * dtvr * iperiod)) * iperiod)
191			PRINT *, 'idissip = ', idissip
192			dtdiss = idissip * dtvr
193			PRINT *, 'dtdiss = ', dtdiss
194
195			END SUBROUTINE inidissip
196
197			end module inidissip_m