/WPS/util/src/calc_ecmwf_p.F
FORTRAN Legacy | 446 lines | 237 code | 104 blank | 105 comment | 58 complexity | f10b4012f397ee3620d75f1454b4c7f6 MD5 | raw file
Possible License(s): AGPL-1.0
- module coefficients
- integer :: n_levels
- real, allocatable, dimension(:) :: a, b
- contains
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: read_coeffs
- !
- ! Notes: Obtain table of coefficients for input by this routine from
- ! http://www.ecmwf.int/products/data/technical/model_levels/index.html
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine read_coeffs()
-
- implicit none
- integer :: i, nlvl, istatus
- open(21,file='ecmwf_coeffs',form='formatted',status='old',iostat=istatus)
-
- n_levels = 0
- if (istatus /= 0) then
- write(6,*) 'ERROR: Error opening ecmwf_coeffs'
- return
- end if
- read(21,*,iostat=istatus) nlvl
- do while (istatus == 0)
- n_levels = n_levels + 1
- read(21,*,iostat=istatus) nlvl
- end do
-
- rewind(21)
- n_levels = n_levels - 1
- allocate(a(0:n_levels))
- allocate(b(0:n_levels))
- write(6,*) ' '
- write(6,*) 'Coefficients for each level:',n_levels
- do i=0,n_levels
- read(21,*,iostat=istatus) nlvl, a(i), b(i)
- write(6,'(i5,5x,f12.6,2x,f12.10)') nlvl, a(i), b(i)
- end do
- write(6,*) ' '
- close(21)
-
- end subroutine read_coeffs
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: cleanup_coeffs
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine cleanup_coeffs()
- implicit none
- n_levels = 0
- if (allocated(a)) deallocate(a)
- if (allocated(b)) deallocate(b)
- end subroutine cleanup_coeffs
- end module coefficients
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! calc_ecmwf_p
- !
- ! The purpose of this program is to compute a 3d pressure field for ECMWF
- ! model-level data sets; the code works in the WPS intermediate file format,
- ! reading a PSFC field from intermediate files, the A and B coefficients
- ! from a text file, ecmwf_coeffs, and writes the pressure data to an
- ! intermediate file.
- !
- ! November 2008
- ! Note: This program now also computes height for each level, this is needed by real
- ! modified by: Daniel van Dijke, MeteoConsult B.V., The Netherlands
- ! Chiel van Heerwaarden, Wageningen University, The Netherlands
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- program calc_ecmwf_p
- use coefficients
- use date_pack
- use gridinfo_module
- use misc_definitions_module
- use module_debug
- use read_met_module
- use stringutil
- use write_met_module
- implicit none
- ! Local variables
- integer :: i, idiff, n_times, t, istatus, fg_idx, counter
- real :: a_full, b_full
- character (len=19) :: valid_date, temp_date
- character (len=128) :: input_name
- real, allocatable, dimension(:,:) :: psfc, hgtsfc, hgtprev, pstart, pend
- real, allocatable, dimension(:,:,:) :: tt, qv, hgt_3ddata
- logical :: is_psfc
- type (met_data) :: ecmwf_data, p_data, rh_data, hgt_data
- !
- ! Setup (read namelist and check on time range)
- !
- call get_namelist_params()
- call set_debug_level(WARN)
- ! Compute number of times that we will process
- call geth_idts(end_date(1), start_date(1), idiff)
- call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=1)
- n_times = idiff / interval_seconds
- ! Check that the interval evenly divides the range of times to process
- call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
- 'In namelist, interval_seconds does not evenly divide '// &
- '(end_date - start_date) for domain %i. Only %i time periods '// &
- 'will be processed.', i1=1, i2=n_times)
- fg_idx = 1
- input_name = fg_name(fg_idx)
- !
- ! Get coefficients for model level pressures
- !
- call read_coeffs()
- !
- ! Loop over all prefixes listed in namelist for fg_name
- !
- do while (input_name /= '*')
- !
- ! Loop over all times to be processed for this domain
- !
- do t=0,n_times
- ! Get the date string for the current time
- call geth_newdate(valid_date, trim(start_date(1)), t*interval_seconds)
- temp_date = ' '
- write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
- ! Initialize the module for reading in the met fields
- call read_met_init(trim(input_name), .false., temp_date(1:13), istatus)
- is_psfc = .false.
- if (istatus == 0) then
- call mprintf(.true.,STDOUT,'Reading from %s at time %s', s1=input_name, s2=temp_date(1:13))
- ! Process all fields and levels from the current file; read_next_met_field()
- ! will return a non-zero status when there are no more fields to be read.
- do while (istatus == 0)
- call read_next_met_field(ecmwf_data, istatus)
- if (istatus == 0) then
- ! Have we found either PSFC or LOGSFP?
- if ((trim(ecmwf_data%field) == 'PSFC' .and. ecmwf_data%xlvl == 200100.) &
- .or. trim(ecmwf_data%field) == 'LOGSFP') then
- p_data = ecmwf_data
- p_data%field = 'PRESSURE '
- p_data%desc = 'Pressure'
- rh_data = ecmwf_data
- rh_data%field = 'RH '
- rh_data%units = '%'
- rh_data%desc = 'Relative humidity'
- if (.not. allocated(psfc)) then
- allocate(psfc(ecmwf_data%nx,ecmwf_data%ny))
- end if
- call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
- s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
- is_psfc = .true.
- if (trim(ecmwf_data%field) == 'LOGSFP') then
- psfc(:,:) = exp(ecmwf_data%slab(:,:))
- else
- psfc(:,:) = ecmwf_data%slab(:,:)
- end if
- !CvH_DvD: - Store geopotential height
- else if (trim(ecmwf_data%field) == 'SOILHGT' .and. ecmwf_data%xlvl == 200100.) then
- hgt_data = ecmwf_data
- hgt_data%field = 'GHT '
- hgt_data%units = 'm'
- hgt_data%desc = 'Height'
- if (.not. allocated(hgtsfc)) then
- allocate(hgtsfc(ecmwf_data%nx,ecmwf_data%ny))
- end if
- call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
- s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
- hgtsfc(:,:) = ecmwf_data%slab(:,:)
- ! Have we found surface geopotential?
- else if (trim(ecmwf_data%field) == 'SOILGEO' .and. ecmwf_data%xlvl == 1.) then
- hgt_data = ecmwf_data
- hgt_data%field = 'GHT '
- hgt_data%units = 'm' ! units on output after conversion to height
- hgt_data%desc = 'Height'
- if (.not. allocated(hgtsfc)) then
- allocate(hgtsfc(ecmwf_data%nx,ecmwf_data%ny))
- end if
- call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
- s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))
-
- hgtsfc(:,:) = ecmwf_data%slab(:,:)/9.81
- ! Have we found temperature?
- else if (trim(ecmwf_data%field) == 'TT') then
- if (.not. allocated(tt)) then
- allocate(tt(ecmwf_data%nx,ecmwf_data%ny,n_levels+1)) ! Extra level is for surface
- end if
- if (nint(ecmwf_data%xlvl) >= 1 .and. &
- nint(ecmwf_data%xlvl) <= n_levels) then
- tt(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
- else if (nint(ecmwf_data%xlvl) == 200100) then
- tt(:,:,n_levels+1) = ecmwf_data%slab
- end if
- ! Have we found specific humidity?
- else if (trim(ecmwf_data%field) == 'SPECHUMD') then
- if (.not. allocated(qv)) then
- allocate(qv(ecmwf_data%nx,ecmwf_data%ny,n_levels+1)) ! Extra level is for surface
- end if
- if (nint(ecmwf_data%xlvl) >= 1 .and. &
- nint(ecmwf_data%xlvl) <= n_levels) then
- qv(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
- else if (nint(ecmwf_data%xlvl) == 200100) then
- qv(:,:,n_levels+1) = ecmwf_data%slab
- end if
- end if
-
- if (associated(ecmwf_data%slab)) deallocate(ecmwf_data%slab)
-
- end if
-
- end do
-
- call read_met_close()
- ! Now write out, for each level, the pressure field
- if (is_psfc) then
- allocate(p_data%slab(p_data%nx,p_data%ny))
- allocate(rh_data%slab(rh_data%nx,rh_data%ny))
- !CvH_DvD: add HGT variable
- allocate(hgt_data%slab(hgt_data%nx,hgt_data%ny))
- call write_met_init(trim(get_path(input_name))//'PRES', .false., temp_date(1:13), istatus)
- if (allocated(tt) .and. allocated(qv)) then
- p_data%xlvl = 200100.
- p_data%slab = psfc
- ! Surface RH should be computed from surface DEWPT by ungrib
- ! rh_data%xlvl = 200100.
- ! call calc_rh(tt(:,:,n_levels+1), qv(:,:,n_levels+1), psfc, rh_data%slab, rh_data%nx, rh_data%ny)
- ! call write_next_met_field(rh_data, istatus)
- call write_next_met_field(p_data, istatus)
- else
- call mprintf(.true.,WARN,'Either TT or SPECHUMD not found. No RH will be computed.')
- end if
-
-
-
- !CvH_DvD: if tt, qv and hgtsfc are available compute hgt
- if (allocated(tt) .and. allocated(qv) .and. allocated(hgtsfc)) then
- if (.not. allocated(hgtprev)) then
- allocate(hgtprev(ecmwf_data%nx,ecmwf_data%ny))
- end if
- if (.not. allocated(pstart)) then
- allocate(pstart(ecmwf_data%nx,ecmwf_data%ny))
- end if
- if (.not. allocated(pend)) then
- allocate(pend(ecmwf_data%nx,ecmwf_data%ny))
- end if
- if (.not. allocated(hgt_3ddata)) then
- allocate(hgt_3ddata(ecmwf_data%nx,ecmwf_data%ny, n_levels))
- end if
-
- ! CvH_DvD: interpolate Q and T if they are not availalbe at all levels
- do i = n_levels, 1, -1
- if (tt(1,1,i) .eq. 0) then ! CvH_DvD: If TT is zero, level is unknown, so moisture is zero too
- if (i .ne. n_levels) then
- counter=1
- ! CvH_DvD: Find first available level
- do while ((tt(1,1,i-counter) .eq. 0) .and. (i-counter .gt. 1))
- counter=counter+1
- end do
- if (tt(1,1,i-counter) .gt. 0) then
- ! CvH_DvD: If TT is zero, we're at the top level which is missing, so use previous level,
- ! mod level will remove it anyway
- tt(:,:,i)=tt(:,:,i+1)+(tt(:,:,i-counter)-tt(:,:,i+1))/(counter+1)
- qv(:,:,i)=qv(:,:,i+1)+(qv(:,:,i-counter)-qv(:,:,i+1))/(counter+1)
- else
- tt(:,:,i)=tt(:,:,i+1)
- qv(:,:,i)=qv(:,:,i+1)
- end if
- else
- write(6,*) 'WARNING First level is missing!' ! CvH_DvD: First level should be there!
- end if
- end if
- end do
-
- ! CvH_DvD: first previous hgt is surface hgt/soilgeo
- hgtprev = hgtsfc
- ! CvH_DvD: - Loop from surface to top, note: 1=top & n_levels=surface
- do i = n_levels, 1, -1
- ! CvH_DvD: compute half level pressure at current half-level
- pend = 0.5 * (a(i-1) + a(i)) + psfc * 0.5 * (b(i-1) + b(i))
- if (i .eq. n_levels) then
- ! CvH_DvD: use Tv not T, use Psfc and T_halflevel_1=AveT lowest level
- hgt_3ddata(:,:,i) = hgtprev + 287.05 * tt(:,:,i) * (1. + 0.61 * qv(:,:,i)) * log(psfc/pend) / 9.81
- else
- ! CvH_DvD: compute half level pressure beneath current half-level
- pstart = 0.5 * (a(i+1) + a(i)) + psfc * 0.5 * (b(i+1) + b(i))
- ! CvH_DvD: use Tv not T, create T at full leve beneath current half level
- hgt_3ddata(:,:,i) = hgtprev + 287.05 * (tt(:,:,i) + tt(:,:,i+1))/2. * &
- (1. + 0.61 * (qv(:,:,i)+qv(:,:,i+1))/2.) * log(pstart/pend) / 9.81
- end if
- hgtprev = hgt_3ddata(:,:,i)
- end do
-
- end if
-
-
- do i = 1, n_levels
- a_full = 0.5 * (a(i-1) + a(i)) ! A and B are dimensioned (0:n_levels)
- b_full = 0.5 * (b(i-1) + b(i))
- p_data%xlvl = real(i)
- p_data%slab = a_full + psfc * b_full
- if (allocated(tt) .and. allocated(qv)) then
- rh_data%xlvl = real(i)
- call calc_rh(tt(:,:,i), qv(:,:,i), p_data%slab, rh_data%slab, rh_data%nx, rh_data%ny)
- call write_next_met_field(rh_data, istatus)
-
- if (allocated(hgtsfc)) then
- ! CvH_DvD: put hgt_3ddata into hgt_data object
- hgt_data%xlvl = real(i)
- hgt_data%slab = hgt_3ddata(:,:,i)
- call write_next_met_field(hgt_data, istatus)
- else
- call mprintf(.true., WARN, &
- 'Either SOILHGT or SOILGEO are required to create 3-d GHT field, which is required '// &
- 'for a correct vertical interpolation in real.')
- end if
- end if
- call write_next_met_field(p_data, istatus)
- end do
- call write_met_close()
-
- deallocate(p_data%slab)
- deallocate(rh_data%slab)
- ! CvH_DvD: deallocate stuff
- deallocate(hgt_data%slab)
- end if
- if (allocated(tt)) deallocate(tt)
- if (allocated(qv)) deallocate(qv)
- if (allocated(psfc)) deallocate(psfc)
- ! CvH_DvD: deallocate stuff
- if (allocated(hgt_3ddata)) deallocate(hgt_3ddata)
- if (allocated(pstart)) deallocate(pstart)
- if (allocated(pend)) deallocate(pend)
- if (allocated(hgtprev)) deallocate(hgtprev)
-
- else
- call mprintf(.true.,ERROR,'Problem opening %s at time %s', s1=input_name, s2=temp_date(1:13))
- end if
- end do
- ! CvH_DvD: only now deallocate hgtsfc,
- ! because oper EC does not have hgtsfc @ fps
- if (allocated(hgtsfc)) deallocate(hgtsfc)
- fg_idx = fg_idx + 1
- input_name = fg_name(fg_idx)
- end do
- call cleanup_coeffs()
- stop
- end program calc_ecmwf_p
- subroutine calc_rh(t, qv, p, rh, nx, ny)
- implicit none
- ! Arguments
- integer, intent(in) :: nx, ny
- real, dimension(nx, ny), intent(in) :: t, qv, p
- real, dimension(nx, ny), intent(out) :: rh
- ! Constants
- real, parameter :: svp1=0.6112
- real, parameter :: svp2=17.67
- real, parameter :: svp3=29.65
- real, parameter :: svpt0=273.15
- real, parameter :: eps = 0.622
- ! Local variables
- integer :: i, j
- real :: es, e
- do j=1,ny
- do i=1,nx
- es=svp1*10.*exp(svp2*(t(i,j)-svpt0)/(t(i,j)-svp3))
- e=qv(i,j)*p(i,j)/100./(eps+qv(i,j)*(1.-eps)) ! qv is specific humidity
- rh(i,j) = 100.0 * e/es
- end do
- end do
- end subroutine calc_rh