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                      end if
 239
 240                      if (nint(ecmwf_data%xlvl) >= 1 .and. &
 241                          nint(ecmwf_data%xlvl) <= n_levels) then
 242                         tt(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
 243                      else if (nint(ecmwf_data%xlvl) == 200100) then
 244                         tt(:,:,n_levels+1) = ecmwf_data%slab
 245                      end if
 246
 247                   ! Have we found specific humidity?
 248                   else if (trim(ecmwf_data%field) == 'SPECHUMD') then
 249
 250                      if (.not. allocated(qv)) then
 251                         allocate(qv(ecmwf_data%nx,ecmwf_data%ny,n_levels+1))  ! Extra level is for surface
 252                      end if
 253
 254                      if (nint(ecmwf_data%xlvl) >= 1 .and. &
 255                          nint(ecmwf_data%xlvl) <= n_levels) then
 256                         qv(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
 257                      else if (nint(ecmwf_data%xlvl) == 200100) then
 258                         qv(:,:,n_levels+1) = ecmwf_data%slab
 259                      end if
 260
 261                   end if
 262
 263                   if (associated(ecmwf_data%slab)) deallocate(ecmwf_data%slab)
 264
 265                end if
 266
 267             end do
 268
 269             call read_met_close()
 270
 271
 272             ! Now write out, for each level, the pressure field
 273             if (is_psfc) then
 274
 275                allocate(p_data%slab(p_data%nx,p_data%ny))
 276                allocate(rh_data%slab(rh_data%nx,rh_data%ny))
 277                !CvH_DvD: add HGT variable
 278                allocate(hgt_data%slab(hgt_data%nx,hgt_data%ny))
 279
 280                call write_met_init(trim(get_path(input_name))//'PRES', .false., temp_date(1:13), istatus)
 281
 282                if (allocated(tt) .and. allocated(qv)) then
 283                   p_data%xlvl  = 200100.
 284                   p_data%slab  = psfc
 285 ! Surface RH should be computed from surface DEWPT by ungrib
 286 !                  rh_data%xlvl = 200100.
 287 !                  call calc_rh(tt(:,:,n_levels+1), qv(:,:,n_levels+1), psfc, rh_data%slab, rh_data%nx, rh_data%ny)
 288 !                  call write_next_met_field(rh_data, istatus)
 289                   call write_next_met_field(p_data, istatus)
 290                else
 291                   call mprintf(.true.,WARN,'Either TT or SPECHUMD not found. No RH will be computed.')
 292                end if
 293
 294
 295
 296                !CvH_DvD: if tt, qv and hgtsfc are available compute hgt
 297                if (allocated(tt) .and. allocated(qv) .and. allocated(hgtsfc)) then
 298
 299                   if (.not. allocated(hgtprev)) then
 300                     allocate(hgtprev(ecmwf_data%nx,ecmwf_data%ny))
 301                   end if
 302                   if (.not. allocated(pstart)) then
 303                     allocate(pstart(ecmwf_data%nx,ecmwf_data%ny))
 304                   end if
 305                   if (.not. allocated(pend)) then
 306                     allocate(pend(ecmwf_data%nx,ecmwf_data%ny))
 307                   end if
 308                   if (.not. allocated(hgt_3ddata)) then
 309                     allocate(hgt_3ddata(ecmwf_data%nx,ecmwf_data%ny, n_levels))
 310                   end if
 311
 312                   ! CvH_DvD: interpolate Q and T if they are not available at all levels
 313                   do i = n_levels, 1, -1
 314                     if (tt(1,1,i) .eq. 0) then  ! CvH_DvD: If TT is zero, level is unknown, so moisture is zero too
 315                       if (i .eq. n_levels) then
 316                         write(6,*) 'WARNING First level is missing!' ! CvH_DvD: First level should be there!
 317                       else if (i .eq. 1) then
 318                         ! DvD: If TT is zero at the top level, so missing --> use previous level,
 319                         ! mod level will remove it anyway
 320                         tt(:,:,i)=tt(:,:,i+1)
 321                         qv(:,:,i)=qv(:,:,i+1)
 322                       else
 323                         counter=1
 324                         ! CvH_DvD: Find first available level
 325                         do while ((i-counter .gt. 1) .and. (tt(1,1,i-counter) .eq. 0))
 326                           counter=counter+1
 327                         end do
 328                         if (tt(1,1,i-counter) .gt. 0) then
 329                           ! DvD: Interpolate tt and qv from next available level
 330                           tt(:,:,i)=tt(:,:,i+1)+(tt(:,:,i-counter)-tt(:,:,i+1))/(counter+1)
 331                           qv(:,:,i)=qv(:,:,i+1)+(qv(:,:,i-counter)-qv(:,:,i+1))/(counter+1)
 332                         else
 333                           ! DvD: No available level found --> should not happen.
 334                           write(6,*) 'WARNING No available level found near level ',i,'!'
 335                         end if
 336                       end if
 337                     end if
 338                   end do
 339
 340                   ! CvH_DvD: first previous hgt is surface hgt/soilgeo
 341                   hgtprev = hgtsfc
 342                   ! CvH_DvD: - Loop from surface to top, note: 1=top & n_levels=surface
 343                   do i = n_levels, 1, -1
 344                       ! CvH_DvD: compute half level pressure at current half-level
 345                       pend   = 0.5 * (a(i-1) + a(i)) + psfc * 0.5 * (b(i-1) + b(i))
 346                       if (i .eq. n_levels) then
 347                         ! CvH_DvD: use Tv not T, use Psfc and T_halflevel_1=AveT lowest level
 348                         hgt_3ddata(:,:,i) = hgtprev + 287.05 * tt(:,:,i) * (1. + 0.61 * qv(:,:,i)) * log(psfc/pend) / 9.81
 349                       else
 350                         ! CvH_DvD: compute half level pressure beneath current half-level
 351                         pstart = 0.5 * (a(i+1) + a(i)) + psfc * 0.5 * (b(i+1) + b(i))
 352                         ! CvH_DvD: use Tv not T, create T at full leve beneath current half level
 353                         hgt_3ddata(:,:,i) = hgtprev + 287.05 * (tt(:,:,i) + tt(:,:,i+1))/2. * &
 354                                                      (1. + 0.61 * (qv(:,:,i)+qv(:,:,i+1))/2.) * log(pstart/pend) / 9.81
 355                       end if
 356                       hgtprev = hgt_3ddata(:,:,i)
 357                   end do
 358
 359                end if
 360
 361
 362                do i = 1, n_levels
 363
 364                   a_full = 0.5 * (a(i-1) + a(i))   ! A and B are dimensioned (0:n_levels)
 365                   b_full = 0.5 * (b(i-1) + b(i))
 366
 367                   p_data%xlvl = real(i)
 368                   p_data%slab = a_full + psfc * b_full
 369
 370                   if (allocated(tt) .and. allocated(qv)) then
 371                      rh_data%xlvl = real(i)
 372                      call calc_rh(tt(:,:,i), qv(:,:,i), p_data%slab, rh_data%slab, rh_data%nx, rh_data%ny)
 373                      call write_next_met_field(rh_data, istatus)
 374
 375                      if (allocated(hgtsfc)) then
 376                         ! CvH_DvD: put hgt_3ddata into hgt_data object
 377                         hgt_data%xlvl = real(i)
 378                         hgt_data%slab = hgt_3ddata(:,:,i)
 379                         call write_next_met_field(hgt_data, istatus)
 380                      else
 381                         call mprintf(.true., WARN, &
 382                                      'Either SOILHGT or SOILGEO are required to create 3-d GHT field, which is required '// &
 383                                      'for a correct vertical interpolation in real.')
 384                      end if
 385                   end if
 386
 387                   call write_next_met_field(p_data, istatus)
 388
 389                end do
 390
 391                call write_met_close()
 392
 393                deallocate(p_data%slab)
 394                deallocate(rh_data%slab)
 395                ! CvH_DvD: deallocate stuff
 396                deallocate(hgt_data%slab)
 397
 398             end if
 399
 400             if (allocated(tt))   deallocate(tt)
 401             if (allocated(qv))   deallocate(qv)
 402             if (allocated(psfc)) deallocate(psfc)
 403             ! CvH_DvD: deallocate stuff
 404             if (allocated(hgt_3ddata)) deallocate(hgt_3ddata)
 405             if (allocated(pstart)) deallocate(pstart)
 406             if (allocated(pend)) deallocate(pend)
 407             if (allocated(hgtprev)) deallocate(hgtprev)
 408
 409
 410          else
 411             call mprintf(.true.,ERROR,'Problem opening %s at time %s', s1=input_name, s2=temp_date(1:13))
 412          end if
 413
 414       end do
 415       ! CvH_DvD: only now deallocate hgtsfc,
 416       ! because oper EC does not have hgtsfc @ fps
 417       if (allocated(hgtsfc)) deallocate(hgtsfc)
 418       fg_idx = fg_idx + 1
 419       input_name = fg_name(fg_idx)
 420
 421    end do
 422
 423    call cleanup_coeffs()
 424
 425    stop
 426
 427 end program calc_ecmwf_p
 428
 429
 430 subroutine calc_rh(t, qv, p, rh, nx, ny)
 431
 432    implicit none
 433
 434    ! Arguments
 435    integer, intent(in) :: nx, ny
 436    real, dimension(nx, ny), intent(in)  :: t, qv, p
 437    real, dimension(nx, ny), intent(out) :: rh
 438
 439    ! Constants
 440    real, parameter :: svp1=0.6112
 441    real, parameter :: svp2=17.67
 442    real, parameter :: svp3=29.65
 443    real, parameter :: svpt0=273.15
 444    real, parameter :: eps = 0.622
 445
 446    ! Local variables
 447    integer :: i, j
 448    real :: es, e
 449
 450    do j=1,ny
 451    do i=1,nx
 452       es=svp1*10.*exp(svp2*(t(i,j)-svpt0)/(t(i,j)-svp3))
 453       e=qv(i,j)*p(i,j)/100./(eps+qv(i,j)*(1.-eps)) ! qv is specific humidity
 454       rh(i,j) = 100.0 * e/es
 455    end do
 456    end do
 457
 458 end subroutine calc_rh