temperature.F90 12.2 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: temperature.F90,v 1.17 2006-03-01 15:54:08 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5 6 7 8 9 10 11
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE:  temperature
!
! !INTERFACE:
   module temperature
!
! !DESCRIPTION:
12 13 14 15 16 17
!  
! In this module, the temperature equation is processed by
! reading in the namelist {\tt temp} and initialising the temperature field
! (this is done in {\tt init\_temperature}), 
! and calculating the advection-diffusion-equation, which includes 
! penetrating short-wave radiation as source term (see {\tt do\_temperature}).
gotm's avatar
gotm committed
18 19
!
! !USES:
20
   use exceptions
gotm's avatar
gotm committed
21 22
   use domain, only: imin,jmin,imax,jmax,H,az
   use domain, only: iimin,jjmin,iimax,jjmax,kmax
23
   use variables_3d, only: T,hn,adv_schemes
24
   use halo_zones, only: update_3d_halo,wait_halo,D_TAG
gotm's avatar
gotm committed
25 26 27 28 29 30 31 32
   IMPLICIT NONE
!
   private
!
! !PUBLIC DATA MEMBERS:
   public init_temperature, do_temperature
!
! !PRIVATE DATA MEMBERS:
kbk's avatar
kbk committed
33 34 35 36
   integer                   :: temp_method=1,temp_format=2
   character(len=PATH_MAX)   :: temp_file="t_and_s.nc"
   character(len=32)         :: temp_name='temp'
   REALTYPE                  :: temp_const=20.
37 38
   integer                   :: temp_hor_adv=1,temp_ver_adv=1
   integer                   :: temp_adv_split=0
kbk's avatar
kbk committed
39
   REALTYPE                  :: temp_AH=-1.
gotm's avatar
gotm committed
40
!
kbk's avatar
kbk committed
41 42 43
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
44 45 46 47 48 49 50 51 52
! !LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
53 54
! !IROUTINE: init_temperature - initialisation of temperature
! \label{sec-init-temperature}
gotm's avatar
gotm committed
55 56 57 58 59
!
! !INTERFACE:
   subroutine init_temperature(adv_method)
!
! !DESCRIPTION:
60 61 62 63 64 65 66 67 68 69 70
!
! Here, the temperature equation is initialised. First, the namelist
! {\tt temp} is read from {\tt getm.inp}. Then, depending on the
! {\tt temp\_method}, the temperature field is read from a
! hotstart file ({\tt temp\_method}=0), initialised with a constant value
! ({\tt temp\_method}=1), initialised and interpolated 
! with horizontally homogeneous
! temperature from a given temperature profile ({\tt temp\_method}=2),
! or read in and interpolated from a 3D netCDF field ({\tt temp\_method}=3).
! Finally, a number of sanity checks are performed for the chosen 
! temperature advection schemes.
gotm's avatar
gotm committed
71 72 73 74 75
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
76
   integer, intent(in)                 :: adv_method
gotm's avatar
gotm committed
77 78 79 80 81 82
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
83
   integer                   :: k,i,j,n
84
   integer, parameter        :: nmax=10000
kbk's avatar
kbk committed
85
   REALTYPE                  :: zlev(nmax),prof(nmax)
86
   integer                   :: temp_field_no=1
kbk's avatar
kbk committed
87 88
   NAMELIST /temp/ &
             temp_method,temp_const,temp_file,              &
89
             temp_format,temp_name,temp_field_no,           &
90
             temp_hor_adv,temp_ver_adv,temp_adv_split,temp_AH
gotm's avatar
gotm committed
91 92 93 94 95 96 97 98 99
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_temperature() # ',Ncall
#endif

kbk's avatar
kbk committed
100
#ifdef NS_FRESHWATER_LENSE_TEST
101
temp_field_no=1
gotm's avatar
gotm committed
102 103
#endif
#ifdef MED_15X15MINS_TEST
104
temp_field_no=1
gotm's avatar
gotm committed
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
#endif

   LEVEL2 'init_temperature()'
   read(NAMLST,temp)

   select case (temp_method)
      case(0)
         LEVEL3 'getting initial fields from hotstart'
      case(1)
         LEVEL3 'setting to constant value'
         forall(i=iimin:iimax,j=jjmin:jjmax, az(i,j) .ne. 0) &
                T(i,j,:) = temp_const
      case(2)
         LEVEL3 'using profile'
         call read_profile(temp_file,nmax,zlev,prof,n)
kbk's avatar
kbk committed
120 121
         call ver_interpol(n,zlev,prof,imin,jmin,imax,jmax,az,H,       &
                           iimin,jjmin,iimax,jjmax,kmax,hn,T)
gotm's avatar
gotm committed
122 123
      case(3)
         LEVEL3 'interpolating from 3D field'
124
         call get_field(temp_file,temp_name,temp_field_no,T)
gotm's avatar
gotm committed
125 126 127 128 129
      case default
         FATAL 'Not valid temp_method specified'
         stop 'init_temperature'
   end select

130 131 132 133 134 135 136 137 138 139
!  Sanity checks for advection specifications
   LEVEL3 'temp_hor_adv=   ',temp_hor_adv
   LEVEL3 'temp_ver_adv=   ',temp_ver_adv
   LEVEL3 'temp_adv_split= ',temp_adv_split
   if(temp_hor_adv .eq. 1) then
      temp_adv_split=-1
      if(temp_ver_adv .ne. 1) then
         LEVEL3 "setting temp_ver_adv to 1 - since temp_hor_adv is 1"
         temp_ver_adv=1
      end if
140
   end if
kbk's avatar
kbk committed
141 142
   LEVEL3 "horizontal: ",trim(adv_schemes(temp_hor_adv))
   LEVEL3 "vertical:   ",trim(adv_schemes(temp_ver_adv))
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190

   select case (temp_adv_split)
      case (-1)
      case (0)
         select case (temp_hor_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "temp_adv_split=0: temp_hor_adv not valid (2-6)")

         end select
         select case (temp_ver_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "temp_adv_split=0: temp_ver_adv not valid (2-6)")
         end select
         LEVEL2 "1D split --> full u, full v, full w"
      case (1)
         select case (temp_hor_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "temp_adv_split=1: temp_hor_adv not valid (2-6)")
         end select
         select case (temp_ver_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "temp_adv_split=1: temp_ver_adv not valid (2-6)")
         end select
         LEVEL2 "1D split --> half u, half v, full w, half v, half u"
      case (2)
         select case (temp_hor_adv)
            case (2,7)
            case default
               call getm_error("init_3d()", &
                    "temp_adv_split=2: temp_hor_adv not valid (2,7)")
         end select
         select case (temp_ver_adv)
            case (2,3,4,5,6)
            case default
               call getm_error("init_3d()", &
                    "temp_adv_split=2: temp_ver_adv not valid (2-6)")
         end select
         LEVEL2 "2D-hor, 1D-vert split --> full uv, full w"
      case default
   end select
191

gotm's avatar
gotm committed
192 193 194 195 196 197 198 199 200 201 202 203 204 205
   call update_3d_halo(T,T,az,iimin,jjmin,iimax,jjmax,kmax,D_TAG)
   call wait_halo(D_TAG)

#ifdef DEBUG
   write(debug,*) 'Leaving init_temperature()'
   write(debug,*)
#endif
   return
   end subroutine init_temperature
!EOC

!-----------------------------------------------------------------------
!BOP
!
206
! !IROUTINE:  do_temperature - temperature equation \label{sec-do-temperature}
gotm's avatar
gotm committed
207 208 209 210 211
!
! !INTERFACE:
   subroutine do_temperature(n)
!
! !DESCRIPTION:
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
!  
! Here, one time step for the temperature equation is performed.
! First, preparations for the call to the advection schemes are
! made, i.e.\ calculating the necessary metric coefficients.
! After the call to the advection schemes, which actually perform
! the advection (and horizontal diffusion) step as an operational
! split step, the solar radiation at the interfaces ({\tt rad(k)}) 
! is calculated
! from given surface radiation ({\tt swr\_loc}) 
! by means of a double exponential
! approach, see equation (\ref{Light}) on page \pageref{Light}). 
! Furthermore, the surface heat flux {\tt sfl\_loc} is given a 
! value.
! The sea surface temperature is limited by the freezing point
! temperature (as a most primitive sea ice model). The next
! step is to set up the tri-diagonal matrix for calculating the
! new temperature by means of a semi-implicit central scheme for the
! vertical diffusion. Source terms which appear on the right hand sides 
! are due to the divergence of the solar radiation at the interfaces.
! The subroutine is completed by solving the tri-diagonal linear
! equation by means of a tri-diagonal solver.
gotm's avatar
gotm committed
233 234 235 236 237 238 239 240 241 242 243 244
!
! !USES:
   use advection_3d, only: do_advection_3d
   use variables_3d, only: dt,cnpar,hn,ho,nuh,uu,vv,ww,hun,hvn,S
   use domain,       only: iimin,iimax,jjmin,jjmax,kmax,az,au,av
   use meteo,        only: swr,shf
   use parameters,   only: rho_0,cp
#if defined(SPHERICAL) || defined(CURVILINEAR)
   use domain, only: dxu,dxv,dyu,dyv,arcd1
#else
   use domain, only: dx,dy,ard1
#endif
245
   use parameters, only: avmolt
gotm's avatar
gotm committed
246 247 248 249 250 251 252 253 254 255
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
  integer, intent(in) :: n
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
256 257 258 259 260 261 262 263 264 265
   integer                   :: i,j,k,rc
   REALTYPE                  :: Res(0:kmax)
   REALTYPE                  :: auxn(1:kmax-1),auxo(1:kmax-1)
   REALTYPE                  :: a1(0:kmax),a2(0:kmax)
   REALTYPE                  :: a3(0:kmax),a4(0:kmax)
   REALTYPE                  :: delxu(I2DFIELD),delxv(I2DFIELD)
   REALTYPE                  :: delyu(I2DFIELD),delyv(I2DFIELD)
   REALTYPE                  :: area_inv(I2DFIELD)
   REALTYPE                  :: swr_loc,shf_loc
   REALTYPE                  :: zz,rad(0:1000),A=0.58,g1=0.35,g2=23.0
gotm's avatar
gotm committed
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_temperature() # ',Ncall
#endif

#if defined(SPHERICAL) || defined(CURVILINEAR)
   delxu=dxu
   delxv=dxv
   delyu=dyu
   delyv=dyv
   area_inv=arcd1
#else
   delxu=dx
   delxv=dx
   delyu=dy
   delyv=dy
   area_inv=ard1
#endif
   call do_advection_3d(dt,T,uu,vv,ww,hun,hvn,ho,hn,   &
kbk's avatar
kbk committed
289
                        delxu,delxv,delyu,delyv,area_inv,az,au,av,     &
290
                        temp_hor_adv,temp_ver_adv,temp_adv_split,temp_AH)
kbk's avatar
kbk committed
291
#ifdef FRESHWATER_LENSE_TEST
292 293 294 295
   T(iimin:iimin+3,jjmin:jjmax,1:kmax)=10.
   T(iimax-3:iimax,jjmin:jjmax,1:kmax)=10.
   T(iimin:iimax,jjmin:jjmin+3,1:kmax)=10.
   T(iimin:iimax,jjmax-3:jjmax,1:kmax)=10.
gotm's avatar
gotm committed
296 297 298 299 300 301 302 303 304 305
#endif

!  Solar radiation and vertical diffusion of temperature

   do j=jjmin,jjmax
      do i=iimin,iimax
         if (az(i,j) .eq. 1) then

! Solar radiation
            swr_loc=swr(i,j)
kbk's avatar
kbk committed
306
            shf_loc=shf(i,j)
307 308
            if (T(i,j,kmax).le.-0.0575*S(i,j,kmax)) then  ! use most primitive 
                                                          ! sea ice model ...
kbk's avatar
kbk committed
309
               shf_loc=max(_ZERO_,shf_loc)
kbk's avatar
kbk committed
310 311 312 313
            end if
            rad(kmax)=(swr_loc+shf_loc)/(rho_0*cp)
            zz = _ZERO_
            do k=kmax-1,0,-1
gotm's avatar
gotm committed
314
               zz=zz+hn(i,j,k+1)
kbk's avatar
kbk committed
315 316
               rad(k)=swr_loc/(rho_0*cp)*(A*exp(-zz/g1)+(1-A)*exp(-zz/g2))
            end do
gotm's avatar
gotm committed
317 318 319 320

            if (kmax.gt.1) then
!     Auxilury terms, old and new time level,
               do k=1,kmax-1
321
                  auxo(k)=2.*(1-cnpar)*dt*(nuh(i,j,k)+avmolt)/ &
kbk's avatar
kbk committed
322
                             (hn(i,j,k+1)+hn(i,j,k))
323
                  auxn(k)=2.*   cnpar *dt*(nuh(i,j,k)+avmolt)/ &
kbk's avatar
kbk committed
324
                             (hn(i,j,k+1)+hn(i,j,k))
kbk's avatar
kbk committed
325
               end do
gotm's avatar
gotm committed
326 327

!        Matrix elements for surface layer
kbk's avatar
kbk committed
328 329 330 331 332
               k=kmax
               a1(k)=-auxn(k-1)
               a2(k)=hn(i,j,k)+auxn(k-1)
               a4(k)=T(i,j,k)*(hn(i,j,k)-auxo(k-1))+T(i,j,k-1)*auxo(k-1)  &
                     +dt*(rad(k)-rad(k-1))
gotm's avatar
gotm committed
333 334 335

!        Matrix elements for inner layers
               do k=2,kmax-1
kbk's avatar
kbk committed
336 337 338 339 340 341 342
                  a3(k)=-auxn(k  )
                  a1(k)=-auxn(k-1)
                  a2(k)=hn(i,j,k)+auxn(k)+auxn(k-1)
                  a4(k)=T(i,j,k+1)*auxo(k)                          &
                       +T(i,j,k  )*(hn(i,j,k)-auxo(k)-auxo(k-1))    &
                       +T(i,j,k-1)*auxo(k-1)                        &
                       +dt*(rad(k)-rad(k-1))
gotm's avatar
gotm committed
343 344 345 346
               end do

!        Matrix elements for bottom layer
               k=1
kbk's avatar
kbk committed
347 348
               a3(k)=-auxn(k  )
               a2(k)=hn(i,j,k)+auxn(k)
gotm's avatar
gotm committed
349
               a4(k)=T(i,j,k+1)*auxo(k)                           &
kbk's avatar
kbk committed
350 351
                    +T(i,j,k  )*(hn(i,j,k)-auxo(k))               &
                    +dt*(rad(k)-rad(k-1))
gotm's avatar
gotm committed
352 353 354

               call getm_tridiagonal(kmax,1,kmax,a1,a2,a3,a4,Res)

kbk's avatar
kbk committed
355 356 357
               do k=1,kmax
                  T(i,j,k)=Res(k)
               end do
gotm's avatar
gotm committed
358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381

            end if
         end if
      end do
   end do

   call update_3d_halo(T,T,az,iimin,jjmin,iimax,jjmax,kmax,D_TAG)
   call wait_halo(D_TAG)

#ifdef DEBUG
   write(debug,*) 'Leaving do_temperature()'
   write(debug,*)
#endif
   return
   end subroutine do_temperature
!EOC

!-----------------------------------------------------------------------

   end module temperature

!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding               !
!-----------------------------------------------------------------------