source: branches/flexpart91_hasod/src_parallel/readwind.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 9 years ago

Added parallel version of Flexpart91

File size: 21.8 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine readwind(indj,n,uuh,vvh,wwh)
23
24  !**********************************************************************
25  !                                                                     *
26  !             TRAJECTORY MODEL SUBROUTINE READWIND                    *
27  !                                                                     *
28  !**********************************************************************
29  !                                                                     *
30  !             AUTHOR:      G. WOTAWA                                  *
31  !             DATE:        1997-08-05                                 *
32  !             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
33  !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
34  !                                 ECMWF grib_api                      *
35  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
36  !                                 ECMWF grib_api                      *
37  !                                                                     *
38  !**********************************************************************
39  !  Changes, Bernd C. Krueger, Feb. 2001:
40  !   Variables tth and qvh (on eta coordinates) in common block
41  !**********************************************************************
42  !                                                                     *
43  ! DESCRIPTION:                                                        *
44  !                                                                     *
45  ! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
46  ! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
47  !                                                                     *
48  ! INPUT:                                                              *
49  ! indj               indicates number of the wind field to be read in *
50  ! n                  temporal index for meteorological fields (1 to 3)*
51  !                                                                     *
52  ! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
53  !                                                                     *
54  ! wfname             File name of data to be read in                  *
55  ! nx,ny,nuvz,nwz     expected field dimensions                        *
56  ! nlev_ec            number of vertical levels ecmwf model            *
57  ! uu,vv,ww           wind fields                                      *
58  ! tt,qv              temperature and specific humidity                *
59  ! ps                 surface pressure                                 *
60  !                                                                     *
61  !**********************************************************************
62
63  use GRIB_API
64  use par_mod
65  use com_mod
66
67  implicit none
68
69  !HSO  parameters for grib_api
70  integer :: ifile
71  integer :: iret
72  integer, dimension(:), allocatable   :: igrib
73  integer :: nfield, ii
74  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
75  integer :: gotGrid
76  !HSO  end
77
78  integer :: indj,i,j,k,n,levdiff2,iumax,iwmax
79  integer :: kz
80  real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
81  real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
82  real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
83
84  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
85
86  ! dimension of isec2 at least (22+n), where n is the number of parallels or
87  ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
88
89  ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
90  ! coordinate parameters
91
92  integer :: isec1(56),isec2(22+nxmax+nymax)
93  real(kind=4) :: zsec4(jpunp)
94  real(kind=4) :: xaux,yaux,xaux0,yaux0
95  real(kind=8) :: xauxin,yauxin
96  real,parameter :: eps=1.e-4
97  real(kind=4),allocatable,dimension(:,:) :: nsss,ewss
98  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
99
100  logical :: hflswitch,strswitch
101
102  !HSO  grib api error messages
103  character(len=24) :: gribErrorMsg = 'Error reading grib file'
104  character(len=20) :: gribFunction = 'readwind'
105
106  !HSO conversion of ECMWF etadot to etadot*dp/deta
107  logical :: etacon=.false.
108  real,parameter :: p00=101325.
109  real :: dak,dbk
110  integer :: stat
111
112  allocate(nsss(0:nxmax-1,0:nymax-1),stat=stat)
113  if (stat.ne.0) print*,'ERROR: could not allocate nsss'
114  allocate(ewss(0:nxmax-1,0:nymax-1),stat=stat)
115  if (stat.ne.0) print*,'ERROR: could not allocate ewss'
116
117  hflswitch=.false.
118  strswitch=.false.
119  levdiff2=nlev_ec-nwz+1
120  iumax=0
121  iwmax=0
122
123  !
124  ! OPENING OF DATA FILE (GRIB CODE)
125  !
126  call grib_open_file(ifile,path(3)(1:length(3)) &
127         //trim(wfname(indj)),'r',iret)
128
129  if (iret.ne.GRIB_SUCCESS) then
130    goto 888   ! ERROR DETECTED
131  endif
132
133  call grib_count_in_file(ifile,nfield)
134  allocate(igrib(nfield),stat=stat)
135  igrib(:) = -1
136
137  do ii = 1,nfield
138    call grib_new_from_file(ifile, igrib(ii), iret)
139  end do
140
141  call grib_close_file(ifile)
142
143  !turn on support for multi fields messages */
144  !call grib_multi_support_on
145
146  gotGrid=0
147
148!$OMP PARALLEL DEFAULT(none) &
149!$OMP SHARED (nfield, igrib, gribFunction, nxfield, ny, nlev_ec, dx, xlon0, ylat0, &
150!$OMP   n, tth, uuh, vvh, iumax, qvh, ps, wwh, iwmax, sd, msl, tcc, u10, v10, tt2, &
151!$OMP   td2, lsprec, convprec, sshf, hflswitch, ssr, ewss, nsss, strswitch, oro,   &
152!$OMP   excessoro, lsm, nymin1) &
153!$OMP PRIVATE(ii, gribVer, iret, isec1, discipl, parCat, parNum, typSurf, valSurf, &
154!$OMP   zsec4, isec2, gribErrorMsg, xauxin, yauxin, xaux, yaux, xaux0,  &
155!$OMP   yaux0, k)  &
156!$OMP REDUCTION(+:gotGrid)
157
158#if (defined STATIC_SCHED)
159!$OMP DO SCHEDULE(static)
160#else
161!$OMP DO SCHEDULE(dynamic, max(1,nfield/100))
162#endif
163
164  fieldloop : do ii = 1,nfield
165
166  !first see if we read GRIB1 or GRIB2
167  call grib_get_int(igrib(ii),'editionNumber',gribVer,iret)
168  call grib_check(iret,gribFunction,gribErrorMsg)
169
170  if (gribVer.eq.1) then ! GRIB Edition 1
171
172    !print*,'GRiB Edition 1'
173    !read the grib2 identifiers
174    call grib_get_int(igrib(ii),'indicatorOfParameter',isec1(6),iret)
175    call grib_check(iret,gribFunction,gribErrorMsg)
176    call grib_get_int(igrib(ii),'level',isec1(8),iret)
177    call grib_check(iret,gribFunction,gribErrorMsg)
178 
179    !change code for etadot to code for omega
180    if (isec1(6).eq.77) then
181      isec1(6)=135
182    endif
183
184  else
185
186    !print*,'GRiB Edition 2'
187    !read the grib2 identifiers
188    call grib_get_int(igrib(ii),'discipline',discipl,iret)
189    call grib_check(iret,gribFunction,gribErrorMsg)
190    call grib_get_int(igrib(ii),'parameterCategory',parCat,iret)
191    call grib_check(iret,gribFunction,gribErrorMsg)
192    call grib_get_int(igrib(ii),'parameterNumber',parNum,iret)
193    call grib_check(iret,gribFunction,gribErrorMsg)
194    call grib_get_int(igrib(ii),'typeOfFirstFixedSurface',typSurf,iret)
195    call grib_check(iret,gribFunction,gribErrorMsg)
196    call grib_get_int(igrib(ii),'level',valSurf,iret)
197    call grib_check(iret,gribFunction,gribErrorMsg)
198 
199    !print*,discipl,parCat,parNum,typSurf,valSurf
200 
201    !convert to grib1 identifiers
202    isec1(6)=-1
203    isec1(7)=-1
204    isec1(8)=-1
205    isec1(8)=valSurf     ! level
206    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
207      isec1(6)=130         ! indicatorOfParameter
208    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
209      isec1(6)=131         ! indicatorOfParameter
210    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
211      isec1(6)=132         ! indicatorOfParameter
212    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
213      isec1(6)=133         ! indicatorOfParameter
214    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
215      isec1(6)=134         ! indicatorOfParameter
216    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
217      isec1(6)=135         ! indicatorOfParameter
218    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
219      isec1(6)=151         ! indicatorOfParameter
220    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
221      isec1(6)=165         ! indicatorOfParameter
222    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
223      isec1(6)=166         ! indicatorOfParameter
224    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
225      isec1(6)=167         ! indicatorOfParameter
226    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
227      isec1(6)=168         ! indicatorOfParameter
228    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
229      isec1(6)=141         ! indicatorOfParameter
230    elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
231      isec1(6)=164         ! indicatorOfParameter
232    elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
233      isec1(6)=142         ! indicatorOfParameter
234    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
235      isec1(6)=143         ! indicatorOfParameter
236    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
237      isec1(6)=146         ! indicatorOfParameter
238    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
239      isec1(6)=176         ! indicatorOfParameter
240    elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
241      isec1(6)=180         ! indicatorOfParameter
242    elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
243      isec1(6)=181         ! indicatorOfParameter
244    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
245      isec1(6)=129         ! indicatorOfParameter
246    elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
247      isec1(6)=160         ! indicatorOfParameter
248    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
249         (typSurf.eq.1)) then ! LSM
250      isec1(6)=172         ! indicatorOfParameter
251    else
252      print*,'***ERROR: undefined GRiB2 message found!',discipl, &
253           parCat,parNum,typSurf
254    endif
255
256  endif
257
258  !HSO  get the size and data of the values array
259  if (isec1(6).ne.-1) then
260    call grib_get_real4_array(igrib(ii),'values',zsec4,iret)
261    call grib_check(iret,gribFunction,gribErrorMsg)
262  endif
263
264  !HSO  get the required fields from section 2 in a gribex compatible manner
265  if (ii.eq.1) then
266    call grib_get_int(igrib(ii),'numberOfPointsAlongAParallel',isec2(2),iret)
267    call grib_check(iret,gribFunction,gribErrorMsg)
268    call grib_get_int(igrib(ii),'numberOfPointsAlongAMeridian',isec2(3),iret)
269    call grib_check(iret,gribFunction,gribErrorMsg)
270    call grib_get_int(igrib(ii),'numberOfVerticalCoordinateValues',isec2(12))
271    call grib_check(iret,gribFunction,gribErrorMsg)
272    ! CHECK GRID SPECIFICATIONS
273    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
274    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
275    if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
276  endif
277
278!$OMP CRITICAL
279  !HSO  get the second part of the grid dimensions only from GRiB1 messages
280  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
281    call grib_get_real8(igrib(ii),'longitudeOfFirstGridPointInDegrees',xauxin,iret)
282    call grib_check(iret,gribFunction,gribErrorMsg)
283    call grib_get_real8(igrib(ii),'latitudeOfLastGridPointInDegrees',yauxin,iret)
284    call grib_check(iret,gribFunction,gribErrorMsg)
285    xaux=real(xauxin)+real(nxshift)*dx
286    yaux=real(yauxin)
287    xaux0=xlon0
288    yaux0=ylat0
289    if(xaux.lt.0.) xaux=xaux+360.
290    if(yaux.lt.0.) yaux=yaux+360.
291    if(xaux0.lt.0.) xaux0=xaux0+360.
292    if(yaux0.lt.0.) yaux0=yaux0+360.
293    if(abs(xaux-xaux0).gt.eps) &
294         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
295    if(abs(yaux-yaux0).gt.eps) &
296         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
297    gotGrid=1
298  endif ! gotGrid
299!$OMP END CRITICAL
300
301! ***  copy data to respective variables  ***
302! *******************************************
303
304  k=isec1(8)
305  select case(isec1(6))
306!! TEMPERATURE
307    case(130)
308      do j=0,nymin1
309        do i=0,nxfield-1
310          tth(i,j,nlev_ec-k+2,n) = zsec4(nxfield*(ny-j-1)+i+1)
311        end do
312      end do
313
314!! U VELOCITY
315    case(131)
316      do j=0,nymin1
317        do i=0,nxfield-1
318          uuh(i,j,nlev_ec-k+2) = zsec4(nxfield*(ny-j-1)+i+1)
319        end do
320      end do
321!$OMP CRITICAL
322      iumax=max(iumax,nlev_ec-k+1)
323!$OMP END CRITICAL
324
325!! V VELOCITY
326    case(132)
327      do j=0,nymin1
328        do i=0,nxfield-1
329          vvh(i,j,nlev_ec-k+2) = zsec4(nxfield*(ny-j-1)+i+1)
330        end do
331      end do
332
333!! SPEC. HUMIDITY
334    case(133)
335      do j=0,nymin1
336        do i=0,nxfield-1
337          qvh(i,j,nlev_ec-k+2,n) = zsec4(nxfield*(ny-j-1)+i+1)
338          if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
339               qvh(i,j,nlev_ec-k+2,n) = 0.
340    !        this is necessary because the gridded data may contain
341    !        spurious negative values
342        end do
343      end do
344
345!! SURF. PRESS.
346    case(134)
347      do j=0,nymin1
348        do i=0,nxfield-1
349          ps(i,j,1,n) = zsec4(nxfield*(ny-j-1)+i+1)
350        end do
351      end do
352
353!! W VELOCITY
354    case(135)
355      do j=0,nymin1
356        do i=0,nxfield-1
357          wwh(i,j,nlev_ec-k+1) = zsec4(nxfield*(ny-j-1)+i+1)
358        end do
359      end do
360!$OMP CRITICAL
361      iwmax=max(iwmax,nlev_ec-k+1)
362!$OMP END CRITICAL
363
364!! SNOW DEPTH
365    case(141)
366      do j=0,nymin1
367        do i=0,nxfield-1
368          sd(i,j,1,n)= zsec4(nxfield*(ny-j-1)+i+1)
369        end do
370      end do
371
372!! SEA LEVEL PRESS.
373    case(151)
374      do j=0,nymin1
375        do i=0,nxfield-1
376          msl(i,j,1,n) = zsec4(nxfield*(ny-j-1)+i+1)
377        end do
378      end do
379
380!! CLOUD COVER
381    case(164)
382      do j=0,nymin1
383        do i=0,nxfield-1
384          tcc(i,j,1,n) =  zsec4(nxfield*(ny-j-1)+i+1)
385        end do
386      end do
387
388!! 10 M U VELOCITY
389    case(165)
390      do j=0,nymin1
391        do i=0,nxfield-1
392          u10(i,j,1,n)= zsec4(nxfield*(ny-j-1)+i+1)
393        end do
394      end do
395
396!! 10 M V VELOCITY
397    case(166)
398      do j=0,nymin1
399        do i=0,nxfield-1
400          v10(i,j,1,n) = zsec4(nxfield*(ny-j-1)+i+1)
401        end do
402      end do
403
404!! 2 M TEMPERATURE
405    case(167)
406      do j=0,nymin1
407        do i=0,nxfield-1
408          tt2(i,j,1,n) = zsec4(nxfield*(ny-j-1)+i+1)
409        end do
410      end do
411
412!! 2 M DEW POINT
413    case(168)
414      do j=0,nymin1
415        do i=0,nxfield-1
416          td2(i,j,1,n) = zsec4(nxfield*(ny-j-1)+i+1)
417        end do
418      end do
419
420!! LARGE SCALE PREC.
421    case(142)
422      do j=0,nymin1
423        do i=0,nxfield-1
424          lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
425          if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
426        end do
427      end do
428
429!! CONVECTIVE PREC.
430    case(143)
431      do j=0,nymin1
432        do i=0,nxfield-1
433          convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
434          if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
435        end do
436      end do
437
438!! SENS. HEAT FLUX
439    case(146)
440      do j=0,nymin1
441        do i=0,nxfield-1
442          sshf(i,j,1,n) = zsec4(nxfield*(ny-j-1)+i+1)
443!$OMP CRITICAL
444          if(zsec4(nxfield*(ny-j-1)+i+1).ne.0.) hflswitch=.true.    ! Heat flux available 
445!$OMP END CRITICAL
446        end do
447      end do
448
449!! SOLAR RADIATION
450    case(176)
451      do j=0,nymin1
452        do i=0,nxfield-1
453          ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
454          if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
455        end do
456      end do
457
458!! EW SURFACE STRESS
459    case(180)
460      do j=0,nymin1
461        do i=0,nxfield-1
462          ewss(i,j) = zsec4(nxfield*(ny-j-1)+i+1)
463!$OMP CRITICAL
464          if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) strswitch=.true.    ! stress available
465!$OMP END CRITICAL
466        end do
467      end do
468
469!! NS SURFACE STRESS
470    case(181)
471     do j=0,nymin1
472       do i=0,nxfield-1
473         nsss(i,j) = zsec4(nxfield*(ny-j-1)+i+1)
474!$OMP CRITICAL
475         if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) strswitch=.true.    ! stress available
476!$OMP END CRITICAL
477       end do
478     end do
479
480!! ECMWF OROGRAPHY
481    case(129)
482      do j=0,nymin1
483        do i=0,nxfield-1
484          oro(i,j) = zsec4(nxfield*(ny-j-1)+i+1)/ga
485        end do
486      end do
487
488!! STANDARD DEVIATION OF OROGRAPHY
489    case(160)
490      do j=0,nymin1
491        do i=0,nxfield-1
492          excessoro(i,j) = zsec4(nxfield*(ny-j-1)+i+1)
493        end do
494      end do
495
496!! ECMWF LAND SEA MASK
497    case(172)
498      do j=0,nymin1
499        do i=0,nxfield-1
500          lsm(i,j) = zsec4(nxfield*(ny-j-1)+i+1)
501        end do
502      end do
503
504  end select
505 
506  call grib_release(igrib(ii))
507
508  end do fieldloop
509!$OMP END DO
510!$OMP END PARALLEL
511
512  !
513  deallocate(igrib)
514  !
515
516  !error message if no fields found with correct first longitude in it
517  if (gotGrid.eq.0) then
518    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
519         'messages'
520    stop
521  endif
522
523  if(levdiff2.eq.0) then
524    iwmax=nlev_ec+1
525    do i=0,nxmin1
526      do j=0,nymin1
527        wwh(i,j,nlev_ec+1)=0.
528      end do
529    end do
530  endif
531
532  ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
533  if(etacon.eqv..true.) then
534    do k=1,nwzmax
535      dak=akm(k+1)-akm(k)
536      dbk=bkm(k+1)-bkm(k)
537      do i=0,nxmin1
538        do j=0,nymin1
539          wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
540          if (k.gt.1) then
541            wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
542          endif
543        end do
544      end do
545    end do
546  endif
547
548  ! For global fields, assign the leftmost data column also to the rightmost
549  ! data column; if required, shift whole grid by nxshift grid points
550  !*************************************************************************
551
552  if (xglobal) then
553    call shift_field_0(ewss,nxfield,ny)
554    call shift_field_0(nsss,nxfield,ny)
555    call shift_field_0(oro,nxfield,ny)
556    call shift_field_0(excessoro,nxfield,ny)
557    call shift_field_0(lsm,nxfield,ny)
558    call shift_field(ps,nxfield,ny,1,1,2,n)
559    call shift_field(sd,nxfield,ny,1,1,2,n)
560    call shift_field(msl,nxfield,ny,1,1,2,n)
561    call shift_field(tcc,nxfield,ny,1,1,2,n)
562    call shift_field(u10,nxfield,ny,1,1,2,n)
563    call shift_field(v10,nxfield,ny,1,1,2,n)
564    call shift_field(tt2,nxfield,ny,1,1,2,n)
565    call shift_field(td2,nxfield,ny,1,1,2,n)
566    call shift_field(lsprec,nxfield,ny,1,1,2,n)
567    call shift_field(convprec,nxfield,ny,1,1,2,n)
568    call shift_field(sshf,nxfield,ny,1,1,2,n)
569    call shift_field(ssr,nxfield,ny,1,1,2,n)
570    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
571    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
572    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
573    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
574    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
575  endif
576
577  do i=0,nxmin1
578    do j=0,nymin1
579      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
580    end do
581  end do
582  deallocate(ewss)
583  deallocate(nsss)
584
585  if ((.not.hflswitch).or.(.not.strswitch)) then
586    write(*,*) 'WARNING: No flux data contained in GRIB file ',wfname(indj)
587
588  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
589  ! As ECMWF has increased the model resolution, such that now the first model
590  ! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
591  ! (3rd model level in FLEXPART) for the profile method
592  !***************************************************************************
593
594    do i=0,nxmin1
595      do j=0,nymin1
596        plev1=akz(3)+bkz(3)*ps(i,j,1,n)
597        pmean=0.5*(ps(i,j,1,n)+plev1)
598        tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
599        fu=-r_air*tv/ga/pmean
600        hlev1=fu*(plev1-ps(i,j,1,n))   ! HEIGTH OF FIRST MODEL LAYER
601        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
602        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
603        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
604             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
605             surfstr(i,j,1,n),sshf(i,j,1,n))
606        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
607        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
608      end do
609    end do
610  endif
611
612
613  ! Assign 10 m wind to model level at eta=1.0 to have one additional model
614  ! level at the ground
615  ! Specific humidity is taken the same as at one level above
616  ! Temperature is taken as 2 m temperature
617  !**************************************************************************
618
619     do i=0,nxmin1
620        do j=0,nymin1
621           uuh(i,j,1)=u10(i,j,1,n)
622           vvh(i,j,1)=v10(i,j,1,n)
623           qvh(i,j,1,n)=qvh(i,j,2,n)
624           tth(i,j,1,n)=tt2(i,j,1,n)
625        end do
626     end do
627
628  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
629  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
630
631  return
632888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
633  write(*,*) ' #### ',wfname(indj),'                    #### '
634  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
635  stop 'Execution terminated'
636999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
637  write(*,*) ' #### ',wfname(indj),'                    #### '
638  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
639  stop 'Execution terminated'
640
641end subroutine readwind
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG