calculate_watercycle.f90 8.32 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

  !*****************************************************************************
  !                                                                            *
  !     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
Sabine's avatar
Sabine committed
37
  integer :: itime,i,j,k,jjjjmmdd,ihmmss,numshortout,numshortall
38
  integer :: ix,jy,ixp,jyp,ind,indz,indzp,il,m,indexh,itage,nage
Sabine's avatar
Sabine committed
39
  integer :: interp_time, ngrid
40
  integer :: partnumber,forparticle ! -1 inititalize q, partnumber real do deposition
Sabine's avatar
Sabine committed
41
  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn
42
  real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2)
Sabine's avatar
Sabine committed
43
  real qv1(2),qvprof(2),qvi,dz,dz1,dz2
Sabine's avatar
Sabine committed
44
  real e_minus_p1(2), e_minus_pi
Sabine's avatar
Sabine committed
45
46
47
48
49
50
51
52
53


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

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

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


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

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

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

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

81
82
  ! Take only valid particles, itime invalid when particle released!
  !*****************************************************************
83
  if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999
Sabine's avatar
Sabine committed
84
85


86
87
  ! Determine age class of the particle - nage is used for the kernel
  !******************************************************************
88
89
90
91
92
93
     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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
   ! Determine which nesting level to be used
   !*****************************************
    ngrid=0
    do k=numbnests,1,-1
      if ((xtra1(i).gt.xln(k)).and.(xtra1(i).lt.xrn(k)).and. &
           (ytra1(i).gt.yln(k)).and.(ytra1(i).lt.yrn(k))) then
        ngrid=k
        goto 23
      endif
    end do
23  continue


! Determine nested grid coordinates
!**********************************

    if (ngrid.gt.0) then
      xtn=(xtra1(i)-xln(ngrid))*xresoln(ngrid)
      ytn=(ytra1(i)-yln(ngrid))*yresoln(ngrid)
      ix=int(xtn)
      jy=int(ytn)
      ddy=ytn-real(jy)
      ddx=xtn-real(ix)
    else
      jy=int(ytra1(i))
      ddx=xtra1(i)-real(ix)
      ddy=ytra1(i)-real(jy)
    endif


Sabine's avatar
Sabine committed
124
125
      ixp=ix+1
      jyp=jy+1
126
      if (jy.eq.numygrid) jyp=numygrid
Sabine's avatar
Sabine committed
127
128
129
130
131
132
133
      rddx=1.-ddx
      rddy=1.-ddy
      p1=rddx*rddy
      p2=ddx*rddy
      p3=rddx*ddy
      p4=ddx*ddy

134
135
  ! Interpolate specific humidity
  !******************************
Sabine's avatar
Sabine committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154

      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)

Sabine's avatar
Sabine committed
155
156
          if (ngrid.eq.0) then
             qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
Sabine's avatar
Sabine committed
157
158
159
               +p2*qv(ixp,jy ,ind,indexh) &
               +p3*qv(ix ,jyp,ind,indexh) &
               +p4*qv(ixp,jyp,ind,indexh)
Sabine's avatar
Sabine committed
160
161
162
163
164
165
           else
             qv1(m)=p1*qvn(ix ,jy ,ind,indexh,ngrid) &
               +p2*qvn(ixp,jy ,ind,indexh,ngrid) &
               +p3*qvn(ix ,jyp,ind,indexh,ngrid) &
               +p4*qvn(ixp,jyp,ind,indexh,ngrid)
           endif
Sabine's avatar
Sabine committed
166
167
168
169
170
171
        end do
        qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt

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

Sabine's avatar
Sabine committed
172
173
174
      do m=1,2
         indexh=memind(m)

Sabine's avatar
Sabine committed
175
176
          if (ngrid.eq.0) then
             e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) &
Sabine's avatar
Sabine committed
177
178
179
                +p2*e_minus_p(ixp,jy ,indexh) &
                +p3*e_minus_p(ix ,jyp,indexh) &
                +p4*e_minus_p(ixp,jyp,indexh)
Sabine's avatar
Sabine committed
180
181
182
183
184
185
          else
             e_minus_p1(m)=p1*e_minus_p_nest(ix ,jy ,indexh,ngrid) &
                +p2*e_minus_p_nest(ixp,jy ,indexh,ngrid) &
                +p3*e_minus_p_nest(ix ,jyp,indexh,ngrid) &
                +p4*e_minus_p_nest(ixp,jyp,indexh,ngrid)
          endif
Sabine's avatar
Sabine committed
186
187
      end do
      e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
Sabine's avatar
Sabine committed
188
!      e_minus_pi=e_minus_p1(memtime(1)) !test how the results change without interpolation - does not change a lot
Sabine's avatar
Sabine committed
189
190


191
192
193
194
!      write(*,*) 'all: ',i,diff,e_minus_p(ix,jy) &
!              ,partnumber,val_q(i),qvi,itime,itra1(i),itage

    if (itage.gt.0) then ! not the initialisation stage, it is assumed a value is already saved
195
196
197
        
  ! Calculate humidity differences 
  !*******************************
198

199
       diff=(qvi-val_q(i))*xmass1(i,1)*ldirect
Sabine's avatar
Sabine committed
200
201
202
203

  ! Determine center of mass position on earth and average height
  !**************************************************************
   
204
205
206
207
       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
Sabine's avatar
Sabine committed
208

209
!     Do this only for the 2nd timestep 
210
!     there has to be a deltaq which is negative and the column has to have precipitation
211
!     e_minus_p units is mm/day .. convert to mm/h
212
213
214

! we cannot be sure 3600 is the second timestep, could also be 1800!

Sabine's avatar
Sabine committed
215
      if (((diff.ge.0).or.((e_minus_pi)/24.gt.-2.0)) .and. (itage.le.3600)) then 
216
217
!            write(*,*) 'terminated: ',i,diff,e_minus_p(ix,jy) &
!              ,partnumber,val_q(i),qvi,itime,itra1(i),itage
218
         itra1(i)=-999999999 
219
      else
220
221
!            write(*,*) 'accounted: ',i,diff,e_minus_p(ix,jy) &
!              ,partnumber,val_q(i),qvi,itime,itra1(i),itage
222
223
224
225
226
         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)
Sabine's avatar
Sabine committed
227
228
            if (nested_output.eq.1)  &
              call drydepokernel_nest(1,diff,centerposx,centerposy,nage,1)
229
230
         else
            call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
Sabine's avatar
Sabine committed
231
232
            if (nested_output.eq.1) & 
              call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1) 
233
         endif
234
      endif
235
    endif 
236

237
238
239
240
241
242
243
244
    ! save the old values 
    !********************
  
     xtra1_q(i)=xtra1(i)
     ytra1_q(i)=ytra1(i)
     val_q(i)=qvi
    
  endif ! for valid trajectories
Sabine's avatar
Sabine committed
245
246
247
248
  end do


end subroutine calculate_watercycle