Skip to content

GitLab

F
flexpart
Commit df96ea65 authored by Ignacio Pisso's avatar Ignacio Pisso
Browse files
Options

resolve merge conflicts with clouds_in_GFS

parents 6d73c4be db91eb78
Hide whitespace changes
Inline Side-by-side
Showing with 136 additions and 24 deletions
src/get_wetscav.f90
View file @ df96ea65
......@@ -297,7 +297,11 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
else !parameterize cloudwater m2/m3
!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
cl=1.6E-6*prec(1)**0.36
! sec test
! cl=1E6*1E-7*prec(1)**0.3 !Sec GFS new
cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new
! cl=2E-7*prec(1)**0.36 !Andreas
! cl=1.6E-6*prec(1)**0.36 !Henrik
endif
!ZHG: Calculate the partition between liquid and water phase water.
......@@ -317,6 +321,7 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
!********
if (dquer(ks).gt.0.) then
S_i= frac_act/cl
! write(*,*) 'Si: ',S_i
! GAS
!****
......@@ -333,7 +338,9 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
else
wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
!SEC wetscav fix
wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
! wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
endif
endif ! positive in-cloud scavenging parameters given in Species file
endif !incloud
......
src/par_mod.f90
View file @ df96ea65
......@@ -147,12 +147,14 @@ module par_mod
! ECMWF
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=359 ! 1.0 degree 92 level
integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0 ! 1.0 degree 138 level
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0 ! 1.0 degree 138 level
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359 ! 1.0 degree 138 level
! integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359 ! 0.5 degree 138 level
! integer,parameter :: nxmax=181,nymax=91,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=0 ! CERA 2.0 degree 92 level
! GFS
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0
integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
integer :: nxshift=0
!*********************************************
......
src/readwind_gfs.f90
View file @ df96ea65
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
......@@ -107,6 +107,7 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
!HSO for grib api error messages
character(len=24) :: gribErrorMsg = 'Error reading grib file'
character(len=20) :: gribFunction = 'readwind_gfs'
character(len=20) :: shortname
hflswitch=.false.
......@@ -119,7 +120,7 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
! OPENING OF DATA FILE (GRIB CODE)
!HSO
5 call grib_open_file(ifile,path(3)(1:length(3)) &
call grib_open_file(ifile,path(3)(1:length(3)) &
//trim(wfname(indj)),'r',iret)
if (iret.ne.GRIB_SUCCESS) then
goto 888 ! ERROR DETECTED
......@@ -161,6 +162,8 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
else ! GRIB Edition 2
!read the grib2 identifiers
call grib_get_string(igrib,'shortName',shortname,iret)
call grib_get_int(igrib,'discipline',discipl,iret)
! call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'parameterCategory',parCat,iret)
......@@ -172,7 +175,8 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
valSurf,iret)
! call grib_check(iret,gribFunction,gribErrorMsg)
! write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
!convert to grib1 identifiers
isec1(6)=-1
isec1(7)=-1
......@@ -213,6 +217,10 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
isec1(6)=34 ! indicatorOfParameter
isec1(7)=105 ! indicatorOfTypeOfLevel
isec1(8)=10
elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
isec1(6)=153 ! indicatorOfParameter
isec1(7)=100 ! indicatorOfTypeOfLevel
isec1(8)=valSurf/100 ! level, convert to hPa
elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
isec1(6)=2 ! indicatorOfParameter
isec1(7)=102 ! indicatorOfTypeOfLevel
......@@ -547,6 +555,13 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
vlev1(i-i181,j)=help
endif
endif
! SEC & IP 12/2018 read GFS clouds
if(isec1(6).eq.153) then !! CLWCR Cloud liquid water content [kg/kg]
clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
readclouds=.true.
sumclouds=.true.
endif
end do
end do
......@@ -674,6 +689,8 @@ subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
! IP & SEC adding GFS Clouds 20181205
call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
endif
do i=0,nxmin1
......
src/verttransform_ecmwf.f90
View file @ df96ea65
......@@ -107,10 +107,11 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
!ZHG 2015 temporary variables for testing
! real :: rcw(0:nxmax-1,0:nymax-1)
! real :: rpc(0:nxmax-1,0:nymax-1)
! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
character(len=60) :: zhgpath='/xnilu_wrk/users/sec/kleinprojekte/hertlfit/'
character(len=60) :: fnameH,fnameI,fnameJ
! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
! CHARACTER(LEN=3) :: aspec
! integer :: virr=0
CHARACTER(LEN=3) :: aspec
integer :: virr=0
!real :: tot_cloud_h
!real :: dbg_height(nzmax)
!ZHG
......@@ -731,8 +732,8 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
! WRITE OUT SOME TEST VARIABLES
!********* TEST ************'**
!teller(:)=0
!virr=virr+1
!WRITE(aspec, '(i3.3)'), virr
virr=virr+1
WRITE(aspec, '(i3.3)'), virr
!if (readclouds) then
!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
......@@ -769,29 +770,34 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
if (1.eq.2) then
fnameH=trim(zhgpath)//trim(aspec)//'tcwc.txt'
fnameI=trim(zhgpath)//trim(aspec)//'prec.txt'
fnameJ=trim(zhgpath)//trim(aspec)//'cloudsh.txt'
write(*,*) 'Writing data to file: ',fnameH
!if (readclouds) then
!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
!else
!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
OPEN(UNIT=115, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
OPEN(UNIT=116, FILE=fnameI,FORM='FORMATTED',STATUS = 'UNKNOWN')
OPEN(UNIT=117, FILE=fnameJ,FORM='FORMATTED',STATUS = 'UNKNOWN')
!endif
!
!do ix=0,nxmin1
do ix=0,nxmin1
!if (readclouds) then
!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
!else
!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1) !integer
!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1) !7.83691406E-02
!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
write(115,*) (ctwc(ix,jy,n),jy=0,nymin1)
write(116,*) (lsprec(ix,jy,1,n)+convprec(ix,jy,1,n),jy=0,nymin1)
write(117,*) (cloudsh(ix,jy,n),jy=0,nymin1)
!endif
!end do
end do
!
!if (readclouds) then
!CLOSE(111)
......@@ -799,9 +805,10 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
!CLOSE(113)
!CLOSE(114)
!else
!CLOSE(115)
!CLOSE(116)
!CLOSE(117)
CLOSE(115)
CLOSE(116)
CLOSE(117)
endif
!endif
!
!END ********* TEST *************** END
......
src/verttransform_gfs.f90
View file @ df96ea65
......@@ -74,7 +74,7 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
integer :: rain_cloud_above,kz_inv
real :: f_qvsat,pressure
real :: rh,lsp,convp
real :: rh,lsp,cloudh_min,convp,prec
real :: rhoh(nuvzmax),pinmconv(nzmax)
real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
real :: xlon,ylat,xlonr,dzdx,dzdy
......@@ -223,6 +223,10 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
vv(ix,jy,1,n)=vvh(ix,jy,llev)
tt(ix,jy,1,n)=tth(ix,jy,llev,n)
qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
! IP & SEC, 201812 add clouds
if (readclouds) then
clwc(ix,jy,1,n)=clwch(ix,jy,llev,n)
endif
pv(ix,jy,1,n)=pvh(ix,jy,llev)
rho(ix,jy,1,n)=rhoh(llev)
pplev(ix,jy,1,n)=akz(llev)
......@@ -230,6 +234,10 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
! IP & SEC, 201812 add clouds
if (readclouds) then
clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
endif
pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
rho(ix,jy,nz,n)=rhoh(nuvz)
pplev(ix,jy,nz,n)=akz(nuvz)
......@@ -241,6 +249,10 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
! IP & SEC, 201812 add clouds
if (readclouds) then
clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
endif
pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
......@@ -257,6 +269,11 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
+tth(ix,jy,kz,n)*dz1)/dz
qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+qvh(ix,jy,kz,n)*dz1)/dz
! IP & SEC, 201812 add clouds
if (readclouds) then
clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2 &
+clwch(ix,jy,kz,n)*dz1)/dz
endif
pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
......@@ -495,6 +512,67 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
endif
!***********************************************************************************
! IP & SEC, 201812 GFS clouds read
if (readclouds) then
! The method is loops all grids vertically and constructs the 3D matrix for clouds
! Cloud top and cloud bottom gid cells are assigned as well as the total column
! cloud water. For precipitating grids, the type and whether it is in or below
! cloud scavenging are assigned with numbers 2-5 (following the old metod).
! Distinction is done for lsp and convp though they are treated the same in regards
! to scavenging. Also clouds that are not precipitating are defined which may be
! to include future cloud processing by non-precipitating-clouds.
!***********************************************************************************
write(*,*) 'Global NCEP fields: using cloud water'
clw(:,:,:,n)=0.0
ctwc(:,:,n)=0.0
clouds(:,:,:,n)=0
! If water/ice are read separately into clwc and ciwc, store sum in clwc
do jy=0,nymin1
do ix=0,nxmin1
lsp=lsprec(ix,jy,1,n)
convp=convprec(ix,jy,1,n)
prec=lsp+convp
! Find clouds in the vertical
do kz=1, nz-1 !go from top to bottom
if (clwc(ix,jy,kz,n).gt.0) then
! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
cloudh_min=min(height(kz+1),height(kz))
endif
end do
! If Precipitation. Define removal type in the vertical
if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
do kz=nz,1,-1 !go Bottom up!
if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
clouds(ix,jy,kz,n)=1 ! is a cloud
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=3 ! lsp in-cloud
else
clouds(ix,jy,kz,n)=2 ! convp in-cloud
endif ! convective or large scale
elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=5 ! lsp dominated washout
else
clouds(ix,jy,kz,n)=4 ! convp dominated washout
endif ! convective or large scale
endif
if (height(kz).ge. 19000) then ! set a max height for removal
clouds(ix,jy,kz,n)=0
endif !clw>0
end do !nz
endif ! precipitation
end do
end do
else
write(*,*) 'Global NCEP fields: using cloud water from Parameterization'
! write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
! create a cloud and rainout/washout field, clouds occur where rh>80%
! total cloudheight is stored at level 0
......@@ -533,6 +611,7 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
end do
end do
end do
endif ! IP & SEC 201812, GFS clouds read
end subroutine verttransform_gfs
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment