source: trunk/LMDZ.COMMON/libf/dyn3d/caladvtrac.F @ 3452

Last change on this file since 3452 was 2126, checked in by emillour, 6 years ago

Common dynamics:
Some work to enforce total tracer mass conservation in the dynamics.
Still to be further studied and validated.
For now these changes are triggered by setting a "force_conserv_tracer"
flag to ".true." in run.def (default is ".false." to not change anything
with respect to previous versions).
When force_conserv_tracer=.true. then:

  1. Rescale tracer mass in caladvtrac after tracer advection computations
  2. Recompute q ratios once atmospheric mass has been updated in integrd

These steps technically ensure total tracer mass conservation but it
might be the tracer advection scheme and/or time-stepping updating
sequence of fields that should be rethought or fixed.
EM

File size: 4.3 KB
Line 
1!
2! $Id: caladvtrac.F 1446 2010-10-22 09:27:25Z emillour $
3!
4c
5c
6            SUBROUTINE caladvtrac(q,pbaru,pbarv ,
7     *                   p ,masse, dq ,  teta,
8     *                   flxw, pk)
9c
10      USE infotrac, ONLY : nqtot
11      USE control_mod, ONLY : iapp_tracvl,planet_type,
12     &                        force_conserv_tracer
13      USE comconst_mod, ONLY: dtvr
14      USE planetary_operations, ONLY: planetary_tracer_amount_from_mass
15      IMPLICIT NONE
16c
17c     Auteurs:   F.Hourdin , P.Le Van, F.Forget, F.Codron 
18c
19c     F.Codron (10/99) : ajout humidite specifique pour eau vapeur
20c=======================================================================
21c
22c       Shema de  Van Leer
23c
24c=======================================================================
25
26
27#include "dimensions.h"
28#include "paramet.h"
29
30c   Arguments:
31c   ----------
32      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm)
33      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot)
34      real :: dq(ip1jmp1,llm,nqtot)
35      REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
36      REAL               :: flxw(ip1jmp1,llm)
37
38c  ..................................................................
39c
40c  .. dq n'est utilise et dimensionne que pour l'eau  vapeur et liqu.
41c
42c  ..................................................................
43c
44c   Local:
45c   ------
46
47      EXTERNAL  advtrac,minmaxq, qminimum
48      INTEGER ij,l, iq, iapptrac
49      REAL finmasse(ip1jmp1,llm), dtvrtrac
50      REAL :: totaltracer_old(nqtot),totaltracer_new(nqtot)
51      REAL :: ratio
52
53! Ehouarn : try to fix tracer conservation issues:
54      if (force_conserv_tracer) then
55        do iq=1,nqtot
56          call planetary_tracer_amount_from_mass(masse,q(:,:,iq),
57     &                                totaltracer_old(iq))
58        enddo
59      endif
60cc
61c
62! Earth-specific stuff for the first 2 tracers (water)
63      if (planet_type.eq."earth") then
64C initialisation
65! CRisi: il faut gérer tous les traceurs si on veut pouvoir faire des
66! isotopes
67!        dq(:,:,1:2)=q(:,:,1:2)
68        dq(:,:,1:nqtot)=q(:,:,1:nqtot)
69       
70c  test des valeurs minmax
71cc        CALL minmaxq(q(1,1,1),1.e33,-1.e33,'Eau vapeur (a) ')
72cc        CALL minmaxq(q(1,1,2),1.e33,-1.e33,'Eau liquide(a) ')
73      endif ! of if (planet_type.eq."earth")
74c   advection
75
76        CALL advtrac( pbaru,pbarv,
77     *       p,  masse,q,iapptrac, teta,
78     .       flxw, pk)
79
80c
81
82      IF( iapptrac.EQ.iapp_tracvl ) THEN
83        if (planet_type.eq."earth") then
84! Earth-specific treatment for the first 2 tracers (water)
85c
86cc          CALL minmaxq(q(1,1,1),1.e33,-1.e33,'Eau vapeur     ')
87cc          CALL minmaxq(q(1,1,2),1.e33,-1.e33,'Eau liquide    ')
88
89cc     ....  Calcul  de deltap  qu'on stocke dans finmasse   ...
90c
91          DO l = 1, llm
92           DO ij = 1, ip1jmp1
93             finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
94           ENDDO
95          ENDDO
96         
97          !write(*,*) 'caladvtrac 87'
98          CALL qminimum( q, nqtot, finmasse )
99          !write(*,*) 'caladvtrac 89'
100
101          CALL SCOPY   ( ip1jmp1*llm, masse, 1, finmasse,       1 )
102          CALL filtreg ( finmasse ,  jjp1,  llm, -2, 2, .TRUE., 1 )
103c
104c   *****  Calcul de dq pour l'eau , pour le passer a la physique ******
105c   ********************************************************************
106c
107          dtvrtrac = iapp_tracvl * dtvr
108c
109           DO iq = 1 , nqtot
110            DO l = 1 , llm
111             DO ij = 1,ip1jmp1
112             dq(ij,l,iq) = ( q(ij,l,iq) - dq(ij,l,iq) ) * finmasse(ij,l)
113     *                               /  dtvrtrac
114             ENDDO
115            ENDDO
116           ENDDO
117c
118        endif ! of if (planet_type.eq."earth")
119       
120        ! Ehouarn : try to fix tracer conservation after tracer advection
121        if (force_conserv_tracer) then
122          do iq=1,nqtot
123            call planetary_tracer_amount_from_mass(masse,q(:,:,iq),
124     &                                  totaltracer_new(iq))
125            ratio=totaltracer_old(iq)/totaltracer_new(iq)
126            q(:,:,iq)=q(:,:,iq)*ratio
127          enddo
128        endif !of if (force_conserv_tracer)
129       
130      ELSE ! i.e. iapptrac.NE.iapp_tracvl
131        if (planet_type.eq."earth") then
132! Earth-specific treatment for the first 2 tracers (water)
133          dq(:,:,1:nqtot)=0.
134        endif ! of if (planet_type.eq."earth")
135      ENDIF ! of IF( iapptrac.EQ.iapp_tracvl )
136
137      END
138
139
Note: See TracBrowser for help on using the repository browser.