Changeset 3522 for palm/trunk
- Timestamp:
- Nov 13, 2018 12:14:36 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
Makefile (modified) (2 diffs)
-
virtual_measurement_mod.f90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r3494 r3522 25 25 # ----------------- 26 26 # $Id$ 27 # Dependencies for virtual measurement module added 28 # 29 # 3494 2018-11-06 14:51:27Z suehring 27 30 # Surface output revised 28 31 # … … 1982 1985 virtual_measurement_mod.o: \ 1983 1986 cpulog_mod.o \ 1987 chemistry_model_mod.o \ 1988 chem_modules.o \ 1984 1989 mod_kinds.o \ 1985 1990 modules.o \ 1986 1991 netcdf_data_input_mod.o \ 1987 netcdf_interface_mod.o 1992 netcdf_interface_mod.o \ 1993 radiation_model_mod.o 1988 1994 wind_turbine_model_mod.o: \ 1989 1995 basic_constants_and_equations_mod.o \ -
palm/trunk/SOURCE/virtual_measurement_mod.f90
r3494 r3522 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Sampling of variables 28 ! 29 ! 3494 2018-11-06 14:51:27Z suehring 27 30 ! Bugfixing 28 31 ! 29 32 ! 3473 2018-10-30 20:50:15Z suehring 30 33 ! Initial revision 31 !32 ! 3472 2018-10-30 20:43:50Z suehring33 34 ! 34 35 ! Authors: 35 36 ! -------- 36 ! @author Matthias Suehring and Klaus Ketelsen 37 ! 38 ! 37 ! @author Matthias Suehring 39 38 ! 40 39 ! Description: … … 49 48 !> which allows for straight-forward comparison of model results with 50 49 !> 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 51 57 !------------------------------------------------------------------------------! 52 58 MODULE virtual_measurement_mod 53 59 54 55 60 USE arrays_3d, & 56 61 ONLY: q, pt, u, v, w, zu, zw 57 62 63 USE chem_modules, & 64 ONLY: nspec 65 66 USE chemistry_model_mod, & 67 ONLY: chem_species 68 58 69 USE control_parameters, & 59 ONLY: dz, message_string, virtual_measurement 70 ONLY: air_chemistry, dz, humidity, neutral, message_string, & 71 virtual_measurement 60 72 61 73 USE cpulog, & … … 66 78 67 79 USE indices, & 68 ONLY: nzb, nzt, nxl, nxr, nys, nyn 80 ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny 69 81 70 82 USE kinds … … 80 92 CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables 81 93 82 INTEGER(iwp) :: ns 83 INTEGER(iwp) :: ntraj !< number of trajectories of a measurement84 INTEGER(iwp) :: nvar !< number of measured variables94 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 85 97 86 98 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 _FillValues89 99 90 100 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction … … 102 112 REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction 103 113 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 109 115 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables 110 116 … … 124 130 CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements 125 131 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 128 135 (/ 'hfls ', & ! surface latent heat flux (W/m2) 129 136 'hfss ', & ! surface sensible heat flux (W/m2) … … 139 146 'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3) 140 147 'ncaa ', & ! number concentation of particles (1/m3) 141 'mfco 2', & ! mole fraction of CO (mol/mol)148 'mfco ', & ! mole fraction of CO (mol/mol) 142 149 'mfco2 ', & ! mole fraction of CO2 (mol/mol) 143 150 'mfch4 ', & ! mole fraction of methane (mol/mol) … … 164 171 'u ', & ! u-component 165 172 '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 flux169 173 'v ', & ! v-component 170 174 'va ', & ! northward wind (is there any difference to v?) 171 'vw ', & ! ? vertical momentum flux - total ?172 'vtheta ', & ! ? horizontal heat flux - total ?173 175 'w ', & ! w-component 174 'wtheta ', & ! ? vertical heat flux - total ?175 176 'rld ', & ! downward longwave radiative flux (W/m2) 176 177 'rlu ', & ! upnward longwave radiative flux (W/m2) … … 307 308 USE arrays_3d, & 308 309 ONLY: zu, zw 309 310 USE control_parameters, &311 ONLY: message_string312 310 313 311 USE grid_variables, & … … 329 327 CHARACTER(LEN=5) :: dum !< dummy string indicate station id 330 328 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 334 330 335 331 INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate … … 338 334 INTEGER(iwp) :: dim_zag !< dimension size of height coordinate 339 335 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-direction341 336 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 step343 337 INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction 344 338 INTEGER(iwp) :: j !< grid index of observation point in x-direction 345 INTEGER(iwp) :: jcov !< index range where observations should be taken in y-direction346 INTEGER(iwp) :: j_prev !< grid index along y for UTM coordinate at previous observation time step347 339 INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction 348 340 INTEGER(iwp) :: k !< grid index of observation point in x-direction 349 INTEGER(iwp) :: k cov !< index range where observations should be taken in z-direction341 INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate 350 342 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 step352 343 INTEGER(iwp) :: ksurf !< topography top index 344 INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate 353 345 INTEGER(iwp) :: l !< running index over all stations 354 346 INTEGER(iwp) :: len_char !< character length of single measured variables without Null character … … 359 351 INTEGER(iwp) :: t !< running index over number of trajectories 360 352 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 361 358 REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM 362 359 REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM … … 370 367 CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm, & 371 368 global_attribute, 'open', '' ) 369 370 ! write(9,*) "num stationi", nvm 371 ! flush(9) 372 372 ! 373 373 !-- ALLOCATE data structure 374 374 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 377 380 ! 378 381 !-- Read station coordinates and further attributes. … … 466 469 vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) ) 467 470 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 470 489 ! 471 490 !-- For the actual measurement ID read the UTM coordinates. Based on these, … … 502 521 !-- trajectory or station 503 522 ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) ) 504 ALLOCATE( vmea(l)%ngp(1:vmea(l)%ntraj) )505 523 ! 506 524 !-- Allocate temporary arrays for UTM and height coordinates. Note, … … 536 554 ENDIF 537 555 ! 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 ! 538 570 !-- Based on UTM coordinates, check if the measurement station or parts 539 571 !-- of the trajectory is on subdomain. This case, setup grid index space 540 572 !-- sample these quantities. 541 ns= 0573 meas_flag = 0 542 574 DO t = 1, vmea(l)%ntraj 543 575 ! … … 553 585 z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n 554 586 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 564 588 ! 565 589 !-- Compute grid indices relative to origin and check if these are … … 567 591 !-- at grid points surrounding the station, hence, check also for 568 592 !-- these grid points. 569 vmea(l)%ngp(t) = 0570 k_prev = -999571 j_prev = -999572 i_prev = -999573 593 DO n = 1, vmea(l)%dim_t(t) 574 594 is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp ) … … 579 599 js >= nys .AND. js <= nyn ) 580 600 ! 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 584 602 IF ( on_pe ) THEN 603 ! 604 !-- Determine vertical index which correspond to the observation 605 !-- height. 585 606 ksurf = get_topography_top_index_ji( js, is, 's' ) 586 607 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 ) 587 617 ENDIF 588 !589 !-- Count the number of observation points in index space on590 !-- subdomain. Only increment if grid indices are different from591 !-- the previous one.592 IF ( on_pe .AND. is /= i_prev .AND. js /= j_prev .AND. &593 ks /= k_prev ) THEN594 ns = ns + 1595 vmea(l)%ngp(t) = vmea(l)%ngp(t) + 1596 ENDIF597 598 !-- Store arrays for next iteration - avoid double counting599 i_prev = is600 j_prev = js601 k_prev = ks602 618 ENDDO 603 619 604 620 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 606 631 ! 607 632 !-- Store number of observation points on subdomain and allocate index 608 633 !-- arrays. 609 634 vmea(l)%ns = ns 635 ns = 0 610 636 611 637 ALLOCATE( vmea(l)%i(1:vmea(l)%ns) ) 612 638 ALLOCATE( vmea(l)%j(1:vmea(l)%ns) ) 613 639 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 614 655 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 663 663 ! 664 664 !-- Allocate array to save the sampled values. … … 677 677 ENDIF 678 678 ENDDO 679 flush(9)679 ! flush(9) 680 680 681 681 ! … … 684 684 CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '', & 685 685 global_attribute, 'close', '' ) 686 ! 687 !-- Dellocate flag array 688 DEALLOCATE( meas_flag ) 689 686 690 END SUBROUTINE vm_init 687 691 … … 694 698 SUBROUTINE vm_sampling 695 699 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 700 711 701 712 IMPLICIT NONE … … 703 714 CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name 704 715 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 711 724 ! 712 725 !-- Loop over all stations. For each possible variable loop over all … … 714 727 DO l = 1, nvm 715 728 ! 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 756 1041 ENDDO 757 1042
Note: See TracChangeset
for help on using the changeset viewer.
