Changeset 6ecb30a in flexpart.git for src/calcpar.f90
- Timestamp:
- Aug 17, 2017, 4:39:17 PM (7 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 5b34509
- Parents:
- 61e07ba
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.