source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn/dissip_p.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 11.9 KB
Line 
1!
2! $Id: $
3!
4      SUBROUTINE dissip_p( vcov,ucov,teta,p, dv,du,dh )
5c
6      USE parallel_lmdz
7      USE write_field_p
8      IMPLICIT NONE
9
10
11c ..  Avec nouveaux operateurs star :  gradiv2 , divgrad2, nxgraro2  ...
12c                                 (  10/01/98  )
13
14c=======================================================================
15c
16c   Auteur:  P. Le Van
17c   -------
18c
19c   Objet:
20c   ------
21c
22c   Dissipation horizontale
23c
24c=======================================================================
25c-----------------------------------------------------------------------
26c   Declarations:
27c   -------------
28
29!-----------------------------------------------------------------------
30!   INCLUDE 'dimensions.h'
31!
32!   dimensions.h contient les dimensions du modele
33!   ndm est tel que iim=2**ndm
34!-----------------------------------------------------------------------
35
36      INTEGER iim,jjm,llm,ndm
37
38      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
39
40!-----------------------------------------------------------------------
41!
42! $Header$
43!
44!
45!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
46!                 veillez  n'utiliser que des ! pour les commentaires
47!                 et  bien positionner les & des lignes de continuation
48!                 (les placer en colonne 6 et en colonne 73)
49!
50!
51!-----------------------------------------------------------------------
52!   INCLUDE 'paramet.h'
53
54      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
55      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
56      INTEGER  ijmllm,mvar
57      INTEGER jcfil,jcfllm
58
59      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
60     &    ,jjp1=jjm+1-1/jjm)
61      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
62      PARAMETER( kftd  = iim/2 -ndm )
63      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
64      PARAMETER( ip1jmi1= ip1jm - iip1 )
65      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
66      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
67      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
68
69!-----------------------------------------------------------------------
70!
71! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
72!
73!-----------------------------------------------------------------------
74! INCLUDE comconst.h
75
76      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
77     &                 iflag_top_bound,mode_top_bound
78      COMMON/comconstr/dtvr,daysec,                                     &
79     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
80     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
81     & ,dissip_pupstart  ,tau_top_bound,                                &
82     & daylen,molmass, ihf
83      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
84
85      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
86      REAL dtvr ! dynamical time step (in s)
87      REAL daysec !length (in s) of a standard day
88      REAL pi    ! something like 3.14159....
89      REAL dtphys ! (s) time step for the physics
90      REAL dtdiss ! (s) time step for the dissipation
91      REAL rad ! (m) radius of the planet
92      REAL r ! Reduced Gas constant r=R/mu
93             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
94      REAL cpp   ! Cp
95      REAL kappa ! kappa=R/Cp
96      REAL cotot
97      REAL unsim ! = 1./iim
98      REAL g ! (m/s2) gravity
99      REAL omeg ! (rad/s) rotation rate of the planet
100! Dissipation factors, for Earth model:
101      REAL dissip_factz,dissip_zref !dissip_deltaz
102! Dissipation factors, for other planets:
103      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
104      REAL dissip_pupstart
105      INTEGER iflag_top_bound,mode_top_bound
106      REAL tau_top_bound
107      REAL daylen ! length of solar day, in 'standard' day length
108      REAL molmass ! (g/mol) molar mass of the atmosphere
109
110      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
111      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
112
113
114!-----------------------------------------------------------------------
115!
116! $Header$
117!
118!CDK comgeom
119      COMMON/comgeom/                                                   &
120     & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),             &
121     & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),                  &
122     & airev(ip1jm),unsaire(ip1jmp1),apoln,apols,                       &
123     & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),                 &
124     & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),              &
125     & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),                &
126     & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),               &
127     & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),           &
128     & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),            &
129     & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),         &
130     & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),           &
131     & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),                           &
132     & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),         &
133     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,                 &
134     & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),    &
135     & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
136
137!
138        REAL                                                            &
139     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,&
140     & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
141     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
142     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,&
143     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
144     & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,&
145     & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
146     & , xprimv
147!
148!
149! $Id: comdissnew.h 1319 2010-02-23 21:29:54Z fairhead $
150!
151!
152!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
153!                 veillez à n'utiliser que des ! pour les commentaires
154!                 et à bien positionner les & des lignes de continuation
155!                 (les placer en colonne 6 et en colonne 73)
156!
157!-----------------------------------------------------------------------
158! INCLUDE 'comdissnew.h'
159
160      COMMON/comdissnew/ lstardis,nitergdiv,nitergrot,niterh,tetagdiv,  &
161     &                   tetagrot,tetatemp,coefdis, vert_prof_dissip
162
163      LOGICAL lstardis
164      INTEGER nitergdiv, nitergrot, niterh
165
166! For the Earth model:
167      integer vert_prof_dissip ! vertical profile of horizontal dissipation
168!     Allowed values:
169!     0: rational fraction, function of pressure
170!     1: tanh of altitude
171
172      REAL     tetagdiv, tetagrot,  tetatemp, coefdis
173
174!
175! ... Les parametres de ce common comdissnew sont  lues par defrun_new
176!              sur le fichier  run.def    ....
177!
178!-----------------------------------------------------------------------
179!
180! $Header$
181!
182!  Attention : ce fichier include est compatible format fixe/format libre
183!                 veillez à n'utiliser que des ! pour les commentaires
184!                 et à bien positionner les & des lignes de continuation
185!                 (les placer en colonne 6 et en colonne 73)
186!-----------------------------------------------------------------------
187! INCLUDE comdissipn.h
188
189      REAL  tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
190!
191      COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,     &
192     &                        cdivu,      crot,         cdivh
193
194!
195!    Les parametres de ce common proviennent des calculs effectues dans
196!             Inidissip  .
197!
198!-----------------------------------------------------------------------
199
200c   Arguments:
201c   ----------
202
203      REAL,INTENT(IN) :: vcov(ip1jm,llm) ! covariant meridional wind
204      REAL,INTENT(IN) :: ucov(ip1jmp1,llm) ! covariant zonal wind
205      REAL,INTENT(IN) :: teta(ip1jmp1,llm) ! potentail temperature
206      REAL,INTENT(IN) :: p(ip1jmp1,llmp1) ! pressure
207      ! tendencies (.../s) on covariant winds and potential temperature
208      REAL,INTENT(OUT) :: dv(ip1jm,llm)
209      REAL,INTENT(OUT) :: du(ip1jmp1,llm)
210      REAL,INTENT(OUT) :: dh(ip1jmp1,llm)
211
212c   Local:
213c   ------
214
215      REAL gdx(ip1jmp1,llm),gdy(ip1jm,llm)
216      REAL grx(ip1jmp1,llm),gry(ip1jm,llm)
217      REAL te1dt(llm),te2dt(llm),te3dt(llm)
218      REAL deltapres(ip1jmp1,llm)
219
220      INTEGER l,ij
221
222      REAL  SSUM
223      integer :: ijb,ije
224c-----------------------------------------------------------------------
225c   initialisations:
226c   ----------------
227
228c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
229      DO l=1,llm
230         te1dt(l) = tetaudiv(l) * dtdiss
231         te2dt(l) = tetaurot(l) * dtdiss
232         te3dt(l) = tetah(l)    * dtdiss
233      ENDDO
234c$OMP END DO NOWAIT
235c      CALL initial0( ijp1llm, du )
236c      CALL initial0( ijmllm , dv )
237c      CALL initial0( ijp1llm, dh )
238     
239      ijb=ij_begin
240      ije=ij_end
241
242c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
243      DO l=1,llm
244        du(ijb:ije,l)=0
245        dh(ijb:ije,l)=0
246      ENDDO
247c$OMP END DO NOWAIT
248     
249      if (pole_sud) ije=ij_end-iip1
250
251c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
252      DO l=1,llm
253        dv(ijb:ije,l)=0
254      ENDDO
255c$OMP END DO NOWAIT
256     
257c-----------------------------------------------------------------------
258c   Calcul de la dissipation:
259c   -------------------------
260
261c   Calcul de la partie   grad  ( div ) :
262c   -------------------------------------
263     
264     
265     
266      IF(lstardis) THEN
267c      IF (.FALSE.) THEN
268         CALL gradiv2_p( llm,ucov,vcov,nitergdiv,gdx,gdy )
269      ELSE
270         CALL gradiv_p ( llm,ucov,vcov,nitergdiv,gdx,gdy )
271      ENDIF
272
273
274      ijb=ij_begin
275      ije=ij_end
276      if (pole_sud) ije=ij_end-iip1
277
278c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
279      DO l=1,llm
280         if (pole_nord) then
281           DO ij = 1, iip1
282              gdx(     ij ,l) = 0.
283           ENDDO
284         endif
285         
286         if (pole_sud) then
287           DO ij = 1, iip1
288              gdx(ij+ip1jm,l) = 0.
289           ENDDO
290         endif
291         
292         if (pole_nord) ijb=ij_begin+iip1
293         DO ij = ijb,ije
294            du(ij,l) = du(ij,l) - te1dt(l) *gdx(ij,l)
295         ENDDO
296
297         if (pole_nord) ijb=ij_begin
298         DO ij = ijb,ije
299            dv(ij,l) = dv(ij,l) - te1dt(l) *gdy(ij,l)
300         ENDDO
301
302       ENDDO
303c$OMP END DO NOWAIT
304c   calcul de la partie   n X grad ( rot ):
305c   ---------------------------------------
306
307      IF(lstardis) THEN
308c      IF (.FALSE.) THEN
309         CALL nxgraro2_p( llm,ucov, vcov, nitergrot,grx,gry )
310      ELSE
311         CALL nxgrarot_p( llm,ucov, vcov, nitergrot,grx,gry )
312      ENDIF
313
314
315
316      ijb=ij_begin
317      ije=ij_end
318      if (pole_sud) ije=ij_end-iip1
319
320c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
321      DO l=1,llm
322         
323         if (pole_nord) then
324           DO ij = 1, iip1
325              grx(ij,l) = 0.
326           ENDDO
327         endif
328         
329         if (pole_nord) ijb=ij_begin+iip1
330         DO ij = ijb,ije
331            du(ij,l) = du(ij,l) - te2dt(l) * grx(ij,l)
332         ENDDO
333         
334         if (pole_nord) ijb=ij_begin
335         DO ij =  ijb, ije
336            dv(ij,l) = dv(ij,l) - te2dt(l) * gry(ij,l)
337         ENDDO
338     
339      ENDDO
340c$OMP END DO NOWAIT
341
342c   calcul de la partie   div ( grad ):
343c   -----------------------------------
344
345       
346      IF(lstardis) THEN
347c      IF (.FALSE.) THEN
348   
349      ijb=ij_begin
350      ije=ij_end
351
352c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
353       DO l = 1, llm
354          DO ij = ijb, ije
355            deltapres(ij,l) = AMAX1( 0.,  p(ij,l) - p(ij,l+1) )
356          ENDDO
357       ENDDO
358c$OMP END DO NOWAIT
359         CALL divgrad2_p( llm,teta, deltapres  ,niterh, gdx )
360      ELSE
361         CALL divgrad_p ( llm,teta, niterh, gdx        )
362      ENDIF
363
364c      call write_field3d_p('gdx2',reshape(gdx,(/iip1,jmp1,llm/)))
365c      stop
366
367      ijb=ij_begin
368      ije=ij_end
369     
370c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
371      DO l = 1,llm
372         DO ij = ijb,ije
373            dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l )
374         ENDDO
375      ENDDO
376c$OMP END DO NOWAIT
377
378      RETURN
379      END
Note: See TracBrowser for help on using the repository browser.