util/src/height_ukmo.F

   1 program height_ukmo
   2
   3    !  This program computes the 3d height field for the UKMO data.  The heights
   4    !  are constant in time, and are a function only of the topography and a few
   5    !  constants.  There are several UKMO data sets available, denoted by the
   6    !  different numbers of vertical levels.  An input file defines the required
   7    !  information for each of these data sets.
   8
   9    use date_pack
  10    use gridinfo_module
  11    use read_met_module
  12    use write_met_module
  13    use misc_definitions_module
  14    use module_debug
  15
  16    implicit none
  17
  18    ! Local variables
  19    integer :: t, istatus
  20    character (len=19) :: valid_date, temp_date, output_date
  21    character (len=128) :: input_name
  22    type (met_data) :: temp_data, soilhgt_data, height_data
  23
  24    integer :: model_levels , first_constant_r_rho_level , k_loop , j_loop , i_loop , temp_count , date_loop
  25    real :: z_top_of_model , etac
  26    real, dimension(:) , allocatable :: eta
  27    real , dimension(:,:) , allocatable :: soil_hold
  28    logical :: still_more_dates_to_process
  29
  30    call get_namelist_params()
  31
  32    call set_debug_level(WARN)
  33
  34    input_name = fg_name(1)
  35
  36    call geth_newdate(valid_date, trim(start_date(1)), 0 )
  37    temp_date = ' '
  38    output_date = ' '
  39    write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
  40    write(output_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
  41
  42    !  Do this twice.  The first time is to find out how many levels of
  43    !  3d fields there are.  Why?  The data from UKMO has *some* of the
  44    !  levels, just not all.  So we can compute the height field, but we
  45    !  are not allowed to just spit out 70 levels.
  46
  47    ! Initialize the module for reading in the met fields
  48    call read_met_init(trim(input_name), .false., temp_date(1:13), istatus)
  49
  50    if (istatus == 0) then
  51       call mprintf(.true.,STDOUT,'Reading from %s at time %s', s1=input_name, s2=temp_date(1:13))
  52
  53       ! Process all fields and levels from the current file; read_next_met_field()
  54       !   will return a non-zero status when there are no more fields to be read.
  55       temp_count = 0
  56       just_temp : do while (istatus == 0)
  57
  58          call read_next_met_field(temp_data, istatus)
  59
  60          if (istatus == 0) then
  61
  62             if (trim(temp_data%field) == 'TT' .and. temp_data%xlvl /= 200100.) then
  63                temp_count = temp_count + 1
  64             else if (trim(temp_data%field) == 'SOILHGT' .and. temp_data%xlvl == 200100.) then
  65                allocate ( soil_hold(temp_data%nx,temp_data%ny) )
  66                do j_loop = 1 , temp_data%ny
  67                   do i_loop = 1 , temp_data%nx
  68                      soil_hold(i_loop,j_loop) = temp_data%slab(i_loop,j_loop)
  69                   end do
  70                end do
  71             end if
  72
  73             deallocate ( temp_data%slab )
  74          end if
  75
  76       end do just_temp
  77       call mprintf(.true.,STDOUT,'Total number of TT levels, excluding 200100 = %i' , i1=temp_count)
  78
  79       call read_met_close()
  80
  81    else
  82       call mprintf(.true.,ERROR,'Problem opening %s at time %s', s1=input_name, s2=temp_date(1:13))
  83    end if
  84
  85    !  Now we are doing this the second time.  We know how many temp levels there are, temp_count,
  86    !  so we will use the same number for height.
  87
  88    ! Initialize the module for reading in the met fields - again!
  89    call read_met_init(trim(input_name), .false., temp_date(1:13), istatus)
  90
  91    if (istatus == 0) then
  92       call mprintf(.true.,STDOUT,'Reading from %s at time %s', s1=input_name, s2=temp_date(1:13))
  93
  94       ! Process all fields and levels from the current file; read_next_met_field()
  95       !   will return a non-zero status when there are no more fields to be read.
  96       all_fields : do while (istatus == 0)
  97
  98
  99          if ( associated ( soilhgt_data%slab ) ) then
 100             deallocate ( soilhgt_data%slab )
 101          end if
 102          call read_next_met_field(soilhgt_data, istatus)
 103
 104          if (istatus == 0) then
 105
 106             if (trim(soilhgt_data%field) == 'SOILHGT  ' .and. soilhgt_data%xlvl == 200100.) then
 107                if (.not. associated(height_data%slab)) then
 108                   allocate(height_data%slab(soilhgt_data%nx,soilhgt_data%ny))
 109                end if
 110                height_data = soilhgt_data
 111                height_data%hdate = '0000-00-00_00:00:00     '
 112                height_data%xfcst = 0.
 113                height_data%field = 'HGT      '
 114                height_data%slab = 0.
 115
 116 !print *,'height_data%nx ,soilhgt_data%nx = ',height_data%nx ,soilhgt_data%nx
 117 !print *,'height_data%ny ,soilhgt_data%ny = ',height_data%ny ,soilhgt_data%ny
 118
 119                if (height_data%nx    /= soilhgt_data%nx .or. &
 120                    height_data%ny    /= soilhgt_data%ny .or. &
 121                    height_data%iproj /= soilhgt_data%iproj) then
 122                   call mprintf(.true.,ERROR,'Mismatch in height field dimensions in file %s', &
 123                                s1=trim(input_name)//':'//temp_date(1:13))
 124                end if
 125
 126                open ( 10 , &
 127                       file = 'util/vertical_grid_70_20m_80km.txt' , &
 128                       form = 'formatted' , &
 129                       access = 'sequential' , &
 130                       status = 'old' , &
 131                       iostat = istatus )
 132
 133                if ( istatus /= 0 ) then
 134                   call mprintf(.true.,ERROR,'Cannot open the UKMO file util/vertical_grid_70_20m_80km.txt')
 135                end if
 136
 137                read (10 , * )
 138                read (10 , * )
 139                read (10 , * )
 140                read (10 , fmt='(30x,i2)' ) model_levels
 141                read (10 , fmt='(30x,i2)' ) first_constant_r_rho_level
 142                read (10 , fmt='(30x,f17.11)' ) z_top_of_model
 143 print *,'model_levels, first_constant_r_rho_level, z_top_of_model = ', model_levels, first_constant_r_rho_level, z_top_of_model
 144                read (10 , * )
 145                read (10 , * )
 146
 147                allocate ( eta ( model_levels ) )
 148                eta = 0.
 149                do k_loop = 1 , first_constant_r_rho_level-1
 150                   read ( 10 , fmt='(65x,f10.7)' ) eta(k_loop)
 151 !print *,k_loop,' eta = ',eta(k_loop)
 152                end do
 153                k_loop = first_constant_r_rho_level
 154                read ( 10 , fmt='(29x,f10.7,26x,f10.7)' ) etac , eta(k_loop)
 155                do k_loop = first_constant_r_rho_level+1 , model_levels
 156                   read ( 10 , fmt='(65x,f10.7)' ) eta(k_loop)
 157 !print *,k_loop,' eta = ',eta(k_loop)
 158                end do
 159
 160                exit all_fields
 161             end if
 162
 163
 164          end if
 165
 166       end do all_fields
 167
 168       call read_met_close()
 169
 170    else
 171       call mprintf(.true.,ERROR,'Problem opening %s at time %s', s1=input_name, s2=temp_date(1:13))
 172    end if
 173
 174    !  Now we have to write the height out for each time.
 175
 176    if (associated(height_data%slab)) then
 177
 178       still_more_dates_to_process = .true.
 179       soilhgt_data%xfcst=  0.
 180       soilhgt_data%hdate=  output_date // '.0000'
 181       height_data%xfcst=  0.
 182       height_data%hdate=  output_date // '.0000'
 183       date_loop = 1
 184
 185       all_dates : do while ( still_more_dates_to_process )
 186
 187 print *,'Generating ',temp_count,' height levels for ',output_date,'.'
 188
 189          call write_met_init('HGT', .false., output_date(1:13), istatus)
 190          soilhgt_data%xlvl = 200100
 191          soilhgt_data%field = 'SOILHGT  '
 192          soilhgt_data%desc = 'Computed Geopotential Height, UKMO diagnostic '
 193
 194          !  Corrections for UKMO
 195
 196 !     soilhgt_data%deltalon=  360. / real(soilhgt_data%nx)
 197
 198          do j_loop = 1 , soilhgt_data%ny
 199             do i_loop = 1 , soilhgt_data%nx
 200                soilhgt_data%slab(i_loop,j_loop) = soil_hold(i_loop,j_loop)
 201             end do
 202          end do
 203
 204          call write_next_met_field(soilhgt_data, istatus)
 205
 206          height_data%version=  soilhgt_data%version
 207          height_data%nx=  soilhgt_data%nx
 208          height_data%ny=  soilhgt_data%ny
 209          height_data%iproj=  soilhgt_data%iproj
 210 !        height_data%xfcst=  soilhgt_data%xfcst
 211 !        height_data%xlvl=  soilhgt_data%xlvl
 212          height_data%startlat=  soilhgt_data%startlat
 213          height_data%startlon=  soilhgt_data%startlon
 214          height_data%starti=  soilhgt_data%starti
 215          height_data%startj=  soilhgt_data%startj
 216          height_data%deltalat=  soilhgt_data%deltalat
 217          height_data%deltalon=  soilhgt_data%deltalon
 218          height_data%dx=  soilhgt_data%dx
 219          height_data%dy=  soilhgt_data%dy
 220          height_data%xlonc=  soilhgt_data%xlonc
 221          height_data%truelat1=  soilhgt_data%truelat1
 222          height_data%truelat2=  soilhgt_data%truelat2
 223          height_data%earth_radius=  soilhgt_data%earth_radius
 224          height_data%is_wind_grid_rel=  soilhgt_data%is_wind_grid_rel
 225          height_data%field=  'HGT      '
 226 !        height_data%hdate=  soilhgt_data%hdate
 227          height_data%units=  soilhgt_data%units
 228          height_data%map_source=  soilhgt_data%map_source
 229          height_data%desc=  soilhgt_data%desc
 230
 231          do k_loop = MIN(model_levels,temp_count) , first_constant_r_rho_level+1 , -1
 232 !        height_data%xlvl = temp_count + 1 - k_loop
 233             height_data%xlvl =                  k_loop
 234             do j_loop = 1 , height_data%ny
 235                do i_loop = 1 , height_data%nx
 236                   height_data%slab(i_loop,j_loop) = eta(k_loop) * z_top_of_model
 237                end do
 238             end do
 239
 240             call write_next_met_field(height_data, istatus)
 241          end do
 242
 243          do k_loop = first_constant_r_rho_level , 1 , -1
 244 !        height_data%xlvl = temp_count + 1 - k_loop
 245             height_data%xlvl =                  k_loop
 246             do j_loop = 1 , height_data%ny
 247                do i_loop = 1 , height_data%nx
 248                   height_data%slab(i_loop,j_loop) = eta(k_loop) * z_top_of_model + &
 249                                           soil_hold(i_loop,j_loop) * ( 1. - eta(k_loop) / etac ) **2
 250                end do
 251             end do
 252
 253             call write_next_met_field(height_data, istatus)
 254          end do
 255
 256          call write_met_close()
 257
 258          call geth_newdate ( output_date , valid_date , interval_seconds * date_loop )
 259          date_loop = date_loop + 1
 260
 261          if ( TRIM(output_date) > TRIM(end_date(1)) ) then
 262             still_more_dates_to_process = .false.
 263          end if
 264
 265       end do all_dates
 266
 267    end if
 268
 269
 270    call mprintf(.true.,STDOUT,' *** Successful completion of program height_ukmo.exe *** ')
 271
 272    stop
 273
 274 end program height_ukmo