calculate_watercycle.f90 10.2 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,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
44
  real e_minus_p1(2), e_minus_pi, threshold
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
  watersynctime=3600*3 ! was loutstep before, but for monthly run not possible
! watersynctime=3600
! watersynctime=loutstep
Sabine's avatar
Sabine committed
53
  threshold = -.5 ! units mm/h
54

Sabine's avatar
Sabine committed
55
  if ((partnumber.eq.1).and.(itime.eq.0)) write(*,*) 'Watersynctime ',watersynctime,lsynctime, ' s, threshold: ',threshold, 'mm/h'
56

Sabine's avatar
Sabine committed
57
58
59
60
61
62
63
  ! Some variables needed for temporal interpolation
  !*************************************************

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

64
  ! if called from timemanager a loop over all particle is necessary
65
66
67
68
  ! if particle is released, the initialization has just to be done 
  ! for the single particle
  !******************************************************************

69
70
  if (partnumber.eq.-1) then
      forparticle=numpart
71
      totalparticle=0
72
       istep=int(itime/watersynctime)
73
!       write(*,*) 'Model step:',istep
74
75
76
77
78
79
80
81
82
83
84
85
86
87
!      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'
88
89
90
  else
      forparticle=1
  endif
Sabine's avatar
Sabine committed
91

92
  ! Loop about all particles
Sabine's avatar
Sabine committed
93
  !*************************
94
95
96
97
98
99
100
  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
101

102
!   write(*,*) 'particle: ',i,xtra1_q(i),xtra1(i),itime,itra1(i),itramem(i),status_q(i)
103
104
  ! Take only valid particles, itime invalid when particle released!
  !*****************************************************************
105
  if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999
Sabine's avatar
Sabine committed
106

107
108
  ! Determine age class of the particle - nage is used for the kernel
  !******************************************************************
109
110
111
112
113
     itage=abs(itra1(i)-itramem(i))
     do nage=1,nageclass
       if (itage.lt.lage(nage)) goto 33
     end do
 33  continue
114
     
115

Sabine's avatar
Sabine committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
   ! 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
140
      ix=int(xtra1(i))
Sabine's avatar
Sabine committed
141
142
143
144
145
      jy=int(ytra1(i))
      ddx=xtra1(i)-real(ix)
      ddy=ytra1(i)-real(jy)
    endif

Sabine's avatar
Sabine committed
146
147
      ixp=ix+1
      jyp=jy+1
148
      if (jy.eq.numygrid) jyp=numygrid
Sabine's avatar
Sabine committed
149
150
151
152
153
154
155
      rddx=1.-ddx
      rddy=1.-ddy
      p1=rddx*rddy
      p2=ddx*rddy
      p3=rddx*ddy
      p4=ddx*ddy

156
157
  ! Interpolate specific humidity - code from partoutput.f90
  !*********************************************************
Sabine's avatar
Sabine committed
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176

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

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

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

Sabine's avatar
Sabine committed
196
197
198
      do m=1,2
         indexh=memind(m)

Sabine's avatar
Sabine committed
199
200
          if (ngrid.eq.0) then
             e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) &
Sabine's avatar
Sabine committed
201
202
203
                +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
204
205
206
207
208
209
          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
210
211
      end do
      e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
212
      eminusp_saved=e_minus_pi
Sabine's avatar
Sabine committed
213

214
215
  ! Calculate humidity differences 
  !*******************************
216

217
      diff=(qvi-val_q(i))*xmass1(i,1)*ldirect
Sabine's avatar
Sabine committed
218
219
220
221

  ! Determine center of mass position on earth and average height
  !**************************************************************
   
222
223
224
225
       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
226

227
!     Do this only for the 2nd timestep 
228
!     there has to be a deltaq which is negative and the column has to have precipitation
229
!     e_minus_p units is mm/day .. convert to mm/h
230

231
232
233
!     e_minus_pi=precip_reference(ix+1,jy-90+2)*(-24.)

!      First timestep and criteria fullfilled, keep it or throw it?
234
      if ( ((diff.ge.0).or.((e_minus_pi/24).gt.threshold)) .and. (status_q(i).eq.1) ) then 
235
!       if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (status_q(i).eq.1) ) then 
236
!       if ( ((diff.ge.0)) .and. (status_q(i).eq.1) ) then 
237
238
239
240
241
!       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
242
243
244
245
246
247
            totalparticle=totalparticle+1

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

248
!            write(*,*) 'accounted: ',i,itime,itra1(i),itage,diff,itramem(i),status_q(i)
249
250
251
252
            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)
253
            else
254
255
256
               call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
               if (nested_output.eq.1) & 
                 call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1) 
257
            endif
258
259
260
261
!            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)
262
263
264
265
266
            xtra1_q(i)=xtra1(i)
            ytra1_q(i)=ytra1(i)
            val_q(i)=qvi
            status_q(i)=2
       endif
267
268

  endif ! at louttimestep interval    
269
270
271
272
273
274
275
276
277

  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

278
  endif ! for valid trajectories - itra1.eq.itime
279

Sabine's avatar
Sabine committed
280
  end do
281
  if ((partnumber.eq.-1).and.(totalparticle.gt.0)) then
282
!      write (*,*) 'calculated watercycle, ',totalparticle,numpart,partnumber,forparticle,itime
283
  endif
Sabine's avatar
Sabine committed
284
285
286


end subroutine calculate_watercycle