Changeset 4853 for LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_gas.F90

Timestamp:

Mar 19, 2024, 3:34:21 PM (2 months ago)

Author:

idelkadi

Message:

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

• version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
• interface routines between lmdz and ecrad grouped in a new "lmdz" directory
• double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
• clean-up in the read.F routine: delete unitules arguments
• modification of default gas model in namelist (default: ECCKD)

File:

• LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_gas.F90 (modified) (5 diffs)

LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_gas.F90

@@ -72 +72 @@
      procedure :: assert_units => assert_units_gas
      procedure :: get        => get_gas
+     procedure :: get_scaling
      procedure :: reverse    => reverse_gas
      procedure :: out_of_physical_bounds
@@ -355 +356 @@
     real(jprb), optional, intent(in)    :: scale_factor
 
-    integer :: ig
+    integer :: jg
 
     ! Scaling factor to convert from old to new
@@ -396 +397 @@
       end if
     else
-      do ig = 1,this%ntype
-        call this%set_units(iunits, igas=this%icode(ig), scale_factor=new_sf)
+      do jg = 1,this%ntype
+        call this%set_units(iunits, igas=this%icode(jg), scale_factor=new_sf)
       end do
     end if
@@ -403 +404 @@
   end subroutine set_units_gas
 
-
+  
+  !---------------------------------------------------------------------
+  ! Return a vector indicating the scaling that one would need to
+  ! apply to each gas in order to obtain the dimension units in
+  ! "iunits" (which can be IVolumeMixingRatio or IMassMixingRatio)
+  subroutine get_scaling(this, iunits, scaling)
+    class(gas_type), intent(in)  :: this
+    integer,         intent(in)  :: iunits
+    real(jprb),      intent(out) :: scaling(NMaxGases)
+    integer :: jg
+    
+    scaling = this%scale_factor
+    do jg = 1,NMaxGases
+      if (iunits == IMassMixingRatio .and. this%iunits(jg) == IVolumeMixingRatio) then
+        scaling(jg) = scaling(jg) * GasMolarMass(jg) / AirMolarMass
+      else if (iunits == IVolumeMixingRatio .and. this%iunits(jg) == IMassMixingRatio) then
+        scaling(jg) = scaling(jg) * AirMolarMass / GasMolarMass(jg)
+      end if
+    end do
+    
+  end subroutine get_scaling
+
+  
   !---------------------------------------------------------------------
   ! Assert that gas mixing ratio units are "iunits", applying to gas
   ! with ID "igas" if present, otherwise to all gases. Otherwise the
-  ! program will exit. Otional argument scale factor specifies any
-  ! subsequent multiplication to apply; for PPMV one would use
-  ! iunits=IVolumeMixingRatio and scale_factor=1.0e6.
-  recursive subroutine assert_units_gas(this, iunits, igas, scale_factor)
+  ! program will exit, except if the optional argument "istatus" is
+  ! provided in which case it will return true if the units are
+  ! correct and false if they are not. Optional argument scale factor
+  ! specifies any subsequent multiplication to apply; for PPMV one
+  ! would use iunits=IVolumeMixingRatio and scale_factor=1.0e6.
+  recursive subroutine assert_units_gas(this, iunits, igas, scale_factor, istatus)
 
     use radiation_io,   only : nulerr, radiation_abort
 
-    class(gas_type),      intent(in) :: this
-    integer,              intent(in) :: iunits
-    integer,    optional, intent(in) :: igas
-    real(jprb), optional, intent(in) :: scale_factor
-
-    integer :: ig
+    class(gas_type),      intent(in)  :: this
+    integer,              intent(in)  :: iunits
+    integer,    optional, intent(in)  :: igas
+    real(jprb), optional, intent(in)  :: scale_factor
+    logical,    optional, intent(out) :: istatus
+
+    integer :: jg
 
     real(jprb) :: sf
@@ -429 +455 @@
     end if
 
+    if (present(istatus)) then
+      istatus = .true.
+    end if
+    
     if (present(igas)) then
       if (this%is_present(igas)) then
         if (iunits /= this%iunits(igas)) then
-          write(nulerr,'(a,a,a)') '*** Error: ', trim(GasName(igas)), &
-               &  ' is not in the required units'
-          call radiation_abort()
+          if (present(istatus)) then
+            istatus = .false.
+          else
+            write(nulerr,'(a,a,a)') '*** Error: ', trim(GasName(igas)), &
+                 &  ' is not in the required units'
+            call radiation_abort()
+          end if
         else if (sf /= this%scale_factor(igas)) then
-          write(nulerr,'(a,a,a,e12.4,a,e12.4)') '*** Error: ', GasName(igas), &
-               &  ' scaling of ', this%scale_factor(igas), &
-               &  ' does not match required ', sf
-          call radiation_abort()
+          if (present(istatus)) then
+            istatus = .false.
+          else
+            write(nulerr,'(a,a,a,e12.4,a,e12.4)') '*** Error: ', GasName(igas), &
+                 &  ' scaling of ', this%scale_factor(igas), &
+                 &  ' does not match required ', sf
+            call radiation_abort()
+          end if
         end if
       end if
     else
-      do ig = 1,this%ntype
-        call this%assert_units(iunits, igas=this%icode(ig), scale_factor=sf)
+      do jg = 1,this%ntype
+        call this%assert_units(iunits, igas=this%icode(jg), scale_factor=sf, istatus=istatus)
       end do
     end if