Simplify klevel/remove kadd

Author: magnusuMET

src/common/bldp.f90

@@ -39,7 +39,7 @@ module bldpML
 subroutine bldp
   USE snapdimML, only: nx,ny,nk
   USE snaptabML, only: cp, g, pmult, pitab
-  USE snapgrdML, only: ahalf, bhalf, vhalf, alevel, blevel, kadd
+  USE snapgrdML, only: ahalf, bhalf, vhalf, alevel, blevel
   USE snapfldML, only: u2, v2, ps2, t2, hbl2, bl2
   USE ftestML, only: ftest
   USE snapdebug, only: iulog
@@ -80,7 +80,7 @@ subroutine bldp
   if(pblbot > p) pblbot=p
   kbltop=2
   kblbot=2
-  nkk=min(nk-2,nk-kadd)
+  nkk=min(nk-2,nk)
   do k=2,nkk
     p=ahalf(k)+bhalf(k)*psurf
     if(p > pbltop) kbltop=k+1

src/common/om2edot.f90

@@ -55,7 +55,7 @@ module om2edotML
 !>                 More than half of the model levels and the upper
 !>                 level are required.
 subroutine om2edot
-  USE snapgrdML, only: ahalf, bhalf, vhalf, kadd, klevel, gparam
+  USE snapgrdML, only: ahalf, bhalf, vhalf, klevel, gparam
   USE snapfldML, only: xm, ym, ps2, u2, v2, w2, field1, field2, field3, field4
   USE snapdimML, only: nx,ny,nk
   USE snapdebug, only: iulog
@@ -67,7 +67,7 @@ subroutine om2edot
 
   write(iulog,*) 'OM2EDOT'
 
-  do k=2,nk-kadd
+  do k=2,nk
 
     deta=vhalf(k-1)-vhalf(k)
 
@@ -85,11 +85,11 @@ subroutine om2edot
 !..check if required model levels are present for subr. edcomp
 
 !..no. of present levels
-  kk=(nk-kadd)-1
+  kk=(nk)-1
 !..assuming that the lower model level is always present above surface level
   km=klevel(2)
 
-  if(kk > km/2 .AND. kadd == 0 .AND. klevel(nk) == 1) then
+  if(kk > km/2 .AND. klevel(nk) == 1) then
 
     d2hx=1./(gparam(7)*2.)
     d2hy=1./(gparam(8)*2.)

src/common/readfield_fi.f90

@@ -72,7 +72,7 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
       hlayer1, hlayer2, bl1, bl2, enspos, precip, t1_abs, t2_abs, &
       field1
     USE snapgrdML, only: alevel, blevel, vlevel, ahalf, bhalf, vhalf, &
-                         gparam, kadd, klevel, ivlevel, imslp, igtype, ivlayer, ivcoor
+                         gparam, klevel, ivlevel, imslp, igtype, ivlayer, ivcoor
     USE snapmetML, only: met_params, xy_wind_units, pressure_units, omega_units, &
                          sigmadot_units, temp_units, requires_precip_deaccumulation
     USE snapdimML, only: nx, ny, nk, output_resolution_factor, hres_field
@@ -215,7 +215,7 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
     end if
 
     ptop = 100.0
-    do k = nk - kadd, 2, -1
+    do k = nk , 2, -1
 
       !..input model level no.
       ilevel = klevel(k)
@@ -269,7 +269,7 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
         end if
       end if
 
-    end do ! k=nk-kadd,2,-1
+    end do ! k=nk,2,-1
 
 !..surface pressure, 10m wind and possibly mean sea level pressure,
 !..precipitation
@@ -333,7 +333,7 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
     if (met_params%temp_is_abs) then
       if (allocated(t2_abs)) t2_abs(:,:,:) = t2
       !..abs.temp. -> pot.temp.
-      do k = 2, nk - kadd
+      do k = 2, nk
         do j = 1, ny
           do i = 1, nx
             p = alevel(k) + blevel(k)*ps2(i,j)
@@ -344,7 +344,7 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
     else
       if (allocated(t2_abs)) then
         ! pot.temp -> abs.temp
-        do k=2,nk-kadd
+        do k=2,nk
           do j = 1, ny
             do i = 1, nx
               p = alevel(k) + blevel(k)*ps2(i,j)
@@ -370,19 +370,6 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
 
 !..no temperature at or near surface (not used, yet)
     t2(:, :, 1) = -999.0
-    if (kadd > 0) then
-      !..levels added at the top
-      dred = 0.5/float(kadd)
-      red = 1.
-      kk = nk - kadd
-      do k = nk - kadd + 1, nk
-        red = red - dred
-        u2(:, :, k) = u2(:, :, kk)
-        v2(:, :, k) = v2(:, :, kk)
-        w2(:, :, k) = w2(:, :, kk)*red
-        t2(:, :, k) = t2(:, :, kk)
-      end do
-    end if
 
     if (backward) then
       ! backward-calculation, switch sign of winds

src/common/readfield_nc.f90

@@ -115,7 +115,7 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
       hlayer1, hlayer2, bl1, bl2, enspos, precip, &
       t1_abs, t2_abs, field1
   USE snapgrdML, only: alevel, blevel, vlevel, ahalf, bhalf, vhalf, &
-      gparam, kadd, klevel, ivlevel, imslp, igtype, ivlayer, ivcoor
+      gparam, klevel, ivlevel, imslp, igtype, ivlayer, ivcoor
   USE snapmetML, only: met_params, requires_precip_deaccumulation
   USE snapdimML, only: nx, ny, nk, output_resolution_factor, hres_field
   USE datetime, only: datetime_t, duration_t
@@ -261,7 +261,7 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
   end if
 
   ptop = 100.0
-  do k=nk-kadd,2,-1
+  do k=nk,2,-1
 
   !..input model level no.
     ilevel=klevel(k)
@@ -318,7 +318,7 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
           start4d, count4d, w2(:,:,k))
     end if
 
-  end do ! k=nk-kadd,2,-1
+  end do ! k=nk,2,-1
 
 
 !..surface pressure, 10m wind and possibly mean sea level pressure,
@@ -401,7 +401,7 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
   if (met_params%temp_is_abs) then
     if (allocated(t2_abs)) t2_abs(:,:,:) = t2
   !..abs.temp. -> pot.temp.
-    do k=2,nk-kadd
+    do k=2,nk
       do j = 1, ny
         do i = 1, nx
           p = alevel(k) + blevel(k)*ps2(i,j)
@@ -412,7 +412,7 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
   else
     if (allocated(t2_abs)) then
       ! pot.temp -> abs.temp
-      do k=2,nk-kadd
+      do k=2,nk
         do j = 1, ny
           do i = 1, nx
             p = alevel(k) + blevel(k)*ps2(i,j)
@@ -438,19 +438,6 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
 
 !..no temperature at or near surface (not used, yet)
   t2(:,:,1) = -999.0
-  if(kadd > 0) then
-  !..levels added at the top
-    dred=0.5/float(kadd)
-    red=1.
-    kk=nk-kadd
-    do k=nk-kadd+1,nk
-      red=red-dred
-      u2(:,:,k) = u2(:,:,kk)
-      v2(:,:,k) = v2(:,:,kk)
-      w2(:,:,k) = w2(:,:,kk)*red
-      t2(:,:,k) = t2(:,:,kk)
-    end do
-  end if
 
   if(backward) then
   ! backward-calculation, switch sign of winds
@@ -916,7 +903,7 @@ end subroutine nfcheckload3d
 
 subroutine compute_vertical_coords(alev, blev, ptop)
   use iso_fortran_env, only: error_unit
-  use snapgrdML, only: alevel, blevel, vlevel, ivcoor, klevel, kadd, &
+  use snapgrdML, only: alevel, blevel, vlevel, ivcoor, klevel, &
                        ahalf, bhalf, vhalf
   use snapmetML, only: met_params
   use snapdimML, only: nk
@@ -929,38 +916,11 @@ subroutine compute_vertical_coords(alev, blev, ptop)
   real :: db, dp, p1, p2
   integer :: k
 
-  do k = 2, nk - kadd
+  do k = 2, nk
     alevel(k) = alev(k)
     blevel(k) = blev(k)
   end do
 
-  if (kadd > 0) then
-    if (ivcoor == 2) then
-      !..sigma levels ... blevel=sigma
-      db = blevel(nk - kadd - 1) - blevel(nk - kadd)
-      db = max(db, blevel(nk - kadd)/float(kadd))
-      do k = nk - kadd + 1, nk
-        blevel(k) = max(blevel(k - 1) - db, 0.)
-      end do
-    elseif (ivcoor == 10) then
-      !..eta (hybrid) levels
-      p1 = alevel(nk - kadd) + blevel(nk - kadd)*1000.
-      p2 = alevel(nk - kadd - 1) + blevel(nk - kadd - 1)*1000.
-      dp = p2 - p1
-      if (p1 - dp*kadd < 10.) dp = (p1 - 10.)/kadd
-      db = blevel(nk - kadd - 1) - blevel(nk - kadd)
-      db = max(db, blevel(nk - kadd)/float(kadd))
-      do k = nk - kadd + 1, nk
-        p1 = p1 - dp
-        blevel(k) = max(blevel(k - 1) - db, 0.)
-        alevel(k) = p1 - blevel(k)*1000.
-      end do
-    else
-      write (error_unit, *) 'PROGRAM ERROR.  ivcoor= ', ivcoor
-      error stop 255
-    end if
-  end if
-
   if (ivcoor == 2) then
     !..sigma levels (norlam)
     do k = 2, nk

src/common/snap.F90

@@ -171,7 +171,7 @@ PROGRAM bsnap
                        ncomp, def_comp, nparnum, &
                        time_profile
   USE snapposML, only: irelpos, nrelpos, release_positions
-  USE snapgrdML, only: modleveldump, ivcoor, kadd, &
+  USE snapgrdML, only: modleveldump, ivcoor, &
                        klevel, imslp, itotcomp, gparam, &
                        igtype, imodlevel, precipitation_in_output
   USE snaptabML, only: tabcon
@@ -232,7 +232,7 @@ PROGRAM bsnap
   integer :: nhrun = 0, nhrel = 0
   logical :: use_random_walk = .true.
   logical :: autodetect_grid_params = .false.
-  integer :: m, np, npl, nlevel = 0, ifltim = 0
+  integer :: m, np, npl, nlevel, ifltim = 0
   logical :: synoptic_output = .false.
   integer :: k, ierror, i, n
   integer :: ih
@@ -459,7 +459,6 @@ PROGRAM bsnap
 
     !..information to log file
     write (iulog, *) 'nx,ny,nk:  ', nx, ny, nk
-    write (iulog, *) 'kadd:      ', kadd
     write (iulog, *) 'klevel:'
     write (iulog, *) (klevel(i), i=1, nk)
     write (iulog, *) 'imslp:     ', imslp
@@ -1584,25 +1583,21 @@ subroutine read_inputfile(snapinput_unit)
         ivcoor = 10
       case ('levels.input')
         !..levels.input=<num_levels, 0,kk,k,k,k,....,1>
-        !..levels.input=<num_levels, 0,kk,k,k,k,....,18,0,0,...>
+        !..levels.input=<num_levels, 0,kk,k,k,k,....,18>
         if (.not. has_value) goto 12
-        if (nlevel /= 0) then
+        if (allocated(klevel)) then
           write (error_unit, *) "re-assigning levels"
           DEALLOCATE(klevel, STAT=AllocateStatus)
         end if
         read (cinput(pname_start:pname_end), *, err=12) nlevel
         nk = nlevel
-        ALLOCATE (klevel(nk), STAT=AllocateStatus)
+        ALLOCATE(klevel(nk), STAT=AllocateStatus)
         IF (AllocateStatus /= 0) ERROR STOP AllocateErrorMessage
-!         ALLOCATE ( ipcount(mdefcomp, nk), STAT = AllocateStatus)
-!         IF (AllocateStatus /= 0) ERROR STOP AllocateErrorMessage
-!         ALLOCATE ( npcount(nk), STAT = AllocateStatus)
-!         IF (AllocateStatus /= 0) ERROR STOP AllocateErrorMessage
 
         read (cinput(pname_start:pname_end), *, err=12) nlevel, (klevel(i), i=1, nlevel)
         if (klevel(1) /= 0 .OR. klevel(2) == 0) goto 12
-        kadd = count(klevel(2:nk) == 0)
-        do i = nk - kadd - 1, 2, -1
+
+        do i = nk - 1, 2, -1
           if (klevel(i) <= klevel(i + 1)) goto 12
         end do
       case ('forecast.hour.min')
@@ -1826,8 +1821,7 @@ subroutine conform_input(ierror)
           endif
         endif
       end if
-      nlevel = size(klevel)
-      nk = nlevel
+      nk = size(klevel)
       write (error_unit, *) "autodetection of grid-param: ", gparam
     end if
 
@@ -1845,7 +1839,7 @@ subroutine conform_input(ierror)
       write (error_unit, *) 'Input model level type (sigma,eta) not specified'
       ierror = 1
     end if
-    if (nlevel == 0) then
+    if (.not.allocated(klevel)) then
       write (error_unit, *) 'Input model levels not specified'
       ierror = 1
     end if

src/common/snapgrdML.f90

@@ -22,7 +22,7 @@ module snapgrdML
 
 !> mapping from internal vertical coordinate to MET vertical coord
 !>
-!> sequence: bottom to top (kk,kk-1,....1)
+!> sequence: bottom to top (0,kk,kk-1,...,1)
 !>
 !> level no.  1 should be 0    (surface fields here)
 !>
@@ -113,10 +113,6 @@ module snapgrdML
 !> * 2=sigma (Norlam)
 !> * 10=eta/hybrid   (Hirlam,...))
     integer, save, public :: ivcoor = 0
-!> levels added at the top (when missing upper model levels)
-!>
-!> (u,v copied up, w reduced, pot.temp. const.)
-    integer, save, public :: kadd
 !> table of level numbers for interpolation
 !>
 !> (key is vlevel*10000)