Changeset 1943 for trunk/LMDZ.TITAN/libf

Timestamp:

Jun 4, 2018, 11:59:42 AM (7 years ago)

Author:

jvatant

Message:

Add XIOS outputs for chemistry including "*_tot" fields (concat GCM+upper atm : 0->1300km )
Only in pseudo-pressure axis for now, pseudo-altitude TBD.
--JVO

Location:

trunk/LMDZ.TITAN/libf/phytitan

Files:

• chem_settings.F90 (modified) (5 diffs)
• comchem_h.F90 (modified) (2 diffs)
• phyetat0_mod.F90 (modified) (1 diff)
• physiq_mod.F90 (modified) (2 diffs)
• xios_output_mod.F90 (modified) (3 diffs)

trunk/LMDZ.TITAN/libf/phytitan/chem_settings.F90

@@ -1 @@
-SUBROUTINE chem_settings(nid,ngrid,indextime)
+SUBROUTINE chem_settings(nid,ngrid,nlayer,indextime)
 
 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -18 +18 @@
 
 USE comchem_h
-
 USE iostart, only: get_field, get_var, inquire_dimension_length
 USE callkeys_mod, only : callchim
+
+USE vertical_layers_mod, only: presnivs, pseudoalt
 
 IMPLICIT NONE
@@ -29 +30 @@
   INTEGER,INTENT(IN)  :: nid            ! Input Netcdf file ID 
   INTEGER,INTENT(IN)  :: ngrid          ! # of horizontal grid points
+  INTEGER,INTENT(IN)  :: nlayer         ! # of vertical layers
   INTEGER,INTENT(IN)  :: indextime      ! position on time axis
 !======================================================================
 ! local variables:
+  REAL :: phi0, phi
+
   INTEGER :: ierr       ! status (returned by NetCDF functions)
   INTEGER :: nvarid     ! ID of NetCDF variable
   INTEGER :: dimid      ! ID of NetCDF dimension
   
-  INTEGER :: iq
+  INTEGER :: ilay, iq
   
   LOGICAL ::  found
@@ -59 +63 @@
   WRITE(*,*) "chem_settings: Upper chemistry pressure grid <preskim> range:", &
                maxval(preskim), minval(preskim) 
+               
+  ! 3. Compute others chemistry grid
+  
+  ! a. Total pressure grid (0->1300km)
+  DO ilay=1,nlayer ! GCM levels
+    preskim_tot(ilay) = presnivs(ilay)
+  ENDDO
+  DO ilay=1,nlaykim_up ! Upper chemistry
+    preskim_tot(ilay+nlayer) = preskim(ilay)
+  ENDDO
+  
+  ! b. Pseudo-altitudes ( TBD - hydrostatic equilibrium or read somewhere ?)          
  
-  ! 3. Inquire ( and load ) upper chemistry composition fields
+  ! 4. Inquire ( and load ) upper chemistry composition fields
  
   CALL get_field("H_up",ykim_up(1,:,:),found,indextime)
@@ -73 +89 @@
     
     DEALLOCATE(ykim_up) ! it will be useless
+    DEALLOCATE(ykim_tot) ! it will be useless
     
   ELSE

trunk/LMDZ.TITAN/libf/phytitan/comchem_h.F90

@@ -78 +78 @@
 
    REAL,SAVE,ALLOCATABLE,DIMENSION(:)     :: preskim  ! Pressure (Pa) of upper chemistry (mid)-layers
-   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: zlay_kim ! Altitude (km) of all chemistry (mid)-layers
+   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: zlaykim  ! Pseudo-altitude (km) of upper chemistry (mid)-layers
    REAL,SAVE,ALLOCATABLE,DIMENSION(:,:,:) :: ykim_up  ! Upper chemistry fields (mol/mol)
-!$OMP_THREADPRIVATE(preskim,zlay_kim,ykim_up)
-
-
+   
+   ! These "_tot" fields are for output only
+   REAL,SAVE,ALLOCATABLE,DIMENSION(:)     :: preskim_tot  ! Pressure (Pa) of total chemistry (mid)-layers
+   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: zlaykim_tot  ! Pseudo-altitude (km) of total chemistry (mid)-layers
+   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:,:) :: ykim_tot     ! Total chemistry fields (mol/mol)
+   
+!$OMP_THREADPRIVATE(preskim,zlaykim,ykim_up)
+!$OMP_THREADPRIVATE(preskim_tot,zlaykim_tot,ykim_tot)
 
    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -116 +121 @@
   
     IF (.NOT.allocated(preskim)) ALLOCATE(preskim(nlaykim_up))
-    IF (.NOT.allocated(zlay_kim)) ALLOCATE(zlay_kim(ngrid,nlaykim_tot))
+    IF (.NOT.allocated(zlaykim)) ALLOCATE(zlaykim(ngrid,nlaykim_up))
     IF (.NOT.allocated(ykim_up)) ALLOCATE(ykim_up(nkim,ngrid,nlaykim_up))
+    
+    IF (.NOT.allocated(preskim_tot)) ALLOCATE(preskim_tot(nlaykim_tot))
+    IF (.NOT.allocated(zlaykim_tot)) ALLOCATE(zlaykim_tot(ngrid,nlaykim_tot))
+    IF (.NOT.allocated(ykim_tot)) ALLOCATE(ykim_tot(nkim,ngrid,nlaykim_tot))
   
   END SUBROUTINE ini_comchem_h

trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90

@@ -274 +274 @@
   ! Call to soil_settings, in order to read upper chemistry
   ! pressure grid as well as composition fields
-  call chem_settings(nid_start,ngrid,indextime)
+  call chem_settings(nid_start,ngrid,nlayer,indextime)
 endif ! of if (startphy_file)
 

trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

@@ -17 +17 @@
       use radcommon_h, only: sigma, glat, grav, BWNV
       use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
-      use comchem_h, only: nkim, cnames
+      use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot
       use comdiurn_h, only: coslat, sinlat, coslon, sinlon
       use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
@@ -1593 +1593 @@
       ! Chemical tracers
       if (callchim) then
+      
+        ! Advected fields
         do iq=1,nkim
-          CALL send_xios_field(cnames(iq),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! mol/mol
+          CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
         enddo
-      endif
+        
+        ! Upper chemistry fields
+        do iq=1,nkim
+          CALL send_xios_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:)) ! mol/mol
+        enddo
+        
+        ! Append fields in ykim_tot for output on the total vertical grid (0->1300km)
+        do iq=1,nkim
+          
+          ! GCM levels
+          do l=1,nlayer
+            ykim_tot(iq,:,l) = zq(:,l,iq+nmicro)/rat_mmol(iq+nmicro)
+          enddo
+          ! Upper levels
+          do l=1,nlaykim_up
+            ykim_tot(iq,:,nlayer+l) = ykim_up(iq,:,l)
+          enddo
+          
+          CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
+          
+        enddo
+        
+      endif ! of 'if callchim'
 
       ! Microphysical tracers

trunk/LMDZ.TITAN/libf/phytitan/xios_output_mod.F90

@@ -35 +35 @@
 #endif
   USE wxios, ONLY: wxios_domain_param, wxios_domain_param_unstructured, wxios_closedef
+  
+  USE comchem_h, ONLY: preskim, preskim_tot, zlaykim, zlaykim_tot
+  
   IMPLICIT NONE
   
@@ -45 +48 @@
   
   
+  INTEGER :: i
   INTEGER :: data_ibegin, data_iend
   TYPE(xios_duration) :: timestep
@@ -60 +64 @@
     CALL xios_set_axis_attr("altitude", n_glo=size(pseudoalt), value=pseudoalt,&
                             unit="km",positive="up")
+    
+    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for preskim"
+    CALL xios_set_axis_attr("preskim", n_glo=size(preskim), value=preskim,&
+                            unit="Pa",positive="down")
+                          
+    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for preskim_tot"
+    CALL xios_set_axis_attr("preskim_tot", n_glo=size(preskim_tot), value=preskim_tot,&
+                            unit="Pa",positive="down")
+                            
+    ! Calculation of pseudo-altitudes is still to be done
+    !IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for zlaykim"
+    !CALL xios_set_axis_attr("zlaykim", n_glo=size(preskim), value=preskim,&
+    !                        unit="km",positive="up")
+    !                      
+    !IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for zlaykim_tot"
+    !CALL xios_set_axis_attr("zlaykim_tot", n_glo=size(zlaykim_tot), value=zlaykim_tot,&
+    !                        unit="km",positive="up")
+    
     
     ! 2. Declare horizontal domain