outgrid_init_nest.f90 7.84 KB
Newer Older
Matthias Langer's avatar
   
Matthias Langer committed
1
2
subroutine outgrid_init_nest

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
!*****************************************************************************
!                                                                            *
!  This routine calculates, for each grid cell of the output nest, the       *
!  volume and the surface area.                                              *
!                                                                            *
!     Author: A. Stohl                                                       *
!                                                                            *
!    30 August 2004                                                          *
!                                                                            *
!*****************************************************************************
!                                                                            *
! Variables:                                                                 *
!                                                                            *
! arean              surface area of all output nest cells                   *
! volumen            volumes of all output nest cells                        *
!                                                                            *
!*****************************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

  use unc_mod
  use outg_mod
  use par_mod
  use com_mod

  implicit none

  integer :: ix,jy,kz,ks,kp,nage,l,iix,jjy,ixp,jyp,i1,j1,j,ngrid
  integer :: stat
  real :: ylat,gridarea,ylatp,ylatm,hzone,cosfactm,cosfactp
  real :: xlon,xl,yl,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn,oroh
  real,parameter :: eps=nxmax/3.e5



36
! gridunc,griduncn        uncertainty of outputted concentrations
Matthias Langer's avatar
   
Matthias Langer committed
37
38
39
40
41
  allocate(griduncn(0:numxgridn-1,0:numygridn-1,numzgrid,maxspec, &
       maxpointspec_act,nclassunc,maxageclass),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'

  if (ldirect.gt.0) then
42
43
44
45
46
47
    allocate(wetgriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
         maxpointspec_act,nclassunc,maxageclass),stat=stat)
    if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
    allocate(drygriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
         maxpointspec_act,nclassunc,maxageclass),stat=stat)
    if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
Matthias Langer's avatar
   
Matthias Langer committed
48
49
  endif

50
51
#ifdef USE_MPIINPLACE
#else
52
53
! Extra field for totals at MPI root process
  if (lroot.and.mpi_mode.gt.0) then
54
55
56
57
58
59
60
61
62
63
! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid_nest(),
! then an aux array is needed for parallel grid reduction
    allocate(griduncn0(0:numxgridn-1,0:numygridn-1,numzgrid,maxspec, &
         maxpointspec_act,nclassunc,maxageclass),stat=stat)
    if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
! allocate a dummy to avoid compilator complaints
  else if (.not.lroot.and.mpi_mode.gt.0) then
    allocate(griduncn0(1,1,1,1,1,1,1),stat=stat)
  end if
#endif
64
65
66
67
68
69
70
71
  if (ldirect.gt.0) then
    if (lroot.and.mpi_mode.gt.0) then
      allocate(wetgriduncn0(0:numxgridn-1,0:numygridn-1,maxspec, &
           maxpointspec_act,nclassunc,maxageclass),stat=stat)
      if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
      allocate(drygriduncn0(0:numxgridn-1,0:numygridn-1,maxspec, &
           maxpointspec_act,nclassunc,maxageclass),stat=stat)
      if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
72
!  endif
73
! allocate a dummy to avoid compilator complaints
74
75
76
77
    else if (.not.lroot.and.mpi_mode.gt.0) then
      allocate(wetgriduncn0(1,1,1,1,1,1),stat=stat)
      allocate(drygriduncn0(1,1,1,1,1,1),stat=stat)
    end if
78
79
80
81
82
  end if

! Compute surface area and volume of each grid cell: area, volume;
! and the areas of the northward and eastward facing walls: areaeast, areanorth
!***********************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
83
84
85
86
87
88
89
90
91

  do jy=0,numygridn-1
    ylat=outlat0n+(real(jy)+0.5)*dyoutn
    ylatp=ylat+0.5*dyoutn
    ylatm=ylat-0.5*dyoutn
    if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
      hzone=dyoutn*r_earth*pi180
    else

92
93
94
! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
!************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

      cosfactp=cos(ylatp*pi180)
      cosfactm=cos(ylatm*pi180)
      if (cosfactp.lt.cosfactm) then
        hzone=sqrt(1-cosfactp**2)- &
             sqrt(1-cosfactm**2)
        hzone=hzone*r_earth
      else
        hzone=sqrt(1-cosfactm**2)- &
             sqrt(1-cosfactp**2)
        hzone=hzone*r_earth
      endif
    endif



111
112
! Surface are of a grid cell at a latitude ylat
!**********************************************
Matthias Langer's avatar
   
Matthias Langer committed
113
114
115
116
117
118

    gridarea=2.*pi*r_earth*hzone*dxoutn/360.

    do ix=0,numxgridn-1
      arean(ix,jy)=gridarea

119
120
! Volume = area x box height
!***************************
Matthias Langer's avatar
   
Matthias Langer committed
121
122
123
124
125
126
127
128
129

      volumen(ix,jy,1)=arean(ix,jy)*outheight(1)
      do kz=2,numzgrid
        volumen(ix,jy,kz)=arean(ix,jy)*(outheight(kz)-outheight(kz-1))
      end do
    end do
  end do


130
131
132
!**************************************************************************
! Determine average height of model topography in nesteed output grid cells
!**************************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
133

134
135
! Loop over all output grid cells
!********************************
Matthias Langer's avatar
   
Matthias Langer committed
136
137
138
139
140

  do jjy=0,numygridn-1
    do iix=0,numxgridn-1
      oroh=0.

141
142
! Take 100 samples of the topography in every grid cell
!******************************************************
Matthias Langer's avatar
   
Matthias Langer committed
143
144
145
146
147
148
149
150

      do j1=1,10
        ylat=outlat0n+(real(jjy)+real(j1)/10.-0.05)*dyoutn
        yl=(ylat-ylat0)/dy
        do i1=1,10
          xlon=outlon0n+(real(iix)+real(i1)/10.-0.05)*dxoutn
          xl=(xlon-xlon0)/dx

151
152
! Determine the nest we are in
!*****************************
Matthias Langer's avatar
   
Matthias Langer committed
153
154
155
156
157
158
159
160
161
162
163

          ngrid=0
          do j=numbnests,1,-1
            if ((xl.gt.xln(j)+eps).and.(xl.lt.xrn(j)-eps).and. &
                 (yl.gt.yln(j)+eps).and.(yl.lt.yrn(j)-eps)) then
              ngrid=j
              goto 43
            endif
          end do
43        continue

164
165
! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
!*****************************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202

          if (ngrid.gt.0) then
            xtn=(xl-xln(ngrid))*xresoln(ngrid)
            ytn=(yl-yln(ngrid))*yresoln(ngrid)
            ix=int(xtn)
            jy=int(ytn)
            ddy=ytn-real(jy)
            ddx=xtn-real(ix)
          else
            ix=int(xl)
            jy=int(yl)
            ddy=yl-real(jy)
            ddx=xl-real(ix)
          endif
          ixp=ix+1
          jyp=jy+1
          rddx=1.-ddx
          rddy=1.-ddy
          p1=rddx*rddy
          p2=ddx*rddy
          p3=rddx*ddy
          p4=ddx*ddy

          if (ngrid.gt.0) then
            oroh=oroh+p1*oron(ix ,jy ,ngrid) &
                 + p2*oron(ixp,jy ,ngrid) &
                 + p3*oron(ix ,jyp,ngrid) &
                 + p4*oron(ixp,jyp,ngrid)
          else
            oroh=oroh+p1*oro(ix ,jy) &
                 + p2*oro(ixp,jy) &
                 + p3*oro(ix ,jyp) &
                 + p4*oro(ixp,jyp)
          endif
        end do
      end do

203
204
! Divide by the number of samples taken
!**************************************
Matthias Langer's avatar
   
Matthias Langer committed
205
206
207
208
209
210
211

      orooutn(iix,jjy)=oroh/100.
    end do
  end do



212
213
214
!*******************************
! Initialization of output grids
!*******************************
Matthias Langer's avatar
   
Matthias Langer committed
215
216

  do kp=1,maxpointspec_act
217
218
219
220
221
222
223
224
225
226
227
228
229
230
    do ks=1,nspec
      do nage=1,nageclass
        do jy=0,numygridn-1
          do ix=0,numxgridn-1
            do l=1,nclassunc
! Deposition fields
              if (ldirect.gt.0) then
                wetgriduncn(ix,jy,ks,kp,l,nage)=0.
                drygriduncn(ix,jy,ks,kp,l,nage)=0.
              endif
! Concentration fields
              do kz=1,numzgrid
                griduncn(ix,jy,kz,ks,kp,l,nage)=0.
              end do
Matthias Langer's avatar
   
Matthias Langer committed
231
232
233
234
235
236
237
238
239
            end do
          end do
        end do
      end do
    end do
  end do


end subroutine outgrid_init_nest