Changeset 2126

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

Files:

• dyn3d/caladvtrac.F (modified) (3 diffs)
• dyn3d/conf_gcm.F90 (modified) (3 diffs)
• dyn3d/integrd.F (modified) (3 diffs)
• dyn3d/planetary_operations.F90 (added)
• dyn3d_common/control_mod.F90 (modified) (1 diff)
• dyn3dpar/advtrac_p.F90 (modified) (3 diffs)
• dyn3dpar/conf_gcm.F90 (modified) (2 diffs)
• dyn3dpar/integrd_p.F (modified) (3 diffs)
• dyn3dpar/planetary_operations_p.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

trunk/LMDZ.COMMON/libf/dyn3d_common/control_mod.F90

@@ -29 +29 @@
   integer,save :: ip_ebil_dyn
   logical,save :: offline
+  logical,save :: force_conserv_tracer ! enforce conservation of tracer mass
   character(len=4),save :: config_inca
   character(len=10),save :: planet_type ! planet type ('earth','mars',...)

trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90

@@ -18 +18 @@
   USE times
   USE infotrac, ONLY: nqtot, iadv
-  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
+  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type, &
+                         force_conserv_tracer
   USE comconst_mod, ONLY: dtvr
+  USE planetary_operations_p, ONLY: planetary_tracer_amount_from_mass_p
   IMPLICIT NONE
   !
@@ -77 +79 @@
   REAL,SAVE :: teta_tmp(ip1jmp1,llm)
   REAL,SAVE :: pk_tmp(ip1jmp1,llm)
+  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_p(masse,q(:,:,iq), &
+                                          totaltracer_old(iq))
+    enddo
+  endif
 
   ijb_u=ij_begin
@@ -467 +479 @@
      endif ! of if (planet_type=="earth")
 
+     ! Ehouarn : try to fix tracer conservation after tracer advection
+     if (force_conserv_tracer) then
+       ijb=ij_begin
+       ije=ij_end
+       do iq=1,nqtot
+         call planetary_tracer_amount_from_mass_p(masse,q(:,:,iq), &
+                                     totaltracer_new(iq))
+         ratio=totaltracer_old(iq)/totaltracer_new(iq)
+         q(ijb:ije,1:llm,iq)=q(ijb:ije,1:llm,iq)*ratio
+       enddo
+     endif !of if (force_conserv_tracer)
+
         !------------------------------------------------------------------
         !   on reinitialise a zero les flux de masse cumules

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

@@ -255 +255 @@
   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
@@ -987 +991 @@
       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/dyn3dpar/integrd_p.F

@@ -7 +7 @@
       USE parallel_lmdz, ONLY: ij_begin, ij_end, pole_nord, pole_sud,
      &                         omp_chunk
-      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
@@ -65 +65 @@
       REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
       REAL massescr( ip1jmp1,llm )
+      REAL :: massratio(ip1jmp1,llm)
 !      REAL finvmasse(ip1jmp1,llm)
       REAL,SAVE :: p(ip1jmp1,llmp1)
@@ -361 +362 @@
       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(ijb:ije,1:llm)=massescr(ijb:ije,1:llm)
+     &                             /masse(ijb:ije,1:llm)
+        do iq=1,nq
+        q(ijb:ije,1:llm,iq)=q(ijb:ije,1:llm,iq)
+     &                        *massratio(ijb:ije,1:llm)
+        enddo
+      endif ! of if (force_conserv_tracer)
 c
 c