Changeset 2435 for LMDZ5/branches/testing/libf/phylmd/cosp/cosp_radar.F90

Timestamp:

Jan 28, 2016, 5:02:13 PM (9 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2396:2434 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/cosp/cosp_radar.F90 (modified) (7 diffs)

LMDZ5/branches/testing/libf/phylmd/cosp/cosp_radar.F90

@@ -26 +26 @@
   USE MOD_COSP_CONSTANTS
   USE MOD_COSP_TYPES
+  USE MOD_COSP_UTILS
   use radar_simulator_types
   use array_lib
@@ -31 +32 @@
   use format_input
   IMPLICIT NONE
-  
+
   INTERFACE
-    subroutine radar_simulator(freq,k2,do_ray,use_gas_abs,use_mie_table,mt, &
-        nhclass,hp,nprof,ngate,nsizes,D,hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix, &
-        rh_matrix,Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe, &
+    subroutine radar_simulator(hp,nprof,ngate,undef, &
+        hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
+        p_matrix,t_matrix,rh_matrix, &
+        Ze_non,Ze_ray,g_to_vol,a_to_vol,dBZe, &
         g_to_vol_in,g_to_vol_out)
-  
-        use m_mrgrnk 
+
+        use m_mrgrnk
         use array_lib
         use math_lib
@@ -44 +46 @@
         use radar_simulator_types
         implicit none
+
         ! ----- INPUTS -----  
-        type(mie), intent(in) :: mt
         type(class_param) :: hp
-        real*8, intent(in) :: freq,k2
-        integer, intent(in) ::  do_ray,use_gas_abs,use_mie_table, &
-            nhclass,nprof,ngate,nsizes
-        real*8, dimension(nsizes), intent(in) :: D
+
+        integer, intent(in) :: nprof,ngate
+
+        real undef
         real*8, dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, &
             t_matrix,rh_matrix
-        real*8, dimension(nhclass,nprof,ngate), intent(in) :: hm_matrix
-        real*8, dimension(nhclass,nprof,ngate), intent(inout) :: re_matrix
+        real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
+        real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
+        real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: Np_matrix
+
         ! ----- OUTPUTS -----
         real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
-            g_atten_to_vol,dBZe,h_atten_to_vol    
+            g_to_vol,dBZe,a_to_vol
         ! ----- OPTIONAL -----
-        real*8, optional, dimension(ngate,nprof) :: &
+        real*8, optional, dimension(nprof,ngate) :: &
             g_to_vol_in,g_to_vol_out
      end subroutine radar_simulator
@@ -73 +77 @@
 
   ! Arguments
-  type(cosp_gridbox),intent(in) :: gbx  ! Gridbox info
+  type(cosp_gridbox),intent(inout) :: gbx  ! Gridbox info
   type(cosp_subgrid),intent(in) :: sgx  ! Subgrid info
   type(cosp_sghydro),intent(in) :: sghydro  ! Subgrid info for hydrometeors
@@ -80 +84 @@
   ! Local variables 
   integer :: & 
-  nsizes                        ! num of discrete drop sizes
+  nsizes            ! num of discrete drop sizes
 
-  real*8 :: &
-  freq, &                       ! radar frequency (GHz)
-  k2                            ! |K|^2, -1=use frequency dependent default
-  
   real*8, dimension(:,:), allocatable :: &
   g_to_vol ! integrated atten due to gases, r>v (dB)
-  
+
   real*8, dimension(:,:), allocatable :: &
-  Ze_non, &                     ! radar reflectivity withOUT attenuation (dBZ)
-  Ze_ray, &                     ! Rayleigh reflectivity (dBZ)
-  h_atten_to_vol, &             ! attenuation by hydromets, radar to vol (dB)
-  g_atten_to_vol, &             ! gaseous atteunation, radar to vol (dB)
-  dBZe, &                       ! effective radar reflectivity factor (dBZ)
-  hgt_matrix, &                 ! height of hydrometeors (km)
+  Ze_non, &         ! radar reflectivity withOUT attenuation (dBZ)
+  Ze_ray, &         ! Rayleigh reflectivity (dBZ)
+  h_atten_to_vol, &     ! attenuation by hydromets, radar to vol (dB)
+  g_atten_to_vol, &     ! gaseous atteunation, radar to vol (dB)
+  dBZe, &           ! effective radar reflectivity factor (dBZ)
+  hgt_matrix, &         ! height of hydrometeors (km)
   t_matrix, &                   !temperature (k)
   p_matrix, &                   !pressure (hPa)
   rh_matrix                     !relative humidity (%)
-  
+
   real*8, dimension(:,:,:), allocatable :: &
-  hm_matrix, &                  ! hydrometeor mixing ratio (g kg^-1)
-  re_matrix
+  hm_matrix, &          ! hydrometeor mixing ratio (g kg^-1)
+  re_matrix, &          ! effective radius (microns).   Optional. 0 ==> use Np_matrix or defaults
+  Np_matrix         ! total number concentration (kg^-1).   Optional 0==> use defaults 
 
   integer, parameter :: one = 1
-  logical :: hgt_reversed
-  integer :: pr,i,j,k,unt
+  ! logical :: hgt_reversed
+  logical :: hgt_descending
+  integer :: pr,i,j,k,unt,ngate
 
 ! ----- main program settings ------
-
-  freq = gbx%radar_freq
-  k2 = gbx%k2
- 
-  !
-  ! note:  intitialization section that was here has been relocated to SUBROUTINE CONSTRUCT_COSP_GRIDBOX by roj, Feb 2008
-  !
-  mt_ttl=gbx%mt_ttl  ! these variables really should be moved into the mt structure rather than kept as global arrays.
-  mt_tti=gbx%mt_tti
 
   ! Inputs to Quickbeam
   allocate(hgt_matrix(gbx%Npoints,gbx%Nlevels),p_matrix(gbx%Npoints,gbx%Nlevels), &
            t_matrix(gbx%Npoints,gbx%Nlevels),rh_matrix(gbx%Npoints,gbx%Nlevels))
-  allocate(hm_matrix(gbx%Nhydro,gbx%Npoints,gbx%Nlevels)) 
+  allocate(hm_matrix(gbx%Nhydro,gbx%Npoints,gbx%Nlevels))
   allocate(re_matrix(gbx%Nhydro,gbx%Npoints,gbx%Nlevels))
+  allocate(Np_matrix(gbx%Nhydro,gbx%Npoints,gbx%Nlevels))
 
   ! Outputs from Quickbeam
@@ -131 +125 @@
   allocate(g_atten_to_vol(gbx%Npoints,gbx%Nlevels))
   allocate(dBZe(gbx%Npoints,gbx%Nlevels))
-  
+
   ! Optional argument. It is computed and returned in the first call to
   ! radar_simulator, and passed as input in the rest
-  allocate(g_to_vol(gbx%Nlevels,gbx%Npoints))
-  
+  allocate(g_to_vol(gbx%Npoints,gbx%Nlevels))
+
+  ! Even if there is no unit conversion, they are needed for type conversion
   p_matrix   = gbx%p/100.0     ! From Pa to hPa
   hgt_matrix = gbx%zlev/1000.0 ! From m to km
-  t_matrix   = gbx%T-273.15    ! From K to C
+  t_matrix   = gbx%T
   rh_matrix  = gbx%q
   re_matrix  = 0.0
-  
-  ! Quickbeam assumes the first row is closest to the radar
-  call order_data(hgt_matrix,hm_matrix,p_matrix,t_matrix, &
-      rh_matrix,gbx%surface_radar,hgt_reversed)
-  
+
+
+  ! set flag denoting position of radar relative to hgt_matrix orientation
+          ngate = size(hgt_matrix,2)
+
+          hgt_descending = hgt_matrix(1,1) > hgt_matrix(1,ngate)
+
+          if ( &
+             (gbx%surface_radar == 1 .and. hgt_descending) .or.  &
+             (gbx%surface_radar == 0 .and. (.not. hgt_descending)) &
+             ) &
+          then
+            gbx%hp%radar_at_layer_one = .false.
+          else
+            gbx%hp%radar_at_layer_one = .true.
+          endif
+
   ! ----- loop over subcolumns -----
   do pr=1,sgx%Ncolumns
+
+      !  NOTE:
       !  atmospheric profiles are the same within the same gridbox
-      !  only hydrometeor profiles will be different
-      if (hgt_reversed) then  
-         do i=1,gbx%Nhydro  
-            hm_matrix(i,:,:) = sghydro%mr_hydro(:,pr,gbx%Nlevels:1:-1,i)*1000.0 ! Units from kg/kg to g/kg
-            if (gbx%use_reff) then
-              re_matrix(i,:,:) = sghydro%Reff(:,pr,gbx%Nlevels:1:-1,i)*1.e6     ! Units from m to micron
-            endif
-         enddo  
-      else  
+      !  only hydrometeor profiles will be different for each subgridbox
+
          do i=1,gbx%Nhydro
             hm_matrix(i,:,:) = sghydro%mr_hydro(:,pr,:,i)*1000.0 ! Units from kg/kg to g/kg
             if (gbx%use_reff) then
               re_matrix(i,:,:) = sghydro%Reff(:,pr,:,i)*1.e6       ! Units from m to micron
+              Np_matrix(i,:,:) = sghydro%Np(:,pr,:,i)              ! Units [#/kg]
             endif
          enddo
-      endif  
 
       !   ----- call radar simulator -----
       if (pr == 1) then ! Compute gaseous attenuation for all profiles
-         j=0
-         if (gbx%Npoints == 53) then
-           unt=10
-           j=1
-         endif
-         if (gbx%Npoints == 153) then
-           unt=11
-           j=101
-         endif
-         call radar_simulator(freq,k2,gbx%do_ray,gbx%use_gas_abs,gbx%use_mie_tables,gbx%mt, &    !  v0.2: mt changed to gbx%mt, roj
-           gbx%Nhydro,gbx%hp,gbx%Npoints,gbx%Nlevels,gbx%nsizes,gbx%D, &                         !  v0.2: hp->gbx%hp, D->gbx%d, nsizes->gbx%nsizes, roj
-           hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix,rh_matrix, &
+         call radar_simulator(gbx%hp,gbx%Npoints,gbx%Nlevels,R_UNDEF, &
+           hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
+           p_matrix,t_matrix,rh_matrix, &
            Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe,g_to_vol_out=g_to_vol)
       else ! Use gaseous atteunuation for pr = 1
-         call radar_simulator(freq,k2,gbx%do_ray,gbx%use_gas_abs,gbx%use_mie_tables,gbx%mt, &
-           gbx%Nhydro,gbx%hp,gbx%Npoints,gbx%Nlevels,gbx%nsizes,gbx%D, &
-           hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix,rh_matrix, &
+         call radar_simulator(gbx%hp,gbx%Npoints,gbx%Nlevels,R_UNDEF, &
+           hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
+           p_matrix,t_matrix,rh_matrix, &
            Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe,g_to_vol_in=g_to_vol)
       endif
-      ! ----- BEGIN output section -----
-      ! spaceborne radar : from TOA to SURFACE
-      if (gbx%surface_radar == 1) then
-        z%Ze_tot(:,pr,:)=dBZe(:,:)
-      else if (gbx%surface_radar == 0) then ! Spaceborne
-        z%Ze_tot(:,pr,:)=dBZe(:,gbx%Nlevels:1:-1)
-      endif
 
+      ! store caluculated dBZe values for later output/processing
+      z%Ze_tot(:,pr,:)=dBZe(:,:)
   enddo !pr
-  
-  ! Change undefined value to one defined in COSP
-  where (z%Ze_tot == -999.0) z%Ze_tot = R_UNDEF
 
   deallocate(hgt_matrix,p_matrix,t_matrix,rh_matrix)
@@ -204 +189 @@
       Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe)
   deallocate(g_to_vol)
- 
-  ! deallocate(mt_ttl,mt_tti)   !v0.2: roj feb 2008 can not be done here,
-                                !these variables now part of gbx structure and dealocated later
-
 END SUBROUTINE COSP_RADAR