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

Last change on this file since 3020 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
RevLine 
[1]1!
[7]2! $Id: caladvtrac.F 1446 2010-10-22 09:27:25Z emillour $
[1]3!
4c
5c
6            SUBROUTINE caladvtrac(q,pbaru,pbarv ,
7     *                   p ,masse, dq ,  teta,
8     *                   flxw, pk)
9c
[66]10      USE infotrac, ONLY : nqtot
[2126]11      USE control_mod, ONLY : iapp_tracvl,planet_type,
12     &                        force_conserv_tracer
[1422]13      USE comconst_mod, ONLY: dtvr
[2126]14      USE planetary_operations, ONLY: planetary_tracer_amount_from_mass
[1]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)
[7]33      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot)
34      real :: dq(ip1jmp1,llm,nqtot)
[1]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
[2126]50      REAL :: totaltracer_old(nqtot),totaltracer_new(nqtot)
51      REAL :: ratio
[1]52
[2126]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
[1]60cc
61c
[7]62! Earth-specific stuff for the first 2 tracers (water)
63      if (planet_type.eq."earth") then
[1]64C initialisation
[1508]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)
[7]69       
[1]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) ')
[7]73      endif ! of if (planet_type.eq."earth")
[1]74c   advection
75
76        CALL advtrac( pbaru,pbarv,
77     *       p,  masse,q,iapptrac, teta,
78     .       flxw, pk)
[7]79
[1]80c
81
[7]82      IF( iapptrac.EQ.iapp_tracvl ) THEN
83        if (planet_type.eq."earth") then
84! Earth-specific treatment for the first 2 tracers (water)
[1]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         
[1508]97          !write(*,*) 'caladvtrac 87'
98          CALL qminimum( q, nqtot, finmasse )
99          !write(*,*) 'caladvtrac 89'
[1]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
[1508]109           DO iq = 1 , nqtot
[1]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
[7]118        endif ! of if (planet_type.eq."earth")
[2126]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
[7]131        if (planet_type.eq."earth") then
132! Earth-specific treatment for the first 2 tracers (water)
[1508]133          dq(:,:,1:nqtot)=0.
[7]134        endif ! of if (planet_type.eq."earth")
135      ENDIF ! of IF( iapptrac.EQ.iapp_tracvl )
[1]136
137      END
138
139
Note: See TracBrowser for help on using the repository browser.