source: trunk/LMDZ.COMMON/libf/dyn3d/planetary_operations.F90 @ 3537

Last change on this file since 3537 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: 3.2 KB
RevLine 
[2126]1!
2! $Id: $
3!
4module planetary_operations
5! module which contains functions to obtain total values over the
6! entire globe (trivial in serial mode, but not in parallel)
7
8implicit none
9
10contains
11
12  subroutine planetary_atm_mass_from_ps(ps,value)
13  implicit none
14  include "dimensions.h"
15  include "paramet.h"
16  include "comgeom.h"
17  real,intent(in) :: ps(ip1jmp1)
18  real,intent(out) :: value
19  integer :: i,j,ij
20
21  ! compute total atmospheric mass over the whole planet
22  ! North Pole
23  value=ps(1)*airesurg(1)*iim
24  do j=2,jjm
25    do i=1,iim
26      ij=i+(j-1)*iip1
27      value=value+ps(ij)*airesurg(ij)
28    enddo
29  enddo
30  ! South pole
31  value=value+ps(ip1jmp1-iim)*airesurg(ip1jmp1-iim)*iim
32 
33  end subroutine planetary_atm_mass_from_ps
34
35  subroutine planetary_atm_mass_from_mass(mass,value)
36  implicit none
37  include "dimensions.h"
38  include "paramet.h"
39  real,intent(in) :: mass(ip1jmp1,llm) ! air mass (kg)
40  real,intent(out) :: value
41  integer :: i,j,ij,l
42  real :: column(ip1jmp1) ! columns amount of air (kg)
43
44  ! compute total atmospheric mass over the whole planet
45  ! 1. build column amout of tracer (kg)
46  column(:)=0
47  do ij=1,ip1jmp1
48    do l=1,llm
49      column(ij)=column(ij)+mass(ij,l)
50    enddo
51  enddo
52 
53  ! 2. Compute total tracer mass over the planet
54  !North pole
55  value=column(1)*iim
56  do j=2,jjm
57    do i=1,iim
58      ij=i+(j-1)*iip1
59      value=value+column(ij)
60    enddo
61  enddo
62  ! South pole
63  value=value+column(ip1jmp1-iim)*iim
64 
65  end subroutine planetary_atm_mass_from_mass
66 
67  subroutine planetary_tracer_amount_from_p(p,q,amount)
68  use comconst_mod, ONLY: g
69  implicit none
70  include "dimensions.h"
71  include "paramet.h"
72  include "comgeom.h"
73  real,intent(in) :: p(ip1jmp1,llmp1) ! pressure at vertical mesh interfaces
74  real,intent(in) :: q(ip1jmp1,llm) ! 3D tracer (kg/kg_air)
75  real,intent(out) :: amount
76  integer :: i,j,ij,l
77  real :: column(ip1jmp1) ! columns amount of tracer (kg.m-2)
78 
79  ! 1. build column amout of tracer (kg.m-2)
80  column(:)=0
81  do ij=1,ip1jmp1
82    do l=1,llm
83      column(ij)=column(ij)+q(ij,l)*(p(ij,l)-p(ij,l+1))/g
84    enddo
85  enddo
86 
87  ! 2. Compute total tracer mass over the planet
88  !North pole
89  amount=column(1)*aire(1)*iim
90  do j=2,jjm
91    do i=1,iim
92      ij=i+(j-1)*iip1
93      amount=amount+column(ij)*aire(ij)
94    enddo
95  enddo
96  ! South pole
97  amount=amount+column(ip1jmp1-iim)*aire(ip1jmp1-iim)*iim
98
99  end subroutine planetary_tracer_amount_from_p
100
101  subroutine planetary_tracer_amount_from_mass(mass,q,amount)
102  implicit none
103  include "dimensions.h"
104  include "paramet.h"
105  real,intent(in) :: mass(ip1jmp1,llm) ! air mass (kg)
106  real,intent(in) :: q(ip1jmp1,llm) ! 3D tracer (kg/kg_air)
107  real,intent(out) :: amount
108  integer :: i,j,ij,l
109  real :: column(ip1jmp1) ! columns amount of tracer (kg)
110 
111  ! 1. build column amout of tracer (kg)
112  column(:)=0
113  do ij=1,ip1jmp1
114    do l=1,llm
115      column(ij)=column(ij)+q(ij,l)*mass(ij,l)
116    enddo
117  enddo
118 
119  ! 2. Compute total tracer mass over the planet
120  !North pole
121  amount=column(1)*iim
122  do j=2,jjm
123    do i=1,iim
124      ij=i+(j-1)*iip1
125      amount=amount+column(ij)
126    enddo
127  enddo
128  ! South pole
129  amount=amount+column(ip1jmp1-iim)*iim
130
131  end subroutine planetary_tracer_amount_from_mass
132
133end module planetary_operations
Note: See TracBrowser for help on using the repository browser.