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