source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn3d_common/tidal_forces.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 10.7 KB
Line 
1      SUBROUTINE tidal_forces (t, du, dv)
2
3      IMPLICIT NONE
4c
5c=======================================================================
6c
7c   Auteur:  B. Charnay  (10/2010)
8c   -------
9c
10c   Objet:
11c   ------
12c
13c   *****************************************************************
14c   ..... calcul du gradient horizontal du potentiel gravitationnel du aux forces de marees causees par Saturne
15c   ..... Formule tiree de Tokano 2002
16c   *****************************************************************
17c          Ces termes sont ajoutes a  d(ucov)/dt et a d(vcov)/dt  ..
18c
19c
20c    du et dv          sont des arguments de sortie pour le s-pg  ....
21c
22c=======================================================================
23c
24!-----------------------------------------------------------------------
25!   INCLUDE 'dimensions.h'
26!
27!   dimensions.h contient les dimensions du modele
28!   ndm est tel que iim=2**ndm
29!-----------------------------------------------------------------------
30
31      INTEGER iim,jjm,llm,ndm
32
33      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
34
35!-----------------------------------------------------------------------
36!
37! $Header$
38!
39!
40!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
41!                 veillez  n'utiliser que des ! pour les commentaires
42!                 et  bien positionner les & des lignes de continuation
43!                 (les placer en colonne 6 et en colonne 73)
44!
45!
46!-----------------------------------------------------------------------
47!   INCLUDE 'paramet.h'
48
49      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
50      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
51      INTEGER  ijmllm,mvar
52      INTEGER jcfil,jcfllm
53
54      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
55     &    ,jjp1=jjm+1-1/jjm)
56      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
57      PARAMETER( kftd  = iim/2 -ndm )
58      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
59      PARAMETER( ip1jmi1= ip1jm - iip1 )
60      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
61      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
62      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
63
64!-----------------------------------------------------------------------
65!
66! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
67!
68!
69! NB: keep items of different kinds in seperate common blocs to avoid
70!     "misaligned commons" issues
71!-----------------------------------------------------------------------
72! INCLUDE 'logic.h'
73
74      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
75     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
76     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
77     &  ,ok_limit,ok_etat0,hybrid                                       &
78     &  ,moyzon_mu,moyzon_ch
79
80      COMMON/logici/ iflag_phys,iflag_trac
81     
82      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
83     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
84     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
85     &  ,ok_limit,ok_etat0
86      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
87                     ! (only used if disvert_type==2)
88      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
89
90      integer iflag_phys,iflag_trac
91!$OMP THREADPRIVATE(/logicl/)
92!$OMP THREADPRIVATE(/logici/)
93!-----------------------------------------------------------------------
94!
95! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
96!
97!-----------------------------------------------------------------------
98!   INCLUDE 'comvert.h'
99
100      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
101     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
102     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
103
104      common/comverti/disvert_type, pressure_exner
105
106      real ap     ! hybrid pressure contribution at interlayers
107      real bp     ! hybrid sigma contribution at interlayer
108      real presnivs ! (reference) pressure at mid-layers
109      real dpres
110      real pa     ! reference pressure (Pa) at which hybrid coordinates
111                  ! become purely pressure
112      real preff  ! reference surface pressure (Pa)
113      real nivsigs
114      real nivsig
115      real aps    ! hybrid pressure contribution at mid-layers
116      real bps    ! hybrid sigma contribution at mid-layers
117      real scaleheight ! atmospheric (reference) scale height (km)
118      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
119                     ! preff and scaleheight
120
121      integer disvert_type ! type of vertical discretization:
122                           ! 1: Earth (default for planet_type==earth),
123                           !     automatic generation
124                           ! 2: Planets (default for planet_type!=earth),
125                           !     using 'z2sig.def' (or 'esasig.def) file
126
127      logical pressure_exner
128!     compute pressure inside layers using Exner function, else use mean
129!     of pressure values at interfaces
130
131 !-----------------------------------------------------------------------
132!
133! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
134!
135!-----------------------------------------------------------------------
136! INCLUDE comconst.h
137
138      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
139     &                 iflag_top_bound,mode_top_bound
140      COMMON/comconstr/dtvr,daysec,                                     &
141     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
142     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
143     & ,dissip_pupstart  ,tau_top_bound,                                &
144     & daylen,molmass, ihf
145      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
146
147      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
148      REAL dtvr ! dynamical time step (in s)
149      REAL daysec !length (in s) of a standard day
150      REAL pi    ! something like 3.14159....
151      REAL dtphys ! (s) time step for the physics
152      REAL dtdiss ! (s) time step for the dissipation
153      REAL rad ! (m) radius of the planet
154      REAL r ! Reduced Gas constant r=R/mu
155             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
156      REAL cpp   ! Cp
157      REAL kappa ! kappa=R/Cp
158      REAL cotot
159      REAL unsim ! = 1./iim
160      REAL g ! (m/s2) gravity
161      REAL omeg ! (rad/s) rotation rate of the planet
162! Dissipation factors, for Earth model:
163      REAL dissip_factz,dissip_zref !dissip_deltaz
164! Dissipation factors, for other planets:
165      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
166      REAL dissip_pupstart
167      INTEGER iflag_top_bound,mode_top_bound
168      REAL tau_top_bound
169      REAL daylen ! length of solar day, in 'standard' day length
170      REAL molmass ! (g/mol) molar mass of the atmosphere
171
172      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
173      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
174
175
176!-----------------------------------------------------------------------
177!
178! $Header$
179!
180!CDK comgeom
181      COMMON/comgeom/                                                   &
182     & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),             &
183     & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),                  &
184     & airev(ip1jm),unsaire(ip1jmp1),apoln,apols,                       &
185     & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),                 &
186     & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),              &
187     & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),                &
188     & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),               &
189     & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),           &
190     & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),            &
191     & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),         &
192     & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),           &
193     & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),                           &
194     & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),         &
195     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,                 &
196     & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),    &
197     & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
198
199!
200        REAL                                                            &
201     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,&
202     & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
203     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
204     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,&
205     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
206     & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,&
207     & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
208     & , xprimv
209!
210!#include "comorbit.h"
211      REAL t        ! jour de l'annee
212      REAL du( ip1jmp1,llm ),  dv( ip1jm,llm )
213
214c     variables locales
215      REAL Vo
216      PARAMETER (Vo=-4.691e-6)
217      INTEGER  l,ij,i,k
218      REAL n                ! 2pi/periode de rotation siderale (en jours)
219      REAL a0               ! angle à l'instant initial entre Titan et le perihelie
220      PARAMETER (a0=0.)     
221
222c     cos et sin de la latitude et longitude, calcules au premiers appel
223      REAL coslonv(ip1jm),sinlonv(ip1jm)
224      REAL sinlatv(ip1jm),coslatv(ip1jm)
225      REAL coslonu(ip1jmp1),sinlonu(ip1jmp1)
226      REAL sinlatu(ip1jmp1),coslatu(ip1jmp1)     
227
228      LOGICAl first     
229
230      SAVE coslonv,coslonu,sinlonu,sinlonv
231      SAVE coslatv,coslatu,sinlatu,sinlatv
232      SAVE first, n
233
234      DATA first /.true./
235
236! Calcul des sin et cos aux points consideres
237
238      IF(first) THEN
239         first=.false.
240         n=2*3.145!*(1+1/673.)
241         do i=1,iip1
242          do k=1,jjm
243            coslonv(i+(k-1)*iip1)=cos(rlonv(i))
244            sinlonv(i+(k-1)*iip1)=sin(rlonv(i))
245            coslatv(i+(k-1)*iip1)=cos(rlatv(k))
246            sinlatv(i+(k-1)*iip1)=sin(rlatv(k))
247          ENDDO
248         ENDDO
249
250
251
252         do i=1,iip1
253          do k=1,jjp1
254            coslonu(i+(k-1)*iip1)=cos(rlonu(i))
255            sinlonu(i+(k-1)*iip1)=sin(rlonu(i))
256            coslatu(i+(k-1)*iip1)=cos(rlatu(k))
257            sinlatu(i+(k-1)*iip1)=sin(rlatu(k))
258          ENDDO
259         ENDDO
260
261
262
263      ENDIF
264
265
266! Tendance du aux forces de maree
267
268      DO l = 1,llm
269
270      DO ij  = 1, ip1jmp1 
271
272       du(ij,l) = cu(ij)*Vo
273     $    *(3*sinlonu(ij)*coslonu(ij)*coslatu(ij)*cos(n*t+a0) 
274     $    -2*coslatu(ij)*(2*coslonu(ij)**2-1)*sin(n*t+a0))       
275      ENDDO
276
277      DO ij  = 1, ip1jm 
278       dv(ij,l) = cv(ij)*Vo
279     $    *(3*sinlatv(ij)*coslatv(ij)*coslonv(ij)**2*cos(n*t+a0) 
280     $    + 4*coslatv(ij)*sinlatv(ij)*sinlonv(ij)*coslonv(ij)
281     $    *sin(n*t+a0))     
282      ENDDO
283
284      ENDDO
285
286
287c
288      RETURN
289      END
Note: See TracBrowser for help on using the repository browser.