meteo.F90 27.1 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: meteo - provides meteorological forcing for \emph{getm}.
!
! !INTERFACE:
   module meteo
!
! !DESCRIPTION:
!  The meteo module provides meteorological forcing for \emph{getm}.
!  The main role of the module is to supply the following 4 fields with
!  sane values - \emph{airp, tausx, tausy, swr} and \emph{shf} - i.e.
!  air pressure [$Pa$], surface stresses [$N/m^2$], short wave radiation
!  [$W/m^2$] and surface heat fluxes [$W/m^2$] on the computational grid.
!  The module provides 3 public functions - \emph{init\_meteo()},
kbk's avatar
kbk committed
17
!  \emph{do\_meteo()} and \emph{clean\_meteo()} - and a number of public
18 19
!  data members to hold the actual meteorological fields. Also included
!  in the module are various constants related to calculating the
kbk's avatar
kbk committed
20
!  meteorological forcing.
gotm's avatar
gotm committed
21 22 23
!  Information about the calculation domain is obtained from the module
!  \emph{domain} and time related information comes from the module
!  \emph{time}.
24 25
!  The meteo module is initialised via a call to \emph{init\_meteo()}
!  that will read a namelist providing all necessary information. Memory
kbk's avatar
kbk committed
26
!  allocation is also done in \emph{init\_meteo()}.
27 28 29 30
!  Obtaining the actual forcing - either by reading from a file or
!  calculating is done via calls to \emph{do\_meteo()}. The actual
!  reading of external data from files is separated completely from
!  \emph{do\_meteo()} and is done in the main time loop via a call to
kbk's avatar
kbk committed
31
!  \emph{do\_input()} where all external file input is handled.
32 33 34 35 36 37 38
!  \emph{meteo} supplies 3 variables which can be used by routines for
!  reading variables. \emph{new\_meteo} is a logical switch which should
!  be set to .true. when new fields have been read. \emph{t\_1} and
!  \emph{t\_2} holds the time (in seconds) since the model run of the two
!  fields surrounding the actual model time - to be used by the temporal
!  interpolation. Finally \emph{clean\_meteo()} should be called when
!  the simulation is over as part of the overall procedure of finalising
kbk's avatar
kbk committed
39
!  the model run.
gotm's avatar
gotm committed
40 41
!
! !SEE ALSO:
42 43
!  short_wave_radiation.F90, solar_zenith_angle.F90, albedo_water.F90, 
!  fluxes.F90, exchange_coefficients.F90
gotm's avatar
gotm committed
44 45 46
!
! !USES:
   use time, only: yearday,secondsofday,timestep
kbk's avatar
kbk committed
47
   use halo_zones, only : H_TAG,update_2d_halo,wait_halo
48
   use domain, only: imin,imax,jmin,jmax,lonc,latc,convc,az
gotm's avatar
gotm committed
49 50 51 52 53
   IMPLICIT NONE
!
   private
!
! !PUBLIC MEMBER FUNCTIONS:
kbk's avatar
kbk committed
54
   public                              :: init_meteo, do_meteo, clean_meteo
gotm's avatar
gotm committed
55 56
!
! !PUBLIC DATA MEMBERS:
kbk's avatar
kbk committed
57 58 59 60
   character(LEN = PATH_MAX), public   :: meteo_file
   logical, public                     :: metforcing=.false.
   logical, public                     :: on_grid=.true.
   logical, public                     :: calc_met=.false.
kbk's avatar
kbk committed
61
   integer, public                     :: met_method
62
   integer, public                     :: albedo_method=1
63 64 65
   integer, public                     :: fwf_method=0
   REALTYPE, public                    :: evap_factor = _ONE_
   REALTYPE, public                    :: precip_factor = _ONE_
66
   REALTYPE, public                    :: w,L,rho_air,qs,qa,ea,es
67 68 69
   REALTYPE, public, dimension(:,:), allocatable, target  :: airp,tausx,tausy
   REALTYPE, public, dimension(:,:), allocatable, target  :: zenith_angle
   REALTYPE, public, dimension(:,:), allocatable, target  :: albedo,swr,shf
70
   REALTYPE, public, dimension(:,:), allocatable, target  :: u10,v10,t2,hum,tcc
Karsten Bolding's avatar
Karsten Bolding committed
71
   REALTYPE, public, dimension(:,:), allocatable, target  :: evap,precip
kbk's avatar
kbk committed
72
   REALTYPE, public                    :: cd_mom,cd_heat,cd_latent
73
   REALTYPE, public                    :: cd_precip = _ZERO_
kbk's avatar
kbk committed
74 75
   REALTYPE, public                    :: t_1=-_ONE_,t_2=-_ONE_
   logical, public                     :: new_meteo=.false.
76 77 78 79
   integer, public, parameter          :: RELATIVE_HUM=1
   integer, public, parameter          :: WET_BULB=2
   integer, public, parameter          :: DEW_POINT=3
   integer, public, parameter          :: SPECIFIC_HUM=4
kbk's avatar
kbk committed
80
   integer, public                     :: hum_method=-1
gotm's avatar
gotm committed
81 82
!
! !DEFINED PARAMETERS:
83 84 85 86
   REALTYPE,public,parameter           :: cpa=1008.d0  !AS that is not exact !
   REALTYPE,public,parameter           :: KELVIN=273.15d0
   REALTYPE,public,parameter           :: emiss=0.97d0
   REALTYPE,public,parameter           :: bolz=5.67d-8
87
!  REALTYPE,public,parameter           :: cpa=1004.67 ! specific heat of dry air- correct
88 89
   REALTYPE,public,parameter           :: cpw=4192.d0   ! specific heat of sea water
   REALTYPE,public,parameter           :: rho_precip = 1000.0d0
gotm's avatar
gotm committed
90 91 92 93 94
!
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
! !LOCAL VARIABLES:
95
   integer                   :: meteo_ramp=0,metfmt=2
kbk's avatar
kbk committed
96
   REALTYPE                  :: tx= _ZERO_ ,ty= _ZERO_
97
   REALTYPE                  :: albedo_const= _ZERO_
kbk's avatar
kbk committed
98
   REALTYPE                  :: swr_const= _ZERO_ ,shf_const= _ZERO_
99
   REALTYPE                  :: evap_const= _ZERO_ ,precip_const= _ZERO_
100 101
   REALTYPE, dimension(:,:), allocatable :: airp_old,airp_new
   REALTYPE, dimension(:,:), allocatable :: tausx_old,tausy_old
gotm's avatar
gotm committed
102
   REALTYPE, dimension(:,:), allocatable :: d_airp,d_tausx,d_tausy
103 104
   REALTYPE, dimension(:,:), allocatable :: tcc_old,tcc_new
   REALTYPE, dimension(:,:), allocatable :: swr_old,shf_old
105
   REALTYPE, dimension(:,:), allocatable :: d_tcc,d_swr,d_shf
106 107
   REALTYPE, dimension(:,:), allocatable :: evap_old,precip_old
   REALTYPE, dimension(:,:), allocatable :: d_evap,d_precip
gotm's avatar
gotm committed
108 109 110 111 112 113 114 115 116 117 118
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_meteo - initialise the \emph{meteo} module.
!
! !INTERFACE:
kbk's avatar
kbk committed
119
   subroutine init_meteo(hotstart)
gotm's avatar
gotm committed
120 121 122 123 124 125 126
   IMPLICIT NONE
!
! !DESCRIPTION:
!  Basically reads the namelist \emph{meteo} from unit NAMLST.  According to
!  the content of the namelist various variables are allocated and initialised.
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
127
   logical, intent(in)                 :: hotstart
gotm's avatar
gotm committed
128 129 130 131 132 133
!
! !REVISION HISTORY:
!
!  See log for module.
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
134
   integer                   :: rc
135 136
   namelist /meteo/ metforcing,on_grid,calc_met,met_method,albedo_method, &
                    fwf_method, &
137
                    meteo_ramp,metfmt,meteo_file, &
138 139
                    tx,ty,albedo_method,albedo_const, &
                    swr_const,shf_const,evap_const,precip_const, &
140
                    precip_factor,evap_factor
gotm's avatar
gotm committed
141 142 143 144 145 146 147 148 149 150 151 152 153
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_meteo() # ',Ncall
#endif
   LEVEL1 'init_meteo'
   read(NAMLST,meteo)

   LEVEL2 'Metforcing=',metforcing
   if (metforcing) then
kbk's avatar
kbk committed
154
      select case (met_method)
gotm's avatar
gotm committed
155 156
            case (1)
               LEVEL2 'Constant forcing is used:'
157 158 159 160
               LEVEL3 'tx     = ',tx
               LEVEL3 'ty     = ',ty
               LEVEL3 'swr    = ',swr_const
               LEVEL3 'shf    = ',shf_const
Knut's avatar
Knut committed
161
               calc_met = .false.
gotm's avatar
gotm committed
162 163 164 165 166 167 168 169 170 171 172 173 174
            case (2)
               if(on_grid) then
                  LEVEL2 'Meteorological fields are on the computational grid'
               else
                  LEVEL2 'Meteorological fields needs to be interpolated'
               end if
               if(calc_met) then
                  LEVEL2 'Stresses and fluxes will be calculated'
               else
                  LEVEL2 'Stresses and fluxes are already calculated'
               end if
         case default
      end select
175

176 177 178 179 180 181 182 183 184 185 186
      LEVEL2 'Albedo method =',albedo_method
      select case (albedo_method)
            case (0)
               LEVEL3 'albedo = ',albedo_const
            case (1)
               LEVEL3 'Albedo according to Payne'
            case (2)
               LEVEL3 'Albedo according to Cogley'
         case default
      end select

187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
      select case (fwf_method)
            case (1)
               LEVEL2 'Constant evaporation/precipitation'
               LEVEL3 'evap   = ',evap_const
               LEVEL3 'precip = ',precip_const
            case (2)
               LEVEL2 'Evaporation/precipitation read from file'
            case (3)
               LEVEL2 'Precipitation read from file'
               LEVEL2 'Evaporation calculated'
            case (4)
               LEVEL2 'No precipitation'
               LEVEL2 'Evaporation calculated'
         case default
      end select

203 204 205
      if (hotstart .and. meteo_ramp .gt. 0) then
         LEVEL2 'hotstart --> meteo_ramp=-1'
         meteo_ramp=-1
kbk's avatar
kbk committed
206
      end if
207 208
      if ( meteo_ramp .gt. 0) then
         LEVEL2 'Forcing will be spun up over ',meteo_ramp,' timesteps'
kbk's avatar
kbk committed
209
      end if
gotm's avatar
gotm committed
210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
   end if

! Allocates memory for the public data members

   allocate(airp(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (airp)'
   airp = _ZERO_

   allocate(tausx(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (tausx)'
   tausx = _ZERO_

   allocate(tausy(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (tausy)'
   tausy = _ZERO_

226 227
   allocate(zenith_angle(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (zenith_angle)'
Knut's avatar
Knut committed
228
   zenith_angle = _ZERO_
229

gotm's avatar
gotm committed
230 231 232 233
   allocate(swr(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (swr)'
   swr = _ZERO_

234 235 236 237
   allocate(albedo(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (albedo)'
   albedo = albedo_const

gotm's avatar
gotm committed
238 239 240 241
   allocate(shf(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (shf)'
   shf = _ZERO_

242 243 244 245 246
   allocate(evap(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (evap)'

   allocate(precip(E2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_meteo: Error allocating memory (precip)'
247 248 249 250 251 252 253
   if (fwf_method .eq. 1) then
      evap = evap_const
      precip = precip_const
   else
      evap = _ZERO_
      precip = _ZERO_
   end if
254

gotm's avatar
gotm committed
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
   if (metforcing) then

      if (calc_met) then
         allocate(u10(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (u10)'
         u10 = _ZERO_

         allocate(v10(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (v10)'
         v10 = _ZERO_

         allocate(t2(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (t2)'
         t2 = _ZERO_

         allocate(hum(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (hum)'
         hum = _ZERO_

274 275 276
         allocate(tcc(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (tcc)'
         tcc = _ZERO_
gotm's avatar
gotm committed
277 278 279 280 281 282 283 284

      else
      end if

      allocate(airp_old(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (airp_old)'
      airp_old = _ZERO_

285 286 287 288
      allocate(airp_new(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (airp_new)'
      airp_new = _ZERO_

gotm's avatar
gotm committed
289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
      allocate(tausx_old(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (tausx_old)'
      tausx_old = _ZERO_

      allocate(tausy_old(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (tausy_old)'
      tausy_old = _ZERO_

      allocate(d_airp(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (d_airp)'
      d_airp = _ZERO_

      allocate(d_tausx(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (d_tausx)'
      d_tausx = _ZERO_

      allocate(d_tausy(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (d_tausy)'
      d_tausy = _ZERO_

309 310 311
      allocate(tcc_old(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (tcc_old)'
      tcc_old = _ZERO_
gotm's avatar
gotm committed
312

313 314 315 316
      allocate(tcc_new(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (tcc_new)'
      tcc_new = _ZERO_

317 318 319 320
      allocate(swr_old(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (swr_old)'
      swr_old = _ZERO_

gotm's avatar
gotm committed
321 322 323 324
      allocate(shf_old(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (shf_old)'
      shf_old = _ZERO_

325 326 327
      allocate(d_tcc(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (d_tcc)'
      d_tcc = _ZERO_
gotm's avatar
gotm committed
328

329 330 331 332
      allocate(d_swr(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (d_swr)'
      d_swr = _ZERO_

gotm's avatar
gotm committed
333 334 335 336
      allocate(d_shf(E2DFIELD),stat=rc)
      if (rc /= 0) stop 'init_meteo: Error allocating memory (d_shf)'
      d_shf = _ZERO_

337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
      if (fwf_method .ge. 2) then
         allocate(evap_old(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (evap_old)'
         evap_old = _ZERO_

         allocate(d_evap(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (d_evap)'
         d_evap = _ZERO_
      end if

      if (fwf_method .eq. 2 .or. fwf_method .eq. 3) then
         allocate(precip_old(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (precip_old)'
         precip_old = _ZERO_

         allocate(d_precip(E2DFIELD),stat=rc)
         if (rc /= 0) stop 'init_meteo: Error allocating memory (d_precip)'
         d_precip = _ZERO_
      end if

gotm's avatar
gotm committed
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372
   end if
#ifdef DEBUG
   write(debug,*) 'Leaving init_meteo()'
   write(debug,*)
#endif
   return
   end subroutine init_meteo
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: do_meteo - update the meteo forcing
!
! !INTERFACE:
   subroutine do_meteo(n,sst)
bjb's avatar
bjb committed
373
!$ use omp_lib
gotm's avatar
gotm committed
374 375 376 377 378 379 380 381 382 383 384 385 386 387
!
! !DESCRIPTION:
!  Should be called once every time step to update the meteorological forcing.
!  \emph{do\_meteo()} is called with two arguments - n is the loop number and
!  the sea surface temperature.
!  The SST is only used in the case where fluxes and stresses are
!  calculated as part of the model simulation.
!  The forcing can be obtained in 3 different way - using constant values,
!  using pre-calculated stresses and heat-fluxes or by calculating the
!  fluxes as part of the model integration. In all 3 cases the following
!  fields are the result \emph{, airp, tausx, tausy, swr} and \emph{shf} - i.e.
!  air pressure, stresses in x and y direction, short wave radiation and
!  surface heat fluxes.
!  The surface heat flux is the sum of the latent and sensible heat + the
388 389
!  net back radiation. Additional if available the surface freshwater fluxes
!  can be set const or read in from the meteo file. The unit in the meteo file
390
!  is assumed to be meter (/day).
391
!  The structure of this routine looks at first glance a bit more complicated
gotm's avatar
gotm committed
392 393 394 395 396
!  than should be necessary. The main reason is we need two fields in order to
!  do any time interpolation - which explains the use of \emph{first}.
!  In addition checks of the logical \emph{new\_meteo} is checked - set by the
!  reading subroutine.
!  Temporal interpolation is done in the principal variables.
Knut's avatar
Knut committed
397 398
!  It is possible to specify a soft start - via the {\tt meteo\_ramp} variable -
!  which is used to calculate a ramp (linearly from 0 to one over {\tt meteo\_ramp}
gotm's avatar
gotm committed
399 400 401
!  time steps).
!  To implement an use a different set of formulae for flux calculations
!  should be a matter of only changing the involved subroutines.
402
!
Knut's avatar
Knut committed
403
! !USES:
404 405
   use getm_timers, only: tic, toc, TIM_METEO
   IMPLICIT NONE
gotm's avatar
gotm committed
406 407
!
! !INPUT/OUTPUT PARAMETERS:
kbk's avatar
kbk committed
408 409
   integer, intent(in)                 :: n
   REALTYPE, optional, intent(inout)   :: sst(I2DFIELD)
gotm's avatar
gotm committed
410 411 412 413 414
!
! !REVISION HISTORY:
!  See module for log.
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
415 416 417
   integer, save             :: k=0
   integer                   :: i,j
   REALTYPE                  :: ramp,hh,t,t_frac
418
   REALTYPE                  :: solar_zenith_angle
kbk's avatar
kbk committed
419
   REALTYPE                  :: short_wave_radiation
420
   REALTYPE                  :: albedo_water
kbk's avatar
kbk committed
421
   REALTYPE                  :: uu,cosconv,vv,sinconv
422
   REALTYPE, parameter       :: pi=3.1415926535897932384626433832795029d0
423
   REALTYPE, parameter       :: deg2rad=pi/180.
kbk's avatar
kbk committed
424 425
   logical,save              :: first=.true.
   logical                   :: have_sst
gotm's avatar
gotm committed
426 427 428 429 430 431 432 433
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_meteo() # ',Ncall
#endif
bjb's avatar
bjb committed
434
   call tic(TIM_METEO)
435
! OMP-NOTE: In this routine some loops, which have to do with read-in of
bjb's avatar
bjb committed
436 437 438 439 440
!    new meteo have not been threaded due to use of heap-allocated scalars
!    in exchange_coefficients() and fluxes().
!    However, for the cases I have tested, calls to short_wave_radiation
!    is by far the most expensive of the present routine.
!       BJB 2009-09-30.
gotm's avatar
gotm committed
441 442 443 444
   if (metforcing) then

      t = n*timestep

445
      if(meteo_ramp .gt. 0 .and. k .lt. meteo_ramp) then
bjb's avatar
bjb committed
446
! BJB-TODO: Replace 1.0 with _ONE_ etc in this file.
447
         ramp = _ONE_*k/meteo_ramp
gotm's avatar
gotm committed
448 449 450 451 452
         k = k + 1
      else
         ramp = _ONE_
      end if

kbk's avatar
kbk committed
453
      select case (met_method)
gotm's avatar
gotm committed
454
         case (1)
455
! BJB-TODO: Why is this called every time step (even after k=meteo_ramp)-
bjb's avatar
bjb committed
456
!    It should all be constant in time after that(?)
gotm's avatar
gotm committed
457 458 459
            airp  =  _ZERO_
            tausx = ramp*tx
            tausy = ramp*ty
bjb's avatar
bjb committed
460
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j, sinconv,cosconv,uu,vv)
gotm's avatar
gotm committed
461
!     Rotation of wind stress due to grid convergence
bjb's avatar
bjb committed
462
!$OMP DO SCHEDULE(RUNTIME)
bjb's avatar
bjb committed
463
             do j=jmin-1,jmax+1
464
               do i=imin-1,imax+1
kbk's avatar
kbk committed
465 466 467
                  if (convc(i,j) .ne. _ZERO_ .and. az(i,j) .gt. 0) then
                     sinconv=sin(-convc(i,j)*deg2rad)
                     cosconv=cos(-convc(i,j)*deg2rad)
kbk's avatar
kbk committed
468 469 470 471 472 473
                     uu=tausx(i,j)
                     vv=tausy(i,j)
                     tausx(i,j)= uu*cosconv+vv*sinconv
                     tausy(i,j)=-uu*sinconv+vv*cosconv
                  end if
               end do
kbk's avatar
kbk committed
474
            end do
bjb's avatar
bjb committed
475 476
!$OMP END DO
!$OMP END PARALLEL
gotm's avatar
gotm committed
477 478
            swr   = swr_const
            shf   = shf_const
479 480 481 482
            if (fwf_method .eq. 1) then
               evap = evap_const
               precip = precip_const
            end if
gotm's avatar
gotm committed
483 484
         case (2)
            if(calc_met) then
485
               have_sst = present(sst)
gotm's avatar
gotm committed
486
               if (new_meteo) then
487 488 489
                  call update_2d_halo(airp,airp,az, &
                                      imin,jmin,imax,jmax,H_TAG)
                  call wait_halo(H_TAG)
gotm's avatar
gotm committed
490
                  if (.not. first) then
491 492
                     airp_old  = airp_new
                     tcc_old   = tcc_new
gotm's avatar
gotm committed
493 494 495
                     tausx_old = tausx
                     tausy_old = tausy
                     shf_old = shf
496 497 498 499 500 501
                     if (fwf_method .ge. 2) then
                        evap_old = evap
                     end if
                     if (fwf_method .eq. 2 .or. fwf_method .eq. 3) then
                        precip_old = precip
                     end if
gotm's avatar
gotm committed
502
                  end if
503 504
                  airp_new = airp
                  tcc_new  = tcc
kbk's avatar
kbk committed
505
                  if (have_sst) then
bjb's avatar
bjb committed
506
! OMP-NOTE: This is an expensive loop, but we cannot thread it as long
507
!    as exchange_coefficients() and fluxes() pass information through
bjb's avatar
bjb committed
508
!    scalars in the meteo module. BJB 2009-09-30.
509 510 511 512 513
                     do j=jmin,jmax
                        do i=imin,imax
                           if (az(i,j) .ge. 1) then
                              call exchange_coefficients( &
                                     u10(i,j),v10(i,j),t2(i,j),airp(i,j), &
kbk's avatar
kbk committed
514
                                     sst(i,j),hum(i,j),hum_method)
515
                              call fluxes(latc(i,j),u10(i,j),v10(i,j),    &
516 517
                                      t2(i,j),tcc(i,j),sst(i,j),precip(i,j), &
                                      shf(i,j),tausx(i,j),tausy(i,j),evap(i,j))
518
                           else
519
! BJB-TODO: This part of the if-block could be omitted, if the entire fields
520
! are zero'ed out during init (unless az(i,j) may change with time).
521 522 523
                              shf(i,j) = _ZERO_
                              tausx(i,j) = _ZERO_
                              tausy(i,j) = _ZERO_
524 525 526 527
                              if (fwf_method .ge. 1) then
                                 evap(i,j) = _ZERO_
                                 precip(i,j) = _ZERO_
                              end if
528 529 530 531
                           end if
                        end do
                     end do
                  else
bjb's avatar
bjb committed
532
! OMP-NOTE: w needs to be a local (stack) variable to thread this loop.
533 534 535
                     do j=jmin,jmax
                        do i=imin,imax
                           if (az(i,j) .ge. 1) then
bjb's avatar
bjb committed
536
! BJB-TODO: Update constants to double.
537 538 539 540 541 542
                              w=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j))
                              tausx(i,j) = 1.25e-3*1.25*w*U10(i,j)
                              tausy(i,j) = 1.25e-3*1.25*w*V10(i,j)
                           end if
                        end do
                     end do
kbk's avatar
kbk committed
543
                  end if
544 545 546 547 548 549
                  call update_2d_halo(tausx,tausx,az, &
                                      imin,jmin,imax,jmax,H_TAG)
                  call wait_halo(H_TAG)
                  call update_2d_halo(tausy,tausy,az, &
                                      imin,jmin,imax,jmax,H_TAG)
                  call wait_halo(H_TAG)
gotm's avatar
gotm committed
550 551 552
                  if (.not. first) then
                     d_tausx = tausx - tausx_old
                     d_tausy = tausy - tausy_old
553
                     d_airp = airp - airp_old
554
                     d_tcc = tcc - tcc_old
gotm's avatar
gotm committed
555
                     d_shf = shf - shf_old
556 557 558 559 560 561
                     if (fwf_method .ge. 2) then
                        d_evap = evap - evap_old
                     end if
                     if (fwf_method .eq. 2 .or. fwf_method .eq. 3) then
                        d_precip = precip - precip_old
                     end if
gotm's avatar
gotm committed
562 563
                  end if
               end if
bjb's avatar
bjb committed
564
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j, t_frac, hh)
gotm's avatar
gotm committed
565
               if (.not. first) then
bjb's avatar
bjb committed
566 567 568 569
                  t_frac = (t-t_1)/(t_2-t_1)
!$OMP DO SCHEDULE(RUNTIME)
                  do j=jmin-HALO,jmax+HALO
!                     do i=imin-HALO,imax+HALO
570 571
                        airp(:,j)  = airp_old(:,j)   + t_frac*d_airp(:,j)
                        tcc(:,j)   = tcc_old(:,j)    + t_frac*d_tcc(:,j)
bjb's avatar
bjb committed
572 573 574 575 576 577 578 579 580 581 582 583
                        shf(:,j)   = shf_old(:,j)    + t_frac*d_shf(:,j)
                        tausx(:,j) = tausx_old(:,j)  + t_frac*d_tausx(:,j)
                        tausy(:,j) = tausy_old(:,j)  + t_frac*d_tausy(:,j)
                        if (fwf_method .ge. 2) then
                           evap(:,j) = evap_old(:,j) + t_frac*d_evap(:,j)
                        end if
                        if (fwf_method .eq. 2 .or. fwf_method .eq. 3) then
                           precip(:,j) = precip_old(:,j) + t_frac*d_precip(:,j)
                        end if
!                     end do
                  end do
!$OMP END DO
gotm's avatar
gotm committed
584
               end if
585
               hh = secondsofday*(_ONE_/3600)
bjb's avatar
bjb committed
586
!$OMP DO SCHEDULE(RUNTIME)
587 588 589
               do j=jmin,jmax
                  do i=imin,imax
                     if (az(i,j) .ge. 1) then
590 591 592 593 594 595 596 597
                        zenith_angle(i,j) = solar_zenith_angle             &
                                (yearday,hh,lonc(i,j),latc(i,j))
                        swr(i,j) = short_wave_radiation(zenith_angle(i,j), &
                                 yearday,lonc(i,j),latc(i,j),tcc(i,j))
                        if (albedo_method .gt. 0) then
                           albedo(i,j) = albedo_water                      &
                                   (albedo_method,zenith_angle(i,j),yearday)
                        end if
598 599 600
                     end if
                  end do
               end do
bjb's avatar
bjb committed
601
!$OMP END DO
bjb's avatar
bjb committed
602
!$OMP END PARALLEL
gotm's avatar
gotm committed
603
            else
604 605
! BJB-TODO: Need to figure out how/if airp&tcc are correctly computed in
!    the following (.not. calc_met):
606
               if (first) then
bjb's avatar
bjb committed
607
! OMP-NOTE: Dont thread simple memory copy:
608 609 610 611
                  tausx_old = tausx
                  tausy_old = tausy
                  swr_old = swr
                  shf_old = shf
612 613 614 615 616 617
                  if (fwf_method .ge. 2) then
                     evap_old = evap
                  end if
                  if (fwf_method .eq. 2 .or. fwf_method .eq. 3) then
                     precip_old = precip
                  end if
618
               end if
bjb's avatar
bjb committed
619
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j, t_frac)
620
               if (new_meteo) then
bjb's avatar
bjb committed
621 622 623 624 625 626
!$OMP DO SCHEDULE(RUNTIME)
                  do j=jmin-HALO,jmax+HALO
                     tausx_old(:,j) = tausx_old(:,j) + d_tausx(:,j)
                     tausy_old(:,j) = tausy_old(:,j) + d_tausy(:,j)
                     swr_old(:,j) = swr_old(:,j) + d_swr(:,j)
                     shf_old(:,j) = shf_old(:,j) + d_shf(:,j)
627

bjb's avatar
bjb committed
628 629 630 631 632 633 634 635 636 637 638 639 640 641
                     d_tausx(:,j) = tausx(:,j) - tausx_old(:,j)
                     d_tausy(:,j) = tausy(:,j) - tausy_old(:,j)
                     d_swr(:,j) = swr(:,j) - swr_old(:,j)
                     d_shf(:,j) = shf(:,j) - shf_old(:,j)
                     if (fwf_method .ge. 2) then
                        evap_old(:,j) = evap_old(:,j) + d_evap(:,j)
                        d_evap(:,j) = evap(:,j) - evap_old(:,j)
                     end if
                     if (fwf_method .eq. 2 .or. fwf_method .eq. 3) then
                        precip_old(:,j) = precip_old(:,j) + d_precip(:,j)
                        d_precip(:,j) = precip(:,j) - precip_old(:,j)
                     end if
                  end do
!$OMP END DO
642 643 644
               end if
               if (.not. first) then
                  t_frac = (t-t_1)/(t_2-t_1)
bjb's avatar
bjb committed
645 646 647 648 649 650 651 652 653 654 655 656 657 658
!$OMP DO SCHEDULE(RUNTIME)
                  do j=jmin-HALO,jmax+HALO
                     tausx(:,j) = tausx_old(:,j) + t_frac*d_tausx(:,j)
                     tausy(:,j) = tausy_old(:,j) + t_frac*d_tausy(:,j)
                     swr(:,j)   = swr_old(:,j)   + t_frac*d_swr(:,j)
                     shf(:,j)   = shf_old(:,j)   + t_frac*d_shf(:,j)
                     if (fwf_method .ge. 2) then
                        evap(:,j) = evap_old(:,j) + t_frac*d_evap(:,j)
                     end if
                     if (fwf_method .eq. 2 .or. fwf_method .eq. 3) then
                        precip(:,j) = precip_old(:,j) + t_frac*d_precip(:,j)
                     end if
                  end do
!$OMP END DO
659
               end if
bjb's avatar
bjb committed
660
!$OMP END PARALLEL
661
            end if
gotm's avatar
gotm committed
662 663 664 665 666 667 668 669
         case default
            FATAL 'A non valid meteo method has been specified.'
            stop 'do_meteo'
      end select

   end if
   first = .false.

bjb's avatar
bjb committed
670
   call toc(TIM_METEO)
gotm's avatar
gotm committed
671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717
#ifdef DEBUG
     write(debug,*) 'Leaving do_meteo()'
     write(debug,*)
#endif
   return
   end subroutine do_meteo
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: clean_meteo() - cleanup the emph{meteo} module.
!
! !INTERFACE:
   subroutine clean_meteo()
   IMPLICIT NONE
!
! !DESCRIPTION:
!  This routine cleans up the \emph{meteo} module.
!
! !REVISION HISTORY:
!  See module for log.
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'clean_meteo() # ',Ncall
#endif

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

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

   end module meteo

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