1 | subroutine sponge(ucov,vcov,h,pext,dt,mode) |
---|
2 | implicit none |
---|
3 | #include "dimensions.h" |
---|
4 | #include "paramet.h" |
---|
5 | #include "comdissip.h" |
---|
6 | #include "comvert.h" |
---|
7 | #include "comgeom2.h" |
---|
8 | #include "sponge.h" |
---|
9 | real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm) |
---|
10 | real h(iip1,jjp1,llm),pext(iip1,jjp1) |
---|
11 | real dt |
---|
12 | real sig_s(llm) !sigma au milieu des couches |
---|
13 | save sig_s |
---|
14 | |
---|
15 | c Local: |
---|
16 | c ------ |
---|
17 | |
---|
18 | REAL vm,um,hm,ptot(jjp1) |
---|
19 | real cst(llm) |
---|
20 | integer mode |
---|
21 | |
---|
22 | INTEGER l,i,j,l0 |
---|
23 | |
---|
24 | real ssum |
---|
25 | |
---|
26 | real echelle,zkm |
---|
27 | logical firstcall |
---|
28 | save firstcall,cst,l0 |
---|
29 | data firstcall/.true./ |
---|
30 | |
---|
31 | if (firstcall) then |
---|
32 | do l=1,llm |
---|
33 | sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2. |
---|
34 | enddo |
---|
35 | |
---|
36 | IF (mode.eq.3) then |
---|
37 | l0=llm-2 |
---|
38 | echelle=10. |
---|
39 | |
---|
40 | PRINT* |
---|
41 | print*,'sponge mode',mode |
---|
42 | print*,'hsponge',hsponge |
---|
43 | print*,'tetasponge n intervient pas' |
---|
44 | print*,'Constantes de dissipations fixees comme les anglais' |
---|
45 | print*,'Coeffs pour la couche en eponge' |
---|
46 | print*,'Z (km) tau' |
---|
47 | cst(llm)=dt/(0.5*88775) |
---|
48 | cst(llm-1)=dt/(88775) |
---|
49 | cst(l0)=dt/(2*88775) |
---|
50 | do l=l0,llm |
---|
51 | zkm=-echelle*log(sig_s(l)) |
---|
52 | print*,zkm,dt/cst(l),cst(l) |
---|
53 | enddo |
---|
54 | firstcall=.false. |
---|
55 | PRINT* |
---|
56 | ELSE |
---|
57 | l0=1 |
---|
58 | echelle=10. |
---|
59 | |
---|
60 | PRINT* |
---|
61 | print*,'sponge mode',mode |
---|
62 | print*,'hsponge tetasponge ',hsponge,tetasponge |
---|
63 | print*,'Coeffs pour la couche en eponge' |
---|
64 | print*,'Z (km) tau' |
---|
65 | do l=1,llm |
---|
66 | zkm=-echelle*log(sig_s(l)) |
---|
67 | cst(l)=.5*(1.+tanh((zkm-hsponge)/echelle)) |
---|
68 | cst(l)= max(tetasponge*1.e-15,cst(l)) |
---|
69 | print*,zkm,1./cst(l)*tetasponge,cst(l)*dt/tetasponge |
---|
70 | cst(l)=cst(l)*dt/tetasponge |
---|
71 | enddo |
---|
72 | firstcall=.false. |
---|
73 | PRINT* |
---|
74 | ENDIF |
---|
75 | endif |
---|
76 | |
---|
77 | c----------------------------------------------------------------------- |
---|
78 | c calcul de la dissipation: |
---|
79 | c ------------------------- |
---|
80 | |
---|
81 | do j=1,jjp1 |
---|
82 | ptot(j)=ssum(iim,pext(1,j),1) |
---|
83 | enddo |
---|
84 | |
---|
85 | c temperature potentielle |
---|
86 | do l=l0,llm |
---|
87 | do j=1,jjp1 |
---|
88 | hm=0. |
---|
89 | do i=1,iim |
---|
90 | hm=hm+h(i,j,l)*pext(i,j) |
---|
91 | enddo |
---|
92 | hm=hm/ptot(j) |
---|
93 | do i=1,iim |
---|
94 | h(i,j,l)=h(i,j,l)-cst(l)*(h(i,j,l)-hm) |
---|
95 | enddo |
---|
96 | h(iip1,j,l)=h(1,j,l) |
---|
97 | enddo |
---|
98 | enddo |
---|
99 | |
---|
100 | c vent zonal |
---|
101 | do l=l0,llm |
---|
102 | do j=2,jjm |
---|
103 | um=0. |
---|
104 | if(mode.ge.1) then |
---|
105 | do i=1,iim |
---|
106 | um=um+0.5*ucov(i,j,l)*(pext(i,j)+pext(i+1,j)) |
---|
107 | s /cu(i,j) |
---|
108 | enddo |
---|
109 | um=um/ptot(j) |
---|
110 | endif |
---|
111 | do i=1,iim |
---|
112 | ucov(i,j,l)=ucov(i,j,l)-cst(l)*(ucov(i,j,l)-um*cu(i,j)) |
---|
113 | enddo |
---|
114 | ucov(iip1,j,l)=ucov(1,j,l) |
---|
115 | enddo |
---|
116 | enddo |
---|
117 | |
---|
118 | c vent meridien |
---|
119 | do l=l0,llm |
---|
120 | do j=1,jjm |
---|
121 | vm=0. |
---|
122 | if(mode.ge.2) then |
---|
123 | do i=1,iim |
---|
124 | vm=vm+vcov(i,j,l)*(pext(i,j)+pext(i,j+1)) |
---|
125 | s /cv(i,j) |
---|
126 | enddo |
---|
127 | vm=vm/(ptot(j)+ptot(j+1)) |
---|
128 | endif |
---|
129 | do i=1,iim |
---|
130 | vcov(i,j,l)=vcov(i,j,l)-cst(l)*(vcov(i,j,l)-vm*cv(i,j)) |
---|
131 | enddo |
---|
132 | vcov(iip1,j,l)=vcov(1,j,l) |
---|
133 | enddo |
---|
134 | enddo |
---|
135 | |
---|
136 | RETURN |
---|
137 | end |
---|