source: trunk/LMDZ.COMMON/libf/dyn3d_common/cpdet_mod.F90 @ 3553

Last change on this file since 3553 was 2135, checked in by slebonnois, 6 years ago

SL, Venus: new keys for flexibility cp0/cp(T) and Held-Suarez type physics
File size: 14.5 KB
Line 
1module cpdet_mod
2
3implicit none
4
5! ADAPTATION OF GCM TO CP(T)
6!======================================================================
7! S. Lebonnois, 10/2010
8!
9! Cp must be computed using cpdet(t) to be valid
10!
11! The Exner function is still pk = RCPD*(play/pref)**RKAPPA
12! (RCPD=cpp, RKAPPA=kappa)
13!
14! One goes from T to teta (potential temperature) using t2tpot(t,teta,pk)
15! One goes from teta to T using tpot2t(teta,t,pk)
16!
17!======================================================================
18
19contains
20
21      SUBROUTINE ini_cpdet
22     
23      USE control_mod, ONLY: cpofT
24      USE comconst_mod, ONLY: nu_venus,t0_venus
25      IMPLICIT none
26!======================================================================
27! Initialization of nu_venus and t0_venus
28!======================================================================
29
30      if (cpofT) then
31          nu_venus=0.35
32          t0_venus=460.
33      else
34          nu_venus=0.
35          t0_venus=0.
36      endif
37
38      return
39      end subroutine ini_cpdet
40
41!======================================================================
42!======================================================================
43
44      FUNCTION cpdet(t)
45
46      USE control_mod, ONLY: cpofT
47      USE comconst_mod, ONLY: cpp,t0_venus,nu_venus
48      IMPLICIT none
49
50! for cpp, nu_venus and t0_venus:
51
52      real,intent(in) :: t
53      real cpdet
54
55      if (cpofT) then
56          cpdet = cpp*(t/t0_venus)**nu_venus
57      else
58          cpdet = cpp
59      endif
60
61      return
62      end function cpdet
63     
64!======================================================================
65!======================================================================
66
67      SUBROUTINE t2tpot(npoints, yt, yteta, ypk)
68!======================================================================
69! Arguments:
70!
71! yt   --------input-R- Temperature
72! yteta-------output-R- Temperature potentielle
73! ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
74!
75!======================================================================
76
77      USE control_mod, ONLY: cpofT
78      USE comconst_mod, ONLY: cpp,t0_venus,nu_venus
79
80      IMPLICIT NONE
81
82      integer,intent(in) :: npoints
83      REAL,intent(in) :: yt(npoints), ypk(npoints)
84      REAL,intent(out) :: yteta(npoints)
85     
86!----------------------
87! ATMOSPHERE PROFONDE
88      real ypklim,ypklim2,mmm0
89      real ratio_mod(npoints)
90      integer k
91! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
92
93      ypklim  = cpp*(  6e6/9.2e6)**0.1914
94      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
95      mmm0 = 43.44
96     
97      DO k = 1, npoints
98        ratio_mod(k) = 1.
99! ATM PROFONDE DESACTIVEE !
100        if (   (1 .EQ. 0)  .AND.(ypk(k).gt.ypklim)) then
101           ratio_mod(k) = mmm0 /                                        &
102     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k)/cpp))                    &
103     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
104           ratio_mod(k) = max(ratio_mod(k),mmm0/(mmm0+0.56))
105        endif
106      ENDDO
107!----------------------
108
109      if (cpofT) then
110          yteta = yt**nu_venus                                          &
111     &            - nu_venus * t0_venus**nu_venus * log(ypk/cpp)
112          yteta = yteta**(1./nu_venus)
113!----------------------
114! ATMOSPHERE PROFONDE
115          yteta = yteta*ratio_mod
116!----------------------
117
118      else
119          yteta = yt * cpp/ypk
120      endif
121
122      return
123      end subroutine t2tpot
124
125!======================================================================
126!======================================================================
127
128      SUBROUTINE tpot2t(npoints,yteta, yt, ypk)
129!======================================================================
130! Arguments:
131!
132! yteta--------input-R- Temperature potentielle
133! yt   -------output-R- Temperature
134! ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
135!
136!======================================================================
137
138      USE control_mod, ONLY: cpofT
139      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
140
141      IMPLICIT NONE
142
143      integer,intent(in) :: npoints
144      REAL,intent(in) :: yteta(npoints), ypk(npoints)
145      REAL,intent(out) :: yt(npoints)
146     
147!----------------------
148! ATMOSPHERE PROFONDE
149      real ypklim,ypklim2,mmm0
150      real ratio_mod(npoints)
151      integer k
152! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
153
154      ypklim  = cpp*(  6e6/9.2e6)**0.1914
155      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
156      mmm0 = 43.44
157
158      DO k = 1, npoints
159        ratio_mod(k) = 1.
160! ATM PROFONDE DESACTIVEE !
161        if (   (1 .EQ. 0)  .AND.(ypk(k).gt.ypklim)) then
162           ratio_mod(k) = mmm0 /                                        &
163     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k)/cpp))                    &
164     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
165           ratio_mod(k) = max(ratio_mod(k),mmm0/(mmm0+0.56))
166        endif
167      ENDDO
168!----------------------
169
170      if (cpofT) then
171
172!----------------------
173! ATMOSPHERE PROFONDE
174!----------------------
175         yt = (yteta/ratio_mod)**nu_venus                               &
176!----------------------
177     &       + nu_venus * t0_venus**nu_venus * log(ypk/cpp)
178         yt = yt**(1./nu_venus)
179      else
180          yt = yteta * ypk/cpp
181      endif
182 
183      return
184      end subroutine tpot2t
185
186!======================================================================
187!======================================================================
188! Routines pour les calculs paralleles
189!======================================================================
190#ifdef CPP_PARA
191!======================================================================
192
193      SUBROUTINE t2tpot_p(nlon,nlev, yt, yteta, ypk)
194! Parallel version of t2tpot, for an arbitrary number of columns
195      USE control_mod, only : cpofT
196      USE parallel_lmdz, only : OMP_CHUNK
197      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
198
199      IMPLICIT none
200
201      integer,intent(in) :: nlon,nlev
202      real,intent(in) :: yt(nlon,nlev)
203      real,intent(out) :: yteta(nlon,nlev)
204      real,intent(in) :: ypk(nlon,nlev)
205! local variable:
206      integer :: l
207     
208!----------------------
209! ATMOSPHERE PROFONDE
210      real ypklim,ypklim2,mmm0
211      real ratio_mod(nlon,nlev)
212      integer k
213! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
214
215      ypklim  = cpp*(  6e6/9.2e6)**0.1914
216      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
217      mmm0 = 43.44
218
219      DO k = 1, nlon
220      DO l = 1, nlev
221        ratio_mod(k,l) = 1.
222! ATM PROFONDE DESACTIVEE !
223        if (   (1 .EQ. 0)  .AND.(ypk(k,l).gt.ypklim)) then
224           ratio_mod(k,l) = mmm0 /                                      &
225     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k,l)/cpp))                  &
226     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
227           ratio_mod(k,l) = max(ratio_mod(k,l),mmm0/(mmm0+0.56))
228        endif
229      ENDDO
230      ENDDO
231!----------------------
232
233      if (cpofT) then
234!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
235        do l=1,nlev
236          yteta(:,l)=yt(:,l)**nu_venus                                  &
237     &                     -nu_venus*t0_venus**nu_venus*                &
238     &                          log(ypk(:,l)/cpp)
239          yteta(:,l)=yteta(:,l)**(1./nu_venus)
240!----------------------
241! ATMOSPHERE PROFONDE
242          yteta(:,l)=yteta(:,l)*ratio_mod(:,l)
243!----------------------
244        enddo
245!$OMP END DO
246      else
247!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
248        do l=1,nlev
249          yteta(:,l)=yt(:,l)*cpp/ypk(:,l)
250        enddo
251!$OMP END DO
252      endif ! of if (cpofT)
253
254      end subroutine t2tpot_p
255
256!======================================================================
257!======================================================================
258
259      SUBROUTINE t2tpot_glo_p(yt, yteta, ypk)
260! Parallel version of t2tpot, over the full dynamics (scalar) grid
261! (more efficient than multiple calls to t2tpot_p() with slices of data)
262      USE parallel_lmdz, only : jj_begin,jj_end,OMP_CHUNK
263      USE control_mod, only : cpofT
264      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
265
266      IMPLICIT none
267! for iip1, jjp1 and llm
268#include "dimensions.h"
269#include "paramet.h"
270
271      real,intent(in) :: yt(iip1,jjp1,llm)
272      real,intent(out) :: yteta(iip1,jjp1,llm)
273      real,intent(in) :: ypk(iip1,jjp1,llm)
274! local variable:
275      integer :: j,l
276      integer :: jjb,jje
277     
278!----------------------
279! ATMOSPHERE PROFONDE
280      real ypklim,ypklim2,mmm0
281      real ratio_mod(iip1,jjp1,llm)
282      integer i
283! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
284
285      ypklim  = cpp*(  6e6/9.2e6)**0.1914
286      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
287      mmm0 = 43.44
288
289      DO i = 1, iip1
290      DO j = 1, jjp1
291      DO l = 1, llm
292        ratio_mod(i,j,l) = 1.
293! ATM PROFONDE DESACTIVEE !
294        if (   (1 .EQ. 0)  .AND.(ypk(i,j,l).gt.ypklim)) then
295           ratio_mod(i,j,l) = mmm0 /                                    &
296     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(i,j,l)/cpp))                &
297     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
298           ratio_mod(i,j,l) = max(ratio_mod(i,j,l),mmm0/(mmm0+0.56))
299        endif
300      ENDDO
301      ENDDO
302      ENDDO
303!----------------------
304
305      jjb=jj_begin
306      jje=jj_end
307
308      if (cpofT) then
309!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
310        do l=1,llm
311          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)**nu_venus                  &
312     &                     -nu_venus*t0_venus**nu_venus*                &
313     &                          log(ypk(:,jjb:jje,l)/cpp)
314          yteta(:,jjb:jje,l)=yteta(:,jjb:jje,l)**(1./nu_venus)
315!----------------------
316! ATMOSPHERE PROFONDE
317          yteta(:,jjb:jje,l)=yteta(:,jjb:jje,l)*ratio_mod(:,jjb:jje,l)
318!----------------------
319        enddo
320!$OMP END DO
321      else
322!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
323        do l=1,llm
324          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)*cpp/ypk(:,jjb:jje,l)
325        enddo
326!$OMP END DO
327      endif ! of if (cpofT)
328
329      end subroutine t2tpot_glo_p
330
331!======================================================================
332!======================================================================
333
334      SUBROUTINE tpot2t_p(nlon,nlev,yteta,yt,ypk)
335! Parallel version of tpot2t, for an arbitrary number of columns
336      USE control_mod, only : cpofT
337      USE parallel_lmdz, only : OMP_CHUNK
338      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
339
340      IMPLICIT none
341
342      integer,intent(in) :: nlon,nlev
343      real,intent(out) :: yt(nlon,nlev)
344      real,intent(in) :: yteta(nlon,nlev)
345      real,intent(in) :: ypk(nlon,nlev)
346
347! local variable:
348      integer :: l
349
350!----------------------
351! ATMOSPHERE PROFONDE
352      real ypklim,ypklim2,mmm0
353      real ratio_mod(nlon,nlev)
354      integer k
355! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
356
357      ypklim  = cpp*(  6e6/9.2e6)**0.1914
358      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
359      mmm0 = 43.44
360
361      DO k = 1, nlon
362      DO l = 1, nlev
363        ratio_mod(k,l) = 1.
364! ATM PROFONDE DESACTIVEE !
365        if (   (1 .EQ. 0)  .AND.(ypk(k,l).gt.ypklim)) then
366           ratio_mod(k,l) = mmm0 /                                      &
367     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k,l)/cpp))                  &
368     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
369           ratio_mod(k,l) = max(ratio_mod(k,l),mmm0/(mmm0+0.56))
370        endif
371      ENDDO
372      ENDDO
373!----------------------
374
375      if (cpofT) then
376!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
377        do l=1,nlev
378!----------------------
379! ATMOSPHERE PROFONDE
380          yt(:,l)=(yteta(:,l)/ratio_mod(:,l))**nu_venus                 &
381!----------------------
382     &                  +nu_venus*t0_venus**nu_venus*                   &
383     &                       log(ypk(:,l)/cpp)
384          yt(:,l)=yt(:,l)**(1./nu_venus)
385        enddo
386!$OMP END DO
387      else
388!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
389        do l=1,nlev
390          yt(:,l)=yteta(:,l)*ypk(:,l)/cpp
391        enddo
392!$OMP END DO
393      endif ! of if (cpofT)
394      end subroutine tpot2t_p
395
396!======================================================================
397!======================================================================
398
399      SUBROUTINE tpot2t_glo_p(yteta,yt,ypk)
400! Parallel version of tpot2t, over the full dynamics (scalar) grid
401! (more efficient than multiple calls to tpot2t_p() with slices of data)
402      USE parallel_lmdz, only : jj_begin,jj_end,OMP_CHUNK
403      USE control_mod, only : cpofT
404      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
405
406      IMPLICIT none
407! for iip1, jjp1 and llm
408#include "dimensions.h"
409#include "paramet.h"
410
411      real,intent(out) :: yt(iip1,jjp1,llm)
412      real,intent(in) :: yteta(iip1,jjp1,llm)
413      real,intent(in) :: ypk(iip1,jjp1,llm)
414! local variable:
415      integer :: j,l
416      integer :: jjb,jje
417     
418!----------------------
419! ATMOSPHERE PROFONDE
420      real ypklim,ypklim2,mmm0
421      real ratio_mod(iip1,jjp1,llm)
422      integer i
423! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
424
425      ypklim  = cpp*(  6e6/9.2e6)**0.1914
426      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
427      mmm0 = 43.44
428
429      DO i = 1, iip1
430      DO j = 1, jjp1
431      DO l = 1, llm
432        ratio_mod(i,j,l) = 1.
433! ATM PROFONDE DESACTIVEE !
434        if (   (1 .EQ. 0)  .AND.(ypk(i,j,l).gt.ypklim)) then
435           ratio_mod(i,j,l) = mmm0 /                                    &
436     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(i,j,l)/cpp))                &
437     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
438           ratio_mod(i,j,l) = max(ratio_mod(i,j,l),mmm0/(mmm0+0.56))
439        endif
440      ENDDO
441      ENDDO
442      ENDDO
443!----------------------
444
445      jjb=jj_begin
446      jje=jj_end
447
448      if (cpofT) then
449!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
450        do l=1,llm
451!----------------------
452! ATMOSPHERE PROFONDE
453          yt(:,jjb:jje,l)=                                              &
454     &          (yteta(:,jjb:jje,l)/ratio_mod(:,jjb:jje,l))**nu_venus   &
455!----------------------
456     &                  +nu_venus*t0_venus**nu_venus*                   &
457     &                       log(ypk(:,jjb:jje,l)/cpp)
458          yt(:,jjb:jje,l)=yt(:,jjb:jje,l)**(1./nu_venus)
459        enddo
460!$OMP END DO
461      else
462!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
463        do l=1,llm
464          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)*ypk(:,jjb:jje,l)/cpp
465        enddo
466!$OMP END DO
467      endif ! of if (cpofT)
468      end subroutine tpot2t_glo_p
469
470!======================================================================
471#endif
472!======================================================================
473! Fin routines specifiques parallele
474!======================================================================
475!======================================================================
476!
477! ATTENTION
478!
479! Si un jour on a besoin, il faudra coder les routines
480!    dt2dtpot / dtpto2dt
481!
482!======================================================================
483!======================================================================
484end module cpdet_mod
Note: See TracBrowser for help on using the repository browser.