Changeset 2126 for trunk/LMDZ.COMMON/libf/dyn3d

Timestamp:

Apr 26, 2019, 11:18:52 AM (6 years ago)

Author:

emillour

Message:

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

Location:

trunk/LMDZ.COMMON/libf/dyn3d

Files:

• caladvtrac.F (modified) (3 diffs)
• conf_gcm.F90 (modified) (3 diffs)
• integrd.F (modified) (3 diffs)
• planetary_operations.F90 (added)

trunk/LMDZ.COMMON/libf/dyn3d/caladvtrac.F

@@ -9 +9 @@
 c
       USE infotrac, ONLY : nqtot
-      USE control_mod, ONLY : iapp_tracvl,planet_type
+      USE control_mod, ONLY : iapp_tracvl,planet_type,
+     &                        force_conserv_tracer
       USE comconst_mod, ONLY: dtvr
- 
+      USE planetary_operations, ONLY: planetary_tracer_amount_from_mass
       IMPLICIT NONE
 c
@@ -47 +48 @@
       INTEGER ij,l, iq, iapptrac
       REAL finmasse(ip1jmp1,llm), dtvrtrac
+      REAL :: totaltracer_old(nqtot),totaltracer_new(nqtot)
+      REAL :: ratio
 
+! Ehouarn : try to fix tracer conservation issues:
+      if (force_conserv_tracer) then
+        do iq=1,nqtot
+          call planetary_tracer_amount_from_mass(masse,q(:,:,iq),
+     &                                totaltracer_old(iq))
+        enddo
+      endif
 cc
 c
@@ -107 +117 @@
 c
         endif ! of if (planet_type.eq."earth")
-      ELSE
+        
+        ! Ehouarn : try to fix tracer conservation after tracer advection
+        if (force_conserv_tracer) then
+          do iq=1,nqtot
+            call planetary_tracer_amount_from_mass(masse,q(:,:,iq),
+     &                                  totaltracer_new(iq))
+            ratio=totaltracer_old(iq)/totaltracer_new(iq)
+            q(:,:,iq)=q(:,:,iq)*ratio
+          enddo
+        endif !of if (force_conserv_tracer)
+        
+      ELSE ! i.e. iapptrac.NE.iapp_tracvl
         if (planet_type.eq."earth") then
 ! Earth-specific treatment for the first 2 tracers (water)

trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F90

@@ -17 +17 @@
                          output_grads_dyn, periodav, planet_type, &
                          raz_date, resetvarc, starttime, timestart, &
-                         ecritstart
+                         ecritstart,force_conserv_tracer
   USE infotrac, ONLY : type_trac
   use assert_m, only: assert
@@ -233 +233 @@
   iapp_tracvl = iperiod
   CALL getin('iapp_tracvl',iapp_tracvl)
+
+! Enforce tracer conservation in integrd & caladvtrac ?
+  force_conserv_tracer = .false. 
+  CALL getin('force_conserv_tracer',force_conserv_tracer)
 
 !Config  Key  = iconser
@@ -949 +953 @@
  write(lunout,*)' iphysiq = ', iphysiq
  write(lunout,*)' iflag_trac = ', iflag_trac
+ write(lunout,*)' iapp_tracvl = ', iapp_tracvl
+ write(lunout,*)' force_conserv_tracer = ', force_conserv_tracer
  write(lunout,*)' clon = ', clon
  write(lunout,*)' clat = ', clat

trunk/LMDZ.COMMON/libf/dyn3d/integrd.F

@@ -7 +7 @@
      &  )
 
-      use control_mod, only : planet_type
+      use control_mod, only : planet_type,force_conserv_tracer
       USE comvert_mod, ONLY: ap,bp
       USE comconst_mod, ONLY: pi
@@ -66 +66 @@
       REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
       REAL massescr( ip1jmp1,llm )
+      REAL :: massratio(ip1jmp1,llm)
 !      REAL finvmasse(ip1jmp1,llm)
       REAL p(ip1jmp1,llmp1)
@@ -242 +243 @@
 
       endif ! of if (planet_type.eq."earth")
+
+      if (force_conserv_tracer) then
+        ! Ehouarn: try to keep total amont of tracers fixed 
+        ! by acounting for mass change in each cell
+        massratio(1:ip1jmp1,1:llm)=massescr(1:ip1jmp1,1:llm)
+     &                             /masse(1:ip1jmp1,1:llm)
+        do iq=1,nq
+        q(1:ip1jmp1,1:llm,iq)=q(1:ip1jmp1,1:llm,iq)
+     &                        *massratio(1:ip1jmp1,1:llm)
+        enddo
+      endif ! of if (force_conserv_tracer)
 c
 c