source: trunk/LMDZ.COMMON/libf/dyn3dpar/addfi_p.F @ 2902

Last change on this file since 2902 was 2242, checked in by jvatant, 5 years ago

+ Update the interface with YAMMS so it now correctly handles the small values of the moments,
requiring dynamics to have a threshold quite low (set to 1D-200) and a sanity check in calmufi
corresponding to this value. Thus it removes 'most' of the unphysical radius obtained in
YAMMS. There are still some, but at least there is no more problem of model stability for the moments
and the code can run.
Still, take care the day you want to calculate opacities from the radii and not the moments.
Although, note that there are some negative values, in the output files, but theses negatives are

harmless, as they are only present in output files, dynamics reseting to epsilon after.

+ Given theses pbs with the radii, also update the optics so that it computes the opacities in
a simpler way, directly for M3, through look-up tables, M3 being a good proxy.
--JVO
File size: 7.2 KB
RevLine 
[1]1!
[7]2! $Id: addfi_p.F 1446 2010-10-22 09:27:25Z emillour $
[1]3!
4      SUBROUTINE addfi_p(pdt, leapf, forward,
5     S          pucov, pvcov, pteta, pq   , pps ,
6     S          pdufi, pdvfi, pdhfi,pdqfi, pdpfi  )
[1019]7      USE parallel_lmdz
[1]8      USE infotrac, ONLY : nqtot
[7]9      USE control_mod, ONLY : planet_type
[1422]10      USE comconst_mod, ONLY: kappa
[1]11      IMPLICIT NONE
12c
13c=======================================================================
14c
15c    Addition of the physical tendencies
16c
17c    Interface :
18c    -----------
19c
20c      Input :
21c      -------
22c      pdt                    time step of integration
23c      leapf                  logical
24c      forward                logical
25c      pucov(ip1jmp1,llm)     first component of the covariant velocity
26c      pvcov(ip1ip1jm,llm)    second component of the covariant velocity
27c      pteta(ip1jmp1,llm)     potential temperature
28c      pts(ip1jmp1,llm)       surface temperature
29c      pdufi(ip1jmp1,llm)     |
30c      pdvfi(ip1jm,llm)       |   respective
[108]31c      pdhfi(ip1jmp1)         |      tendencies  (unit/s)
[1]32c      pdtsfi(ip1jmp1)        |
33c
34c      Output :
35c      --------
36c      pucov
37c      pvcov
38c      ph
39c      pts
40c
41c
42c=======================================================================
43c
44c-----------------------------------------------------------------------
45c
46c    0.  Declarations :
47c    ------------------
48c
49#include "dimensions.h"
50#include "paramet.h"
51#include "comgeom.h"
52c
53c    Arguments :
54c    -----------
55c
[1189]56      REAL,INTENT(IN) :: pdt ! time step for the integration (s)
[1]57c
[1189]58      REAL,INTENT(INOUT) :: pvcov(ip1jm,llm) ! covariant meridional wind
59      REAL,INTENT(INOUT) :: pucov(ip1jmp1,llm) ! covariant zonal wind
60      REAL,INTENT(INOUT) :: pteta(ip1jmp1,llm) ! potential temperature
61      REAL,INTENT(INOUT) :: pq(ip1jmp1,llm,nqtot) ! tracers
62      REAL,INTENT(INOUT) :: pps(ip1jmp1) ! surface pressure (Pa)
63c respective tendencies (.../s) to add
64      REAL,INTENT(IN) :: pdvfi(ip1jm,llm)
65      REAL,INTENT(IN) :: pdufi(ip1jmp1,llm)
66      REAL,INTENT(IN) :: pdqfi(ip1jmp1,llm,nqtot)
67      REAL,INTENT(IN) :: pdhfi(ip1jmp1,llm)
68      REAL,INTENT(IN) :: pdpfi(ip1jmp1)
[1]69c
[1189]70      LOGICAL,INTENT(IN) :: leapf,forward ! not used
[1]71c
72c
73c    Local variables :
74c    -----------------
75c
76      REAL xpn(iim),xps(iim),tpn,tps
77      INTEGER j,k,iq,ij
[1189]78      REAL,PARAMETER :: qtestw = 1.0e-15
79      REAL,PARAMETER :: qtestt = 1.0e-40
[2242]80      REAL,PARAMETER :: qtestt2 = 1.0D-200
[1]81
82      REAL SSUM
83      EXTERNAL SSUM
84     
85      INTEGER :: ijb,ije
86c
87c-----------------------------------------------------------------------
88     
89      ijb=ij_begin
90      ije=ij_end
91     
92c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
93      DO k = 1,llm
94         DO j = ijb,ije
95            pteta(j,k)= pteta(j,k) + pdhfi(j,k) * pdt
96         ENDDO
97      ENDDO
[1348]98c$OMP END DO
[1]99
100      if (pole_nord) then
101c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
102        DO  k    = 1, llm
103         DO  ij   = 1, iim
104           xpn(ij) = aire(   ij   ) * pteta(  ij    ,k)
105         ENDDO
106         tpn      = SSUM(iim,xpn,1)/ apoln
107
108         DO ij   = 1, iip1
109           pteta(   ij   ,k)  = tpn
110         ENDDO
111       ENDDO
[1348]112c$OMP END DO
[1]113      endif
114
115      if (pole_sud) then
116c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
117        DO  k    = 1, llm
118         DO  ij   = 1, iim
119           xps(ij) = aire(ij+ip1jm) * pteta(ij+ip1jm,k)
120         ENDDO
121         tps      = SSUM(iim,xps,1)/ apols
122
123         DO ij   = 1, iip1
124           pteta(ij+ip1jm,k)  = tps
125         ENDDO
126       ENDDO
[1348]127c$OMP END DO
[1]128      endif
129c
[1238]130!***********************
131! Correction on teta due to surface pressure changes
132c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
133      DO k = 1,llm
134        DO j = ijb,ije
135           pteta(j,k)= pteta(j,k)*(1+pdpfi(j)*pdt/pps(j))**kappa
136        ENDDO
137      ENDDO
[1348]138c$OMP END DO
[1238]139!***********************
[1]140
141      ijb=ij_begin
142      ije=ij_end
143      if (pole_nord) ijb=ij_begin+iip1
144      if (pole_sud)  ije=ij_end-iip1
145
146c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
147      DO k = 1,llm
148         DO j = ijb,ije
149            pucov(j,k)= pucov(j,k) + pdufi(j,k) * pdt
150         ENDDO
151      ENDDO
[1348]152c$OMP END DO
[1]153
154      if (pole_nord) ijb=ij_begin
155
156c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
157      DO k = 1,llm
158         DO j = ijb,ije
159            pvcov(j,k)= pvcov(j,k) + pdvfi(j,k) * pdt
160         ENDDO
161      ENDDO
[1348]162c$OMP END DO
[1]163
164c
165      if (pole_sud)  ije=ij_end
166c$OMP MASTER
167      DO j = ijb,ije
168         pps(j) = pps(j) + pdpfi(j) * pdt
169      ENDDO
170c$OMP END MASTER
171 
[7]172      if (planet_type=="earth") then
173      ! earth case, special treatment for first 2 tracers (water)
174       DO iq = 1, 2
[1]175c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
176         DO k = 1,llm
177            DO j = ijb,ije
178               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
179               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestw )
180            ENDDO
181         ENDDO
[1348]182c$OMP END DO
[7]183       ENDDO
[1]184
[7]185       DO iq = 3, nqtot
[1]186c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
187         DO k = 1,llm
188            DO j = ijb,ije
189               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
190               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt )
191            ENDDO
192         ENDDO
[1348]193c$OMP END DO
[7]194       ENDDO
[2242]195      else if (planet_type=="titan") then
196      ! Titan : needs to be able to deal with very low values of tracers
197      ! values for microphysics -> threshold at 1D-200 - JVO 20
198       DO iq = 1, nqtot
199c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
200         DO k = 1,llm
201            DO j = ijb,ije
202               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
203               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt2 )
204            ENDDO
205         ENDDO
206c$OMP END DO
207       ENDDO
[7]208      else
209      ! general case, treat all tracers equally)
210       DO iq = 1, nqtot
211c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
212         DO k = 1,llm
213            DO j = ijb,ije
214               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
215               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt )
216            ENDDO
217         ENDDO
[1348]218c$OMP END DO
[7]219       ENDDO
220      endif ! of if (planet_type=="earth")
[1]221
222c$OMP MASTER
223      if (pole_nord) then
224     
225        DO  ij   = 1, iim
226          xpn(ij) = aire(   ij   ) * pps(  ij     )
227        ENDDO
228
229        tpn      = SSUM(iim,xpn,1)/apoln
230
231        DO ij   = 1, iip1
232          pps (   ij     )  = tpn
233        ENDDO
234     
235      endif
236
237      if (pole_sud) then
238     
239        DO  ij   = 1, iim
240          xps(ij) = aire(ij+ip1jm) * pps(ij+ip1jm )
241        ENDDO
242
243        tps      = SSUM(iim,xps,1)/apols
244
245        DO ij   = 1, iip1
246          pps ( ij+ip1jm )  = tps
247        ENDDO
248     
249      endif
250c$OMP END MASTER
251
252      if (pole_nord) then
253        DO iq = 1, nqtot
254c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
255          DO  k    = 1, llm
256            DO  ij   = 1, iim
257              xpn(ij) = aire(   ij   ) * pq(  ij    ,k,iq)
258            ENDDO
259            tpn      = SSUM(iim,xpn,1)/apoln
260 
261            DO ij   = 1, iip1
262              pq (   ij   ,k,iq)  = tpn
263            ENDDO
264          ENDDO
[1348]265c$OMP END DO
[1]266        ENDDO
267      endif
268     
269      if (pole_sud) then
270        DO iq = 1, nqtot
271c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
272          DO  k    = 1, llm
273            DO  ij   = 1, iim
274              xps(ij) = aire(ij+ip1jm) * pq(ij+ip1jm,k,iq)
275            ENDDO
276            tps      = SSUM(iim,xps,1)/apols
277 
278            DO ij   = 1, iip1
279              pq (ij+ip1jm,k,iq)  = tps
280            ENDDO
281          ENDDO
[1348]282c$OMP END DO
[1]283        ENDDO
284      endif
285     
286     
287      RETURN
288      END
Note: See TracBrowser for help on using the repository browser.