Changeset bb579a9 in flexpart.git
- Timestamp:
- Sep 13, 2017, 11:55:16 AM (7 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 08a38b5
- Parents:
- aa8c34a (diff), dd6a4aa (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Files:
-
- 8 added
- 7 deleted
- 28 edited
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
src/FLEXPART.f90
rc8fc724 rd8eed02 33 33 ! * 34 34 !***************************************************************************** 35 ! Changes: * 36 ! Unified ECMWF and GFS builds * 37 ! Marian Harustak, 12.5.2017 * 38 ! - Added detection of metdata format using gributils routines * 39 ! - Distinguished calls to ecmwf/gfs gridcheck versions based on * 40 ! detected metdata format * 41 ! - Passed metdata format down to timemanager * 42 !***************************************************************************** 35 43 ! * 36 44 ! Variables: * … … 46 54 use netcdf_output_mod, only: writeheader_netcdf 47 55 use random_mod, only: gasdev1 56 use class_gribfile 48 57 49 58 implicit none … … 52 61 integer :: idummy = -320 53 62 character(len=256) :: inline_options !pathfile, flexversion, arg2 63 integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN 64 integer :: detectformat 65 54 66 55 67 … … 70 82 ! FLEXPART version string 71 83 flexversion_major = '10' ! Major version number, also used for species file names 72 flexversion='Version '//trim(flexversion_major)//'. 1beta (2016-11-02)'84 flexversion='Version '//trim(flexversion_major)//'.2beta (2017-08-01)' 73 85 verbosity=0 74 86 … … 172 184 call readavailable 173 185 186 ! Detect metdata format 187 !********************** 188 189 metdata_format = detectformat() 190 191 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 192 print *,'ECMWF metdata detected' 193 elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then 194 print *,'NCEP metdata detected' 195 else 196 print *,'Unknown metdata format' 197 return 198 endif 199 200 201 174 202 ! If nested wind fields are used, allocate arrays 175 203 !************************************************ … … 188 216 endif 189 217 190 call gridcheck 218 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 219 call gridcheck_ecmwf 220 else 221 call gridcheck_gfs 222 end if 191 223 192 224 if (verbosity.gt.1) then … … 411 443 endif 412 444 413 call timemanager 445 call timemanager(metdata_format) 414 446 415 447 ! NIK 16.02.2005 -
src/FLEXPART_MPI.f90
r5184a7c rd8eed02 33 33 ! * 34 34 !***************************************************************************** 35 ! Changes: * 36 ! Unified ECMWF and GFS builds * 37 ! Marian Harustak, 12.5.2017 * 38 ! - Added detection of metdata format using gributils routines * 39 ! - Distinguished calls to ecmwf/gfs gridcheck versions based on * 40 ! detected metdata format * 41 ! - Passed metdata format down to timemanager * 42 !***************************************************************************** 35 43 ! * 36 44 ! Variables: * … … 47 55 use netcdf_output_mod, only: writeheader_netcdf 48 56 use random_mod, only: gasdev1 57 use class_gribfile 49 58 50 59 implicit none … … 53 62 integer :: idummy = -320 54 63 character(len=256) :: inline_options !pathfile, flexversion, arg2 64 integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN 65 integer :: detectformat 66 55 67 56 68 … … 80 92 ! FLEXPART version string 81 93 flexversion_major = '10' ! Major version number, also used for species file names 82 flexversion='Ver. '//trim(flexversion_major)//'. 1beta MPI (2016-11-02)'94 flexversion='Ver. '//trim(flexversion_major)//'.2beta MPI (2017-08-01)' 83 95 verbosity=0 84 96 … … 197 209 call readavailable 198 210 211 ! Detect metdata format 212 !********************** 213 214 metdata_format = detectformat() 215 216 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 217 print *,'ECMWF metdata detected' 218 elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then 219 print *,'NCEP metdata detected' 220 else 221 print *,'Unknown metdata format' 222 return 223 endif 224 225 226 199 227 ! If nested wind fields are used, allocate arrays 200 228 !************************************************ … … 213 241 endif 214 242 215 call gridcheck 243 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 244 call gridcheck_ecmwf 245 else 246 call gridcheck_gfs 247 end if 248 216 249 217 250 if (verbosity.gt.1 .and. lroot) then … … 456 489 457 490 458 call timemanager 491 call timemanager(metdata_format) 459 492 460 493 -
src/calcmatrix.f90
re200b7a rd8eed02 20 20 !********************************************************************** 21 21 22 subroutine calcmatrix(lconv,delt,cbmf )22 subroutine calcmatrix(lconv,delt,cbmf,metdata_format) 23 23 ! o i o 24 24 !***************************************************************************** … … 30 30 ! Petra Seibert, Bernd C. Krueger, 2000-2001 * 31 31 ! * 32 !***************************************************************************** 33 ! Changes: * 32 34 ! changed by C. Forster, November 2003 - February 2004 * 33 35 ! array fmassfrac(nconvlevmax,nconvlevmax) represents * 34 36 ! the convective redistribution matrix for the particles * 35 37 ! * 38 ! Unified ECMWF and GFS builds * 39 ! Marian Harustak, 12.5.2017 * 40 ! - Merged calcmatrix and calcmatrix_gfs into one routine using if-then * 41 ! for meteo-type dependent code * 42 !***************************************************************************** 43 ! * 36 44 ! lconv indicates whether there is convection in this cell, or not * 37 45 ! delt time step for convection [s] * 38 46 ! cbmf cloud base mass flux * 47 ! metdata_format format of metdata (ecmwf/gfs) * 39 48 ! * 40 49 !***************************************************************************** … … 43 52 use com_mod 44 53 use conv_mod 54 use class_gribfile 45 55 46 56 implicit none 47 57 48 58 real :: rlevmass,summe 59 integer :: metdata_format 49 60 50 61 integer :: iflag, k, kk, kuvz … … 77 88 do kuvz = 2,nuvz 78 89 k = kuvz-1 90 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 79 91 pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv) 80 92 phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv) 93 else 94 phconv(kuvz) = 0.5*(pconv(kuvz)+pconv(k)) 95 endif 81 96 dpr(k) = phconv(k) - phconv(kuvz) 82 97 qsconv(k) = f_qvsat( pconv(k), tconv(k) ) … … 85 100 do kk=1,nconvlev 86 101 fmassfrac(k,kk)=0. 87 end do88 end do102 end do 103 end do 89 104 90 105 -
src/calcpar.f90
re200b7a r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine calcpar(n,uuh,vvh,pvh )22 subroutine calcpar(n,uuh,vvh,pvh,metdata_format) 23 23 ! i i i o 24 24 !***************************************************************************** … … 37 37 ! new variables in call to richardson * 38 38 ! * 39 !***************************************************************************** 40 ! Changes, Bernd C. Krueger, Feb. 2001: 41 ! Variables tth and qvh (on eta coordinates) in common block 39 ! CHANGE 17/11/2005 Caroline Forster NCEP GFS version * 40 ! * 41 ! Changes, Bernd C. Krueger, Feb. 2001: * 42 ! Variables tth and qvh (on eta coordinates) in common block * 43 ! * 44 ! Unified ECMWF and GFS builds * 45 ! Marian Harustak, 12.5.2017 * 46 ! - Merged calcpar and calcpar_gfs into one routine using if-then * 47 ! for meteo-type dependent code * 48 !***************************************************************************** 49 42 50 !***************************************************************************** 43 51 ! * 44 52 ! Variables: * 45 53 ! n temporal index for meteorological fields (1 to 3) * 54 ! uuh * 55 ! vvh * 56 ! pvh * 57 ! metdata_format format of metdata (ecmwf/gfs) * 46 58 ! * 47 59 ! Constants: * … … 56 68 use par_mod 57 69 use com_mod 70 use class_gribfile 58 71 59 72 implicit none 60 73 61 integer :: n,ix,jy,i,kz,lz,kzmin 74 integer :: metdata_format 75 integer :: n,ix,jy,i,kz,lz,kzmin,llev,loop_start 62 76 real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus 63 77 real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat 64 real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax) 78 real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy 65 79 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) 66 80 real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax) … … 111 125 !*********************************************** 112 126 127 if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then 128 ! NCEP version: find first level above ground 129 llev = 0 130 do i=1,nuvz 131 if (ps(ix,jy,1,n).lt.akz(i)) llev=i 132 end do 133 llev = llev+1 134 if (llev.gt.nuvz) llev = nuvz-1 135 ! NCEP version 136 137 ! calculate inverse Obukhov length scale with tth(llev) 113 138 ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), & 114 tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm) 139 tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format) 140 else 141 llev=0 142 ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), & 143 tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format) 144 end if 145 115 146 if (ol.ne.0.) then 116 147 oli(ix,jy,1,n)=1./ol … … 130 161 end do 131 162 163 if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then 164 ! NCEP version hmix has been read in in readwind.f, is therefore not calculated here 132 165 call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, & 133 166 ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), & 134 td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus) 167 td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format) 168 else 169 call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, & 170 ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), & 171 td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus,metdata_format) 172 end if 135 173 136 174 if(lsubgrid.eq.1) then … … 173 211 !****************************************************** 174 212 175 ! 1) Calculate altitudes of ECMWFmodel levels176 !*************************************** ******213 ! 1) Calculate altitudes of model levels 214 !*************************************** 177 215 178 216 tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ & … … 180 218 pold=ps(ix,jy,1,n) 181 219 zold=0. 182 do kz=2,nuvz 220 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 221 loop_start=2 222 else 223 loop_start=llev 224 end if 225 do kz=loop_start,nuvz 183 226 pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers 184 227 tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n)) … … 199 242 !************************************************************************ 200 243 201 do kz=1,nuvz 244 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 245 loop_start=1 246 else 247 loop_start=llev 248 end if 249 250 do kz=loop_start,nuvz 202 251 if (zlev(kz).ge.altmin) then 203 252 kzmin=kz -
src/calcpar_nests.f90
rdb712a8 r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine calcpar_nests(n,uuhn,vvhn,pvhn )22 subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format) 23 23 ! i i i o 24 24 !***************************************************************************** … … 39 39 ! * 40 40 !***************************************************************************** 41 ! Changes, Bernd C. Krueger, Feb. 2001: 42 ! Variables tth and qvh (on eta coordinates) in common block 41 ! Changes, Bernd C. Krueger, Feb. 2001: * 42 ! Variables tth and qvh (on eta coordinates) in common block * 43 ! * 44 ! Unified ECMWF and GFS builds * 45 ! Marian Harustak, 12.5.2017 * 46 ! - Added passing of metdata_format as it was needed by called routines * 43 47 !***************************************************************************** 44 48 ! * 45 49 ! Variables: * 46 50 ! n temporal index for meteorological fields (1 to 3) * 51 ! metdata_format format of metdata (ecmwf/gfs) * 47 52 ! * 48 53 ! Constants: * … … 60 65 implicit none 61 66 67 integer :: metdata_format 62 68 integer :: n,ix,jy,i,l,kz,lz,kzmin 63 real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus 69 real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus,dummyakzllev 64 70 real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat 65 71 real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax) … … 110 116 ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), & 111 117 td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), & 112 sshfn(ix,jy,1,n,l),akm,bkm )118 sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev,metdata_format) 113 119 if (ol.ne.0.) then 114 120 olin(ix,jy,1,n,l)=1./ol … … 131 137 qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), & 132 138 tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), & 133 wstarn(ix,jy,1,n,l),hmixplus )139 wstarn(ix,jy,1,n,l),hmixplus,metdata_format) 134 140 135 141 if(lsubgrid.eq.1) then -
src/convmix.f90
r8a65cb0 r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine convmix(itime )22 subroutine convmix(itime,metdata_format) 23 23 ! i 24 24 !************************************************************** … … 31 31 !CHANGES by A. Stohl: 32 32 ! various run-time optimizations - February 2005 33 ! CHANGES by C. Forster, November 2005, NCEP GFS version 34 ! in the ECMWF version convection is calculated on the 35 ! original eta-levels 36 ! in the GFS version convection is calculated on the 37 ! FLEXPART levels 38 ! 39 ! Unified ECMWF and GFS builds 40 ! Marian Harustak, 12.5.2017 41 ! - Merged convmix and convmix_gfs into one routine using if-then 42 ! for meteo-type dependent code 33 43 !************************************************************** 34 44 … … 37 47 use com_mod 38 48 use conv_mod 49 use class_gribfile 39 50 40 51 implicit none … … 44 55 integer :: jy, kpart, ktop, ngrid,kz 45 56 integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests) 57 integer :: metdata_format 46 58 47 59 ! itime [s] current time … … 104 116 105 117 ngrid=0 118 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 106 119 do j=numbnests,1,-1 107 120 if ( x.gt.xln(j)+eps .and. x.lt.xrn(j)-eps .and. & … … 111 124 endif 112 125 end do 126 else 127 do j=numbnests,1,-1 128 if ( x.gt.xln(j) .and. x.lt.xrn(j) .and. & 129 y.gt.yln(j) .and. y.lt.yrn(j) ) then 130 ngrid=j 131 goto 23 132 endif 133 end do 134 endif 113 135 23 continue 114 136 … … 167 189 td2conv=(td2(ix,jy,1,mind1)*dt2+td2(ix,jy,1,mind2)*dt1)*dtt 168 190 !!$ do kz=1,nconvlev+1 !old 191 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 169 192 do kz=1,nuvz-1 !bugfix 170 193 tconv(kz)=(tth(ix,jy,kz+1,mind1)*dt2+ & … … 173 196 qvh(ix,jy,kz+1,mind2)*dt1)*dtt 174 197 end do 198 else 199 do kz=1,nuvz-1 !bugfix 200 pconv(kz)=(pplev(ix,jy,kz,mind1)*dt2+ & 201 pplev(ix,jy,kz,mind2)*dt1)*dtt 202 tconv(kz)=(tt(ix,jy,kz,mind1)*dt2+ & 203 tt(ix,jy,kz,mind2)*dt1)*dtt 204 qconv(kz)=(qv(ix,jy,kz,mind1)*dt2+ & 205 qv(ix,jy,kz,mind2)*dt1)*dtt 206 end do 207 end if 175 208 176 209 ! Calculate translocation matrix 177 call calcmatrix(lconv,delt,cbaseflux(ix,jy) )210 call calcmatrix(lconv,delt,cbaseflux(ix,jy),metdata_format) 178 211 igrold = igr 179 212 ktop = 0 … … 252 285 ! calculate translocation matrix 253 286 !******************************* 254 call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest) )287 call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest),metdata_format) 255 288 igrold = igr 256 289 ktop = 0 -
src/getfields.f90
rd1a8707 r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine getfields(itime,nstop )22 subroutine getfields(itime,nstop,metdata_format) 23 23 ! i o 24 24 !***************************************************************************** … … 38 38 ! * 39 39 !***************************************************************************** 40 ! Changes, Bernd C. Krueger, Feb. 2001: 41 ! Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block. 42 ! Function of nstop extended. 40 ! Changes, Bernd C. Krueger, Feb. 2001: * 41 ! Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block. * 42 ! Function of nstop extended. * 43 ! * 44 ! Unified ECMWF and GFS builds * 45 ! Marian Harustak, 12.5.2017 * 46 ! - Added passing of metdata_format as it was needed by called routines * 43 47 !***************************************************************************** 44 48 ! * … … 58 62 ! tt(0:nxmax,0:nymax,nuvzmax,2) temperature [K] * 59 63 ! ps(0:nxmax,0:nymax,2) surface pressure [Pa] * 64 ! metdata_format format of metdata (ecmwf/gfs) * 60 65 ! * 61 66 ! Constants: * … … 67 72 use par_mod 68 73 use com_mod 74 use class_gribfile 69 75 70 76 implicit none 71 77 72 78 integer :: indj,itime,nstop,memaux 79 integer :: metdata_format 73 80 74 81 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) … … 125 132 do indj=indmin,numbwf-1 126 133 if (ldirect*wftime(indj+1).gt.ldirect*itime) then 127 call readwind(indj+1,memind(2),uuh,vvh,wwh) 134 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 135 call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh) 136 else 137 call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh) 138 end if 128 139 call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) 129 call calcpar(memind(2),uuh,vvh,pvh) 130 call calcpar_nests(memind(2),uuhn,vvhn,pvhn) 131 call verttransform(memind(2),uuh,vvh,wwh,pvh) 140 call calcpar(memind(2),uuh,vvh,pvh,metdata_format) 141 call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format) 142 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 143 call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh) 144 else 145 call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh) 146 end if 132 147 call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) 133 148 memtime(2)=wftime(indj+1) … … 152 167 (ldirect*wftime(indj+1).gt.ldirect*itime)) then 153 168 memind(1)=1 154 call readwind(indj,memind(1),uuh,vvh,wwh) 169 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 170 call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh) 171 else 172 call readwind_gfs(indj,memind(1),uuh,vvh,wwh) 173 end if 155 174 call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn) 156 call calcpar(memind(1),uuh,vvh,pvh) 157 call calcpar_nests(memind(1),uuhn,vvhn,pvhn) 158 call verttransform(memind(1),uuh,vvh,wwh,pvh) 175 call calcpar(memind(1),uuh,vvh,pvh,metdata_format) 176 call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format) 177 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 178 call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh) 179 else 180 call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh) 181 end if 159 182 call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn) 160 183 memtime(1)=wftime(indj) 161 184 memind(2)=2 162 call readwind(indj+1,memind(2),uuh,vvh,wwh) 185 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 186 call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh) 187 else 188 call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh) 189 end if 163 190 call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) 164 call calcpar(memind(2),uuh,vvh,pvh) 165 call calcpar_nests(memind(2),uuhn,vvhn,pvhn) 166 call verttransform(memind(2),uuh,vvh,wwh,pvh) 191 call calcpar(memind(2),uuh,vvh,pvh,metdata_format) 192 call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format) 193 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 194 call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh) 195 else 196 call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh) 197 end if 167 198 call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) 168 199 memtime(2)=wftime(indj+1) -
src/getfields_mpi.f90
r78e62dc r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine getfields(itime,nstop,memstat )22 subroutine getfields(itime,nstop,memstat,metdata_format) 23 23 ! i o o 24 24 !***************************************************************************** … … 58 58 ! memstat=0: no new fields to be read 59 59 ! 60 ! Unified ECMWF and GFS builds 61 ! Marian Harustak, 12.5.2017 62 ! - Added passing of metdata_format as it was needed by called routines 63 ! 60 64 !***************************************************************************** 61 65 ! * … … 77 81 ! tt(0:nxmax,0:nymax,nuvzmax,numwfmem) temperature [K] * 78 82 ! ps(0:nxmax,0:nymax,numwfmem) surface pressure [Pa] * 83 ! metdata_format format of metdata (ecmwf/gfs) * 79 84 ! * 80 85 ! Constants: * … … 87 92 use com_mod 88 93 use mpi_mod, only: lmpreader,lmp_use_reader,lmp_sync 94 use class_gribfile 89 95 90 96 implicit none 91 97 98 integer :: metdata_format 92 99 integer :: indj,itime,nstop,memaux,mindread 93 100 ! mindread: index of where to read 3rd field … … 204 211 if (ldirect*wftime(indj+1).gt.ldirect*itime) then 205 212 if (lmpreader.or..not. lmp_use_reader) then 206 call readwind(indj+1,mindread,uuh,vvh,wwh) 213 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 214 call readwind_ecmwf(indj+1,mindread,uuh,vvh,wwh) 215 else 216 call readwind_gfs(indj+1,mindread,uuh,vvh,wwh) 217 end if 207 218 call readwind_nests(indj+1,mindread,uuhn,vvhn,wwhn) 208 call calcpar(mindread,uuh,vvh,pvh) 209 call calcpar_nests(mindread,uuhn,vvhn,pvhn) 210 call verttransform(mindread,uuh,vvh,wwh,pvh) 219 call calcpar(mindread,uuh,vvh,pvh,metdata_format) 220 call calcpar_nests(mindread,uuhn,vvhn,pvhn,metdata_format) 221 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 222 call verttransform_ecmwf(mindread,uuh,vvh,wwh,pvh) 223 else 224 call verttransform_gfs(mindread,uuh,vvh,wwh,pvh) 225 end if 211 226 call verttransform_nests(mindread,uuhn,vvhn,wwhn,pvhn) 212 227 end if … … 230 245 memind(1)=1 231 246 if (lmpreader.or..not.lmp_use_reader) then 232 call readwind(indj,memind(1),uuh,vvh,wwh) 247 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 248 call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh) 249 else 250 call readwind_gfs(indj,memind(1),uuh,vvh,wwh) 251 end if 233 252 call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn) 234 call calcpar(memind(1),uuh,vvh,pvh) 235 call calcpar_nests(memind(1),uuhn,vvhn,pvhn) 236 call verttransform(memind(1),uuh,vvh,wwh,pvh) 253 call calcpar(memind(1),uuh,vvh,pvh,metdata_format) 254 call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format) 255 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 256 call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh) 257 else 258 call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh) 259 end if 237 260 call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn) 238 261 end if … … 240 263 memind(2)=2 241 264 if (lmpreader.or..not.lmp_use_reader) then 242 call readwind(indj+1,memind(2),uuh,vvh,wwh) 265 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 266 call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh) 267 else 268 call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh) 269 end if 243 270 call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) 244 call calcpar(memind(2),uuh,vvh,pvh) 245 call calcpar_nests(memind(2),uuhn,vvhn,pvhn) 246 call verttransform(memind(2),uuh,vvh,wwh,pvh) 271 call calcpar(memind(2),uuh,vvh,pvh,metdata_format) 272 call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format) 273 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 274 call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh) 275 else 276 call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh) 277 end if 247 278 call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) 248 279 end if -
src/gridcheck_ecmwf.f90
rdb712a8 r61e07ba 20 20 !********************************************************************** 21 21 22 subroutine gridcheck 22 subroutine gridcheck_ecmwf 23 23 24 24 !********************************************************************** … … 37 37 ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * 38 38 ! ECMWF grib_api * 39 ! * 40 ! Unified ECMWF and GFS builds * 41 ! Marian Harustak, 12.5.2017 * 42 ! - Renamed from gridcheck to gridcheck_ecmwf * 39 43 ! * 40 44 !********************************************************************** … … 117 121 ! OPENING OF DATA FILE (GRIB CODE) 118 122 ! 119 5 120 123 5 call grib_open_file(ifile,path(3)(1:length(3)) & 124 //trim(wfname(ifn)),'r',iret) 121 125 if (iret.ne.GRIB_SUCCESS) then 122 126 goto 999 ! ERROR DETECTED … … 127 131 gotGrid=0 128 132 ifield=0 129 10 133 10 ifield=ifield+1 130 134 131 135 ! … … 145 149 if (gribVer.eq.1) then ! GRIB Edition 1 146 150 147 !print*,'GRiB Edition 1'148 !read the grib2 identifiers149 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)150 call grib_check(iret,gribFunction,gribErrorMsg)151 call grib_get_int(igrib,'level',isec1(8),iret)152 call grib_check(iret,gribFunction,gribErrorMsg)153 154 !change code for etadot to code for omega155 if (isec1(6).eq.77) then156 isec1(6)=135157 endif158 159 !print*,isec1(6),isec1(8)151 !print*,'GRiB Edition 1' 152 !read the grib2 identifiers 153 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 154 call grib_check(iret,gribFunction,gribErrorMsg) 155 call grib_get_int(igrib,'level',isec1(8),iret) 156 call grib_check(iret,gribFunction,gribErrorMsg) 157 158 !change code for etadot to code for omega 159 if (isec1(6).eq.77) then 160 isec1(6)=135 161 endif 162 163 !print*,isec1(6),isec1(8) 160 164 161 165 else 162 166 163 !print*,'GRiB Edition 2' 164 !read the grib2 identifiers 165 call grib_get_int(igrib,'discipline',discipl,iret) 166 call grib_check(iret,gribFunction,gribErrorMsg) 167 call grib_get_int(igrib,'parameterCategory',parCat,iret) 168 call grib_check(iret,gribFunction,gribErrorMsg) 169 call grib_get_int(igrib,'parameterNumber',parNum,iret) 170 call grib_check(iret,gribFunction,gribErrorMsg) 171 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 172 call grib_check(iret,gribFunction,gribErrorMsg) 173 call grib_get_int(igrib,'level',valSurf,iret) 174 call grib_check(iret,gribFunction,gribErrorMsg) 175 call grib_get_int(igrib,'paramId',parId,iret) 176 call grib_check(iret,gribFunction,gribErrorMsg) 177 178 !print*,discipl,parCat,parNum,typSurf,valSurf 179 180 !convert to grib1 identifiers 181 isec1(6)=-1 182 isec1(7)=-1 183 isec1(8)=-1 184 isec1(8)=valSurf ! level 185 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 186 isec1(6)=130 ! indicatorOfParameter 187 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 188 isec1(6)=131 ! indicatorOfParameter 189 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 190 isec1(6)=132 ! indicatorOfParameter 191 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 192 isec1(6)=133 ! indicatorOfParameter 193 !ZHG FOR CLOUDS FROM GRIB 194 elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc 195 isec1(6)=246 ! indicatorOfParameter 196 elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc 197 isec1(6)=247 ! indicatorOfParameter 198 !ZHG end 199 ! ESO qc(=clwc+ciwc) 200 elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc 201 isec1(6)=201031 ! indicatorOfParameter 202 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 203 isec1(6)=134 ! indicatorOfParameter 204 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 205 isec1(6)=135 ! indicatorOfParameter 206 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot 207 isec1(6)=135 ! indicatorOfParameter 208 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 209 isec1(6)=151 ! indicatorOfParameter 210 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 211 isec1(6)=165 ! indicatorOfParameter 212 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 213 isec1(6)=166 ! indicatorOfParameter 214 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 215 isec1(6)=167 ! indicatorOfParameter 216 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 217 isec1(6)=168 ! indicatorOfParameter 218 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 219 isec1(6)=141 ! indicatorOfParameter 220 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 221 isec1(6)=164 ! indicatorOfParameter 222 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 223 isec1(6)=142 ! indicatorOfParameter 224 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 225 isec1(6)=143 ! indicatorOfParameter 226 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 227 isec1(6)=146 ! indicatorOfParameter 228 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 229 isec1(6)=176 ! indicatorOfParameter 230 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 231 isec1(6)=180 ! indicatorOfParameter 232 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 233 isec1(6)=181 ! indicatorOfParameter 234 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 235 isec1(6)=129 ! indicatorOfParameter 236 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 237 isec1(6)=160 ! indicatorOfParameter 238 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 239 (typSurf.eq.1)) then ! LSM 240 isec1(6)=172 ! indicatorOfParameter 241 else 242 print*,'***ERROR: undefined GRiB2 message found!',discipl, & 243 parCat,parNum,typSurf 244 endif 245 if(parId .ne. isec1(6) .and. parId .ne. 77) then 246 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 247 ! stop 248 endif 249 250 endif 167 !print*,'GRiB Edition 2' 168 !read the grib2 identifiers 169 call grib_get_int(igrib,'discipline',discipl,iret) 170 call grib_check(iret,gribFunction,gribErrorMsg) 171 call grib_get_int(igrib,'parameterCategory',parCat,iret) 172 call grib_check(iret,gribFunction,gribErrorMsg) 173 call grib_get_int(igrib,'parameterNumber',parNum,iret) 174 call grib_check(iret,gribFunction,gribErrorMsg) 175 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 176 call grib_check(iret,gribFunction,gribErrorMsg) 177 call grib_get_int(igrib,'level',valSurf,iret) 178 call grib_check(iret,gribFunction,gribErrorMsg) 179 call grib_get_int(igrib,'paramId',parId,iret) 180 call grib_check(iret,gribFunction,gribErrorMsg) 181 182 !print*,discipl,parCat,parNum,typSurf,valSurf 183 184 !convert to grib1 identifiers 185 isec1(6)=-1 186 isec1(7)=-1 187 isec1(8)=-1 188 isec1(8)=valSurf ! level 189 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 190 isec1(6)=130 ! indicatorOfParameter 191 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 192 isec1(6)=131 ! indicatorOfParameter 193 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 194 isec1(6)=132 ! indicatorOfParameter 195 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 196 isec1(6)=133 ! indicatorOfParameter 197 !ZHG FOR CLOUDS FROM GRIB 198 elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc 199 isec1(6)=246 ! indicatorOfParameter 200 elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc 201 isec1(6)=247 ! indicatorOfParameter 202 !ZHG end 203 ! ESO qc(=clwc+ciwc) 204 elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc 205 isec1(6)=201031 ! indicatorOfParameter 206 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 207 isec1(6)=134 ! indicatorOfParameter 208 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 209 isec1(6)=135 ! indicatorOfParameter 210 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot 211 isec1(6)=135 ! indicatorOfParameter 212 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 213 isec1(6)=151 ! indicatorOfParameter 214 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 215 isec1(6)=165 ! indicatorOfParameter 216 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 217 isec1(6)=166 ! indicatorOfParameter 218 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 219 isec1(6)=167 ! indicatorOfParameter 220 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 221 isec1(6)=168 ! indicatorOfParameter 222 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 223 isec1(6)=141 ! indicatorOfParameter 224 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 225 isec1(6)=164 ! indicatorOfParameter 226 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 227 isec1(6)=142 ! indicatorOfParameter 228 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 229 isec1(6)=143 ! indicatorOfParameter 230 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 231 isec1(6)=146 ! indicatorOfParameter 232 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 233 isec1(6)=176 ! indicatorOfParameter 234 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 235 isec1(6)=180 ! indicatorOfParameter 236 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 237 isec1(6)=181 ! indicatorOfParameter 238 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 239 isec1(6)=129 ! indicatorOfParameter 240 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 241 isec1(6)=160 ! indicatorOfParameter 242 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 243 (typSurf.eq.1)) then ! LSM 244 isec1(6)=172 ! indicatorOfParameter 245 else 246 print*,'***ERROR: undefined GRiB2 message found!',discipl, & 247 parCat,parNum,typSurf 248 endif 249 if(parId .ne. isec1(6) .and. parId .ne. 77) then 250 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 251 ! stop 252 endif 253 254 endif 255 256 CALL grib_get_int(igrib,'numberOfPointsAlongAParallel', & 257 isec2(2),iret) 258 ! nx=isec2(2) 259 ! WRITE(*,*) nx,nxmax 260 IF (isec2(2).GT.nxmax) THEN 261 WRITE(*,*) 'FLEXPART error: Too many grid points in x direction.' 262 WRITE(*,*) 'Reduce resolution of wind fields.' 263 WRITE(*,*) 'Or change parameter settings in file ecmwf_mod.' 264 WRITE(*,*) nx,nxmax 265 ! STOP 266 ENDIF 251 267 252 268 !get the size and data of the values array … … 258 274 if (ifield.eq.1) then 259 275 260 !HSO get the required fields from section 2 in a gribex compatible manner276 !HSO get the required fields from section 2 in a gribex compatible manner 261 277 call grib_get_int(igrib,'numberOfPointsAlongAParallel', & 262 278 isec2(2),iret) … … 272 288 call grib_check(iret,gribFunction,gribErrorMsg) 273 289 274 ! get the size and data of the vertical coordinate array290 ! get the size and data of the vertical coordinate array 275 291 call grib_get_real4_array(igrib,'pv',zsec2,iret) 276 292 call grib_check(iret,gribFunction,gribErrorMsg) … … 279 295 ny=isec2(3) 280 296 nlev_ec=isec2(12)/2-1 281 297 endif 282 298 283 299 !HSO get the second part of the grid dimensions only from GRiB1 messages … … 308 324 dyconst=180./(dy*r_earth*pi) 309 325 gotGrid=1 310 ! Check whether fields are global311 ! If they contain the poles, specify polar stereographic map312 ! projections using the stlmbr- and stcm2p-calls313 !***********************************************************326 ! Check whether fields are global 327 ! If they contain the poles, specify polar stereographic map 328 ! projections using the stlmbr- and stcm2p-calls 329 !*********************************************************** 314 330 315 331 xauxa=abs(xaux2+dx-360.-xaux1) … … 332 348 if (xglobal.and.xauxa.lt.0.001) then 333 349 sglobal=.true. ! field contains south pole 334 ! Enhance the map scale by factor 3 (*2=6) compared to north-south335 ! map scale350 ! Enhance the map scale by factor 3 (*2=6) compared to north-south 351 ! map scale 336 352 sizesouth=6.*(switchsouth+90.)/dy 337 353 call stlmbr(southpolemap,-90.,0.) … … 346 362 if (xglobal.and.xauxa.lt.0.001) then 347 363 nglobal=.true. ! field contains north pole 348 ! Enhance the map scale by factor 3 (*2=6) compared to north-south349 ! map scale364 ! Enhance the map scale by factor 3 (*2=6) compared to north-south 365 ! map scale 350 366 sizenorth=6.*(90.-switchnorth)/dy 351 367 call stlmbr(northpolemap,90.,0.) … … 363 379 364 380 if (nx.gt.nxmax) then 365 write(*,*) 'FLEXPART error: Too many grid points in x direction.'381 write(*,*) 'FLEXPART error: Too many grid points in x direction.' 366 382 write(*,*) 'Reduce resolution of wind fields.' 367 383 write(*,*) 'Or change parameter settings in file par_mod.' … … 371 387 372 388 if (ny.gt.nymax) then 373 write(*,*) 'FLEXPART error: Too many grid points in y direction.'389 write(*,*) 'FLEXPART error: Too many grid points in y direction.' 374 390 write(*,*) 'Reduce resolution of wind fields.' 375 391 write(*,*) 'Or change parameter settings in file par_mod.' … … 410 426 ! 411 427 412 30 428 30 call grib_close_file(ifile) 413 429 414 430 !error message if no fields found with correct first longitude in it … … 469 485 470 486 numskip=nlev_ec-nuvz ! number of ecmwf model layers not used 471 487 ! by trajectory model 472 488 !do 8940 i=1,244 473 489 ! write (*,*) 'zsec2:',i,ifield,zsec2(i),numskip … … 482 498 k=nlev_ec+1+numskip+i 483 499 akm(nwz-i+1)=zsec2(j) 484 ! write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)500 ! write (*,*) 'ifield:',ifield,k,j,zsec2(10+j) 485 501 bkm(nwz-i+1)=zsec2(k) 486 502 end do … … 498 514 bkz(1)=1.0 499 515 do i=1,nuvz 500 501 516 akz(i+1)=0.5*(akm(i+1)+akm(i)) 517 bkz(i+1)=0.5*(bkm(i+1)+bkm(i)) 502 518 end do 503 519 nuvz=nuvz+1 … … 539 555 if (pint.lt.5000.) goto 96 540 556 end do 541 96 557 96 nconvlev=i 542 558 if (nconvlev.gt.nconvlevmax-1) then 543 559 nconvlev=nconvlevmax-1 … … 548 564 return 549 565 550 999 566 999 write(*,*) 551 567 write(*,*) ' ###########################################'// & 552 568 '###### ' … … 568 584 endif 569 585 570 end subroutine gridcheck 571 586 end subroutine gridcheck_ecmwf 587 -
src/gridcheck_gfs.f90
r4c64400 rd8eed02 20 20 !********************************************************************** 21 21 22 subroutine gridcheck 22 subroutine gridcheck_gfs 23 23 24 24 !********************************************************************** … … 38 38 ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * 39 39 ! ECMWF grib_api * 40 ! * 41 ! Unified ECMWF and GFS builds * 42 ! Marian Harustak, 12.5.2017 * 43 ! - Renamed routine from gridcheck to gridcheck_gfs * 40 44 ! * 41 45 !********************************************************************** … … 227 231 yaux2in,iret) 228 232 call grib_check(iret,gribFunction,gribErrorMsg) 233 234 ! Fix for flexpart.eu ticket #48 235 if (xaux2in.lt.0) xaux2in = 359.0 229 236 230 237 xaux1=xaux1in … … 537 544 endif 538 545 539 end subroutine gridcheck 546 end subroutine gridcheck_gfs -
src/makefile
r6985a98 rd8eed02 12 12 # 'make -j ecmwf gcc=4.9', 13 13 # also set environment variable LD_LIBRARY_PATH to point to compiler libraries 14 # 15 # Makefile was modified to produce unified executable for both ECMWF and GFS meteo data formats 16 # gributils were included to detect format of meteo data 17 # 18 # Cpp directives USE_MPIINPLACE were added to three source files. The effect of these directives 19 # are to enable the MPI_IN_PLACE option only if compiled with a -DUSE_MPIINPLACE directive. 20 # Otherwise, a safer option (which requires the allocation of another array) is used by default. 21 # In makefile added the -x f95-cpp-input flag for compiling of cpp directives. 14 22 # 15 23 # USAGE 16 # Compile serial FLEXPART (ECMWF)17 # make [-j] ecmwf18 # 19 # Compile parallel FLEXPART (ECMWF)20 # make [-j] ecmwf-mpi24 # Compile serial FLEXPART 25 # make [-j] serial 26 # 27 # Compile parallel FLEXPART 28 # make [-j] mpi 21 29 # 22 # Compile for debugging parallel FLEXPART (ECMWF) 23 # make [-j] ecmwf-mpi-dbg 24 # 25 # Compile serial FLEXPART (GFS) 26 # make [-j] gfs 27 # 28 # Compile parallel FLEXPART (GFS) 29 # make [-j] gfs-mpi 30 # Compile for debugging parallel FLEXPART 31 # make [-j] mpi-dbg 30 32 # 31 33 ################################################################################ 32 34 33 35 ## PROGRAMS 34 FLEXPART-ECMWF-MPI = FP_ecmwf_MPI 35 FLEXPART-ECMWF-MPI-DBG = DBG_FP_ecmwf_MPI 36 FLEXPART-ECMWF = FP_ecmwf_gfortran 37 FLEXPART-GFS = FP_gfs_gfortran 38 FLEXPART-GFS-MPI = FP_gfs_MPI 36 # Unified executable names 37 # The same executable is used for both ECMWF and GFS metdata 38 39 # Parallel processing executable 40 FLEXPART-MPI = FLEXPART_MPI 41 42 # Parallel processing executable with debugging info 43 FLEXPART-MPI-DBG = DBG_FLEXPART_MPI 44 45 # Serial processing executable 46 FLEXPART-SERIAL = FLEXPART 47 39 48 40 49 ifeq ($(gcc), 4.9) … … 61 70 62 71 72 # path to gributils used to detect meteodata format 73 VPATH = gributils/ 74 75 63 76 ## OPTIMIZATION LEVEL 64 77 O_LEV = 2 # [0,1,2,3,g,s,fast] … … 68 81 LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp 69 82 70 FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g - m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) #-Warray-bounds -fcheck=all # -march=native71 72 DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 - m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all73 74 LDFLAGS = $(FFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)83 FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) #-Warray-bounds -fcheck=all # -march=native 84 85 DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all 86 87 LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2) 75 88 LDDEBUG = $(DBGFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2) 76 89 … … 82 95 xmass_mod.o flux_mod.o \ 83 96 point_mod.o outg_mod.o \ 84 mean_mod.o random_mod.o 97 mean_mod.o random_mod.o \ 98 class_gribfile_mod.o 85 99 86 100 MPI_MODOBJS = \ … … 99 113 redist.o \ 100 114 concoutput_surf.o concoutput_surf_nest.o \ 101 getfields.o 115 getfields.o \ 116 readwind_ecmwf.o 102 117 103 118 ## For MPI version … … 112 127 redist_mpi.o \ 113 128 concoutput_surf_mpi.o concoutput_surf_nest_mpi.o \ 114 getfields_mpi.o 115 116 ### WINDFIELDS 117 ## For ECMWF (serial) version: 118 OBJECTS_ECMWF = \ 119 calcpar.o readwind.o \ 120 richardson.o verttransform.o \ 121 obukhov.o gridcheck.o \ 122 convmix.o calcmatrix.o \ 123 ecmwf_mod.o 124 125 126 ## For ECMWF MPI version: 127 OBJECTS_ECMWF_MPI = \ 128 gridcheck.o readwind_mpi.o \ 129 calcpar.o \ 130 richardson.o verttransform.o \ 131 obukhov.o \ 132 convmix.o calcmatrix.o \ 133 ecmwf_mod.o 134 135 ## For GFS (serial) version: 136 OBJECTS_GFS = \ 137 calcpar_gfs.o readwind_gfs.o \ 138 richardson_gfs.o verttransform_gfs.o \ 139 obukhov_gfs.o gridcheck_gfs.o \ 140 convmix_gfs.o calcmatrix_gfs.o \ 141 gfs_mod.o 129 getfields_mpi.o \ 130 readwind_ecmwf_mpi.o 142 131 143 132 OBJECTS = \ … … 154 143 ew.o readreleases.o \ 155 144 readdepo.o get_vdep_prob.o \ 156 get_wetscav.o \145 get_wetscav.o readwind_gfs.o \ 157 146 psim.o outgrid_init.o \ 158 outgrid_init_nest.o \147 outgrid_init_nest.o calcmatrix.o \ 159 148 photo_O1D.o readlanduse.o \ 160 149 interpol_wind.o readoutgrid.o \ 161 150 interpol_all.o readpaths.o \ 162 getrb.o \163 getrc.o \151 getrb.o obukhov.o \ 152 getrc.o convmix.o \ 164 153 getvdep.o readspecies.o \ 165 interpol_misslev.o \166 scalev.o \167 pbl_profile.o readOHfield.o \168 juldate.o \154 interpol_misslev.o richardson.o \ 155 scalev.o verttransform_ecmwf.o \ 156 pbl_profile.o readOHfield.o \ 157 juldate.o verttransform_gfs.o \ 169 158 interpol_vdep.o interpol_rain.o \ 170 159 hanna.o wetdepokernel.o \ 171 160 calcpar.o wetdepo.o \ 172 161 hanna_short.o windalign.o \ 173 hanna1.o 174 162 hanna1.o gridcheck_ecmwf.o \ 163 gridcheck_gfs.o gridcheck_nests.o \ 175 164 readwind_nests.o calcpar_nests.o \ 176 165 verttransform_nests.o interpol_all_nests.o \ 177 166 interpol_wind_nests.o interpol_misslev_nests.o \ 178 167 interpol_vdep_nests.o interpol_rain_nests.o \ 179 readageclasses.o \168 readageclasses.o detectformat.o \ 180 169 calcfluxes.o fluxoutput.o \ 181 170 qvsat.o skplin.o \ … … 201 190 %.o: %.mod 202 191 203 ecmwf: $(FLEXPART-ECMWF) 204 ecmwf: FC := $(F90) 205 206 ecmwf-mpi: $(FLEXPART-ECMWF-MPI) 207 ecmwf-mpi: FC := $(MPIF90) 208 209 ecmwf-mpi-dbg: $(FLEXPART-ECMWF-MPI-DBG) 210 ecmwf-mpi-dbg: FFLAGS := $(DBGFLAGS) 211 ecmwf-mpi-dbg: LDFLAGS:= $(LDDEBUG) 212 ecmwf-mpi-dbg: FC := $(MPIF90) 213 214 gfs: $(FLEXPART-GFS) 215 gfs: FC := $(F90) 216 217 gfs-mpi: $(FLEXPART-GFS-MPI) 218 gfs-mpi: FC := $(MPIF90) 219 220 #all: $(FLEXPART-ECMWF) 221 #all: $(FLEXPART-ECMWF-MPI) 222 223 ## This allows for switching between ECMWF/GFS without editing source code. 224 wind_mod = ecmwf_mod.o # default wind 225 # ifeq ($(MAKECMDGOALS),ecmwf) 226 # wind_mod = ecmwf_mod.o 227 # endif 228 # ifeq ($(MAKECMDGOALS),ecmwf-mpi) 229 # wind_mod = ecmwf_mod.o 230 # endif 231 # ifeq ($(MAKECMDGOALS),ecmwf-mpi-dbg) 232 # wind_mod = ecmwf_mod.o 233 # endif 234 235 ifeq ($(MAKECMDGOALS),gfs) 236 wind_mod = gfs_mod.o 237 endif 238 ifeq ($(MAKECMDGOALS),gfs-mpi) 239 wind_mod = gfs_mod.o 240 endif 241 242 $(FLEXPART-ECMWF): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_ECMWF) 243 +$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_ECMWF) $(LDFLAGS) 244 245 $(FLEXPART-ECMWF-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) $(OBJECTS_ECMWF_MPI) 192 # serial executable 193 serial: $(FLEXPART-SERIAL) 194 serial: FC := $(F90) 195 196 # parallel processing executable 197 mpi: $(FLEXPART-MPI) 198 mpi: FC := $(MPIF90) 199 200 # parallel processing with debugging info 201 mpi-dbg: $(FLEXPART-MPI-DBG) 202 mpi-dbg: FFLAGS := $(DBGFLAGS) 203 mpi-dbg: LDFLAGS:= $(LDDEBUG) 204 mpi-dbg: FC := $(MPIF90) 205 206 $(FLEXPART-SERIAL): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) 207 +$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(LDFLAGS) 208 209 $(FLEXPART-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) 246 210 +$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \ 247 $(OBJECTS_ECMWF_MPI) $(LDFLAGS) 248 # +$(FC) -o $@ *.o $(LDFLAGS) 249 250 $(FLEXPART-ECMWF-MPI-DBG): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \ 251 $(OBJECTS_ECMWF_MPI) 211 $(LDFLAGS) 212 213 $(FLEXPART-MPI-DBG): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) 252 214 +$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \ 253 $(OBJECTS_ECMWF_MPI) $(LDFLAGS) 254 255 $(FLEXPART-GFS): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_GFS) 256 +$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_GFS) $(LDFLAGS) 257 258 $(FLEXPART-GFS-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) $(OBJECTS_GFS) 259 +$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \ 260 $(OBJECTS_GFS) $(LDFLAGS) 215 $(LDFLAGS) 261 216 262 217 %.o: %.f90 … … 267 222 268 223 cleanall: 269 \rm -f *.o *.mod $(FLEXPART- ECMWF-MPI) $(FLEXPART-ECMWF-MPI-DBG) $(FLEXPART-ECMWF) \270 $(FLEXPART-GFS-MPI) $(FLEXPART-GFS) 224 \rm -f *.o *.mod $(FLEXPART-MPI) $(FLEXPART-MPI-DBG) $(FLEXPART-SERIAL) 225 271 226 272 227 .SUFFIXES = $(SUFFIXES) .f90 … … 282 237 random_mod.o 283 238 calcfluxes.o: com_mod.o flux_mod.o outg_mod.o par_mod.o 284 calcmatrix.o: com_mod.o conv_mod.o par_mod.o 285 calcmatrix_gfs.o: com_mod.o conv_mod.o par_mod.o 286 calcpar.o: com_mod.o par_mod.o 287 calcpar_gfs.o: com_mod.o par_mod.o 239 calcmatrix.o: com_mod.o conv_mod.o par_mod.o class_gribfile_mod.o 240 calcpar.o: com_mod.o par_mod.o class_gribfile_mod.o 288 241 calcpar_nests.o: com_mod.o par_mod.o 289 242 calcpv.o: com_mod.o par_mod.o … … 311 264 conv_mod.o: par_mod.o 312 265 convect43c.o: conv_mod.o par_mod.o 313 convmix.o: com_mod.o conv_mod.o flux_mod.o par_mod.o 314 convmix_gfs.o: com_mod.o conv_mod.o par_mod.o 266 convmix.o: com_mod.o conv_mod.o flux_mod.o par_mod.o class_gribfile_mod.o 315 267 coordtrafo.o: com_mod.o par_mod.o point_mod.o 268 detectformat.o: com_mod.o par_mod.o class_gribfile_mod.o 316 269 distance.o: par_mod.o 317 270 distance2.o: par_mod.o … … 319 272 drydepokernel_nest.o: com_mod.o par_mod.o unc_mod.o 320 273 erf.o: par_mod.o 321 FLEXPART.o: com_mod.o conv_mod.o par_mod.o point_mod.o random_mod.o netcdf_output_mod.o 274 FLEXPART.o: com_mod.o conv_mod.o par_mod.o point_mod.o random_mod.o netcdf_output_mod.o class_gribfile_mod.o 322 275 FLEXPART_MPI.o: com_mod.o conv_mod.o mpi_mod.o par_mod.o point_mod.o \ 323 random_mod.o netcdf_output_mod.o 276 random_mod.o netcdf_output_mod.o class_gribfile_mod.o 324 277 fluxoutput.o: com_mod.o flux_mod.o outg_mod.o par_mod.o 325 278 get_settling.o: com_mod.o par_mod.o 326 getfields.o: com_mod.o par_mod.o 327 getfields_mpi.o: com_mod.o par_mod.o mpi_mod.o 279 getfields.o: com_mod.o par_mod.o class_gribfile_mod.o 280 getfields_mpi.o: com_mod.o par_mod.o mpi_mod.o class_gribfile_mod.o 328 281 gethourlyOH.o: com_mod.o oh_mod.o par_mod.o 329 282 getrb.o: par_mod.o … … 331 284 getvdep.o: com_mod.o par_mod.o 332 285 getvdep_nests.o: com_mod.o par_mod.o 333 gridcheck .o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o286 gridcheck_ecmwf.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o 334 287 gridcheck_emos.o: com_mod.o conv_mod.o par_mod.o 335 288 gridcheck_fnl.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o … … 366 319 mpi_mod.o: com_mod.o par_mod.o unc_mod.o 367 320 netcdf_output_mod.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o 368 obukhov.o: par_mod.o 369 obukhov_gfs.o: par_mod.o 321 obukhov.o: par_mod.o class_gribfile_mod.o 370 322 ohreaction.o: com_mod.o oh_mod.o par_mod.o 371 323 openouttraj.o: com_mod.o par_mod.o point_mod.o … … 374 326 outgrid_init.o: com_mod.o flux_mod.o oh_mod.o outg_mod.o par_mod.o unc_mod.o 375 327 outgrid_init_nest.o: com_mod.o outg_mod.o par_mod.o unc_mod.o 376 par_mod.o : $(wind_mod)377 328 part0.o: par_mod.o 378 329 partdep.o: par_mod.o … … 402 353 readreleases.o: com_mod.o par_mod.o point_mod.o xmass_mod.o 403 354 readspecies.o: com_mod.o par_mod.o 404 readwind .o: com_mod.o par_mod.o355 readwind_ecmwf.o: com_mod.o par_mod.o 405 356 readwind_emos.o: com_mod.o par_mod.o 406 357 readwind_gfs.o: com_mod.o par_mod.o 407 358 readwind_gfs_emos.o: com_mod.o par_mod.o 408 readwind_ mpi.o: com_mod.o mpi_mod.o par_mod.o359 readwind_ecmwf_mpi.o: com_mod.o mpi_mod.o par_mod.o 409 360 readwind_nests.o: com_mod.o par_mod.o 410 361 readwind_nests_emos.o: com_mod.o par_mod.o … … 415 366 releaseparticles_mpi.o: com_mod.o mpi_mod.o par_mod.o point_mod.o \ 416 367 random_mod.o xmass_mod.o 417 richardson.o: par_mod.o 418 richardson_gfs.o: par_mod.o 368 richardson.o: par_mod.o class_gribfile_mod.o 419 369 scalev.o: par_mod.o 420 370 shift_field.o: par_mod.o … … 425 375 par_mod.o point_mod.o unc_mod.o xmass_mod.o netcdf_output_mod.o 426 376 unc_mod.o: par_mod.o 427 verttransform .o: cmapf_mod.o com_mod.o par_mod.o377 verttransform_ecmwf.o: cmapf_mod.o com_mod.o par_mod.o 428 378 verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o 429 379 verttransform_nests.o: com_mod.o par_mod.o -
src/mpi_mod.f90
raa8c34a rbb579a9 59 59 ! * 60 60 !***************************************************************************** 61 ! * 62 ! Modification by DJM, 2017-05-09 - added #ifdef USE_MPIINPLACE cpp * 63 ! directive to mpif_tm_reduce_grid() to insure that MPI_IN_PLACE is * 64 ! used in the MPI_Reduce() only if specifically compiled with that * 65 ! directive. * 66 ! * 67 !***************************************************************************** 61 68 62 69 use mpi … … 2425 2432 2426 2433 ! 2) Using in-place reduction 2434 2435 !!!-------------------------------------------------------------------- 2436 !!! DJM - 2017-05-09 change - MPI_IN_PLACE option for MPI_Reduce() causes 2437 !!! severe numerical problems in some cases. I'm guessing there are memory 2438 !!! overrun issues in this complex code, but have so far failed to identify 2439 !!! a specific problem. And/or, when one searches the Internet for this 2440 !!! issue, there is "some" hint that the implementation may be buggy. 2441 !!! 2442 !!! At this point, with the following CPP USE_MPIINPLACE directive, the 2443 !!! default behaviour will be to not use the MPI_IN_PLACE option. 2444 !!! Users will have to compile with -DUSE_MPIINPLACE if they want that option. 2445 !!! Introduction of the CPP directives also requires that the code be compiled 2446 !!! with the "-x f95-cpp-input" option. 2447 !!! 2448 !!! Modification of this section requires the addition of array gridunc0, which 2449 !!! requires an allocation in outgrid_init.f90 and initial declaration in 2450 !!! unc_mod.f90. 2451 !!!-------------------------------------------------------------------- 2452 2453 #ifdef USE_MPIINPLACE 2454 2427 2455 if (lroot) then 2428 2456 call MPI_Reduce(MPI_IN_PLACE, gridunc, grid_size3d, mp_sp, MPI_SUM, id_root, & … … 2433 2461 & mp_comm_used, mp_ierr) 2434 2462 end if 2463 2464 #else 2465 2466 call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, & 2467 & mp_comm_used, mp_ierr) 2468 if (lroot) gridunc = gridunc0 2469 2470 #endif 2435 2471 2436 2472 if ((WETDEP).and.(ldirect.gt.0)) then -
src/obukhov.f90
re200b7a r6ecb30a 20 20 !********************************************************************** 21 21 22 real function obukhov(ps,tsurf,tdsurf,tlev,ustar,hf,akm,bkm )22 real function obukhov(ps,tsurf,tdsurf,tlev,ustar,hf,akm,bkm,plev,metdata_format) 23 23 24 24 !******************************************************************** … … 27 27 ! Date: 1994-06-27 * 28 28 ! * 29 ! Update: A. Stohl, 2000-09-25, avoid division by zero by*30 ! setting ustar to minimum value*29 ! This program calculates Obukhov scale height from surface * 30 ! meteorological data and sensible heat flux. * 31 31 ! * 32 32 !******************************************************************** 33 33 ! * 34 ! This program calculates Obukhov scale height from surface * 35 ! meteorological data and sensible heat flux. * 34 ! Update: A. Stohl, 2000-09-25, avoid division by zero by * 35 ! setting ustar to minimum value * 36 ! CHANGE: 17/11/2005 Caroline Forster NCEP GFS version * 37 ! * 38 ! Unified ECMWF and GFS builds * 39 ! Marian Harustak, 12.5.2017 * 40 ! - Merged obukhov and obukhov_gfs into one routine using * 41 ! if-then for meteo-type dependent code * 36 42 ! * 37 43 !******************************************************************** … … 47 53 ! akm ECMWF vertical discretization parameter * 48 54 ! bkm ECMWF vertical discretization parameter * 55 ! plev * 56 ! metdata_format format of metdata (ecmwf/gfs) * 49 57 ! * 50 58 !******************************************************************** 51 59 52 60 use par_mod 61 use class_gribfile 53 62 54 63 implicit none 55 64 65 integer :: metdata_format 56 66 real :: akm(nwzmax),bkm(nwzmax) 57 67 real :: ps,tsurf,tdsurf,tlev,ustar,hf,e,ew,tv,rhoa,plev … … 62 72 tv=tsurf*(1.+0.378*e/ps) ! virtual temperature 63 73 rhoa=ps/(r_air*tv) ! air density 74 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 64 75 ak1=(akm(1)+akm(2))/2. 65 76 bk1=(bkm(1)+bkm(2))/2. 66 77 plev=ak1+bk1*ps ! Pressure level 1 78 end if 67 79 theta=tlev*(100000./plev)**(r_air/cpa) ! potential temperature 68 80 if (ustar.le.0.) ustar=1.e-8 -
src/outgrid_init.f90
r6b22af9 r6ecb30a 19 19 ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * 20 20 !********************************************************************** 21 ! * 22 ! DJM - 2017-05-09 - added #ifdef USE_MPIINPLACE cpp directive to * 23 ! enable allocation of a gridunc0 array if required by MPI code in * 24 ! mpi_mod.f90 * 25 ! * 26 !********************************************************************** 27 21 28 22 29 subroutine outgrid_init … … 215 222 ! Extra field for totals at MPI root process 216 223 if (lroot.and.mpi_mode.gt.0) then 217 ! allocate(gridunc0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, & 218 ! maxpointspec_act,nclassunc,maxageclass),stat=stat) 219 ! if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc0' 224 225 #ifdef USE_MPIINPLACE 226 #else 227 ! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid(), 228 ! then an aux array is needed for parallel grid reduction 229 allocate(gridunc0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, & 230 maxpointspec_act,nclassunc,maxageclass),stat=stat) 231 if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc0' 232 #endif 220 233 if (ldirect.gt.0) then 221 234 allocate(wetgridunc0(0:numxgrid-1,0:numygrid-1,maxspec, & -
src/par_mod.f90
r6985a98 rd8eed02 37 37 38 38 module par_mod 39 40 !************************************************************************41 ! wind_mod: is gfs_mod.f90 for target gfs, ecmwf_mod.f90 for target ecmwf42 !************************************************************************43 use wind_mod44 39 45 40 implicit none … … 146 141 !********************************************* 147 142 148 ! Moved to ecmwf_mod.f90 (for ECMWF) / gfs_mod.f90 (for GFS) 143 ! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new 144 integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 145 146 ! integer,parameter :: nxmax=181,nymax=91,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 147 148 ! INTEGER,PARAMETER :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !NCEP data 149 150 ! integer,parameter :: nxshift=359 ! for ECMWF 151 integer,parameter :: nxshift=0 ! for GFS 152 153 !********************************************* 154 ! Maximum dimensions of the nested input grids 155 !********************************************* 156 157 integer,parameter :: maxnests=1,nxmaxn=451,nymaxn=226 158 159 ! nxmax,nymax maximum dimension of wind fields in x and y 160 ! direction, respectively 161 ! nuvzmax,nwzmax maximum dimension of (u,v) and (w) wind fields in z 162 ! direction (for fields on eta levels) 163 ! nzmax maximum dimension of wind fields in z direction 164 ! for the transformed Cartesian coordinates 165 ! nxshift for global grids (in x), the grid can be shifted by 166 ! nxshift grid points, in order to accomodate nested 167 ! grids, and output grids overlapping the domain "boundary" 168 ! nxshift must not be negative; "normal" setting would be 0 169 149 170 150 171 integer,parameter :: nconvlevmax = nuvzmax-1 … … 269 290 integer,parameter :: icmv=-9999 270 291 271 ! Parameters for testing272 !*******************************************273 ! integer :: verbosity=0274 292 275 293 end module par_mod -
src/readwind_gfs.f90
r8a65cb0 r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine readwind (indj,n,uuh,vvh,wwh)22 subroutine readwind_gfs(indj,n,uuh,vvh,wwh) 23 23 24 24 !*********************************************************************** … … 38 38 !* CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * 39 39 !* ECMWF grib_api * 40 ! * 41 ! Unified ECMWF and GFS builds * 42 ! Marian Harustak, 12.5.2017 * 43 ! - Renamed routine from readwind to readwind_gfs * 40 44 !* * 41 45 !*********************************************************************** … … 716 720 stop 'Execution terminated' 717 721 718 end subroutine readwind 722 end subroutine readwind_gfs -
src/richardson.f90
rc2162ce r6ecb30a 21 21 22 22 subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, & 23 akz,bkz,hf,tt2,td2,h,wst,hmixplus )23 akz,bkz,hf,tt2,td2,h,wst,hmixplus,metdata_format) 24 24 ! i i i i i i i 25 25 ! i i i i i o o o … … 41 41 ! Meteor. 81, 245-269. * 42 42 ! * 43 !**************************************************************************** 44 ! * 43 45 ! Update: 1999-02-01 by G. Wotawa * 44 46 ! * … … 46 48 ! instead of first model level. * 47 49 ! New input variables tt2, td2 introduced. * 50 ! * 51 ! CHANGE: 17/11/2005 Caroline Forster NCEP GFS version * 52 ! * 53 ! Unified ECMWF and GFS builds * 54 ! Marian Harustak, 12.5.2017 * 55 ! - Merged richardson and richardson_gfs into one routine using * 56 ! if-then for meteo-type dependent code * 48 57 ! * 49 58 !**************************************************************************** … … 55 64 ! tv virtual temperature * 56 65 ! wst convective velocity scale * 66 ! metdata_format format of metdata (ecmwf/gfs) * 57 67 ! * 58 68 ! Constants: * … … 62 72 63 73 use par_mod 74 use class_gribfile 64 75 65 76 implicit none 66 77 67 integer :: i,k,nuvz,iter 78 integer :: metdata_format 79 integer :: i,k,nuvz,iter,llev,loop_start 68 80 real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri 69 81 real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew … … 77 89 iter=0 78 90 91 if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then 92 ! NCEP version: find first model level above ground 93 !************************************************** 94 95 llev = 0 96 do i=1,nuvz 97 if (psurf.lt.akz(i)) llev=i 98 end do 99 llev = llev+1 100 ! sec llev should not be 1! 101 if (llev.eq.1) llev = 2 102 if (llev.gt.nuvz) llev = nuvz-1 103 ! NCEP version 104 end if 105 106 79 107 ! Compute virtual temperature and virtual potential temperature at 80 108 ! reference level (2 m) … … 95 123 ! Integrate z up to one level above zt 96 124 !************************************* 97 98 do k=2,nuvz 125 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 126 loop_start=2 127 else 128 loop_start=llev 129 end if 130 do k=loop_start,nuvz 99 131 pint=akz(k)+bkz(k)*psurf ! pressure on model layers 100 132 tv=ttlev(k)*(1.+0.608*qvlev(k)) -
src/timemanager.f90
ra76d954 r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine timemanager 22 subroutine timemanager(metdata_format) 23 23 24 24 !***************************************************************************** … … 50 50 ! For compatibility with MPI version, * 51 51 ! variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod * 52 ! Unified ECMWF and GFS builds * 53 ! Marian Harustak, 12.5.2017 * 54 ! - Added passing of metdata_format as it was needed by called routines * 52 55 !***************************************************************************** 53 56 ! * … … 83 86 ! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) = * 84 87 ! spatial positions of trajectories * 88 ! metdata_format format of metdata (ecmwf/gfs) * 85 89 ! * 86 90 ! Constants: * … … 102 106 implicit none 103 107 108 integer :: metdata_format 104 109 integer :: j,ks,kp,l,n,itime=0,nstop,nstop1 105 110 ! integer :: ksp … … 194 199 write (*,*) 'timemanager> call convmix -- backward' 195 200 endif 196 call convmix(itime )201 call convmix(itime,metdata_format) 197 202 if (verbosity.gt.1) then 198 203 !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max) … … 207 212 write (*,*) 'timemanager> call getfields' 208 213 endif 209 call getfields(itime,nstop1 )214 call getfields(itime,nstop1,metdata_format) 210 215 if (verbosity.gt.1) then 211 216 CALL SYSTEM_CLOCK(count_clock) … … 266 271 write (*,*) 'timemanager> call convmix -- forward' 267 272 endif 268 call convmix(itime )273 call convmix(itime,metdata_format) 269 274 endif 270 275 -
src/timemanager_mpi.f90
rb5127f9 rd8eed02 20 20 !********************************************************************** 21 21 22 subroutine timemanager 22 subroutine timemanager(metdata_format) 23 23 24 24 !***************************************************************************** … … 50 50 ! MPI version * 51 51 ! Variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod * 52 ! Unified ECMWF and GFS builds * 53 ! Marian Harustak, 12.5.2017 * 54 ! - Added passing of metdata_format as it was needed by called routines * 52 55 !***************************************************************************** 53 56 ! * … … 83 86 ! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) = * 84 87 ! spatial positions of trajectories * 88 ! metdata_format format of metdata (ecmwf/gfs) * 85 89 ! * 86 90 ! Constants: * … … 103 107 implicit none 104 108 109 integer :: metdata_format 105 110 logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete 106 111 integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0 … … 205 210 endif 206 211 ! readwind process skips this step 207 if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime )212 if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime,metdata_format) 208 213 endif 209 214 … … 218 223 if (mp_measure_time) call mpif_mtime('getfields',0) 219 224 220 call getfields(itime,nstop1,memstat )225 call getfields(itime,nstop1,memstat,metdata_format) 221 226 222 227 if (mp_measure_time) call mpif_mtime('getfields',1) … … 346 351 write (*,*) 'timemanager> call convmix -- forward' 347 352 endif 348 call convmix(itime )353 call convmix(itime,metdata_format) 349 354 endif 350 355 -
src/unc_mod.f90
r4c64400 r6ecb30a 19 19 ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * 20 20 !********************************************************************** 21 ! * 22 ! DJM - 2017-05-09 - added #ifdef USE_MPIINPLACE cpp directive to * 23 ! enable declaration of a gridunc0 array if required by MPI code in * 24 ! mpi_mod.f90 * 25 ! * 26 !********************************************************************** 21 27 22 28 module unc_mod … … 27 33 28 34 real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc 35 #ifdef USE_MPIINPLACE 36 #else 37 ! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid(), 38 ! then an aux array is needed for parallel grid reduction 39 real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc0 40 #endif 29 41 real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn 30 42 real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc -
src/verttransform_gfs.f90
r4fbe7a5 r6ecb30a 20 20 !********************************************************************** 21 21 22 subroutine verttransform (n,uuh,vvh,wwh,pvh)22 subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh) 23 23 ! i i i i i 24 24 !***************************************************************************** … … 47 47 ! Changes, Bernd C. Krueger, Feb. 2001: 48 48 ! Variables tth and qvh (on eta coordinates) from common block 49 ! 50 ! Unified ECMWF and GFS builds 51 ! Marian Harustak, 12.5.2017 52 ! - Renamed routine from verttransform to verttransform_gfs 53 ! 49 54 !***************************************************************************** 50 55 ! * … … 531 536 532 537 533 end subroutine verttransform 538 end subroutine verttransform_gfs -
src/writeheader.f90
rf13406c r6ecb30a 35 35 !***************************************************************************** 36 36 ! * 37 ! Modified to remove TRIM around the output of flexversion so that * 38 ! it will be a constant length (defined in com_mod.f90) in output header * 39 ! * 40 ! Don Morton, Boreal Scientific Computing * 41 ! 07 May 2017 * 42 ! * 43 !***************************************************************************** 44 ! * 37 45 ! Variables: * 38 46 ! * … … 67 75 68 76 if (ldirect.eq.1) then 69 write(unitheader) ibdate,ibtime, trim(flexversion)77 write(unitheader) ibdate,ibtime, flexversion 70 78 else 71 write(unitheader) iedate,ietime, trim(flexversion)79 write(unitheader) iedate,ietime, flexversion 72 80 endif 73 81 -
src/writeheader_nest.f90
r4fbe7a5 r6ecb30a 35 35 !***************************************************************************** 36 36 ! * 37 ! Modified to remove TRIM around the output of flexversion so that * 38 ! it will be a constant length (defined in com_mod.f90) in output header * 39 ! * 40 ! Don Morton, Boreal Scientific Computing * 41 ! 07 May 2017 * 42 ! * 43 !***************************************************************************** 44 ! * 37 45 ! Variables: * 38 46 ! * … … 67 75 68 76 if (ldirect.eq.1) then 69 write(unitheader) ibdate,ibtime, trim(flexversion)77 write(unitheader) ibdate,ibtime, flexversion 70 78 else 71 write(unitheader) iedate,ietime, trim(flexversion)79 write(unitheader) iedate,ietime, flexversion 72 80 endif 73 81 -
src/writeheader_nest_surf.f90
rf13406c r6ecb30a 35 35 !***************************************************************************** 36 36 ! * 37 ! Modified to remove TRIM around the output of flexversion so that * 38 ! it will be a constant length (defined in com_mod.f90) in output header * 39 ! * 40 ! Don Morton, Boreal Scientific Computing * 41 ! 07 May 2017 * 42 ! * 43 !***************************************************************************** 44 ! * 37 45 ! Variables: * 38 46 ! * … … 67 75 68 76 if (ldirect.eq.1) then 69 write(unitheader) ibdate,ibtime, trim(flexversion)77 write(unitheader) ibdate,ibtime,flexversion 70 78 else 71 write(unitheader) iedate,ietime, trim(flexversion)79 write(unitheader) iedate,ietime,flexversion 72 80 endif 73 81 -
src/writeheader_surf.f90
rf13406c r6ecb30a 35 35 !***************************************************************************** 36 36 ! * 37 ! Modified to remove TRIM around the output of flexversion so that * 38 ! it will be a constant length (defined in com_mod.f90) in output header * 39 ! * 40 ! Don Morton, Boreal Scientific Computing * 41 ! 07 May 2017 * 42 ! * 43 !***************************************************************************** 44 ! * 37 45 ! Variables: * 38 46 ! * … … 67 75 68 76 if (ldirect.eq.1) then 69 write(unitheader) ibdate,ibtime, trim(flexversion)77 write(unitheader) ibdate,ibtime, flexversion 70 78 else 71 write(unitheader) iedate,ietime, trim(flexversion)79 write(unitheader) iedate,ietime, flexversion 72 80 endif 73 81 -
tests/NILU/basic_tests
rca350ba rdd6a4aa 1 1 #WORKSPACE=/xnilu_wrk/flexbuild 2 2 #WORKSPACE=/home/ignacio/repos/ 3 FP_exec=$WORKSPACE/src/F P_ecmwf_gfortran3 FP_exec=$WORKSPACE/src/FLEXPART 4 4 path_flextest=$WORKSPACE/flextest/ 5 5 declare -a test_names=('1' 'HelloWorld' 'Fwd1' 'Fwd2' 'Bwd1' 'Volc' '2') -
tests/NILU/run_tests
rca350ba rdd6a4aa 9 9 #FP_exec=/home/ignacio/repos/flexpart/src/FP_ecmwf_gfortran 10 10 #FP_exec=/xnilu_wrk/flexbuild/tests/NILU/FP_ecmwf_gfortran 11 FP_exec=$WORKSPACE/src/FP_ecmwf_gfortran 11 #FP_exec=$WORKSPACE/src/FP_ecmwf_gfortran 12 # CTBTO unified executable 13 FP_exec=$WORKSPACE/src/FLEXPART 12 14 #path_flextest=/home/ignacio/repos/flextest/ 13 15 #path_flextest=/xnilu_wrk/flexbuild/flextest/ -
src/coordtrafo.f90
re200b7a r1d072c0 46 46 47 47 integer :: i,j,k 48 real :: yrspc ! small real number relative to x 48 49 49 50 if (numpoint.eq.0) goto 30 … … 64 65 ! CHECK IF RELEASE POINTS ARE WITHIN DOMAIN 65 66 !****************************************** 66 67 68 yrspc = spacing(real(nymin1,kind=sp)) 69 67 70 do i=1,numpoint 68 71 if (sglobal.and.(ypoint1(i).lt.1.e-6)) ypoint1(i)=1.e-6 69 if (nglobal.and.(ypoint2(i).gt.real(nymin1 )-1.e-5)) &70 ypoint2(i)=real(nymin1 )-1.e-571 if ((ypoint1(i).lt.1.e-6).or.(ypoint1(i).ge.real(nymin1)-1.e-6) &72 .or.(ypoint2(i).lt.1.e-6).or.(ypoint2(i).ge.real(nymin1 )-1.e-6) &72 if (nglobal.and.(ypoint2(i).gt.real(nymin1,kind=dp)-1.e-5)) & 73 ypoint2(i)=real(nymin1,kind=dp)-10*yrspc 74 if ((ypoint1(i).lt.1.e-6).or.(ypoint1(i).ge.real(nymin1,kind=dp)-1.e-6) & 75 .or.(ypoint2(i).lt.1.e-6).or.(ypoint2(i).ge.real(nymin1,kind=dp)-yrspc) & 73 76 .or.((.not.xglobal).and.((xpoint1(i).lt.1.e-6).or. & 74 (xpoint1(i).ge.real(nxmin1 )-1.e-6).or.(xpoint2(i).lt.1.e-6).or. &75 (xpoint2(i).ge.real(nxmin1 )-1.e-6)))) then77 (xpoint1(i).ge.real(nxmin1,kind=dp)-1.e-6).or.(xpoint2(i).lt.1.e-6).or. & 78 (xpoint2(i).ge.real(nxmin1,kind=dp)-1.e-6)))) then 76 79 write(*,*) ' NOTICE: RELEASE POINT OUT OF DOMAIN DETECTED.' 77 80 write(*,*) ' IT IS REMOVED NOW ... ' 78 if (i. ge.1000) then81 if (i.le.1000) then 79 82 write(*,*) ' COMMENT: ',compoint(i) 80 83 else -
src/readspecies.f90
rc8fc724 raa8c34a 67 67 character(len=16) :: pspecies 68 68 real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer 69 real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, p spec_ass, pkao69 real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pkao 70 70 real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero 71 integer :: readerror 71 integer :: readerror, pspec_ass 72 72 73 73 ! declare namelist
Note: See TracChangeset
for help on using the changeset viewer.