Changeset 2558 for trunk/LMDZ.MARS/util/simu_MCS.F90

Timestamp:

Sep 7, 2021, 9:57:37 AM (3 years ago)

Author:

abierjon

Message:

Mars GCM:
simu_MCS.F90 : now computes the column-integrated dust optical depths of the MCS and GCM files
(between levels with no missing values in both files, which vary with space and time) and outputs
the ratio tau_GCM/tau_MCS to be used to normalize the GCM dust opacity profiles when comparing to MCSake simu_MCS.F90 compatible with DYNAMICO lon-lat output files, with a check on the latitude array order when doing the interpolation

AB

File:

• trunk/LMDZ.MARS/util/simu_MCS.F90 (modified) (13 diffs)

trunk/LMDZ.MARS/util/simu_MCS.F90

@@ -28 +28 @@
 !                           # dso(/dsodust/qdust)+rho --->   ddust,ndust (dust opacity/km)
 !                           # h2o_ice+rho             --->   dwice,nwice (water ice opacity/km)
+!
+! For dust especially, if a dust opacity variable has been created in the output file, the program also
+! computes the ratio of the Column-integrated Dust Optical Depths from the binned GCM file and the
+! MCS file, that can then be used to normalize the dust opacity profiles when doing comparisons.
 !
 ! Eventually, the program outputs a netcdf file, filled with the GCM binned data for dayside
@@ -91 +95 @@
 !    3.4 Write the data in the netcdf output file
 ! ..var loop ends]
-!  4. END OF THE DAY/NIGHT LOOP
+!  4. CDOD RATIO (tau_GCM/tau_MCS)
+!    4.1 Preliminary checks
+!    4.2 tau_ratio computation
+!    4.3 Write tau_ratio in the netcdf output file
+!  5. END OF THE DAY/NIGHT LOOP
 !..day/night loop ends]
-!  5. CLOSE THE FILES
+!  6. CLOSE THE FILES
 !
 !  Subroutines
@@ -159 +167 @@
 integer :: nbdim                                       ! nb of dimensions of a variable
 real,dimension(:,:,:,:),allocatable :: GCM_var,OBS_var ! the 4D variable extracted from GCM & OBS files
-real,dimension(:,:,:,:),allocatable :: rho             ! atmospheric density for opacities' computation
+real,dimension(:,:,:,:),allocatable :: GCM_rho         ! atmospheric density for opacities' computation
 character(len=64) :: long_name,units,comment           ! netcdf attributes of the variable
 
@@ -199 +207 @@
 integer :: m       ! sol binning loops iteration index
 integer :: v,vnot  ! variable loops indices
+
+!------------------------
+! CDOD ratio (tau_GCM/tau_MCS)
+logical :: tau_ratio_ok ! true if it is possible to compute the CDOD ratio
+
+character(len=64)::out_dztauname,OBS_dztauname,OBS_tempname ! input variable names in files
+integer :: out_dztauid,OBS_dztauid,OBS_tempid               ! input variable ID in files
+
+real,dimension(:,:,:,:),allocatable :: out_dztau,OBS_dztau ! the 4D dust opacities from outfile and obsfile
+real,dimension(:,:,:,:),allocatable :: OBS_temp            ! temperature from obsfile
+
+real,dimension(:),allocatable :: OBS_plev ! pressure levels encompassing each Observer pressure levels
+
+real :: OBS_rho              ! temporary OBS atm density computed from OBS temp and pressure
+real,parameter :: r_atm=191. ! Mars atmosphere specific gas constant
+real,parameter :: g=3.72     ! Mars gravity
+real :: out_tau,OBS_tau      ! integrated dust columns from GCM (outfile) and OBS
+
+real,dimension(:,:,:),allocatable :: tau_ratio ! ratio out_tau/OBS_tau
 
 
@@ -655 +682 @@
   else
     ! Length allocation for each dimension of the 4D variable
-    allocate(rho(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
+    allocate(GCM_rho(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
 
     ! Load datasets
     if (.not.is_stats) then
-      status=NF90_GET_VAR(gcmfid,GCMvarid,rho)
+      status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_rho)
       error_text="Failed to load rho"
       call status_check(status,error_text)
     else
       ! if it is a stats file, we load only the first sol, and then copy it to all the other sols
-      status=NF90_GET_VAR(gcmfid,GCMvarid,rho(:,:,:,1:GCMstatstimelen))
+      status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_rho(:,:,:,1:GCMstatstimelen))
       error_text="Failed to load rho"
       call status_check(status,error_text)
@@ -671 +698 @@
       do l=(GCMstatstimelen+1),GCMtimelen
         if (modulo(l,GCMstatstimelen).ne.0) then
-          rho(:,:,:,l) = rho(:,:,:,modulo(l,GCMstatstimelen))
+          GCM_rho(:,:,:,l) = GCM_rho(:,:,:,modulo(l,GCMstatstimelen))
         else ! if l is a multiple of GCMstatstimelen, since the index modulo(l,GCMstatstimelen)=0
              ! doesn't exist, we make a special case
-          rho(:,:,:,l) = rho(:,:,:,GCMstatstimelen)
+          GCM_rho(:,:,:,l) = GCM_rho(:,:,:,GCMstatstimelen)
         endif
       enddo
@@ -1024 +1051 @@
     GCMvarname = gcm_vars(v)
     
-    ! Detect the dust and wice opacities special cases
+    ! Detect the dust and wice opacities special cases         
     if (trim(GCMvarname).eq."dust") then
+      ! Default : dso > dsodust > dustq
       if (dustok1) then ! "dso" is detected in gcmfile
         GCMvarname="dso"
@@ -1038 +1066 @@
       write(*,*) "Computing dust opacity..."
     endif
+    
     if (trim(GCMvarname).eq."wice") then ! "h2o_ice" detected in gcmfile
       GCMvarname="h2o_ice"
@@ -1142 +1171 @@
             if (GCM_var(i,j,k,l).ne.GCMmiss_val) then
               ! Multiply by rho to have opacity [1/km]
-              GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * rho(i,j,k,l) *1000.
+              GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * GCM_rho(i,j,k,l) *1000.
 
               if (long_name(index(long_name,"(")+1:index(long_name,")")-1).eq."TES") then
@@ -1168 +1197 @@
               ! assuming a rho_dust=3000 kg/m3 and an effective radius reff=1.06 microns
               ! + the opacity is here in 1/km and the MMR in kg/kg
-              GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * rho(i,j,k,l) / 0.012 * 1000 
+              GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * GCM_rho(i,j,k,l) / 0.012 * 1000 
             endif
           enddo
@@ -1197 +1226 @@
              ! We assume a rho_wice=920 kg/m3, an effective radius reff=3um,
              ! an extinction coefficient Qext(wvl,reff)=1.54471
-             GCM_var(i,j,k,l) = 750*1.54471* GCM_var(i,j,k,l) * rho(i,j,k,l) / (920*3e-6)
+             GCM_var(i,j,k,l) = 750*1.54471* GCM_var(i,j,k,l) * GCM_rho(i,j,k,l) / (920*3e-6)
            endif
          enddo
@@ -1512 +1541 @@
     deallocate(outvar)
   ENDDO VAR ! end loop on variables
+
+
+!=================================================================================== 
+! 4. CDOD RATIO (tau_GCM/tau_MCS)
+!===================================================================================
+!===============================================================================  
+! 4.1 Preliminary checks
+!===============================================================================
+  if (OBSalttype.eq.'p') then
+    tau_ratio_ok = .true.
+  else
+    tau_ratio_ok = .false.
+  endif
+  
+  IF (dayornight.EQ."dayside") THEN
+    out_dztauname = "d"
+    OBS_dztauname = "ddust"
+    OBS_tempname = "dtemp"
+  ELSE ! i.e. dayornight="nightside"
+    out_dztauname = "n"
+    OBS_dztauname = "ndust"
+    OBS_tempname = "ntemp"
+  ENDIF
+  
+  ! Check that the output file contains dust opacity
+  ! NB: a dust opacity in outfile requires the dust
+  ! opacity variable to be present in the obsfile,
+  ! so we don't have to check it out
+  if (tau_ratio_ok) then
+    out_dztauname = out_dztauname(1:1)//"dust"
+    status=nf90_inq_varid(outfid,trim(out_dztauname),out_dztauid)
+    if (status.ne.nf90_noerr) then
+      ! if no "d/ndust" in outfile, look for "d/nopadust"
+      out_dztauname = out_dztauname(1:1)//"opadust"
+      status=nf90_inq_varid(outfid,trim(out_dztauname),out_dztauid)
+      if (status.ne.nf90_noerr) then
+        tau_ratio_ok = .false.
+        write(*,*)"No dust opacity in the outfile, we'll pass the computation of the CDOD ratio"
+      endif
+    endif
+  endif
+  
+  ! Check that the obsfile contains temperature (for rho computation)
+  if (tau_ratio_ok) then
+    status=nf90_inq_varid(obsfid,trim(OBS_tempname),OBS_tempid)
+    if (status.ne.nf90_noerr) then
+      tau_ratio_ok = .false.
+      write(*,*)"No temperature in the outfile, we'll pass the computation of the CDOD ratio"
+    endif
+  endif
+
+
+!===============================================================================  
+! 4.2 tau_ratio computation
+!===============================================================================
+  if (tau_ratio_ok) then
+    write(*,*)"Computing the Column Dust Optical Depths ratio (tau_ratio = tau_GCM/tau_MCS)..."
+    write(*,*)"(it is recommended to use it to normalize the GCM dust profiles when comparing with MCS dust profiles)"
+    
+    ! Creation of the OBS_plev table
+    IF (dayornight.EQ."dayside") THEN ! only do it once
+      allocate(OBS_plev(OBSaltlen+1))
+      DO k=2,OBSaltlen
+        ! mid-altitude pressure between 2 OBS pressure layers (OBSalt)
+        ! NB: OBS_plev(k) > OBSalt(k) > OBS_plev(k+1)
+        OBS_plev(k) = sqrt(OBSalt(k-1)*OBSalt(k))
+      ENDDO
+      OBS_plev(1)=2000.
+      OBS_plev(OBSaltlen+1)=0
+      write(*,*)"OBS_plev = ",OBS_plev
+    ENDIF
+    
+    ! Arrays allocation
+    allocate(out_dztau(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
+    allocate(OBS_dztau(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
+    allocate(OBS_temp(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
+    
+    allocate(tau_ratio(OBSlonlen,OBSlatlen,OBSLslen))
+    tau_ratio(:,:,:) = 0
+    
+    ! Load GCM opacity
+    status=NF90_GET_VAR(outfid,out_dztauid,out_dztau)
+    error_text="Failed to load outfile dust opacity"
+    call status_check(status,error_text)
+    
+    ! Load OBS opacity
+    status=nf90_inq_varid(obsfid,trim(OBS_dztauname),OBS_dztauid)
+    status=NF90_GET_VAR(OBSfid,OBS_dztauid,OBS_dztau)
+    error_text="Failed to load obsfile dust opacity"
+    call status_check(status,error_text)
+    
+    ! Load OBS temperature
+    status=NF90_GET_VAR(OBSfid,OBS_tempid,OBS_temp)
+    error_text="Failed to load obsfile temperature"
+    call status_check(status,error_text)
+    
+    do l=1,OBSLslen
+     do j=1,OBSlatlen
+      do i=1,OBSlonlen
+        ! initialize both columns to 0
+        out_tau = 0
+        OBS_tau = 0
+        do k=1,OBSlatlen
+          ! Compute optical depths where both GCM and OBS have non-missing values
+          ! (especially above the highest (in altitude) of the lowest (in altitude) valid values in the GCM and Observer profile)
+          IF ((out_dztau(i,j,k,l).ne.OBSmiss_val).and.(OBS_dztau(i,j,k,l).ne.OBSmiss_val).and.(OBS_temp(i,j,k,l).ne.OBSmiss_val)) then
+            OBS_rho = OBSalt(k) / (r_atm*OBS_temp(i,j,k,l))
+            
+            OBS_tau = OBS_tau + OBS_dztau(i,j,k,l)*1.e-3 / OBS_rho /g *(OBS_plev(k)-OBS_plev(k+1))
+            
+            out_tau = out_tau + out_dztau(i,j,k,l)*1.e-3 / OBS_rho /g *(OBS_plev(k)-OBS_plev(k+1))
+            
+          ENDIF
+        enddo !alt
+        
+        ! Compute the Column Dust Optical Depth ratio
+        if ((out_tau.eq.0).or.(OBS_tau.eq.0)) then
+          tau_ratio(i,j,l) = LTmiss_val
+        else
+          tau_ratio(i,j,l) = out_tau/OBS_tau
+        endif
+        
+      enddo !lon
+     enddo !lat
+    enddo !Ls
+
+
+!===============================================================================  
+! 4.3 Write tau_ratio in the netcdf output file
+!===============================================================================
+    ! Switch to netcdf define mode
+    status=nf90_redef(outfid)
+    error_text="Error: could not switch to define mode in the outfile"
+    call status_check(status,error_text)
+
+    ! Definition of the variable
+    SELECT CASE (dayornight)
+    CASE ("dayside")
+      outvarname = "dtau_ratio"
+    CASE ("nightside")
+      outvarname = "ntau_ratio"
+    END SELECT
+    status=NF90_DEF_VAR(outfid,trim(outvarname),nf90_float,LTshape,outvarid)
+    error_text="Error: could not define the variable "//trim(outvarname)//" in the outfile"
+    call status_check(status,error_text)
+
+    ! Write the attributes
+    long_name = "Integrated optical depths ratio tau_GCM/tau_MCS"
+    SELECT CASE (dayornight)
+    CASE ("dayside")
+      long_name = trim(long_name)//" - day side"
+      status=nf90_put_att(outfid,outvarid,"long_name",long_name)
+    CASE ("nightside")
+      long_name = trim(long_name)//" - night side"
+      status=nf90_put_att(outfid,outvarid,"long_name",long_name)
+    END SELECT
+    status=nf90_put_att(outfid,outvarid,"units","/")
+    status=nf90_put_att(outfid,outvarid,"_FillValue",LTmiss_val)
+    comment = "Reference numbin: dust+temp"
+    status=nf90_put_att(outfid,outvarid,"comment",comment)
+
+    write(*,*)trim(outvarname)," has been created in the outfile"
+    write(*,'("  with missing_value attribute : ",1pe12.5)')LTmiss_val
+    write(*,*)""
+
+    ! End the netcdf define mode (and thus enter the "data writing" mode)
+    status=nf90_enddef(outfid)
+    error_text="Error: could not close the define mode of the outfile"
+    call status_check(status,error_text)
+
+    ! Write data
+    status = nf90_put_var(outfid, outvarid, tau_ratio)
+    error_text="Error: could not write "//trim(outvarname)//" data in the outfile"
+    call status_check(status,error_text)
+
+    ! Deallocations 
+    deallocate(out_dztau)
+    deallocate(OBS_dztau)
+    deallocate(OBS_temp)
+    deallocate(tau_ratio)
+
+  endif !tau_ratio_ok
+  
   
 !===================================================================================  
-! 4. END OF THE DAY/NIGHT LOOP
+! 5. END OF THE DAY/NIGHT LOOP
 !===================================================================================
   IF (dayornight.EQ."dayside") THEN
@@ -1534 +1746 @@
 
 !=================================================================================== 
-! 5. CLOSE THE FILES
+! 6. CLOSE THE FILES
 !===================================================================================
 status = nf90_close(gcmfid)