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

Last change on this file since 1592 was 1591, checked in by slebonnois, 9 years ago

SL: Mise a jour de la haute atmosphere, du transfert radiatif (solaire=Rainer Haus; IR: reglages de juin 2016), et implementation des variations de temperature potentielle dans la basse atmosphere (variation de la masse molaire moyenne)
File size: 14.4 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: planet_type
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 (planet_type.eq."venus") 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: planet_type
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 (planet_type.eq."venus") 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: planet_type
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        if (ypk(k).gt.ypklim) then
100           ratio_mod(k) = mmm0 /                                        &
101     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k)/cpp))                    &
102     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
103           ratio_mod(k) = max(ratio_mod(k),mmm0/(mmm0+0.56))
104        endif
105      ENDDO
106!----------------------
107
108      if (planet_type.eq."venus") then
109          yteta = yt**nu_venus                                          &
110     &            - nu_venus * t0_venus**nu_venus * log(ypk/cpp)
111          yteta = yteta**(1./nu_venus)
112!----------------------
113! ATMOSPHERE PROFONDE
114          yteta = yteta*ratio_mod
115!----------------------
116
117      else
118          yteta = yt * cpp/ypk
119      endif
120
121      return
122      end subroutine t2tpot
123
124!======================================================================
125!======================================================================
126
127      SUBROUTINE tpot2t(npoints,yteta, yt, ypk)
128!======================================================================
129! Arguments:
130!
131! yteta--------input-R- Temperature potentielle
132! yt   -------output-R- Temperature
133! ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
134!
135!======================================================================
136
137      USE control_mod, ONLY: planet_type
138      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
139
140      IMPLICIT NONE
141
142      integer,intent(in) :: npoints
143      REAL,intent(in) :: yteta(npoints), ypk(npoints)
144      REAL,intent(out) :: yt(npoints)
145     
146!----------------------
147! ATMOSPHERE PROFONDE
148      real ypklim,ypklim2,mmm0
149      real ratio_mod(npoints)
150      integer k
151! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
152
153      ypklim  = cpp*(  6e6/9.2e6)**0.1914
154      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
155      mmm0 = 43.44
156
157      DO k = 1, npoints
158        ratio_mod(k) = 1.
159        if (ypk(k).gt.ypklim) then
160           ratio_mod(k) = mmm0 /                                        &
161     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k)/cpp))                    &
162     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
163           ratio_mod(k) = max(ratio_mod(k),mmm0/(mmm0+0.56))
164        endif
165      ENDDO
166!----------------------
167
168      if (planet_type.eq."venus") then
169
170!----------------------
171! ATMOSPHERE PROFONDE
172!----------------------
173         yt = (yteta/ratio_mod)**nu_venus                               &
174!----------------------
175     &       + nu_venus * t0_venus**nu_venus * log(ypk/cpp)
176         yt = yt**(1./nu_venus)
177      else
178          yt = yteta * ypk/cpp
179      endif
180 
181      return
182      end subroutine tpot2t
183
184!======================================================================
185!======================================================================
186! Routines pour les calculs paralleles
187!======================================================================
188#ifdef CPP_PARA
189!======================================================================
190
191      SUBROUTINE t2tpot_p(nlon,nlev, yt, yteta, ypk)
192! Parallel version of t2tpot, for an arbitrary number of columns
193      USE control_mod, only : planet_type
194      USE parallel_lmdz, only : OMP_CHUNK
195      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
196
197      IMPLICIT none
198
199      integer,intent(in) :: nlon,nlev
200      real,intent(in) :: yt(nlon,nlev)
201      real,intent(out) :: yteta(nlon,nlev)
202      real,intent(in) :: ypk(nlon,nlev)
203! local variable:
204      integer :: l
205     
206!----------------------
207! ATMOSPHERE PROFONDE
208      real ypklim,ypklim2,mmm0
209      real ratio_mod(nlon,nlev)
210      integer k
211! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
212
213      ypklim  = cpp*(  6e6/9.2e6)**0.1914
214      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
215      mmm0 = 43.44
216
217      DO k = 1, nlon
218      DO l = 1, nlev
219        ratio_mod(k,l) = 1.
220        if (ypk(k,l).gt.ypklim) then
221           ratio_mod(k,l) = mmm0 /                                      &
222     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k,l)/cpp))                  &
223     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
224           ratio_mod(k,l) = max(ratio_mod(k,l),mmm0/(mmm0+0.56))
225        endif
226      ENDDO
227      ENDDO
228!----------------------
229
230      if (planet_type.eq."venus") then
231!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
232        do l=1,nlev
233          yteta(:,l)=yt(:,l)**nu_venus                                  &
234     &                     -nu_venus*t0_venus**nu_venus*                &
235     &                          log(ypk(:,l)/cpp)
236          yteta(:,l)=yteta(:,l)**(1./nu_venus)
237!----------------------
238! ATMOSPHERE PROFONDE
239          yteta(:,l)=yteta(:,l)*ratio_mod(:,l)
240!----------------------
241        enddo
242!$OMP END DO
243      else
244!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
245        do l=1,nlev
246          yteta(:,l)=yt(:,l)*cpp/ypk(:,l)
247        enddo
248!$OMP END DO
249      endif ! of if (planet_type.eq."venus")
250
251      end subroutine t2tpot_p
252
253!======================================================================
254!======================================================================
255
256      SUBROUTINE t2tpot_glo_p(yt, yteta, ypk)
257! Parallel version of t2tpot, over the full dynamics (scalar) grid
258! (more efficient than multiple calls to t2tpot_p() with slices of data)
259      USE parallel_lmdz, only : jj_begin,jj_end,OMP_CHUNK
260      USE control_mod, only : planet_type
261      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
262
263      IMPLICIT none
264! for iip1, jjp1 and llm
265#include "dimensions.h"
266#include "paramet.h"
267
268      real,intent(in) :: yt(iip1,jjp1,llm)
269      real,intent(out) :: yteta(iip1,jjp1,llm)
270      real,intent(in) :: ypk(iip1,jjp1,llm)
271! local variable:
272      integer :: j,l
273      integer :: jjb,jje
274     
275!----------------------
276! ATMOSPHERE PROFONDE
277      real ypklim,ypklim2,mmm0
278      real ratio_mod(iip1,jjp1,llm)
279      integer i
280! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
281
282      ypklim  = cpp*(  6e6/9.2e6)**0.1914
283      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
284      mmm0 = 43.44
285
286      DO i = 1, iip1
287      DO j = 1, jjp1
288      DO l = 1, llm
289        ratio_mod(i,j,l) = 1.
290        if (ypk(i,j,l).gt.ypklim) then
291           ratio_mod(i,j,l) = mmm0 /                                    &
292     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(i,j,l)/cpp))                &
293     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
294           ratio_mod(i,j,l) = max(ratio_mod(i,j,l),mmm0/(mmm0+0.56))
295        endif
296      ENDDO
297      ENDDO
298      ENDDO
299!----------------------
300
301      jjb=jj_begin
302      jje=jj_end
303
304      if (planet_type.eq."venus") then
305!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
306        do l=1,llm
307          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)**nu_venus                  &
308     &                     -nu_venus*t0_venus**nu_venus*                &
309     &                          log(ypk(:,jjb:jje,l)/cpp)
310          yteta(:,jjb:jje,l)=yteta(:,jjb:jje,l)**(1./nu_venus)
311!----------------------
312! ATMOSPHERE PROFONDE
313          yteta(:,jjb:jje,l)=yteta(:,jjb:jje,l)*ratio_mod(:,jjb:jje,l)
314!----------------------
315        enddo
316!$OMP END DO
317      else
318!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
319        do l=1,llm
320          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)*cpp/ypk(:,jjb:jje,l)
321        enddo
322!$OMP END DO
323      endif ! of if (planet_type.eq."venus")
324
325      end subroutine t2tpot_glo_p
326
327!======================================================================
328!======================================================================
329
330      SUBROUTINE tpot2t_p(nlon,nlev,yteta,yt,ypk)
331! Parallel version of tpot2t, for an arbitrary number of columns
332      USE control_mod, only : planet_type
333      USE parallel_lmdz, only : OMP_CHUNK
334      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
335
336      IMPLICIT none
337
338      integer,intent(in) :: nlon,nlev
339      real,intent(out) :: yt(nlon,nlev)
340      real,intent(in) :: yteta(nlon,nlev)
341      real,intent(in) :: ypk(nlon,nlev)
342
343! local variable:
344      integer :: l
345
346!----------------------
347! ATMOSPHERE PROFONDE
348      real ypklim,ypklim2,mmm0
349      real ratio_mod(nlon,nlev)
350      integer k
351! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
352
353      ypklim  = cpp*(  6e6/9.2e6)**0.1914
354      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
355      mmm0 = 43.44
356
357      DO k = 1, nlon
358      DO l = 1, nlev
359        ratio_mod(k,l) = 1.
360        if (ypk(k,l).gt.ypklim) then
361           ratio_mod(k,l) = mmm0 /                                      &
362     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(k,l)/cpp))                  &
363     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
364           ratio_mod(k,l) = max(ratio_mod(k,l),mmm0/(mmm0+0.56))
365        endif
366      ENDDO
367      ENDDO
368!----------------------
369
370      if (planet_type.eq."venus") then
371!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
372        do l=1,nlev
373!----------------------
374! ATMOSPHERE PROFONDE
375          yt(:,l)=(yteta(:,l)/ratio_mod(:,l))**nu_venus                 &
376!----------------------
377     &                  +nu_venus*t0_venus**nu_venus*                   &
378     &                       log(ypk(:,l)/cpp)
379          yt(:,l)=yt(:,l)**(1./nu_venus)
380        enddo
381!$OMP END DO
382      else
383!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
384        do l=1,nlev
385          yt(:,l)=yteta(:,l)*ypk(:,l)/cpp
386        enddo
387!$OMP END DO
388      endif ! of if (planet_type.eq."venus")
389      end subroutine tpot2t_p
390
391!======================================================================
392!======================================================================
393
394      SUBROUTINE tpot2t_glo_p(yteta,yt,ypk)
395! Parallel version of tpot2t, over the full dynamics (scalar) grid
396! (more efficient than multiple calls to tpot2t_p() with slices of data)
397      USE parallel_lmdz, only : jj_begin,jj_end,OMP_CHUNK
398      USE control_mod, only : planet_type
399      USE comconst_mod, ONLY: cpp,nu_venus,t0_venus
400
401      IMPLICIT none
402! for iip1, jjp1 and llm
403#include "dimensions.h"
404#include "paramet.h"
405
406      real,intent(out) :: yt(iip1,jjp1,llm)
407      real,intent(in) :: yteta(iip1,jjp1,llm)
408      real,intent(in) :: ypk(iip1,jjp1,llm)
409! local variable:
410      integer :: j,l
411      integer :: jjb,jje
412     
413!----------------------
414! ATMOSPHERE PROFONDE
415      real ypklim,ypklim2,mmm0
416      real ratio_mod(iip1,jjp1,llm)
417      integer i
418! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
419
420      ypklim  = cpp*(  6e6/9.2e6)**0.1914
421      ypklim2 = cpp*(8.9e6/9.2e6)**0.1914
422      mmm0 = 43.44
423
424      DO i = 1, iip1
425      DO j = 1, jjp1
426      DO l = 1, llm
427        ratio_mod(i,j,l) = 1.
428        if (ypk(i,j,l).gt.ypklim) then
429           ratio_mod(i,j,l) = mmm0 /                                    &
430     &  (mmm0+0.56*(log(ypklim/cpp)-log(ypk(i,j,l)/cpp))                &
431     &            /(log(ypklim/cpp)-log(ypklim2/cpp)))
432           ratio_mod(i,j,l) = max(ratio_mod(i,j,l),mmm0/(mmm0+0.56))
433        endif
434      ENDDO
435      ENDDO
436      ENDDO
437!----------------------
438
439      jjb=jj_begin
440      jje=jj_end
441
442      if (planet_type.eq."venus") then
443!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
444        do l=1,llm
445!----------------------
446! ATMOSPHERE PROFONDE
447          yt(:,jjb:jje,l)=                                              &
448     &          (yteta(:,jjb:jje,l)/ratio_mod(:,jjb:jje,l))**nu_venus   &
449!----------------------
450     &                  +nu_venus*t0_venus**nu_venus*                   &
451     &                       log(ypk(:,jjb:jje,l)/cpp)
452          yt(:,jjb:jje,l)=yt(:,jjb:jje,l)**(1./nu_venus)
453        enddo
454!$OMP END DO
455      else
456!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
457        do l=1,llm
458          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)*ypk(:,jjb:jje,l)/cpp
459        enddo
460!$OMP END DO
461      endif ! of if (planet_type.eq."venus")
462      end subroutine tpot2t_glo_p
463
464!======================================================================
465#endif
466!======================================================================
467! Fin routines specifiques parallele
468!======================================================================
469!======================================================================
470!
471! ATTENTION
472!
473! Si un jour on a besoin, il faudra coder les routines
474!    dt2dtpot / dtpto2dt
475!
476!======================================================================
477!======================================================================
478end module cpdet_mod
Note: See TracBrowser for help on using the repository browser.