calculate_watercycle.f90 7.24 KB
Newer Older
Sabine's avatar
Sabine committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

22
subroutine calculate_watercycle(partnumber,itime)
Sabine's avatar
Sabine committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

  !*****************************************************************************
  !                                                                            *
  !     calculate difference in q and dump it to the dry/wetdepofield          *
  !     VAR: real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q        * 
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod

  implicit none

  real(kind=dp) :: jul
  integer :: itime,i,j,jjjjmmdd,ihmmss,numshortout,numshortall
38
  integer :: ix,jy,ixp,jyp,ind,indz,indzp,il,m,indexh,itage,nage
39
  integer :: interp_time
40
41
  integer :: partnumber,forparticle ! -1 inititalize q, partnumber real do deposition
  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,topo
42
  real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2), prec
Sabine's avatar
Sabine committed
43
44
45
46
47
48
49
50
51
52
  real qv1(2),qvprof(2),qvi,dz,dz1,dz2


  ! Some variables needed for temporal interpolation
  !*************************************************

  dt1=real(itime-memtime(1))
  dt2=real(memtime(2)-itime)
  dtt=1./(dt1+dt2)

53
54
55
56
57
  indexh=memind(2)
  if (abs(memtime(1)-itime).lt.abs(memtime(2)-itime)) &
      indexh=memind(1)


58
  ! if called from timemanager a loop over all particle is necessary
59
60
61
62
  ! if particle is released, the initialization has just to be done 
  ! for the single particle
  !******************************************************************

63
64
  if (partnumber.eq.-1) then
      forparticle=numpart
65
      write (*,*) 'calculate watercycle, ',numpart,partnumber,forparticle,itime
66
67
68
  else
      forparticle=1
  endif
Sabine's avatar
Sabine committed
69

70
 ! Loop about all particles
Sabine's avatar
Sabine committed
71
  !*************************
72
73
74
75
76
77
  do j=1,forparticle

   if (partnumber.eq.-1) then
       i=j ! loop over all particle
   else
       i=partnumber ! just for one particle
78
       firsttimestep(i)=.true.     
79
   endif
Sabine's avatar
Sabine committed
80

81
82
  ! Take only valid particles, itime invalid when particle released!
  !*****************************************************************
Sabine's avatar
Sabine committed
83

84
   if ((itra1(i).eq.itime).or.(forparticle.eq.1)) then
Sabine's avatar
Sabine committed
85
86
87
88
89

  !*****************************************************************************
  ! Interpolate several variables (PV, specific humidity, etc.) to particle position
  !*****************************************************************************

90
91
92
93
94
95
96
97
! Determine age class of the particle - nage is used for the kernel
!******************************************************************
     itage=abs(itra1(i)-itramem(i))
     do nage=1,nageclass
       if (itage.lt.lage(nage)) goto 33
     end do
 33  continue
 
Sabine's avatar
Sabine committed
98
99
100
101
      ix=xtra1(i)
      jy=ytra1(i)
      ixp=ix+1
      jyp=jy+1
102
      if (jy.eq.numygrid) jyp=numygrid
Sabine's avatar
Sabine committed
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
      ddx=xtra1(i)-real(ix)
      ddy=ytra1(i)-real(jy)
      rddx=1.-ddx
      rddy=1.-ddy
      p1=rddx*rddy
      p2=ddx*rddy
      p3=rddx*ddy
      p4=ddx*ddy

  ! Topography
  !***********

      topo=p1*oro(ix ,jy) &
           + p2*oro(ixp,jy) &
           + p3*oro(ix ,jyp) &
           + p4*oro(ixp,jyp)

  ! Potential vorticity, specific humidity, temperature, and density
  !*****************************************************************

      do il=2,nz
        if (height(il).gt.ztra1(i)) then
          indz=il-1
          indzp=il
          goto 6
        endif
      end do
6     continue

      dz1=ztra1(i)-height(indz)
      dz2=height(indzp)-ztra1(i)
      dz=1./(dz1+dz2)


      do ind=indz,indzp
        do m=1,2
          indexh=memind(m)

  ! Specific humidity
          qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
               +p2*qv(ixp,jy ,ind,indexh) &
               +p3*qv(ix ,jyp,ind,indexh) &
               +p4*qv(ixp,jyp,ind,indexh)

        end do
        qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt

      end do
      qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz

153
154
155
156
157
158
159


      prec=p1*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) &
          +p2*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) &
          +p3*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) &
          +p4*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) 

160
      if (partnumber.eq.-1) then ! not the initialisation stage
Sabine's avatar
Sabine committed
161
162
163
164

     ! Calculate humidity differences 
     !*******************************

165
         diff=(qvi-val_q(i))*xmass1(i,1)*ldirect
Sabine's avatar
Sabine committed
166
167
168
169

  ! Determine center of mass position on earth and average height
  !**************************************************************
   
170
171
172
173
174
        xl(1)=xlon0+xtra1_q(i)*dx
        xl(2)=xlon0+xtra1(i)*dx
        yl(1)=ylat0+ytra1_q(i)*dy
        yl(2)=ylat0+ytra1(i)*dy
     endif
Sabine's avatar
Sabine committed
175

176
  ! save the old values - for particle after release - particlenumber = -1
Sabine's avatar
Sabine committed
177
178
  !********************

179
180
181
     xtra1_q(i)=xtra1(i)
     ytra1_q(i)=ytra1(i)
     val_q(i)=qvi
Sabine's avatar
Sabine committed
182
      
183
184
    if (partnumber.eq.-1) then ! not the initialisation stage
        
185
186
!     Do this only for the 2nd timestep i
!     there has to be a deltaq which is negative and the column has to have precipitation
187
188
189
!     e_minus_p units is mm/day .. convert to mm/h
      if ( ((diff.ge.0).or.(e_minus_p(ix,jy)/24.gt.-2.0)) .and. firsttimestep(i)) then 
!      if ( ((diff.ge.0)) .and. firsttimestep(i)) then 
190
         itra1(i)=-999999999 
191
192
193
194
195
196
197
198
199
200
201
202
203
!        write(*,*) 'Trajectory terminated: ',i,diff,e_minus_p(ix,jy),firsttimestep(i)
      else
!        write(*,*) 'Trajectory accounted: ',i,diff,e_minus_p(ix,jy),firsttimestep(i),&
!                          val_q(i),qvi,xmass1(i,1)
         prec_q(i)=prec
         call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
         centerposx=(centerposx_real-xlon0)/dx 
         centerposy=(centerposy_real-ylat0)/dy
         if (diff.gt.0) then
            call drydepokernel(1,diff,centerposx,centerposy,nage,1)
         else
            call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
         endif
204
      endif
205
      firsttimestep(i)=.false.
206
    endif
207

208
  endif
Sabine's avatar
Sabine committed
209
210
211
212
  end do


end subroutine calculate_watercycle