Ignore:
Timestamp:
Nov 13, 2018 12:14:36 PM (5 years ago)
Author:
suehring
Message:

Sampling of variables extended in virtual measurement module

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/virtual_measurement_mod.f90

    r3494 r3522  
    2525! -----------------
    2626! $Id$
     27! Sampling of variables
     28!
     29! 3494 2018-11-06 14:51:27Z suehring
    2730! Bugfixing
    2831!
    2932! 3473 2018-10-30 20:50:15Z suehring
    3033! Initial revision
    31 !
    32 ! 3472 2018-10-30 20:43:50Z suehring
    3334!
    3435! Authors:
    3536! --------
    36 ! @author Matthias Suehring and Klaus Ketelsen
    37 !
    38 !
     37! @author Matthias Suehring
    3938!
    4039! Description:
     
    4948!> which allows for straight-forward comparison of model results with
    5049!> observations.
     50!>
     51!> @todo list_of_allowed variables needs careful checking
     52!> @todo output (binary or NetCDF) needs to be implemented
     53!> @todo clean-up anything from current test modus
     54!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
     55!>       the (UC)2 standard
     56!> @note Fluxes are not processed
    5157!------------------------------------------------------------------------------!
    5258 MODULE virtual_measurement_mod
    5359
    54 
    5560    USE arrays_3d,                                                             &
    5661        ONLY:  q, pt, u, v, w, zu, zw
    5762
     63    USE chem_modules,                                                          &
     64        ONLY:  nspec
     65
     66    USE chemistry_model_mod,                                                   &
     67        ONLY:  chem_species
     68       
    5869    USE control_parameters,                                                    &
    59         ONLY:  dz, message_string, virtual_measurement
     70        ONLY:  air_chemistry, dz, humidity, neutral, message_string,           &
     71               virtual_measurement
    6072
    6173    USE cpulog,                                                                &
     
    6678
    6779    USE indices,                                                               &
    68         ONLY:  nzb, nzt, nxl, nxr, nys, nyn
     80        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny
    6981
    7082    USE kinds
     
    8092       CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE ::  measured_vars_name !< name of the measured variables
    8193   
    82        INTEGER(iwp) ::  ns   !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
    83        INTEGER(iwp) ::  ntraj !< number of trajectories of a measurement
    84        INTEGER(iwp) ::  nvar  !< number of measured variables
     94       INTEGER(iwp) ::  ns = 0 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
     95       INTEGER(iwp) ::  ntraj  !< number of trajectories of a measurement
     96       INTEGER(iwp) ::  nvar   !< number of measured variables
    8597       
    8698       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
    87        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp   !< number of grid points where observations for a site took place,
    88                                                          !<individual for each trajectory or station that are no _FillValues
    8999       
    90100       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i  !< grid index for measurement position in x-direction
     
    102112       REAL(wp) ::  origin_x_obs                      !< origin of the observation in UTM coordiates in x-direction
    103113       REAL(wp) ::  origin_y_obs                      !< origin of the observation in UTM coordiates in y-direction
    104        
    105        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  xmea   !< measurement x-position in absolute UTM coordinates 
    106        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ymea   !< measurement y-position in absolute UTM coordinates 
    107        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zmea   !< measurement z-position in height above ground level
    108        
     114             
    109115       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars  !< measured variables
    110116       
     
    124130    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
    125131    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
    126    
    127     CHARACTER(LEN=10), DIMENSION(1:53), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
     132!
     133!-- MS: List requires careful revision!
     134    CHARACTER(LEN=10), DIMENSION(1:47), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
    128135       (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
    129136          'hfss      ',  & ! surface sensible heat flux (W/m2)
     
    139146          'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
    140147          'ncaa      ',  & ! number concentation of particles (1/m3)
    141           'mfco2     ',  & ! mole fraction of CO (mol/mol)
     148          'mfco      ',  & ! mole fraction of CO (mol/mol)
    142149          'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
    143150          'mfch4     ',  & ! mole fraction of methane (mol/mol)
     
    164171          'u         ',  & ! u-component
    165172          'ua        ',  & ! eastward wind (is there any difference to u?)
    166           'uw        ',  & ! ? vertical momentum flux - total ?
    167           'utheta    ',  & ! ? horizontal heat flux - total ?
    168           'uv        ',  & ! upward-northward horizontal momentum flux
    169173          'v         ',  & ! v-component
    170174          'va        ',  & ! northward wind (is there any difference to v?)
    171           'vw        ',  & ! ? vertical momentum flux - total ?
    172           'vtheta    ',  & ! ? horizontal heat flux - total ?
    173175          'w         ',  & ! w-component
    174           'wtheta    ',  & ! ? vertical heat flux - total ?
    175176          'rld       ',  & ! downward longwave radiative flux (W/m2)
    176177          'rlu       ',  & ! upnward longwave radiative flux (W/m2)
     
    307308    USE arrays_3d,                                                             &
    308309        ONLY:  zu, zw
    309  
    310     USE control_parameters,                                                    &
    311         ONLY:  message_string
    312310       
    313311    USE grid_variables,                                                        &
     
    329327    CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
    330328    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
    331     CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed
    332    
    333     LOGICAL ::  on_pe !< flag indicating that the respective measurement coordinate is on subdomain
     329    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed   
    334330   
    335331    INTEGER(iwp) ::  dim_eutm  !< dimension size of UTM easting coordinate
     
    338334    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
    339335    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
    340     INTEGER(iwp) ::  icov      !< index range where observations should be taken in x-direction
    341336    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
    342     INTEGER(iwp) ::  i_prev    !< grid index along x for UTM coordinate at previous observation time step
    343337    INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
    344338    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
    345     INTEGER(iwp) ::  jcov      !< index range where observations should be taken in y-direction
    346     INTEGER(iwp) ::  j_prev    !< grid index along y for UTM coordinate at previous observation time step
    347339    INTEGER(iwp) ::  js        !< grid index of real observation point of the respective station in y-direction
    348340    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
    349     INTEGER(iwp) ::  kcov      !< index range where observations should be taken in z-direction
     341    INTEGER(iwp) ::  kl        !< lower vertical index of surrounding grid points of an observation coordinate
    350342    INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
    351     INTEGER(iwp) ::  k_prev    !< grid index along z for UTM coordinate at previous observation time step
    352343    INTEGER(iwp) ::  ksurf     !< topography top index
     344    INTEGER(iwp) ::  ku        !< upper vertical index of surrounding grid points of an observation coordinate
    353345    INTEGER(iwp) ::  l         !< running index over all stations
    354346    INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
     
    359351    INTEGER(iwp) ::  t         !< running index over number of trajectories
    360352   
     353    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
     354   
     355    LOGICAL ::  chem_include !< flag indicating that chemical species is considered in modelled mechanism
     356    LOGICAL ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
     357   
    361358    REAL(wp)     ::  fill_eutm !< _FillValue for coordinate array E_UTM
    362359    REAL(wp)     ::  fill_nutm !< _FillValue for coordinate array N_UTM
     
    370367    CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm,   &
    371368                                global_attribute, 'open', '' )
     369                               
     370!     write(9,*) "num stationi", nvm
     371!     flush(9)
    372372!
    373373!-- ALLOCATE data structure
    374374    ALLOCATE( vmea(1:nvm) )
    375    
    376 !    print*, "nvm", nvm
     375!
     376!-- Allocate flag array
     377    ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     378    meas_flag = 0
     379 
    377380!
    378381!-- Read station coordinates and further attributes.
     
    466469          vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
    467470       ENDDO
    468        
    469 !        print*, "numvars", vmea(l)%nvar, vmea(l)%measured_vars_name(1:vmea(l)%nvar)
     471!
     472!--    In case of chemistry, check if species is considered in the modelled
     473!--    chemistry mechanism.
     474       IF ( air_chemistry )  THEN
     475          DO  ll = 1, vmea(l)%nvar
     476             chem_include = .FALSE.
     477             DO  n = 1, nspec
     478                IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
     479                     TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
     480             ENDDO
     481             IF ( .NOT. chem_include )  THEN
     482                message_string = TRIM( vmea(l)%measured_vars_name(ll) ) //     &
     483                                 ' is not considered in the modelled '  //     &
     484                                 'chemistry mechanism'
     485                CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
     486             ENDIF
     487          ENDDO
     488       ENDIF
    470489!
    471490!--    For the actual measurement ID read the UTM coordinates. Based on these,
     
    502521!--       trajectory or station
    503522          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
    504           ALLOCATE( vmea(l)%ngp(1:vmea(l)%ntraj)   )
    505523!
    506524!--       Allocate temporary arrays for UTM and height coordinates. Note,
     
    536554          ENDIF
    537555!
     556!-- For testing:
     557          e_utm = e_utm - e_utm(1,1)
     558          n_utm = n_utm - n_utm(1,1)
     559         
     560!             
     561!--       First, compute relative x- and y-coordinates with respect to the
     562!--       lower-left origin of the model domain, which is the difference
     563!--       betwen UTM coordinates. 
     564!           e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t))            &
     565!                                       - init_model%origin_x
     566!           n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t))            &
     567!                                       - init_model%origin_y
     568
     569!
    538570!--       Based on UTM coordinates, check if the measurement station or parts
    539571!--       of the trajectory is on subdomain. This case, setup grid index space
    540572!--       sample these quantities.
    541           ns = 0
     573          meas_flag = 0
    542574          DO  t = 1, vmea(l)%ntraj
    543575!
     
    553585                     z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
    554586             ENDDO
    555 !             
    556 !--          First, compute relative x- and y-coordinates with respect to the
    557 !--          lower-left origin of the model domain, which is the difference
    558 !--          betwen UTM coordinates.
    559            
    560              e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t))         &
    561                                          - init_model%origin_x
    562              n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t))         &
    563                                          - init_model%origin_y
     587
    564588!
    565589!--          Compute grid indices relative to origin and check if these are
     
    567591!--          at grid points surrounding the station, hence, check also for
    568592!--          these grid points.
    569              vmea(l)%ngp(t) = 0
    570              k_prev = -999
    571              j_prev = -999
    572              i_prev = -999
    573593             DO  n = 1, vmea(l)%dim_t(t)
    574594                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
     
    579599                          js >= nys  .AND.  js <= nyn )
    580600!
    581 !--             If the measurement is on subdomain, determine vertical index
    582 !--             which refers to the observation height above ground level.
    583                 ks = k_prev
     601!--             Check if observation coordinate is on subdomain
    584602                IF ( on_pe )  THEN
     603!
     604!--                Determine vertical index which correspond to the observation
     605!--                height.
    585606                   ksurf = get_topography_top_index_ji( js, is, 's' )
    586607                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
     608!
     609!--                Set mask array at the observation coordinates. Also, flag the
     610!--                surrounding coordinate points, but first check whether the
     611!--                surrounding coordinate points are on the subdomain.
     612                   kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 > ksurf )
     613                   ku = MERGE( ks+1, ks, ks+1 <= nzt+1 )
     614                 
     615                   meas_flag(kl:ku,js-1:js+1,is-1:is+1) =                      &
     616                               IBSET( meas_flag(kl:ku,js-1:js+1,is-1:is+1), 0 )
    587617                ENDIF
    588 !
    589 !--             Count the number of observation points in index space on
    590 !--             subdomain. Only increment if grid indices are different from
    591 !--             the previous one.
    592                 IF ( on_pe  .AND.  is /= i_prev  .AND.  js /= j_prev  .AND.    &
    593                                    ks /= k_prev )  THEN
    594                    ns             = ns             + 1
    595                    vmea(l)%ngp(t) = vmea(l)%ngp(t) + 1
    596                 ENDIF
    597 
    598 !--             Store arrays for next iteration - avoid double counting
    599                 i_prev = is
    600                 j_prev = js
    601                 k_prev = ks
    602618             ENDDO
    603619             
    604620          ENDDO
    605 
     621!
     622!--       Based on the flag array count the number of observation points.
     623          ns = 0
     624          DO  is = nxl-1, nxr+1
     625             DO  js = nys-1, nyn+1
     626                DO  ks = nzb, nzt+1
     627                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(ks,js,is), 0 ) )
     628                ENDDO
     629             ENDDO
     630          ENDDO
    606631!
    607632!--       Store number of observation points on subdomain and allocate index
    608633!--       arrays.
    609634          vmea(l)%ns = ns
     635          ns = 0
    610636         
    611637          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
    612638          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
    613639          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
     640!
     641!--       Based on the flag array store the grid indices which correspond to
     642!--       the observation coordinates.
     643          DO  is = nxl-1, nxr+1
     644             DO  js = nys-1, nyn+1
     645                DO  ks = nzb, nzt+1
     646                   IF ( BTEST( meas_flag(ks,js,is), 0 ) )  THEN
     647                      ns = ns + 1
     648                      vmea(l)%i(ns) = is
     649                      vmea(l)%j(ns) = js
     650                      vmea(l)%k(ns) = ks
     651                   ENDIF
     652                ENDDO
     653             ENDDO
     654          ENDDO
    614655         
    615 !           print*, "Num ns: ", vmea(l)%ns, "per traj", vmea(l)%ngp(:)
    616 !
    617 !--       Repeat the prior loop and save the grid indices relevant for
    618 !--       sampling.
    619           ns = 0
    620           DO  t = 1, vmea(l)%ntraj
    621 !
    622 !--          Compute grid indices relative to origin and check if these are
    623 !--          on the subdomain. Note, virtual measurements will be taken also
    624 !--          at grid points surrounding the station, hence, check also for
    625 !--          these grid points.
    626              k_prev = -999
    627              j_prev = -999
    628              i_prev = -999
    629              DO  n = 1, vmea(l)%dim_t(t)
    630                 is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
    631                 js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
    632 !
    633 !--             Is the observation point on subdomain?
    634                 on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
    635                           js >= nys  .AND.  js <= nyn )
    636 !
    637 !--             If the measurement is on subdomain, determine vertical index
    638 !--             which refers to the observation height above ground level.
    639                 ks = k_prev
    640                 IF ( on_pe )  THEN
    641                    ksurf = get_topography_top_index_ji( js, is, 's' )
    642                    ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
    643                 ENDIF
    644 !
    645 !--             Count the number of observation points in index space on
    646 !--             subdomain. Only increment if grid indices are different from
    647 !--             the previous one.
    648                 IF ( on_pe  .AND.  is /= i_prev  .AND.  js /= j_prev  .AND.    &
    649                                    ks /= k_prev )  THEN
    650                    ns             = ns + 1
    651                    vmea(l)%i(ns)  = is
    652                    vmea(l)%j(ns)  = js
    653                    vmea(l)%k(ns)  = ks
    654                 ENDIF
    655 !
    656 !--             Store arrays for next iteration - avoid double counting
    657                 i_prev = is
    658                 j_prev = js
    659                 k_prev = ks
    660              ENDDO
    661              
    662           ENDDO
     656!           write(9,*) l, "NS: ", ns
     657!           IF ( ns > 0 )   THEN
     658!              write(9,*) "i ", vmea(l)%i
     659!              write(9,*) "j ", vmea(l)%j
     660!              write(9,*) "k ", vmea(l)%k
     661!              write(9,*)
     662!           ENDIF
    663663!
    664664!--       Allocate array to save the sampled values.
     
    677677       ENDIF
    678678    ENDDO
    679     flush(9)
     679!     flush(9)
    680680   
    681681!
     
    684684    CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '',              &
    685685                                global_attribute, 'close', '' )
     686!                               
     687!-- Dellocate flag array
     688    DEALLOCATE( meas_flag )
     689       
    686690  END SUBROUTINE vm_init
    687691 
     
    694698  SUBROUTINE vm_sampling
    695699
    696     USE arrays_3d !,                                                             &
    697 !         ONLY:  pt
    698 
    699     USE surface_mod
     700    USE arrays_3d,                                                             &
     701        ONLY:  exner, pt, q, u, v, w
     702
     703    USE basic_constants_and_equations_mod,                                     &
     704        ONLY:  pi
     705   
     706    USE radiation_model_mod,                                                   &
     707        ONLY:  radiation 
     708
     709    USE surface_mod,                                                           &
     710        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
    700711   
    701712     IMPLICIT NONE
     
    703714     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
    704715     
    705      INTEGER(iwp) ::  l !<
    706      INTEGER(iwp) ::  m !<
    707      INTEGER(iwp) ::  var !<
    708      
    709      INTEGER(iwp) ::  mm, j, i
    710 
     716     INTEGER(iwp) ::  i  !< grid index in x-direction
     717     INTEGER(iwp) ::  j  !< grid index in y-direction
     718     INTEGER(iwp) ::  k  !< grid index in z-direction
     719     INTEGER(iwp) ::  l  !< running index over the number of stations
     720     INTEGER(iwp) ::  m  !< running index over all virtual observation coordinates
     721     INTEGER(iwp) ::  mm !< index of surface element which corresponds to the virtual observation coordinate
     722     INTEGER(iwp) ::  n  !< running index over all measured variables at a station
     723     INTEGER(iwp) ::  nn !< running index over the number of chemcal species
    711724!
    712725!--  Loop over all stations. For each possible variable loop over all
     
    714727     DO  l = 1, nvm
    715728!
    716 !--     Loop over all measured variables. Please note, for the moment
    717 !--     the same indices for scalar and velocity components are used.
    718 !--     ToDo: Revise this later.
    719         DO  m = 1, vmea(l)%ns
    720            j = vmea(l)%j(m)
    721            i = vmea(l)%i(m)
    722 !           
    723 !            IF ( i >= nxl  .AND.  i <= nxr  .AND.                               &
    724 !                 j >= nys  .AND.  j <= nyn )  THEN
    725 !               IF ( surf_def_h(0)%start_index(j,i) <= &
    726 !                    surf_def_h(0)%end_index(j,i) )  THEN
    727 !                  mm = surf_def_h(0)%end_index(j,i)
    728 !               
    729 !                  surf_def_h(0)%pt_surface(mm) = -99.0
    730 !               ENDIF
    731 !            ENDIF
    732 !         ENDDO
    733 !         DO  var = 1, vmea(l)%nvar
    734 !            trimvar = TRIM( vmea(l)%measured_vars_name(var) )
    735 !           
    736 !            IF ( TRIM( trimvar ) == 'theta' )  THEN           
    737 !               DO  m = 1, vmea(l)%ns
    738 !                  vmea(l)%measured_vars(var,m) = pt(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
    739 !               ENDDO
    740 !            ENDIF           
    741 !            IF ( TRIM( trimvar ) == 'w' )  THEN           
    742 !               DO  m = 1, vmea(l)%ns
    743 !                  vmea(l)%measured_vars(var,m) = w(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
    744 !               ENDDO
    745 !            ENDIF
    746 !            IF ( TRIM( trimvar ) == 'v' )  THEN           
    747 !               DO  m = 1, vmea(l)%ns
    748 !                  vmea(l)%measured_vars(var,m) = v(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
    749 !               ENDDO
    750 !            ENDIF
    751 !            IF ( TRIM( trimvar ) == 'u' )  THEN           
    752 !               DO  m = 1, vmea(l)%ns
    753 !                  vmea(l)%measured_vars(var,m) = u(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
    754 !               ENDDO
    755 !            ENDIF
     729!--     Loop over all measured variables at this station. Please note,
     730!--     velocity components are interpolated onto scalar grid. 
     731        DO  n = 1, vmea(l)%nvar
     732       
     733           SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
     734           
     735              CASE ( 'theta' )
     736                 IF ( .NOT. neutral )  THEN
     737                    DO  m = 1, vmea(l)%ns
     738                       k = vmea(l)%k(m)
     739                       j = vmea(l)%j(m)
     740                       i = vmea(l)%i(m)
     741                       vmea(l)%measured_vars(n,m) = pt(k,j,i)
     742                    ENDDO
     743                 ENDIF
     744                 
     745              CASE ( 'ta', 't_va' )
     746                 IF ( .NOT. neutral )  THEN
     747                    DO  m = 1, vmea(l)%ns
     748                       k = vmea(l)%k(m)
     749                       j = vmea(l)%j(m)
     750                       i = vmea(l)%i(m)
     751                       vmea(l)%measured_vars(n,m) = pt(k,j,i) * exner( k )
     752                    ENDDO
     753                 ENDIF
     754                 
     755              CASE ( 'hus', 'haa' )
     756                 IF ( humidity )  THEN
     757                    DO  m = 1, vmea(l)%ns
     758                       k = vmea(l)%k(m)
     759                       j = vmea(l)%j(m)
     760                       i = vmea(l)%i(m)
     761                       vmea(l)%measured_vars(n,m) = q(k,j,i)
     762                    ENDDO
     763                 ENDIF
     764                 
     765              CASE ( 'u', 'ua' )
     766                 DO  m = 1, vmea(l)%ns
     767                    k = vmea(l)%k(m)
     768                    j = vmea(l)%j(m)
     769                    i = vmea(l)%i(m)
     770                    vmea(l)%measured_vars(n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
     771                 ENDDO
     772                 
     773              CASE ( 'v', 'va' )
     774                 DO  m = 1, vmea(l)%ns
     775                    k = vmea(l)%k(m)
     776                    j = vmea(l)%j(m)
     777                    i = vmea(l)%i(m)
     778                    vmea(l)%measured_vars(n,m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
     779                 ENDDO
     780                 
     781              CASE ( 'w' )
     782                 DO  m = 1, vmea(l)%ns
     783                    k = vmea(l)%k(m)
     784                    j = vmea(l)%j(m)
     785                    i = vmea(l)%i(m)
     786                    vmea(l)%measured_vars(n,m) = w(k,j,i) !0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     787                 ENDDO
     788                 
     789              CASE ( 'wspeed' )
     790                 DO  m = 1, vmea(l)%ns
     791                    k = vmea(l)%k(m)
     792                    j = vmea(l)%j(m)
     793                    i = vmea(l)%i(m)
     794                    vmea(l)%measured_vars(n,m) = SQRT(                         &
     795                                   ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
     796                                   ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2   &
     797                                                     )
     798                 ENDDO
     799                 
     800              CASE ( 'wdir' )
     801                 DO  m = 1, vmea(l)%ns
     802                    k = vmea(l)%k(m)
     803                    j = vmea(l)%j(m)
     804                    i = vmea(l)%i(m)
     805                   
     806                    vmea(l)%measured_vars(n,m) = ATAN2(                        &
     807                                       - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ),   &
     808                                       - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )    &
     809                                                      ) * 180.0_wp / pi
     810                 ENDDO
     811!
     812!--           MS: list of variables needs extension.
     813              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mcco', 'mcco2', 'mcbcda',  &
     814                     'ncaa', 'mfco2', 'mfco', 'mfch4', 'mfnh3', 'mfno',        &
     815                     'mfno2', 'mfso2', 'mfh20', 'tr03' )
     816                 IF ( air_chemistry )  THEN                 
     817                    DO  nn = 1, nspec                   
     818                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
     819                            TRIM( chem_species(nn)%name ) )  THEN                           
     820                          DO  m = 1, vmea(l)%ns             
     821                             k = vmea(l)%k(m)
     822                             j = vmea(l)%j(m)
     823                             i = vmea(l)%i(m)                   
     824                             vmea(l)%measured_vars(n,m) =                      &
     825                                                   chem_species(nn)%conc(k,j,i)
     826                          ENDDO
     827                       ENDIF
     828                    ENDDO
     829                 ENDIF
     830                 
     831              CASE ( 'us' )
     832                 DO  m = 1, vmea(l)%ns
     833!
     834!--                 Surface data is only available on inner subdomains, not
     835!--                 on ghost points. Hence, limit the indices.
     836                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     837                    j = MERGE( j           , nyn, j            > nyn )
     838                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     839                    i = MERGE( i           , nxr, i            > nxr )
     840                   
     841                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
     842                             surf_def_h(0)%end_index(j,i)
     843                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%us(mm)
     844                    ENDDO
     845                    DO  mm = surf_lsm_h%start_index(j,i),                      &
     846                             surf_lsm_h%end_index(j,i)
     847                       vmea(l)%measured_vars(n,m) = surf_lsm_h%us(mm)
     848                    ENDDO
     849                    DO  mm = surf_usm_h%start_index(j,i),                      &
     850                             surf_usm_h%end_index(j,i)
     851                       vmea(l)%measured_vars(n,m) = surf_usm_h%us(mm)
     852                    ENDDO
     853                 ENDDO
     854                 
     855              CASE ( 'ts' )
     856                 DO  m = 1, vmea(l)%ns
     857!
     858!--                 Surface data is only available on inner subdomains, not
     859!--                 on ghost points. Hence, limit the indices.
     860                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     861                    j = MERGE( j           , nyn, j            > nyn )
     862                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     863                    i = MERGE( i           , nxr, i            > nxr )
     864                   
     865                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
     866                             surf_def_h(0)%end_index(j,i)
     867                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%ts(mm)
     868                    ENDDO
     869                    DO  mm = surf_lsm_h%start_index(j,i),                      &
     870                             surf_lsm_h%end_index(j,i)
     871                       vmea(l)%measured_vars(n,m) = surf_lsm_h%ts(mm)
     872                    ENDDO
     873                    DO  mm = surf_usm_h%start_index(j,i),                      &
     874                             surf_usm_h%end_index(j,i)
     875                       vmea(l)%measured_vars(n,m) = surf_usm_h%ts(mm)
     876                    ENDDO
     877                 ENDDO
     878                 
     879              CASE ( 'hfls' )
     880                 DO  m = 1, vmea(l)%ns
     881!
     882!--                 Surface data is only available on inner subdomains, not
     883!--                 on ghost points. Hence, limit the indices.
     884                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     885                    j = MERGE( j           , nyn, j            > nyn )
     886                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     887                    i = MERGE( i           , nxr, i            > nxr )
     888                   
     889                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
     890                             surf_def_h(0)%end_index(j,i)
     891                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%qsws(mm)
     892                    ENDDO
     893                    DO  mm = surf_lsm_h%start_index(j,i),                      &
     894                             surf_lsm_h%end_index(j,i)
     895                       vmea(l)%measured_vars(n,m) = surf_lsm_h%qsws(mm)
     896                    ENDDO
     897                    DO  mm = surf_usm_h%start_index(j,i),                      &
     898                             surf_usm_h%end_index(j,i)
     899                       vmea(l)%measured_vars(n,m) = surf_usm_h%qsws(mm)
     900                    ENDDO
     901                 ENDDO
     902                 
     903              CASE ( 'hfss' )
     904                 DO  m = 1, vmea(l)%ns
     905!
     906!--                 Surface data is only available on inner subdomains, not
     907!--                 on ghost points. Hence, limit the indices.
     908                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     909                    j = MERGE( j           , nyn, j            > nyn )
     910                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     911                    i = MERGE( i           , nxr, i            > nxr )
     912                   
     913                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
     914                             surf_def_h(0)%end_index(j,i)
     915                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%shf(mm)
     916                    ENDDO
     917                    DO  mm = surf_lsm_h%start_index(j,i),                      &
     918                             surf_lsm_h%end_index(j,i)
     919                       vmea(l)%measured_vars(n,m) = surf_lsm_h%shf(mm)
     920                    ENDDO
     921                    DO  mm = surf_usm_h%start_index(j,i),                      &
     922                             surf_usm_h%end_index(j,i)
     923                       vmea(l)%measured_vars(n,m) = surf_usm_h%shf(mm)
     924                    ENDDO
     925                 ENDDO
     926                 
     927              CASE ( 'rnds' )
     928                 IF ( radiation )  THEN
     929                    DO  m = 1, vmea(l)%ns
     930!
     931!--                    Surface data is only available on inner subdomains, not
     932!--                    on ghost points. Hence, limit the indices.
     933                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     934                       j = MERGE( j           , nyn, j            > nyn )
     935                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     936                       i = MERGE( i           , nxr, i            > nxr )
     937                   
     938                       DO  mm = surf_lsm_h%start_index(j,i),                   &
     939                                surf_lsm_h%end_index(j,i)
     940                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_net(mm)
     941                       ENDDO
     942                       DO  mm = surf_usm_h%start_index(j,i),                   &
     943                                surf_usm_h%end_index(j,i)
     944                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_net(mm)
     945                       ENDDO
     946                    ENDDO
     947                 ENDIF
     948                 
     949              CASE ( 'rsus', 'rsu' )
     950                 IF ( radiation )  THEN
     951                    DO  m = 1, vmea(l)%ns
     952!
     953!--                    Surface data is only available on inner subdomains, not
     954!--                    on ghost points. Hence, limit the indices.
     955                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     956                       j = MERGE( j           , nyn, j            > nyn )
     957                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     958                       i = MERGE( i           , nxr, i            > nxr )
     959                   
     960                       DO  mm = surf_lsm_h%start_index(j,i),                   &
     961                                surf_lsm_h%end_index(j,i)
     962                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_out(mm)
     963                       ENDDO
     964                       DO  mm = surf_usm_h%start_index(j,i),                   &
     965                                surf_usm_h%end_index(j,i)
     966                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_out(mm)
     967                       ENDDO
     968                    ENDDO
     969                 ENDIF
     970                 
     971              CASE ( 'rsds', 'rsd' )
     972                 IF ( radiation )  THEN
     973                    DO  m = 1, vmea(l)%ns
     974!
     975!--                    Surface data is only available on inner subdomains, not
     976!--                    on ghost points. Hence, limit the indices.
     977                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     978                       j = MERGE( j           , nyn, j            > nyn )
     979                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     980                       i = MERGE( i           , nxr, i            > nxr )
     981                   
     982                       DO  mm = surf_lsm_h%start_index(j,i),                   &
     983                                surf_lsm_h%end_index(j,i)
     984                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_in(mm)
     985                       ENDDO
     986                       DO  mm = surf_usm_h%start_index(j,i),                   &
     987                                surf_usm_h%end_index(j,i)
     988                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_in(mm)
     989                       ENDDO
     990                    ENDDO
     991                 ENDIF
     992                 
     993              CASE ( 'rlus', 'rlu' )
     994                 IF ( radiation )  THEN
     995                    DO  m = 1, vmea(l)%ns
     996!
     997!--                    Surface data is only available on inner subdomains, not
     998!--                    on ghost points. Hence, limit the indices.
     999                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     1000                       j = MERGE( j           , nyn, j            > nyn )
     1001                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     1002                       i = MERGE( i           , nxr, i            > nxr )
     1003                   
     1004                       DO  mm = surf_lsm_h%start_index(j,i),                   &
     1005                                surf_lsm_h%end_index(j,i)
     1006                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_out(mm)
     1007                       ENDDO
     1008                       DO  mm = surf_usm_h%start_index(j,i),                   &
     1009                                surf_usm_h%end_index(j,i)
     1010                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_out(mm)
     1011                       ENDDO
     1012                    ENDDO
     1013                 ENDIF
     1014                 
     1015              CASE ( 'rlds', 'rld' )
     1016                 IF ( radiation )  THEN
     1017                    DO  m = 1, vmea(l)%ns
     1018!
     1019!--                    Surface data is only available on inner subdomains, not
     1020!--                    on ghost points. Hence, limit the indices.
     1021                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
     1022                       j = MERGE( j           , nyn, j            > nyn )
     1023                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
     1024                       i = MERGE( i           , nxr, i            > nxr )
     1025                   
     1026                       DO  mm = surf_lsm_h%start_index(j,i),                   &
     1027                                surf_lsm_h%end_index(j,i)
     1028                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_in(mm)
     1029                       ENDDO
     1030                       DO  mm = surf_usm_h%start_index(j,i),                   &
     1031                                surf_usm_h%end_index(j,i)
     1032                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_in(mm)
     1033                       ENDDO
     1034                    ENDDO
     1035                 ENDIF
     1036!
     1037!--           More will follow ...
     1038                 
     1039           END SELECT
     1040
    7561041        ENDDO
    7571042
Note: See TracChangeset for help on using the changeset viewer.