calculate_watercycle.f90 9.85 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,watersynctime)
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,totalparticle ! -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
45
46
  
  real  precip_reference(360,180),eminusp_saved
47
  integer ixr,jyr,idumx,idumy,istep, watersynctime
48

Sabine's avatar
Sabine committed
49
50
51
52
53
54
55
56

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

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

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

62
63
  if (partnumber.eq.-1) then
      forparticle=numpart
64
      totalparticle=0
65
       istep=int(itime/watersynctime)
66
       write(*,*) 'Model step:',istep
67
68
69
70
71
72
73
74
75
76
77
78
79
80
!      open(43,file='precip_ref_input.dat',status='old')
!      write(*,*) 'Reading eminusp from file'
!      do 75 jyr=1,90
!          do 75 ixr=1,360
!75          read(43,'(2i4,e15.3)') idumx,idumy,precip_reference(ixr,jyr)
!      close(43)
!     open(43,file='precip_ref_plus_eminusp.dat')
!     write(*,*) 'Writing eminusp from file'
!     do 76 jyr=1,90
!         do 76 ixr=1,360
!6           write(43,'(2i4,3e15.3)') ixr,jyr,precip_reference(ixr,jyr), &
!                  e_minus_p(ixr,jyr,1), e_minus_p(ixr,jyr+90,1)
!     close(43)
!     write(*,*) 'Written eminusp from file'
81
82
83
  else
      forparticle=1
  endif
Sabine's avatar
Sabine committed
84

85
  ! Loop about all particles
Sabine's avatar
Sabine committed
86
  !*************************
87
88
89
90
91
92
93
  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
94

95
!   write(*,*) 'particle: ',i,xtra1_q(i),xtra1(i),itime,itra1(i),itramem(i),status_q(i)
96
97
  ! Take only valid particles, itime invalid when particle released!
  !*****************************************************************
98
  if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999
Sabine's avatar
Sabine committed
99

100
101
  ! Determine age class of the particle - nage is used for the kernel
  !******************************************************************
102
103
104
105
106
     itage=abs(itra1(i)-itramem(i))
     do nage=1,nageclass
       if (itage.lt.lage(nage)) goto 33
     end do
 33  continue
107
     
108

Sabine's avatar
Sabine committed
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
   ! 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
Sabine's avatar
Sabine committed
133
      ix=int(xtra1(i))
Sabine's avatar
Sabine committed
134
135
136
137
138
      jy=int(ytra1(i))
      ddx=xtra1(i)-real(ix)
      ddy=ytra1(i)-real(jy)
    endif

Sabine's avatar
Sabine committed
139
140
      ixp=ix+1
      jyp=jy+1
141
      if (jy.eq.numygrid) jyp=numygrid
Sabine's avatar
Sabine committed
142
143
144
145
146
147
148
      rddx=1.-ddx
      rddy=1.-ddy
      p1=rddx*rddy
      p2=ddx*rddy
      p3=rddx*ddy
      p4=ddx*ddy

149
150
  ! Interpolate specific humidity - code from partoutput.f90
  !*********************************************************
Sabine's avatar
Sabine committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169

      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
170
171
          if (ngrid.eq.0) then
             qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
Sabine's avatar
Sabine committed
172
173
174
               +p2*qv(ixp,jy ,ind,indexh) &
               +p3*qv(ix ,jyp,ind,indexh) &
               +p4*qv(ixp,jyp,ind,indexh)
Sabine's avatar
Sabine committed
175
176
177
178
179
180
           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
181
182
183
184
185
186
        end do
        qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt

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

187
      if ( (mod(itage,watersynctime).eq.0).and.(itage.gt.0) ) then  !trajectory is loutstep after simulation start, not a new particle
188

Sabine's avatar
Sabine committed
189
190
191
      do m=1,2
         indexh=memind(m)

Sabine's avatar
Sabine committed
192
193
          if (ngrid.eq.0) then
             e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) &
Sabine's avatar
Sabine committed
194
195
196
                +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
197
198
199
200
201
202
          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
203
204
      end do
      e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
205
      eminusp_saved=e_minus_pi
Sabine's avatar
Sabine committed
206

207
208
  ! Calculate humidity differences 
  !*******************************
209

210
      diff=(qvi-val_q(i))*xmass1(i,1)*ldirect
Sabine's avatar
Sabine committed
211
212
213
214

  ! Determine center of mass position on earth and average height
  !**************************************************************
   
215
216
217
218
       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
219

220
!     Do this only for the 2nd timestep 
221
!     there has to be a deltaq which is negative and the column has to have precipitation
222
!     e_minus_p units is mm/day .. convert to mm/h
223

224
225
226
!     e_minus_pi=precip_reference(ix+1,jy-90+2)*(-24.)

!      First timestep and criteria fullfilled, keep it or throw it?
227
228
       if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (status_q(i).eq.1) ) then 
!       if ( ((diff.ge.0)) .and. (status_q(i).eq.1) ) then 
229
230
231
232
233
!       if  (status_q(i).eq.-100) then  ! never! all trajectories.
!          write(*,*) 'terminated: ',i,diff,e_minus_pi,itage,itramem(i),status_q(i),qvi,val_q(i),xmass1(i,1)
          status_q(i)=-9 ! this particle is invalid
          itra1(i)=-999999999
       else
234
235
236
237
238
239
            totalparticle=totalparticle+1

            call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
            centerposx=(centerposx_real-xlon0)/dx 
            centerposy=(centerposy_real-ylat0)/dy

240
!            write(*,*) 'accounted: ',i,itime,itra1(i),itage,diff,itramem(i),status_q(i)
241
242
243
244
            if (diff.gt.0) then
               call drydepokernel(1,diff,centerposx,centerposy,nage,1)
               if (nested_output.eq.1)  &
                 call drydepokernel_nest(1,diff,centerposx,centerposy,nage,1)
245
            else
246
247
248
               call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
               if (nested_output.eq.1) & 
                 call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1) 
249
            endif
250
251
252
253
!            if (istep.lt.100) then
!              waterfieldp(istep,ix+1,jy+1)=waterfieldp(istep,ix+1,jy+1)+xmass1(i,1)
!            endif
!            waterfielde(ix+1,jy+1)=waterfielde(ix+1,jy+1)+xmass1(i,1)
254
255
256
257
258
            xtra1_q(i)=xtra1(i)
            ytra1_q(i)=ytra1(i)
            val_q(i)=qvi
            status_q(i)=2
       endif
259
260

  endif ! at louttimestep interval    
261
262
263
264
265
266
267
268
269

  if (status_q(i).eq.-1) then ! first time, a value is saved 
!        write(*,*) 'saved: ',i,xtra1(i),ytra1(i),qvi,itime,itage,itramem(i),status_q(i)
        xtra1_q(i)=xtra1(i)
        ytra1_q(i)=ytra1(i)
        val_q(i)=qvi
        status_q(i)=1
  endif

270
  endif ! for valid trajectories - itra1.eq.itime
271

Sabine's avatar
Sabine committed
272
  end do
273
  if ((partnumber.eq.-1).and.(totalparticle.gt.0)) then
274
275
      write (*,*) 'calculated watercycle, ',totalparticle,numpart,partnumber,forparticle,itime
  endif
Sabine's avatar
Sabine committed
276
277
278


end subroutine calculate_watercycle