util/src/calc_ecmwf_p.F

   1 module coefficients
   2
   3    integer :: n_levels
   4    real, allocatable, dimension(:) :: a, b
   5
   6    contains
   7
   8    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   9    ! Name: read_coeffs
  10    !
  11    ! Notes: Obtain table of coefficients for input by this routine from the link
  12    !        below that corresponds to the correct number of levels:
  13    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/16-model-levels
  14    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/19-model-levels
  15    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/31-model-levels
  16    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/40-model-levels
  17    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/50-model-levels
  18    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/60-model-levels
  19    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/62-model-levels
  20    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/91-model-levels
  21    !        http://www.ecmwf.int/en/forecasts/documentation-and-support/137-model-levels
  22    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  23    subroutine read_coeffs()
  24
  25       implicit none
  26
  27       integer :: i, nlvl, istatus
  28
  29       open(21,file='ecmwf_coeffs',form='formatted',status='old',iostat=istatus)
  30
  31       n_levels = 0
  32
  33       if (istatus /= 0) then
  34          write(6,*) 'ERROR: Error opening ecmwf_coeffs'
  35          return
  36       end if
  37
  38       read(21,*,iostat=istatus) nlvl
  39       do while (istatus == 0)
  40          n_levels = n_levels + 1
  41          read(21,*,iostat=istatus) nlvl
  42       end do
  43
  44       rewind(21)
  45
  46       n_levels = n_levels - 1
  47
  48       allocate(a(0:n_levels))
  49       allocate(b(0:n_levels))
  50
  51       write(6,*) ' '
  52       write(6,*) 'Coefficients for each level:',n_levels
  53       do i=0,n_levels
  54          read(21,*,iostat=istatus) nlvl, a(i), b(i)
  55          write(6,'(i5,5x,f12.6,2x,f12.10)') nlvl, a(i), b(i)
  56       end do
  57       write(6,*) ' '
  58
  59       close(21)
  60
  61    end subroutine read_coeffs
  62
  63
  64    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  65    ! Name: cleanup_coeffs
  66    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  67    subroutine cleanup_coeffs()
  68
  69       implicit none
  70
  71       n_levels = 0
  72       if (allocated(a)) deallocate(a)
  73       if (allocated(b)) deallocate(b)
  74
  75    end subroutine cleanup_coeffs
  76
  77 end module coefficients
  78
  79
  80 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  81 ! calc_ecmwf_p
  82 !
  83 ! The purpose of this program is to compute a 3d pressure field for ECMWF
  84 !    model-level data sets; the code works in the WPS intermediate file format,
  85 !    reading a PSFC field from intermediate files, the A and B coefficients
  86 !    from a text file, ecmwf_coeffs, and writes the pressure data to an
  87 !    intermediate file.
  88 !
  89 ! November 2008
  90 ! Note: This program now also computes height for each level, this is needed by real
  91 !       modified by: Daniel van Dijke, MeteoConsult B.V., The Netherlands
  92 !                    Chiel van Heerwaarden, Wageningen University, The Netherlands
  93 !
  94 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  95 program calc_ecmwf_p
  96
  97    use coefficients
  98    use date_pack
  99    use gridinfo_module
 100    use misc_definitions_module
 101    use module_debug
 102    use read_met_module
 103    use stringutil
 104    use write_met_module
 105
 106    implicit none
 107
 108    ! Local variables
 109    integer :: i, idiff, n_times, t, istatus, fg_idx, counter
 110    real :: a_full, b_full
 111    character (len=19) :: valid_date, temp_date
 112    character (len=128) :: input_name
 113    real, allocatable, dimension(:,:) :: psfc, hgtsfc, hgtprev, pstart, pend
 114    real, allocatable, dimension(:,:,:) :: tt, qv, hgt_3ddata
 115    logical :: is_psfc
 116    type (met_data) :: ecmwf_data, p_data, rh_data, hgt_data
 117
 118
 119    !
 120    ! Setup (read namelist and check on time range)
 121    !
 122    call get_namelist_params()
 123
 124    call set_debug_level(WARN)
 125
 126    ! Compute number of times that we will process
 127    call geth_idts(end_date(1), start_date(1), idiff)
 128    call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=1)
 129
 130    n_times = idiff / interval_seconds
 131
 132    ! Check that the interval evenly divides the range of times to process
 133    call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
 134                 'In namelist, interval_seconds does not evenly divide '// &
 135                 '(end_date - start_date) for domain %i. Only %i time periods '// &
 136                 'will be processed.', i1=1, i2=n_times)
 137
 138    fg_idx = 1
 139
 140    input_name = fg_name(fg_idx)
 141
 142    !
 143    ! Get coefficients for model level pressures
 144    !
 145    call read_coeffs()
 146
 147
 148    !
 149    ! Loop over all prefixes listed in namelist for fg_name
 150    !
 151    do while (input_name /= '*')
 152
 153       !
 154       ! Loop over all times to be processed for this domain
 155       !
 156       do t=0,n_times
 157
 158          ! Get the date string for the current time
 159          call geth_newdate(valid_date, trim(start_date(1)), t*interval_seconds)
 160          temp_date = ' '
 161          write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
 162
 163          ! Initialize the module for reading in the met fields
 164          call read_met_init(trim(input_name), .false., temp_date(1:13), istatus)
 165
 166          is_psfc = .false.
 167
 168          if (istatus == 0) then
 169             call mprintf(.true.,STDOUT,'Reading from %s at time %s', s1=input_name, s2=temp_date(1:13))
 170
 171             ! Process all fields and levels from the current file; read_next_met_field()
 172             !   will return a non-zero status when there are no more fields to be read.
 173             do while (istatus == 0)
 174
 175                call read_next_met_field(ecmwf_data, istatus)
 176
 177                if (istatus == 0) then
 178
 179                   ! Have we found either PSFC or LOGSFP?
 180                   if ((trim(ecmwf_data%field) == 'PSFC' .and. ecmwf_data%xlvl == 200100.) &
 181                       .or. trim(ecmwf_data%field) == 'LOGSFP') then
 182                      p_data = ecmwf_data
 183                      p_data%field = 'PRESSURE '
 184                      p_data%desc  = 'Pressure'
 185
 186                      rh_data = ecmwf_data
 187                      rh_data%field = 'RH       '
 188                      rh_data%units = '%'
 189                      rh_data%desc  = 'Relative humidity'
 190
 191                      if (.not. allocated(psfc)) then
 192                         allocate(psfc(ecmwf_data%nx,ecmwf_data%ny))
 193                      end if
 194                      call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
 195                                   s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
 196
 197                      is_psfc   = .true.
 198                      if (trim(ecmwf_data%field) == 'LOGSFP') then
 199                         psfc(:,:) = exp(ecmwf_data%slab(:,:))
 200                      else
 201                         psfc(:,:) = ecmwf_data%slab(:,:)
 202                      end if
 203
 204                   !CvH_DvD: - Store geopotential height
 205                   else if (trim(ecmwf_data%field) == 'SOILHGT' .and. ecmwf_data%xlvl == 200100.) then
 206                      hgt_data = ecmwf_data
 207                      hgt_data%field = 'GHT      '
 208                      hgt_data%units = 'm'
 209                      hgt_data%desc  = 'Height'
 210                      if (.not. allocated(hgtsfc)) then
 211                         allocate(hgtsfc(ecmwf_data%nx,ecmwf_data%ny))
 212                      end if
 213                      call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
 214                                   s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
 215
 216                      hgtsfc(:,:) = ecmwf_data%slab(:,:)
 217
 218                   ! Have we found surface geopotential?
 219                   else if (trim(ecmwf_data%field) == 'SOILGEO' .and. ecmwf_data%xlvl == 1.) then
 220                      hgt_data = ecmwf_data
 221                      hgt_data%field = 'GHT      '
 222                      hgt_data%units = 'm'      ! units on output after conversion to height
 223                      hgt_data%desc  = 'Height'
 224                      if (.not. allocated(hgtsfc)) then
 225                         allocate(hgtsfc(ecmwf_data%nx,ecmwf_data%ny))
 226                      end if
 227                      call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
 228                                   s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
 229
 230                      hgtsfc(:,:) = ecmwf_data%slab(:,:)/9.81
 231
 232
 233                   ! Have we found temperature?
 234                   else if (trim(ecmwf_data%field) == 'TT') then
 235
 236                      if (.not. allocated(tt)) then
 237                         allocate(tt(ecmwf_data%nx,ecmwf_data%ny,n_levels+1))  ! Extra level is for surface
 238                         tt(:,:,:) = 0.0
 239                      end if
 240
 241                      if (nint(ecmwf_data%xlvl) >= 1 .and. &
 242                          nint(ecmwf_data%xlvl) <= n_levels) then
 243                         tt(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
 244                      else if (nint(ecmwf_data%xlvl) == 200100) then
 245                         tt(:,:,n_levels+1) = ecmwf_data%slab
 246                      end if
 247
 248                   ! Have we found specific humidity?
 249                   else if (trim(ecmwf_data%field) == 'SPECHUMD') then
 250
 251                      if (.not. allocated(qv)) then
 252                         allocate(qv(ecmwf_data%nx,ecmwf_data%ny,n_levels+1))  ! Extra level is for surface
 253                         qv(:,:,:) = 0.0
 254                      end if
 255
 256                      if (nint(ecmwf_data%xlvl) >= 1 .and. &
 257                          nint(ecmwf_data%xlvl) <= n_levels) then
 258                         qv(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
 259                      else if (nint(ecmwf_data%xlvl) == 200100) then
 260                         qv(:,:,n_levels+1) = ecmwf_data%slab
 261                      end if
 262
 263                   end if
 264
 265                   if (associated(ecmwf_data%slab)) deallocate(ecmwf_data%slab)
 266
 267                end if
 268
 269             end do
 270
 271             call read_met_close()
 272
 273
 274             ! Now write out, for each level, the pressure field
 275             if (is_psfc) then
 276
 277                allocate(p_data%slab(p_data%nx,p_data%ny))
 278                allocate(rh_data%slab(rh_data%nx,rh_data%ny))
 279                !CvH_DvD: add HGT variable
 280                allocate(hgt_data%slab(hgt_data%nx,hgt_data%ny))
 281
 282                call write_met_init(trim(get_path(input_name))//'PRES', .false., temp_date(1:13), istatus)
 283
 284                if (allocated(tt) .and. allocated(qv)) then
 285                   p_data%xlvl  = 200100.
 286                   p_data%slab  = psfc
 287 ! Surface RH should be computed from surface DEWPT by ungrib
 288 !                  rh_data%xlvl = 200100.
 289 !                  call calc_rh(tt(:,:,n_levels+1), qv(:,:,n_levels+1), psfc, rh_data%slab, rh_data%nx, rh_data%ny)
 290 !                  call write_next_met_field(rh_data, istatus)
 291                   call write_next_met_field(p_data, istatus)
 292                else
 293                   call mprintf(.true.,WARN,'Either TT or SPECHUMD not found. No RH will be computed.')
 294                end if
 295
 296
 297
 298                !CvH_DvD: if tt, qv and hgtsfc are available compute hgt
 299                if (allocated(tt) .and. allocated(qv) .and. allocated(hgtsfc)) then
 300
 301                   if (.not. allocated(hgtprev)) then
 302                     allocate(hgtprev(ecmwf_data%nx,ecmwf_data%ny))
 303                   end if
 304                   if (.not. allocated(pstart)) then
 305                     allocate(pstart(ecmwf_data%nx,ecmwf_data%ny))
 306                   end if
 307                   if (.not. allocated(pend)) then
 308                     allocate(pend(ecmwf_data%nx,ecmwf_data%ny))
 309                   end if
 310                   if (.not. allocated(hgt_3ddata)) then
 311                     allocate(hgt_3ddata(ecmwf_data%nx,ecmwf_data%ny, n_levels))
 312                   end if
 313
 314                   ! CvH_DvD: interpolate Q and T if they are not available at all levels
 315                   do i = n_levels, 1, -1
 316                     if (tt(1,1,i) .eq. 0) then  ! CvH_DvD: If TT is zero, level is unknown, so moisture is zero too
 317                       if (i .eq. n_levels) then
 318                         write(6,*) 'WARNING First level is missing!' ! CvH_DvD: First level should be there!
 319                       else if (i .eq. 1) then
 320                         ! DvD: If TT is zero at the top level, so missing --> use previous level,
 321                         ! mod level will remove it anyway
 322                         tt(:,:,i)=tt(:,:,i+1)
 323                         qv(:,:,i)=qv(:,:,i+1)
 324                       else
 325                         counter=1
 326                         ! CvH_DvD: Find first available level
 327                         do while ((i-counter .gt. 1) .and. (tt(1,1,i-counter) .eq. 0))
 328                           counter=counter+1
 329                         end do
 330                         if (tt(1,1,i-counter) .gt. 0) then
 331                           ! DvD: Interpolate tt and qv from next available level
 332                           tt(:,:,i)=tt(:,:,i+1)+(tt(:,:,i-counter)-tt(:,:,i+1))/(counter+1)
 333                           qv(:,:,i)=qv(:,:,i+1)+(qv(:,:,i-counter)-qv(:,:,i+1))/(counter+1)
 334                         else
 335                           ! DvD: No available level found --> should not happen.
 336                           write(6,*) 'WARNING No available level found near level ',i,'!'
 337                         end if
 338                       end if
 339                     end if
 340                   end do
 341
 342                   ! CvH_DvD: first previous hgt is surface hgt/soilgeo
 343                   hgtprev = hgtsfc
 344                   ! CvH_DvD: - Loop from surface to top, note: 1=top & n_levels=surface
 345                   do i = n_levels, 1, -1
 346                       ! CvH_DvD: compute half level pressure at current half-level
 347                       pend   = 0.5 * (a(i-1) + a(i)) + psfc * 0.5 * (b(i-1) + b(i))
 348                       if (i .eq. n_levels) then
 349                         ! CvH_DvD: use Tv not T, use Psfc and T_halflevel_1=AveT lowest level
 350                         hgt_3ddata(:,:,i) = hgtprev + 287.05 * tt(:,:,i) * (1. + 0.61 * qv(:,:,i)) * log(psfc/pend) / 9.81
 351                       else
 352                         ! CvH_DvD: compute half level pressure beneath current half-level
 353                         pstart = 0.5 * (a(i+1) + a(i)) + psfc * 0.5 * (b(i+1) + b(i))
 354                         ! CvH_DvD: use Tv not T, create T at full leve beneath current half level
 355                         hgt_3ddata(:,:,i) = hgtprev + 287.05 * (tt(:,:,i) + tt(:,:,i+1))/2. * &
 356                                                      (1. + 0.61 * (qv(:,:,i)+qv(:,:,i+1))/2.) * log(pstart/pend) / 9.81
 357                       end if
 358                       hgtprev = hgt_3ddata(:,:,i)
 359                   end do
 360
 361                end if
 362
 363
 364                do i = 1, n_levels
 365
 366                   a_full = 0.5 * (a(i-1) + a(i))   ! A and B are dimensioned (0:n_levels)
 367                   b_full = 0.5 * (b(i-1) + b(i))
 368
 369                   p_data%xlvl = real(i)
 370                   p_data%slab = a_full + psfc * b_full
 371
 372                   if (allocated(tt) .and. allocated(qv)) then
 373                      rh_data%xlvl = real(i)
 374                      call calc_rh(tt(:,:,i), qv(:,:,i), p_data%slab, rh_data%slab, rh_data%nx, rh_data%ny)
 375                      call write_next_met_field(rh_data, istatus)
 376
 377                      if (allocated(hgtsfc)) then
 378                         ! CvH_DvD: put hgt_3ddata into hgt_data object
 379                         hgt_data%xlvl = real(i)
 380                         hgt_data%slab = hgt_3ddata(:,:,i)
 381                         call write_next_met_field(hgt_data, istatus)
 382                      else
 383                         call mprintf(.true., WARN, &
 384                                      'Either SOILHGT or SOILGEO are required to create 3-d GHT field, which is required '// &
 385                                      'for a correct vertical interpolation in real.')
 386                      end if
 387                   end if
 388
 389                   call write_next_met_field(p_data, istatus)
 390
 391                end do
 392
 393                call write_met_close()
 394
 395                deallocate(p_data%slab)
 396                deallocate(rh_data%slab)
 397                ! CvH_DvD: deallocate stuff
 398                deallocate(hgt_data%slab)
 399
 400             end if
 401
 402             if (allocated(tt))   deallocate(tt)
 403             if (allocated(qv))   deallocate(qv)
 404             if (allocated(psfc)) deallocate(psfc)
 405             ! CvH_DvD: deallocate stuff
 406             if (allocated(hgt_3ddata)) deallocate(hgt_3ddata)
 407             if (allocated(pstart)) deallocate(pstart)
 408             if (allocated(pend)) deallocate(pend)
 409             if (allocated(hgtprev)) deallocate(hgtprev)
 410
 411
 412          else
 413             call mprintf(.true.,ERROR,'Problem opening %s at time %s', s1=input_name, s2=temp_date(1:13))
 414          end if
 415
 416       end do
 417       ! CvH_DvD: only now deallocate hgtsfc,
 418       ! because oper EC does not have hgtsfc @ fps
 419       if (allocated(hgtsfc)) deallocate(hgtsfc)
 420       fg_idx = fg_idx + 1
 421       input_name = fg_name(fg_idx)
 422
 423    end do
 424
 425    call cleanup_coeffs()
 426
 427    stop
 428
 429 end program calc_ecmwf_p
 430
 431
 432 subroutine calc_rh(t, qv, p, rh, nx, ny)
 433
 434    implicit none
 435
 436    ! Arguments
 437    integer, intent(in) :: nx, ny
 438    real, dimension(nx, ny), intent(in)  :: t, qv, p
 439    real, dimension(nx, ny), intent(out) :: rh
 440
 441    ! Constants
 442    real, parameter :: svp1=0.6112
 443    real, parameter :: svp2=17.67
 444    real, parameter :: svp3=29.65
 445    real, parameter :: svpt0=273.15
 446    real, parameter :: eps = 0.622
 447
 448    ! Local variables
 449    integer :: i, j
 450    real :: es, e
 451
 452    do j=1,ny
 453    do i=1,nx
 454       es=svp1*10.*exp(svp2*(t(i,j)-svpt0)/(t(i,j)-svp3))
 455       e=qv(i,j)*p(i,j)/100./(eps+qv(i,j)*(1.-eps)) ! qv is specific humidity
 456       rh(i,j) = 100.0 * e/es
 457    end do
 458    end do
 459
 460 end subroutine calc_rh