Changeset db712a8 in flexpart.git
- Timestamp:
- Mar 3, 2016, 12:34:56 PM (8 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 38b7917
- Parents:
- b0434e1
- Location:
- src
- Files:
-
- 24 edited
Legend:
- Unmodified
- Added
- Removed
-
src/FLEXPART.f90
rb0434e1 rdb712a8 315 315 endif 316 316 317 !! testing !!318 ! open(999,file=trim(path(1))//'OH_FIELDS/jscalar_50N.txt',action='write',status='new')319 ! open(998,file=trim(path(1))//'OH_FIELDS/jscalar_50S.txt',action='write',status='new')320 321 322 317 ! Write basic information on the simulation to a file "header" 323 318 ! and open files that are to be kept open throughout the simulation -
src/FLEXPART_MPI.f90
rb0434e1 rdb712a8 348 348 endif 349 349 350 !! testing !!351 ! open(999,file=trim(path(1))//'OH_FIELDS/jscalar_50N.txt',action='write',status='new')352 ! open(998,file=trim(path(1))//'OH_FIELDS/jscalar_50S.txt',action='write',status='new')353 354 355 350 ! Write basic information on the simulation to a file "header" 356 351 ! and open files that are to be kept open throughout the simulation -
src/calcpar_nests.f90
re200b7a rdb712a8 103 103 ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), & 104 104 td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l)) 105 if (ustarn(ix,jy,1,n,l).le.1.e-8) ustarn(ix,jy,1,n,l)=1.e-8 105 106 106 107 ! 2) Calculation of inverse Obukhov length scale … … 135 136 subsceff=min(excessoron(ix,jy,l),hmixplus) 136 137 else 137 subsceff=0 138 subsceff=0.0 138 139 endif 139 140 ! … … 149 150 150 151 if (DRYDEP) then 151 z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga 152 z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga 152 ! z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga 153 ! z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga 154 z0(7)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga 153 155 154 156 ! Calculate relative humidity at surface … … 228 230 end do 229 231 230 231 call calcpv_nests(l,n,uuhn,vvhn,pvhn) 232 ! Calculation of potential vorticity on 3-d grid 233 !*********************************************** 234 235 call calcpv_nests(l,n,uuhn,vvhn,pvhn) 232 236 233 237 end do -
src/calcpv.f90
r8a65cb0 rdb712a8 70 70 enddo 71 71 72 ! :DBG: halts with SIGFPE if compiling with -ffpe-trap73 ! where(ppml.le.0) ppml=10000.074 72 ! ppmk(:,:,1:nuvz)=(100000./ppml(:,:,1:nuvz))**kappa 75 ppmk =(100000./ppml)**kappa73 ppmk(0:nxmin1,0:nymin1,1:nuvz)=(100000./ppml(0:nxmin1,0:nymin1,1:nuvz))**kappa 76 74 77 75 do jy=0,nymin1 -
src/calcpv_nests.f90
r4fbe7a5 rdb712a8 68 68 enddo 69 69 enddo 70 ppmk=(100000./ppml)**kappa 70 ! ppmk=(100000./ppml)**kappa 71 ppmk(0:nxn(l)-1,0:nyn(l)-1,1:nuvz)=(100000./ppml(0:nxn(l)-1,0:nyn(l)-1,1:nuvz))**kappa 71 72 72 73 do jy=0,nyn(l)-1 -
src/com_mod.f90
rb0434e1 rdb712a8 134 134 logical :: readclouds 135 135 !ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed) 136 logical :: sumclouds 136 logical :: sumclouds 137 138 logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest 137 139 138 140 … … 370 372 !ZHG Sep 2015 371 373 real :: icloud_stats(0:nxmax-1,0:nymax-1,5,numwfmem) 372 real :: clh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)373 374 real :: clw4(0:nxmax-1,0:nymax-1,numwfmem) ! eso: =icloud_stats(:,:,4,:) 374 375 … … 487 488 488 489 real,allocatable,dimension(:,:,:,:,:) :: uun, vvn, wwn, ttn, qvn, pvn,& 489 & rhon, drhodzn, tthn, qvhn 490 integer,allocatable,dimension(:,:,:,:) :: cloudsnh 490 & rhon, drhodzn, tthn, qvhn, clwcn, ciwcn, clwn, clwchn, ciwchn 491 real,allocatable,dimension(:,:,:,:) :: clw4n 492 integer,allocatable,dimension(:,:,:,:) :: cloudshn 491 493 integer(kind=1),allocatable,dimension(:,:,:,:,:) :: cloudsn 492 494 … … 788 790 allocate(qvn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 789 791 allocate(pvn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 790 allocate(cloudsn(0:nxmaxn-1,0:nymaxn-1,0:nzmax,numwfmem,numbnests)) 791 allocate(cloudsnh(0:nxmaxn-1,0:nymaxn-1,numwfmem,numbnests)) 792 allocate(clwcn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 793 allocate(ciwcn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 794 allocate(clwn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 795 796 allocate(cloudsn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 797 allocate(cloudshn(0:nxmaxn-1,0:nymaxn-1,numwfmem,numbnests)) 792 798 allocate(rhon(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 793 799 allocate(drhodzn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests)) 794 800 allocate(tthn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests)) 795 801 allocate(qvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests)) 796 802 allocate(clwchn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests)) 803 allocate(ciwchn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests)) 804 allocate(clw4n(0:nxmax-1,0:nymax-1,numwfmem,numbnests)) 805 806 ! clw4n(:,:,:,:)=0. 807 clwcn(:,:,:,:,:)=0. 808 ciwcn(:,:,:,:,:)=0. 809 clwchn(:,:,:,:,:)=0. 810 ciwchn(:,:,:,:,:)=0. 811 ! clwn(:,:,:,:,:)=0. 812 813 ! cloudsn(:,:,:,:,:)=0 814 ! cloudshn(:,:,:,:)=0 797 815 798 816 end subroutine com_mod_allocate_nests -
src/concoutput.f90
r6a678e3 rdb712a8 626 626 !************************* 627 627 628 do ks=1,nspec 629 do kp=1,maxpointspec_act 630 do i=1,numreceptor 631 creceptor(i,ks)=0. 632 end do 633 do jy=0,numygrid-1 634 do ix=0,numxgrid-1 635 do l=1,nclassunc 636 do nage=1,nageclass 637 do kz=1,numzgrid 638 gridunc(ix,jy,kz,ks,kp,l,nage)=0. 639 end do 640 end do 641 end do 642 end do 643 end do 644 end do 645 end do 628 ! do ks=1,nspec 629 ! do kp=1,maxpointspec_act 630 ! do i=1,numreceptor 631 ! creceptor(i,ks)=0. 632 ! end do 633 ! do jy=0,numygrid-1 634 ! do ix=0,numxgrid-1 635 ! do l=1,nclassunc 636 ! do nage=1,nageclass 637 ! do kz=1,numzgrid 638 ! gridunc(ix,jy,kz,ks,kp,l,nage)=0. 639 ! end do 640 ! end do 641 ! end do 642 ! end do 643 ! end do 644 ! end do 645 ! end do 646 creceptor(:,:)=0. 647 gridunc(:,:,:,:,:,:,:)=0. 646 648 647 649 -
src/concoutput_nest.f90
r6a678e3 rdb712a8 545 545 !************************* 546 546 547 do ks=1,nspec 548 do kp=1,maxpointspec_act 549 do i=1,numreceptor 550 creceptor(i,ks)=0. 551 end do 552 do jy=0,numygridn-1 553 do ix=0,numxgridn-1 554 do l=1,nclassunc 555 do nage=1,nageclass 556 do kz=1,numzgrid 557 griduncn(ix,jy,kz,ks,kp,l,nage)=0. 558 end do 559 end do 560 end do 561 end do 562 end do 563 end do 564 end do 547 ! do ks=1,nspec 548 ! do kp=1,maxpointspec_act 549 ! do i=1,numreceptor 550 ! creceptor(i,ks)=0. 551 ! end do 552 ! do jy=0,numygridn-1 553 ! do ix=0,numxgridn-1 554 ! do l=1,nclassunc 555 ! do nage=1,nageclass 556 ! do kz=1,numzgrid 557 ! griduncn(ix,jy,kz,ks,kp,l,nage)=0. 558 ! end do 559 ! end do 560 ! end do 561 ! end do 562 ! end do 563 ! end do 564 ! end do 565 creceptor(:,:)=0. 566 griduncn(:,:,:,:,:,:,:)=0. 565 567 566 568 -
src/ecmwf_mod.f90
rb0434e1 rdb712a8 41 41 !********************************************* 42 42 43 integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 43 ! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new 44 integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 44 45 integer,parameter :: nxshift=359 46 ! integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !test BUG 47 ! integer,parameter :: nxshift=0 45 48 46 49 ! … … 49 52 !********************************************* 50 53 51 integer,parameter :: maxnests=1,nxmaxn=361,nymaxn=181 !ECMWF 54 integer,parameter :: maxnests=1,nxmaxn=361,nymaxn=181 55 ! integer,parameter :: maxnests=1,nxmaxn=86,nymaxn=31 52 56 53 57 -
src/gridcheck.f90
rb0434e1 rdb712a8 194 194 elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc 195 195 isec1(6)=246 ! indicatorOfParameter 196 ! readclouds=.true.197 ! sumclouds=.false.198 196 elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc 199 197 isec1(6)=247 ! indicatorOfParameter … … 202 200 elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc 203 201 isec1(6)=201031 ! indicatorOfParameter 204 ! readclouds=.true.205 ! sumclouds=.true.206 202 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 207 203 isec1(6)=134 ! indicatorOfParameter … … 366 362 endif ! gotGrid 367 363 364 if (nx.gt.nxmax) then 365 write(*,*) 'FLEXPART error: Too many grid points in x direction.' 366 write(*,*) 'Reduce resolution of wind fields.' 367 write(*,*) 'Or change parameter settings in file par_mod.' 368 write(*,*) nx,nxmax 369 stop 370 endif 371 372 if (ny.gt.nymax) then 373 write(*,*) 'FLEXPART error: Too many grid points in y direction.' 374 write(*,*) 'Reduce resolution of wind fields.' 375 write(*,*) 'Or change parameter settings in file par_mod.' 376 write(*,*) ny,nymax 377 stop 378 endif 379 368 380 k=isec1(8) 369 381 if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1) … … 410 422 nwz =iwmax 411 423 if(nuvz.eq.nlev_ec) nwz=nlev_ec+1 412 413 if (nx.gt.nxmax) then414 write(*,*) 'FLEXPART error: Too many grid points in x direction.'415 write(*,*) 'Reduce resolution of wind fields.'416 write(*,*) 'Or change parameter settings in file par_mod.'417 write(*,*) nx,nxmax418 stop419 endif420 421 if (ny.gt.nymax) then422 write(*,*) 'FLEXPART error: Too many grid points in y direction.'423 write(*,*) 'Reduce resolution of wind fields.'424 write(*,*) 'Or change parameter settings in file par_mod.'425 write(*,*) ny,nymax426 stop427 endif428 424 429 425 if (nuvz+1.gt.nuvzmax) then -
src/gridcheck_nests.f90
rb0434e1 rdb712a8 172 172 elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc 173 173 isec1(6)=246 ! indicatorOfParameter 174 ! readclouds=.true.175 ! sumclouds=.false.176 174 elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc 177 175 isec1(6)=247 ! indicatorOfParameter … … 180 178 elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc 181 179 isec1(6)=201031 ! indicatorOfParameter 182 ! readclouds=.true.183 ! sumclouds=.true.184 180 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 185 181 isec1(6)=134 ! indicatorOfParameter 186 182 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 187 183 isec1(6)=135 ! indicatorOfParameter 188 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridche k.f90189 isec1(6)=135 ! indicatorOfParameter ! ! " " " " " " " " " " " " " " " " " " " " " " " " " " "184 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridcheck.f90 185 isec1(6)=135 ! indicatorOfParameter ! 190 186 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 191 187 isec1(6)=151 ! indicatorOfParameter … … 258 254 endif ! ifield 259 255 256 if (nxn(l).gt.nxmaxn) then 257 write(*,*) 'FLEXPART error: Too many grid points in x direction.' 258 write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)' 259 write(*,*) 'for nesting level ',l 260 write(*,*) 'Or change parameter settings in file par_mod.' 261 write(*,*) nxn(l),nxmaxn 262 stop 263 endif 264 265 if (nyn(l).gt.nymaxn) then 266 write(*,*) 'FLEXPART error: Too many grid points in y direction.' 267 write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)' 268 write(*,*) 'for nesting level ',l 269 write(*,*) 'Or change parameter settings in file par_mod.' 270 write(*,*) nyn(l),nymaxn 271 stop 272 endif 273 260 274 !HSO get the second part of the grid dimensions only from GRiB1 messages 261 275 if (isec1(6) .eq. 167 .and. (gotGrib.eq.0)) then !added by mc to make it consistent with new gridchek.f90 note that gotGrid must be changed in gotGrib!! … … 333 347 if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1 334 348 335 if (nxn(l).gt.nxmaxn) then336 write(*,*) 'FLEXPART error: Too many grid points in x direction.'337 write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'338 write(*,*) 'for nesting level ',l339 write(*,*) 'Or change parameter settings in file par_mod.'340 write(*,*) nxn(l),nxmaxn341 stop342 endif343 344 if (nyn(l).gt.nymaxn) then345 write(*,*) 'FLEXPART error: Too many grid points in y direction.'346 write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'347 write(*,*) 'for nesting level ',l348 write(*,*) 'Or change parameter settings in file par_mod.'349 write(*,*) nyn(l),nymaxn350 stop351 endif352 353 349 if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax)) then 354 350 write(*,*) 'FLEXPART error: Nested wind fields have too many'// & -
src/initialize.f90
r8a65cb0 rdb712a8 99 99 jyp=jy+1 100 100 101 h=max(hmix(ix ,jy 101 h=max(hmix(ix ,jy,1,memind(1)), & 102 102 hmix(ixp,jy ,1,memind(1)), & 103 103 hmix(ix ,jyp,1,memind(1)), & … … 157 157 wp=rannumb(nrand+2) 158 158 if (.not.turbswitch) then ! modified by mc 159 wp=wp*sigw 160 else if (cblflag.eq.1) then ! modified by mc 161 if(-h/ol.gt.5) then 162 !if (ol.lt.0.) then 163 !if (ol.gt.0.) then !by mc : only for test correct is lt.0 164 call initialize_cbl_vel(idummy,zt,ust,wst,h,sigw,wp,ol) 165 else 159 166 wp=wp*sigw 160 else if (cblflag.eq.1) then ! modified by mc 161 if(-h/ol.gt.5) then 162 !if (ol.lt.0.) then 163 !if (ol.gt.0.) then !by mc : only for test correct is lt.0 164 call initialize_cbl_vel(idummy,zt,ust,wst,h,sigw,wp,ol) 165 else 166 wp=wp*sigw 167 end if 167 end if 168 168 end if 169 169 -
src/interpol_rain_nests.f90
r8a65cb0 rdb712a8 88 88 !*********************************************************** 89 89 90 if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) & 91 xt=real(nxn(ngrid)-1)-0.0001 92 if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) & 93 yt=real(nyn(ngrid)-1)-0.0001 90 ! if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) & 91 ! xt=real(nxn(ngrid)-1)-0.0001 92 ! if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) & 93 ! yt=real(nyn(ngrid)-1)-0.0001 94 95 ! ESO make it consistent with interpol_rain 96 if (xt.ge.(real(nxn(ngrid)-1))) xt=real(nxn(ngrid)-1)-0.00001 97 if (yt.ge.(real(nyn(ngrid)-1))) yt=real(nyn(ngrid)-1)-0.00001 94 98 95 99 … … 105 109 ix=int(xt) 106 110 jy=int(yt) 111 107 112 ixp=ix+1 108 113 jyp=jy+1 -
src/makefile
rb0434e1 rdb712a8 388 388 readdepo.o: com_mod.o par_mod.o 389 389 readlanduse.o: com_mod.o par_mod.o 390 readlanduse_int1.o: com_mod.o par_mod.o390 #readlanduse_int1.o: com_mod.o par_mod.o 391 391 readOHfield.o: com_mod.o oh_mod.o par_mod.o 392 392 readoutgrid.o: com_mod.o outg_mod.o par_mod.o -
src/mpi_mod.f90
r6a678e3 rdb712a8 860 860 ! step 861 861 ! 862 ! TODO863 ! GFS version864 862 ! 865 863 !******************************************************************************* … … 922 920 923 921 ! The non-reader processes need to know if cloud water were read. 924 ! TODO: only at first step or always?925 922 call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr) 926 923 if (mp_ierr /= 0) goto 600 … … 1017 1014 if (mp_ierr /= 0) goto 600 1018 1015 1019 call MPI_Bcast(z0,numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)1020 if (mp_ierr /= 0) goto 6001021 1022 1016 if (mp_measure_time) call mpif_mtime('commtime',1) 1023 1017 … … 1045 1039 ! step 1046 1040 ! 1047 ! TODO1048 ! Transfer cloud water/ice if and when available for nested1049 1041 ! 1050 1042 !*********************************************************************** … … 1060 1052 1061 1053 ! Common array sizes used for communications 1062 integer :: d3_size1 = nxmaxn*nymaxn*nzmax *maxnests1063 integer :: d3_size2 = nxmaxn*nymaxn*nuvzmax *maxnests1064 integer :: d2_size1 = nxmaxn*nymaxn *maxnests1065 integer :: d2_size2 = nxmaxn*nymaxn*maxspec *maxnests1066 integer :: d2_size3 = nxmaxn*nymaxn *maxnests1054 integer :: d3_size1 = nxmaxn*nymaxn*nzmax 1055 integer :: d3_size2 = nxmaxn*nymaxn*nuvzmax 1056 integer :: d2_size1 = nxmaxn*nymaxn 1057 integer :: d2_size2 = nxmaxn*nymaxn*maxspec 1058 integer :: d2_size3 = nxmaxn*nymaxn 1067 1059 1068 1060 integer :: d3s1,d3s2,d2s1,d2s2 … … 1106 1098 !********************************************************************** 1107 1099 1100 ! The non-reader processes need to know if cloud water were read. 1101 call MPI_Bcast(readclouds_nest,maxnests,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr) 1102 if (mp_ierr /= 0) goto 600 1103 1108 1104 ! Static fields/variables sent only at startup 1109 1105 if (first_call) then … … 1120 1116 1121 1117 ! MPI prefers contiguous arrays for sending (else a buffer is created), 1122 ! hence the loop 1123 1118 ! hence the loop over nests 1119 !********************************************************************** 1124 1120 do i=1, numbnests 1125 1121 ! 3D fields … … 1145 1141 if (mp_ierr /= 0) goto 600 1146 1142 call MPI_Bcast(cloudsn(:,:,:,li:ui,i),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr) 1147 if (mp_ierr /= 0) goto 600 1143 if (mp_ierr /= 0) goto 600 1144 1145 ! cloud water/ice: 1146 if (readclouds_nest(i)) then 1147 ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr) 1148 ! if (mp_ierr /= 0) goto 600 1149 call MPI_Bcast(clw4n(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr) 1150 if (mp_ierr /= 0) goto 600 1151 end if 1148 1152 1149 1153 ! 2D fields 1150 call MPI_Bcast(clouds nh(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)1154 call MPI_Bcast(cloudshn(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr) 1151 1155 if (mp_ierr /= 0) goto 600 1152 1156 call MPI_Bcast(vdepn(:,:,:,li:ui,i),d2s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr) … … 1206 1210 ! mind -- index where to place new fields 1207 1211 ! 1208 ! TODO1209 1212 ! 1210 1213 !******************************************************************************* … … 1392 1395 ! memstat -- input, used to resolve windfield index being received 1393 1396 ! 1394 ! TODO1395 1397 ! 1396 1398 !******************************************************************************* … … 2019 2021 ! In the implementation with 3 fields, the processes may have posted 2020 2022 ! MPI_Irecv requests that should be cancelled here 2021 !! TODO:2022 2023 ! if (.not.lmp_sync) then 2023 2024 ! r=mp_pid*nvar_async … … 2105 2106 ! 2106 2107 ! 2107 ! TODO2108 2108 ! 2109 2109 !******************************************************************************* … … 2132 2132 uu(:,:,:,li:ui) = 10.0 2133 2133 vv(:,:,:,li:ui) = 0.0 2134 uupol(:,:,:,li:ui) = 10.0 ! TODO check if ok2134 uupol(:,:,:,li:ui) = 10.0 2135 2135 vvpol(:,:,:,li:ui)=0.0 2136 2136 ww(:,:,:,li:ui)=0. … … 2164 2164 tropopause(:,:,:,li:ui)=10000. 2165 2165 oli(:,:,:,li:ui)=0.01 2166 z0=1.02167 2166 2168 2167 end subroutine set_fields_synthetic -
src/netcdf_output_mod.f90
r6a678e3 rdb712a8 79 79 implicit none 80 80 81 private 82 83 public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,& 84 &concoutput_nest_netcdf,concoutput_surf_netcdf 85 81 86 ! include 'netcdf.inc' 82 87 … … 97 102 real,parameter :: eps=nxmax/3.e5 98 103 99 private:: writemetadata, output_units, nf90_err104 ! private:: writemetadata, output_units, nf90_err 100 105 101 106 ! switch output of release point information on/off … … 628 633 629 634 ! volume 630 call nf90_err(nf90_put_var(ncid, volID, volume(:,:,:))) 635 if (lnest) then 636 call nf90_err(nf90_put_var(ncid, volID, volumen(:,:,:))) 637 else 638 call nf90_err(nf90_put_var(ncid, volID, volume(:,:,:))) 639 end if 631 640 632 641 ! area 633 call nf90_err(nf90_put_var(ncid, areaID, area(:,:))) 642 if (lnest) then 643 call nf90_err(nf90_put_var(ncid, areaID, arean(:,:))) 644 else 645 call nf90_err(nf90_put_var(ncid, areaID, area(:,:))) 646 end if 634 647 635 648 if (write_releases.eqv..true.) then -
src/par_mod.f90
rb0434e1 rdb712a8 139 139 !integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 140 140 !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF 141 !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26142 141 !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64 143 !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58144 142 145 143 ! integer,parameter :: nxshift=359 ! for ECMWF … … 216 214 !************************************************** 217 215 218 integer,parameter :: maxpart=400000 00216 integer,parameter :: maxpart=400000 219 217 integer,parameter :: maxspec=6 220 integer,parameter :: minmass=0.0 !0.0001218 real,parameter :: minmass=0.0 !0.0001 221 219 222 220 ! maxpart Maximum number of particles -
src/readavailable.f90
r5f9d14a rdb712a8 283 283 stop 284 284 285 999 write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '285 999 write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE #### ' 286 286 write(*,'(a)') ' '//path(4)(1:length(4)) 287 287 write(*,*) ' #### CANNOT BE OPENED #### ' -
src/readwind_nests.f90
rb0434e1 rdb712a8 347 347 ! -add check for if one of clwc/ciwc missing (error), 348 348 ! also if all 3 cw fields present, use qc and disregard the others 349 ! -use same flags readclouds/sumclouds as in mother grid? this assumes350 ! that both the nested and mother grids contain CW in same format351 349 if(isec1(6).eq.246) then !! CLWC Cloud liquid water content [kg/kg] 352 clwch (i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)353 ! readclouds=.true.354 ! sumclouds=.false.350 clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1) 351 readclouds_nest(l)=.true. 352 sumclouds_nest(l)=.false. 355 353 endif 356 354 if(isec1(6).eq.247) then !! CIWC Cloud ice water content 357 ciwch (i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)355 ciwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1) 358 356 endif 359 357 !ZHG end 360 358 !ESO read qc (=clwc+ciwc) 361 359 if(isec1(6).eq.201031) then !! QC Cloud liquid water content [kg/kg] 362 clwch (i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)363 ! readclouds=.true.364 ! sumclouds=.true.360 clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1) 361 readclouds_nest(l)=.true. 362 sumclouds_nest(l)=.true. 365 363 endif 366 364 -
src/timemanager_mpi.f90
r6a678e3 rdb712a8 237 237 if (lmp_sync.and.lmp_use_reader.and.memstat.gt.0) then 238 238 call mpif_gf_send_vars(memstat) 239 call mpif_gf_send_vars_nest(memstat)239 if (numbnests>0) call mpif_gf_send_vars_nest(memstat) 240 240 ! Version 2 (lmp_sync=.false., see below) is also used whenever 2 new fields are 241 241 ! read (as at first time step), in which case async send/recv is impossible. 242 242 else if (.not.lmp_sync.and.lmp_use_reader.and.memstat.ge.32) then 243 243 call mpif_gf_send_vars(memstat) 244 call mpif_gf_send_vars_nest(memstat)244 if (numbnests>0) call mpif_gf_send_vars_nest(memstat) 245 245 end if 246 246 … … 260 260 ! Issued at start of each new field period. 261 261 if (memstat.ne.0.and.memstat.lt.32.and.lmp_use_reader) then 262 ! TODO: z0(7) changes with time, so should be dimension (numclass,2) to263 ! allow transfer of the future value in the background264 call MPI_Bcast(z0,numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)265 262 call mpif_gf_request 266 263 end if -
src/verttransform.f90
rb0434e1 rdb712a8 67 67 use com_mod 68 68 use cmapf_mod, only: cc2gll 69 ! eso TODO:70 ! only used for timing of CPU measurement. Remove this (and calls to mpif_mtime below)71 ! as this routine should not be dependent on MPI72 69 ! use mpi_mod 73 ! :TODO74 70 75 71 implicit none 76 72 77 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym 78 integer :: rain_cloud_above(0:nxmax-1,0:nymax-1),kz_inv,idx(0:nxmax-1,0:nymax-1) 79 real :: f_qvsat,pressure,cosf(0:nymax-1) 80 real :: rh,lsp,convp,tim,tim2,rhmin,precmin,prec 81 real :: uvzlev(0:nxmax-1,0:nymax-1,nuvzmax),rhoh(0:nxmax-1,0:nymax-1,nuvzmax) 82 real :: pinmconv(0:nxmax-1,0:nymax-1,nzmax) 83 real :: ew,pint(0:nxmax-1,0:nymax-1),tv(0:nxmax-1,0:nymax-1) 84 real :: tvold(0:nxmax-1,0:nymax-1),pold(0:nxmax-1,0:nymax-1),dz1,dz2,dz,ui,vi 73 real,intent(in),dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: uuh,vvh,pvh 74 real,intent(in),dimension(0:nxmax-1,0:nymax-1,nwzmax) :: wwh 75 76 real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev 77 real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv 78 real,dimension(0:nxmax-1,0:nymax-1) :: tvold,pold,pint,tv 79 real,dimension(0:nymax-1) :: cosf 80 81 integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx 82 83 integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv 84 real :: f_qvsat,pressure,rh,lsp,convp,prec 85 real :: ew,dz1,dz2,dz 85 86 real :: xlon,ylat,xlonr,dzdx,dzdy 86 87 real :: dzdx1,dzdx2,dzdy1,dzdy2 87 88 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy 88 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)89 real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)90 real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)91 real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)92 real :: wzlev(0:nxmax-1,0:nymax-1,nwzmax)93 89 real,parameter :: const=r_air/ga 94 ! integer:: ihr, imin, isec, i100th,ihr2, imin2, isec2, i100th2 95 parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics 90 real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics 96 91 97 92 logical :: init = .true. 98 93 99 94 !ZHG SEP 2014 tests 100 integer :: cloud_ver,cloud_min, cloud_max101 integer ::teller(5), convpteller=0, lspteller=0102 real :: cloud_col_wat, cloud_water95 ! integer :: cloud_ver,cloud_min, cloud_max 96 ! integer ::teller(5), convpteller=0, lspteller=0 97 ! real :: cloud_col_wat, cloud_water 103 98 !ZHG 2015 temporary variables for testing 104 real :: rcw(0:nxmax-1,0:nymax-1)105 real :: rpc(0:nxmax-1,0:nymax-1)99 ! real :: rcw(0:nxmax-1,0:nymax-1) 100 ! real :: rpc(0:nxmax-1,0:nymax-1) 106 101 ! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/' 107 102 ! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH 108 CHARACTER(LEN=3) :: aspec109 integer :: virr=0103 ! CHARACTER(LEN=3) :: aspec 104 ! integer :: virr=0 110 105 real :: tot_cloud_h 106 real :: dbg_height(nzmax) 111 107 !ZHG 112 108 … … 143 139 144 140 145 tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew (td2(ixm,jym,1,n))/ &141 tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew*(td2(ixm,jym,1,n))/ & 146 142 ps(ixm,jym,1,n)) 147 143 pold(ixm,jym)=ps(ixm,jym,1,n) … … 182 178 init=.false. 183 179 180 dbg_height = height 181 184 182 endif 185 183 … … 191 189 do jy=0,nymin1 192 190 do ix=0,nxmin1 193 tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew (td2(ix,jy,1,n))/ &191 tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew*(td2(ix,jy,1,n))/ & 194 192 ps(ix,jy,1,n)) 195 193 enddo … … 250 248 qv(:,:,1,n)=qvh(:,:,1,n) 251 249 !hg adding the cloud water 252 clwc(:,:,1,n)=clwch(:,:,1,n) 253 if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n) 250 if (readclouds) then 251 clwc(:,:,1,n)=clwch(:,:,1,n) 252 if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n) 253 end if 254 254 !hg 255 255 pv(:,:,1,n)=pvh(:,:,1) 256 256 rho(:,:,1,n)=rhoh(:,:,1) 257 257 258 uu(:,:,nz,n)=uuh(:,:,nuvz) 258 259 vv(:,:,nz,n)=vvh(:,:,nuvz) 259 260 tt(:,:,nz,n)=tth(:,:,nuvz,n) 260 261 qv(:,:,nz,n)=qvh(:,:,nuvz,n) 261 262 262 !hg adding the cloud water 263 clwc(:,:,nz,n)=clwch(:,:,nuvz,n) 264 if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n) 263 if (readclouds) then 264 clwc(:,:,nz,n)=clwch(:,:,nuvz,n) 265 if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n) 266 end if 265 267 !hg 266 268 pv(:,:,nz,n)=pvh(:,:,nuvz) … … 279 281 qv(ix,jy,iz,n)=qv(ix,jy,nz,n) 280 282 !hg adding the cloud water 281 clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n) 282 if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n) 283 if (readclouds) then 284 clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n) 285 if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n) 286 end if 283 287 !hg 284 288 pv(ix,jy,iz,n)=pv(ix,jy,nz,n) … … 309 313 +qvh(ix,jy,kz,n)*dz1)/dz 310 314 !hg adding the cloud water 311 clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz 312 if (.not.sumclouds) & 313 &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz 315 if (readclouds) then 316 clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz 317 if (.not.sumclouds) & 318 &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz 319 end if 314 320 !hg 315 321 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz … … 601 607 icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz)) ! Cloud BOT height stats [m] 602 608 !ZHG 2015 extra for testing 603 clh(ix,jy,kz,n)=height(kz+1)-height(kz)609 ! clh(ix,jy,kz,n)=height(kz+1)-height(kz) 604 610 icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m] 605 611 icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz)) ! Cloud TOP height [m] -
src/verttransform_nests.f90
rb0434e1 rdb712a8 21 21 22 22 subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn) 23 ! i i i i i 24 !***************************************************************************** 25 ! * 26 ! This subroutine transforms temperature, dew point temperature and * 27 ! wind components from eta to meter coordinates. * 28 ! The vertical wind component is transformed from Pa/s to m/s using * 29 ! the conversion factor pinmconv. * 30 ! In addition, this routine calculates vertical density gradients * 31 ! needed for the parameterization of the turbulent velocities. * 32 ! It is similar to verttransform, but makes the transformations for * 33 ! the nested grids. * 34 ! * 35 ! Author: A. Stohl, G. Wotawa * 36 ! * 37 ! 12 August 1996 * 38 ! Update: 16 January 1998 * 39 ! * 40 ! Major update: 17 February 1999 * 41 ! by G. Wotawa * 42 ! * 43 ! - Vertical levels for u, v and w are put together * 44 ! - Slope correction for vertical velocity: Modification of calculation * 45 ! procedure * 46 ! * 47 !***************************************************************************** 48 ! Changes, Bernd C. Krueger, Feb. 2001: (marked "C-cv") 49 ! Variables tthn and qvhn (on eta coordinates) from common block 50 !***************************************************************************** 51 ! Sabine Eckhardt, March 2007 52 ! add the variable cloud for use with scavenging - descr. in com_mod 53 !***************************************************************************** 54 ! * 55 ! Variables: * 56 ! nxn,nyn,nuvz,nwz field dimensions in x,y and z direction * 57 ! uun wind components in x-direction [m/s] * 58 ! vvn wind components in y-direction [m/s] * 59 ! wwn wind components in z-direction [deltaeta/s]* 60 ! ttn temperature [K] * 61 ! pvn potential vorticity (pvu) * 62 ! psn surface pressure [Pa] * 63 ! * 64 !***************************************************************************** 23 ! i i i i i 24 !***************************************************************************** 25 ! * 26 ! This subroutine transforms temperature, dew point temperature and * 27 ! wind components from eta to meter coordinates. * 28 ! The vertical wind component is transformed from Pa/s to m/s using * 29 ! the conversion factor pinmconv. * 30 ! In addition, this routine calculates vertical density gradients * 31 ! needed for the parameterization of the turbulent velocities. * 32 ! It is similar to verttransform, but makes the transformations for * 33 ! the nested grids. * 34 ! * 35 ! Author: A. Stohl, G. Wotawa * 36 ! * 37 ! 12 August 1996 * 38 ! Update: 16 January 1998 * 39 ! * 40 ! Major update: 17 February 1999 * 41 ! by G. Wotawa * 42 ! * 43 ! - Vertical levels for u, v and w are put together * 44 ! - Slope correction for vertical velocity: Modification of calculation * 45 ! procedure * 46 ! * 47 !***************************************************************************** 48 ! Changes, Bernd C. Krueger, Feb. 2001: (marked "C-cv") 49 ! Variables tthn and qvhn (on eta coordinates) from common block 50 !***************************************************************************** 51 ! Sabine Eckhardt, March 2007 52 ! add the variable cloud for use with scavenging - descr. in com_mod 53 !***************************************************************************** 54 ! ESO, 2016 55 ! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than 56 ! the actual field dimensions 57 !***************************************************************************** 58 ! * 59 ! Variables: * 60 ! nxn,nyn,nuvz,nwz field dimensions in x,y and z direction * 61 ! uun wind components in x-direction [m/s] * 62 ! vvn wind components in y-direction [m/s] * 63 ! wwn wind components in z-direction [deltaeta/s]* 64 ! ttn temperature [K] * 65 ! pvn potential vorticity (pvu) * 66 ! psn surface pressure [Pa] * 67 ! * 68 !***************************************************************************** 65 69 66 70 use par_mod … … 69 73 implicit none 70 74 71 integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp 72 integer :: rain_cloud_above,kz_inv 73 real :: f_qvsat,pressure,rh,lsp,convp 74 real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax) 75 real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax) 76 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf 75 real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: uuhn,vvhn,pvhn 76 real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn 77 78 real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev 79 real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv 80 real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv 81 real,dimension(0:nymaxn-1) :: cosf 82 83 integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx 84 85 integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv 86 real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec 87 88 real :: ew,dz1,dz2,dz 77 89 real :: dzdx,dzdy 78 90 real :: dzdx1,dzdx2,dzdy1,dzdy2 79 real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)80 real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)81 real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)82 real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)83 91 real,parameter :: const=r_air/ga 84 85 86 ! Loop over all nests 87 !******************** 92 real :: tot_cloud_h 93 integer :: nxm1, nym1 94 95 ! real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics 96 97 ! Loop over all nests 98 !******************** 88 99 89 100 do l=1,numbnests 90 91 ! Loop over the whole grid 92 !************************* 93 94 do jy=0,nyn(l)-1 95 do ix=0,nxn(l)-1 96 97 tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ & 98 psn(ix,jy,1,n,l)) 99 pold=psn(ix,jy,1,n,l) 100 uvwzlev(ix,jy,1)=0. 101 wzlev(1)=0. 102 rhoh(1)=pold/(r_air*tvold) 103 104 105 ! Compute heights of eta levels 106 !****************************** 107 108 do kz=2,nuvz 109 pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l) 110 tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l)) 111 rhoh(kz)=pint/(r_air*tv) 112 113 if (abs(tv-tvold).gt.0.2) then 114 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* & 115 (tv-tvold)/log(tv/tvold) 116 else 117 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv 118 endif 119 120 tvold=tv 121 pold=pint 101 nxm1=nxn(l)-1 102 nym1=nyn(l)-1 103 104 ! Loop over the whole grid 105 !************************* 106 107 do jy=0,nyn(l)-1 108 do ix=0,nxn(l)-1 109 tvold(ix,jy)=tt2n(ix,jy,1,n,l)*(1.+0.378*ew*(td2n(ix,jy,1,n,l))/ & 110 psn(ix,jy,1,n,l)) 122 111 end do 123 124 125 do kz=2,nwz-1 126 wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2. 127 end do 128 wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1) 129 130 ! pinmconv=(h2-h1)/(p2-p1) 131 132 pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ & 133 ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- & 134 (aknew(1)+bknew(1)*psn(ix,jy,1,n,l))) 135 do kz=2,nz-1 136 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ & 137 ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- & 138 (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l))) 139 end do 140 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ & 141 ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- & 142 (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l))) 143 144 145 ! Levels, where u,v,t and q are given 146 !************************************ 147 148 uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l) 149 vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l) 150 ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l) 151 qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l) 152 pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l) 153 rhon(ix,jy,1,n,l)=rhoh(1) 154 uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l) 155 vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l) 156 ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l) 157 qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l) 158 pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l) 159 rhon(ix,jy,nz,n,l)=rhoh(nuvz) 160 kmin=2 161 do iz=2,nz-1 162 do kz=kmin,nuvz 163 if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 112 end do 113 114 pold(0:nxm1,0:nym1)=psn(0:nxm1,0:nym1,1,n,l) 115 uvzlev(0:nxm1,0:nym1,1)=0. 116 wzlev(0:nxm1,0:nym1,1)=0. 117 rhohn(0:nxm1,0:nym1,1)=pold(0:nxm1,0:nym1)/(r_air*tvold(0:nxm1,0:nym1)) 118 119 ! Compute heights of eta levels 120 !****************************** 121 122 do kz=2,nuvz 123 pint(0:nxm1,0:nym1)=akz(kz)+bkz(kz)*psn(0:nxm1,0:nym1,1,n,l) 124 tv(0:nxm1,0:nym1)=tthn(0:nxm1,0:nym1,kz,n,l)*(1.+0.608*qvhn(0:nxm1,0:nym1,kz,n,l)) 125 rhohn(0:nxm1,0:nym1,kz)=pint(0:nxm1,0:nym1)/(r_air*tv(0:nxm1,0:nym1)) 126 127 where (abs(tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1)).gt.0.2) 128 uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*& 129 &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))* & 130 (tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1))/& 131 &log(tv(0:nxm1,0:nym1)/tvold(0:nxm1,0:nym1)) 132 elsewhere 133 uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*& 134 &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))*tv(0:nxm1,0:nym1) 135 endwhere 136 137 tvold(0:nxm1,0:nym1)=tv(0:nxm1,0:nym1) 138 pold(0:nxm1,0:nym1)=pint(0:nxm1,0:nym1) 139 140 end do 141 142 do kz=2,nwz-1 143 wzlev(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)+uvzlev(0:nxm1,0:nym1,kz))/2. 144 end do 145 wzlev(0:nxm1,0:nym1,nwz)=wzlev(0:nxm1,0:nym1,nwz-1)+ & 146 uvzlev(0:nxm1,0:nym1,nuvz)-uvzlev(0:nxm1,0:nym1,nuvz-1) 147 148 149 pinmconv(0:nxm1,0:nym1,1)=(uvzlev(0:nxm1,0:nym1,2))/ & 150 ((aknew(2)+bknew(2)*psn(0:nxm1,0:nym1,1,n,l))- & 151 (aknew(1)+bknew(1)*psn(0:nxm1,0:nym1,1,n,l))) 152 do kz=2,nz-1 153 pinmconv(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)-uvzlev(0:nxm1,0:nym1,kz-1))/ & 154 ((aknew(kz+1)+bknew(kz+1)*psn(0:nxm1,0:nym1,1,n,l))- & 155 (aknew(kz-1)+bknew(kz-1)*psn(0:nxm1,0:nym1,1,n,l))) 156 end do 157 pinmconv(0:nxm1,0:nym1,nz)=(uvzlev(0:nxm1,0:nym1,nz)-uvzlev(0:nxm1,0:nym1,nz-1))/ & 158 ((aknew(nz)+bknew(nz)*psn(0:nxm1,0:nym1,1,n,l))- & 159 (aknew(nz-1)+bknew(nz-1)*psn(0:nxm1,0:nym1,1,n,l))) 160 161 ! Levels, where u,v,t and q are given 162 !************************************ 163 164 uun(0:nxm1,0:nym1,1,n,l)=uuhn(0:nxm1,0:nym1,1,l) 165 vvn(0:nxm1,0:nym1,1,n,l)=vvhn(0:nxm1,0:nym1,1,l) 166 ttn(0:nxm1,0:nym1,1,n,l)=tthn(0:nxm1,0:nym1,1,n,l) 167 qvn(0:nxm1,0:nym1,1,n,l)=qvhn(0:nxm1,0:nym1,1,n,l) 168 if (readclouds_nest(l)) then 169 clwcn(0:nxm1,0:nym1,1,n,l)=clwchn(0:nxm1,0:nym1,1,n,l) 170 if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,1,n,l)=ciwchn(0:nxm1,0:nym1,1,n,l) 171 end if 172 pvn(0:nxm1,0:nym1,1,n,l)=pvhn(0:nxm1,0:nym1,1,l) 173 rhon(0:nxm1,0:nym1,1,n,l)=rhohn(0:nxm1,0:nym1,1) 174 175 uun(0:nxm1,0:nym1,nz,n,l)=uuhn(0:nxm1,0:nym1,nuvz,l) 176 vvn(0:nxm1,0:nym1,nz,n,l)=vvhn(0:nxm1,0:nym1,nuvz,l) 177 ttn(0:nxm1,0:nym1,nz,n,l)=tthn(0:nxm1,0:nym1,nuvz,n,l) 178 qvn(0:nxm1,0:nym1,nz,n,l)=qvhn(0:nxm1,0:nym1,nuvz,n,l) 179 if (readclouds_nest(l)) then 180 clwcn(0:nxm1,0:nym1,nz,n,l)=clwchn(0:nxm1,0:nym1,nuvz,n,l) 181 if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,nz,n,l)=ciwchn(0:nxm1,0:nym1,nuvz,n,l) 182 end if 183 pvn(0:nxm1,0:nym1,nz,n,l)=pvhn(0:nxm1,0:nym1,nuvz,l) 184 rhon(0:nxm1,0:nym1,nz,n,l)=rhohn(0:nxm1,0:nym1,nuvz) 185 186 187 kmin=2 188 idx=kmin 189 do iz=2,nz-1 190 do jy=0,nyn(l)-1 191 do ix=0,nxn(l)-1 192 if(height(iz).gt.uvzlev(ix,jy,nuvz)) then 164 193 uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l) 165 194 vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l) 166 195 ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l) 167 196 qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l) 197 !hg adding the cloud water 198 if (readclouds_nest(l)) then 199 clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l) 200 if (.not.sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l) 201 end if 202 !hg 168 203 pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l) 169 204 rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l) 170 goto 30 205 else 206 innuvz: do kz=idx(ix,jy),nuvz 207 if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. & 208 (height(iz).le.uvzlev(ix,jy,kz))) then 209 idx(ix,jy)=kz 210 exit innuvz 211 endif 212 enddo innuvz 171 213 endif 172 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 173 (height(iz).le.uvwzlev(ix,jy,kz))) then 174 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 175 dz2=uvwzlev(ix,jy,kz)-height(iz) 176 dz=dz1+dz2 177 uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ & 178 uuhn(ix,jy,kz,l)*dz1)/dz 179 vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ & 180 vvhn(ix,jy,kz,l)*dz1)/dz 181 ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ & 182 tthn(ix,jy,kz,n,l)*dz1)/dz 183 qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ & 184 qvhn(ix,jy,kz,n,l)*dz1)/dz 185 pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ & 186 pvhn(ix,jy,kz,l)*dz1)/dz 187 rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz 188 kmin=kz 189 goto 30 214 enddo 215 enddo 216 do jy=0,nyn(l)-1 217 do ix=0,nxn(l)-1 218 if(height(iz).le.uvzlev(ix,jy,nuvz)) then 219 kz=idx(ix,jy) 220 dz1=height(iz)-uvzlev(ix,jy,kz-1) 221 dz2=uvzlev(ix,jy,kz)-height(iz) 222 dz=dz1+dz2 223 uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+uuhn(ix,jy,kz,l)*dz1)/dz 224 vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+vvhn(ix,jy,kz,l)*dz1)/dz 225 ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2 & 226 +tthn(ix,jy,kz,n,l)*dz1)/dz 227 qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2 & 228 +qvhn(ix,jy,kz,n,l)*dz1)/dz 229 !hg adding the cloud water 230 if (readclouds_nest(l)) then 231 clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2+clwchn(ix,jy,kz,n,l)*dz1)/dz 232 if (.not.sumclouds_nest(l)) & 233 &ciwcn(ix,jy,iz,n,l)=(ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz 234 end if 235 !hg 236 pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+pvhn(ix,jy,kz,l)*dz1)/dz 237 rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)*dz2+rhohn(ix,jy,kz)*dz1)/dz 190 238 endif 239 enddo 240 enddo 241 enddo 242 243 244 245 ! do kz=kmin,nuvz 246 ! if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 247 ! uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l) 248 ! vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l) 249 ! ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l) 250 ! qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l) 251 ! pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l) 252 ! rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l) 253 ! goto 30 254 ! endif 255 ! if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 256 ! (height(iz).le.uvwzlev(ix,jy,kz))) then 257 ! dz1=height(iz)-uvwzlev(ix,jy,kz-1) 258 ! dz2=uvwzlev(ix,jy,kz)-height(iz) 259 ! dz=dz1+dz2 260 ! uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ & 261 ! uuhn(ix,jy,kz,l)*dz1)/dz 262 ! vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ & 263 ! vvhn(ix,jy,kz,l)*dz1)/dz 264 ! ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ & 265 ! tthn(ix,jy,kz,n,l)*dz1)/dz 266 ! qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ & 267 ! qvhn(ix,jy,kz,n,l)*dz1)/dz 268 ! pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ & 269 ! pvhn(ix,jy,kz,l)*dz1)/dz 270 ! rhon(ix,jy,iz,n,l)=(rhohn(kz-1)*dz2+rhohn(kz)*dz1)/dz 271 ! kmin=kz 272 ! goto 30 273 ! endif 274 ! end do 275 ! 30 continue 276 ! end do 277 278 279 ! Levels, where w is given 280 !************************* 281 282 wwn(0:nxm1,0:nym1,1,n,l)=wwhn(0:nxm1,0:nym1,1,l)*pinmconv(0:nxm1,0:nym1,1) 283 wwn(0:nxm1,0:nym1,nz,n,l)=wwhn(0:nxm1,0:nym1,nwz,l)*pinmconv(0:nxm1,0:nym1,nz) 284 kmin=2 285 idx=kmin 286 do iz=2,nz 287 do jy=0,nym1 288 do ix=0,nxm1 289 inn: do kz=idx(ix,jy),nwz 290 if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. & 291 height(iz).le.wzlev(ix,jy,kz)) then 292 idx(ix,jy)=kz 293 exit inn 294 endif 295 enddo inn 296 enddo 297 enddo 298 do jy=0,nym1 299 do ix=0,nxm1 300 kz=idx(ix,jy) 301 dz1=height(iz)-wzlev(ix,jy,kz-1) 302 dz2=wzlev(ix,jy,kz)-height(iz) 303 dz=dz1+dz2 304 wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(ix,jy,kz-1)*dz2 & 305 +wwhn(ix,jy,kz,l)*pinmconv(ix,jy,kz)*dz1)/dz 306 enddo 307 enddo 308 enddo 309 310 ! wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1) 311 ! wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz) 312 ! kmin=2 313 ! do iz=2,nz 314 ! do kz=kmin,nwz 315 ! if ((height(iz).gt.wzlev(kz-1)).and. & 316 ! (height(iz).le.wzlev(kz))) then 317 ! dz1=height(iz)-wzlev(kz-1) 318 ! dz2=wzlev(kz)-height(iz) 319 ! dz=dz1+dz2 320 ! wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 & 321 ! +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz 322 ! kmin=kz 323 ! goto 40 324 ! endif 325 ! end do 326 ! 40 continue 327 ! end do 328 329 ! Compute density gradients at intermediate levels 330 !************************************************* 331 332 drhodzn(0:nxm1,0:nym1,1,n,l)=(rhon(0:nxm1,0:nym1,2,n,l)-rhon(0:nxm1,0:nym1,1,n,l))/ & 333 (height(2)-height(1)) 334 do kz=2,nz-1 335 drhodzn(0:nxm1,0:nym1,kz,n,l)=(rhon(0:nxm1,0:nym1,kz+1,n,l)-rhon(0:nxm1,0:nym1,kz-1,n,l))/ & 336 (height(kz+1)-height(kz-1)) 337 end do 338 drhodzn(0:nxm1,0:nym1,nz,n,l)=drhodzn(0:nxm1,0:nym1,nz-1,n,l) 339 340 341 ! drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ & 342 ! (height(2)-height(1)) 343 ! do kz=2,nz-1 344 ! drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- & 345 ! rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1)) 346 ! end do 347 ! drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l) 348 349 350 351 ! end do 352 ! end do 353 354 355 !**************************************************************** 356 ! Compute slope of eta levels in windward direction and resulting 357 ! vertical wind correction 358 !**************************************************************** 359 360 do jy=1,nyn(l)-2 361 ! cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180) 362 cosf(jy)=1./cos((real(jy)*dyn(l)+ylat0n(l))*pi180) 363 end do 364 365 kmin=2 366 idx=kmin 367 do iz=2,nz-1 368 do jy=1,nyn(l)-2 369 do ix=1,nxn(l)-2 370 371 inneta: do kz=idx(ix,jy),nz 372 if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. & 373 (height(iz).le.uvzlev(ix,jy,kz))) then 374 idx(ix,jy)=kz 375 exit inneta 376 endif 377 enddo inneta 378 enddo 379 enddo 380 381 do jy=1,nyn(l)-2 382 do ix=1,nxn(l)-2 383 kz=idx(ix,jy) 384 dz1=height(iz)-uvzlev(ix,jy,kz-1) 385 dz2=uvzlev(ix,jy,kz)-height(iz) 386 dz=dz1+dz2 387 ix1=ix-1 388 jy1=jy-1 389 ixp=ix+1 390 jyp=jy+1 391 392 dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2. 393 dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2. 394 dzdx=(dzdx1*dz2+dzdx2*dz1)/dz 395 396 dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2. 397 dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2. 398 dzdy=(dzdy1*dz2+dzdy2*dz1)/dz 399 400 wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy)+& 401 &dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)) 402 191 403 end do 192 30 continue 404 193 405 end do 194 195 196 ! Levels, where w is given 197 !************************* 198 199 wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1) 200 wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz) 201 kmin=2 202 do iz=2,nz 203 do kz=kmin,nwz 204 if ((height(iz).gt.wzlev(kz-1)).and. & 205 (height(iz).le.wzlev(kz))) then 206 dz1=height(iz)-wzlev(kz-1) 207 dz2=wzlev(kz)-height(iz) 208 dz=dz1+dz2 209 wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 & 210 +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz 211 kmin=kz 212 goto 40 213 endif 406 end do 407 408 409 ! do jy=1,nyn(l)-2 410 ! do ix=1,nxn(l)-2 411 ! kmin=2 412 ! do iz=2,nz-1 413 414 ! ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf(jy) 415 ! vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l) 416 417 ! do kz=kmin,nz 418 ! if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 419 ! (height(iz).le.uvwzlev(ix,jy,kz))) then 420 ! dz1=height(iz)-uvwzlev(ix,jy,kz-1) 421 ! dz2=uvwzlev(ix,jy,kz)-height(iz) 422 ! dz=dz1+dz2 423 ! kl=kz-1 424 ! klp=kz 425 ! kmin=kz 426 ! goto 47 427 ! endif 428 ! end do 429 430 ! 47 ix1=ix-1 431 ! jy1=jy-1 432 ! ixp=ix+1 433 ! jyp=jy+1 434 435 ! dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2. 436 ! dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2. 437 ! dzdx=(dzdx1*dz2+dzdx2*dz1)/dz 438 439 ! dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2. 440 ! dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2. 441 ! dzdy=(dzdy1*dz2+dzdy2*dz1)/dz 442 443 ! wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi) 444 445 ! end do 446 447 ! end do 448 ! end do 449 450 451 !write (*,*) 'initializing nested cloudsn, n:',n 452 ! create a cloud and rainout/washout field, cloudsn occur where rh>80% 453 454 455 !*********************************************************************************** 456 if (readclouds_nest(l)) then !HG METHOD 457 ! The method is loops all grids vertically and constructs the 3D matrix for clouds 458 ! Cloud top and cloud bottom gid cells are assigned as well as the total column 459 ! cloud water. For precipitating grids, the type and whether it is in or below 460 ! cloud scavenging are assigned with numbers 2-5 (following the old metod). 461 ! Distinction is done for lsp and convp though they are treated the same in regards 462 ! to scavenging. Also clouds that are not precipitating are defined which may be 463 ! to include future cloud processing by non-precipitating-clouds. 464 !*********************************************************************************** 465 write(*,*) 'Nested ECMWF fields: using cloud water' 466 clwn(0:nxm1,0:nym1,:,n,l)=0.0 467 ! icloud_stats(0:nxm1,0:nym1,:,n)=0.0 468 clw4n(0:nxm1,0:nym1,n,l)=0.0 469 cloudsn(0:nxm1,0:nym1,:,n,l)=0 470 ! If water/ice are read separately into clwc and ciwc, store sum in clwcn 471 if (.not.sumclouds_nest(l)) then 472 clwcn(0:nxm1,0:nym1,:,n,l) = clwcn(0:nxm1,0:nym1,:,n,l) + ciwcn(0:nxm1,0:nym1,:,n,l) 473 end if 474 do jy=0,nym1 475 do ix=0,nxm1 476 lsp=lsprecn(ix,jy,1,n,l) 477 convp=convprecn(ix,jy,1,n,l) 478 prec=lsp+convp 479 tot_cloud_h=0 480 ! Find clouds in the vertical 481 do kz=1, nz-1 !go from top to bottom 482 if (clwcn(ix,jy,kz,n,l).gt.0) then 483 ! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3 484 clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l))*(height(kz+1)-height(kz)) 485 tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 486 clw4n(ix,jy,n,l) = clw4n(ix,jy,n,l)+clwn(ix,jy,kz,n,l) 487 ! icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n) ! Column cloud water [m3/m3] 488 ! icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz)) ! Cloud BOT height stats [m] 489 cloudh_min=min(height(kz+1),height(kz)) 490 !ZHG 2015 extra for testing 491 ! clh(ix,jy,kz,n)=height(kz+1)-height(kz) 492 ! icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m] 493 ! icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz)) ! Cloud TOP height [m] 494 !ZHG 495 endif 496 end do 497 498 ! If Precipitation. Define removal type in the vertical 499 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 500 501 do kz=nz,1,-1 !go Bottom up! 502 if (clwn(ix,jy,kz,n,l).gt.0.0) then ! is in cloud 503 cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1) 504 cloudsn(ix,jy,kz,n,l)=1 ! is a cloud 505 if (lsp.ge.convp) then 506 cloudsn(ix,jy,kz,n,l)=3 ! lsp in-cloud 507 else 508 cloudsn(ix,jy,kz,n,l)=2 ! convp in-cloud 509 endif ! convective or large scale 510 else if((clwn(ix,jy,kz,n,l).le.0.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud 511 if (lsp.ge.convp) then 512 cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout 513 else 514 cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout 515 endif ! convective or large scale 516 endif 517 518 if (height(kz).ge. 19000) then ! set a max height for removal 519 cloudsn(ix,jy,kz,n,l)=0 520 endif !clw>0 521 end do !nz 522 endif ! precipitation 214 523 end do 215 40 continue216 524 end do 217 218 ! Compute density gradients at intermediate levels 219 !************************************************* 220 221 drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ & 222 (height(2)-height(1)) 223 do kz=2,nz-1 224 drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- & 225 rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1)) 525 !******************************************************************** 526 else ! old method: 527 !******************************************************************** 528 write(*,*) 'Nested fields: using cloud water from Parameterization' 529 do jy=0,nyn(l)-1 530 do ix=0,nxn(l)-1 531 rain_cloud_above(ix,jy)=0 532 lsp=lsprecn(ix,jy,1,n,l) 533 convp=convprecn(ix,jy,1,n,l) 534 cloudshn(ix,jy,n,l)=0 535 do kz_inv=1,nz-1 536 kz=nz-kz_inv+1 537 pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l) 538 rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l)) 539 cloudsn(ix,jy,kz,n,l)=0 540 if (rh.gt.0.8) then ! in cloud 541 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then 542 rain_cloud_above(ix,jy)=1 543 cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ & 544 height(kz)-height(kz-1) 545 if (lsp.ge.convp) then 546 cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout 547 else 548 cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout 549 endif 550 else ! no precipitation 551 cloudsn(ix,jy,kz,n,l)=1 ! cloud 552 endif 553 else ! no cloud 554 if (rain_cloud_above(ix,jy).eq.1) then ! scavenging 555 if (lsp.ge.convp) then 556 cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout 557 else 558 cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout 559 endif 560 endif 561 endif 562 end do 563 end do 226 564 end do 227 drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l) 228 229 end do 230 end do 231 232 233 !**************************************************************** 234 ! Compute slope of eta levels in windward direction and resulting 235 ! vertical wind correction 236 !**************************************************************** 237 238 do jy=1,nyn(l)-2 239 cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180) 240 do ix=1,nxn(l)-2 241 242 kmin=2 243 do iz=2,nz-1 244 245 ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf 246 vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l) 247 248 do kz=kmin,nz 249 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 250 (height(iz).le.uvwzlev(ix,jy,kz))) then 251 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 252 dz2=uvwzlev(ix,jy,kz)-height(iz) 253 dz=dz1+dz2 254 kl=kz-1 255 klp=kz 256 kmin=kz 257 goto 47 258 endif 259 end do 260 261 47 ix1=ix-1 262 jy1=jy-1 263 ixp=ix+1 264 jyp=jy+1 265 266 dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2. 267 dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2. 268 dzdx=(dzdx1*dz2+dzdx2*dz1)/dz 269 270 dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2. 271 dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2. 272 dzdy=(dzdy1*dz2+dzdy2*dz1)/dz 273 274 wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi) 275 276 end do 277 278 end do 279 end do 280 281 282 !write (*,*) 'initializing nested cloudsn, n:',n 283 ! create a cloud and rainout/washout field, cloudsn occur where rh>80% 284 write(*,*) 'Nested fields: using cloud water from Parameterization' 285 do jy=0,nyn(l)-1 286 do ix=0,nxn(l)-1 287 rain_cloud_above=0 288 lsp=lsprecn(ix,jy,1,n,l) 289 convp=convprecn(ix,jy,1,n,l) 290 cloudsnh(ix,jy,n,l)=0 291 do kz_inv=1,nz-1 292 kz=nz-kz_inv+1 293 pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l) 294 rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l)) 295 cloudsn(ix,jy,kz,n,l)=0 296 if (rh.gt.0.8) then ! in cloud 297 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then 298 rain_cloud_above=1 299 cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ & 300 height(kz)-height(kz-1) 301 if (lsp.ge.convp) then 302 cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout 303 else 304 cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout 305 endif 306 else ! no precipitation 307 cloudsn(ix,jy,kz,n,l)=1 ! cloud 308 endif 309 else ! no cloud 310 if (rain_cloud_above.eq.1) then ! scavenging 311 if (lsp.ge.convp) then 312 cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout 313 else 314 cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout 315 endif 316 endif 317 endif 318 end do 319 end do 320 end do 321 322 end do 565 end if 566 567 end do ! end loop over nests 323 568 324 569 end subroutine verttransform_nests -
src/wetdepo.f90
r41d8574 rdb712a8 21 21 22 22 subroutine wetdepo(itime,ltsample,loutnext) 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 23 ! i i i 24 !***************************************************************************** 25 ! * 26 ! Calculation of wet deposition using the concept of scavenging coefficients.* 27 ! For lack of detailed information, washout and rainout are jointly treated. * 28 ! It is assumed that precipitation does not occur uniformly within the whole * 29 ! grid cell, but that only a fraction of the grid cell experiences rainfall. * 30 ! This fraction is parameterized from total cloud cover and rates of large * 31 ! scale and convective precipitation. * 32 ! * 33 ! Author: A. Stohl * 34 ! * 35 ! 1 December 1996 * 36 ! * 37 ! Correction by Petra Seibert, Sept 2002: * 38 ! use centred precipitation data for integration * 39 ! Code may not be correct for decay of deposition! * 40 ! * 41 !***************************************************************************** 42 ! * 43 ! Variables: * 44 ! cc [0-1] total cloud cover * 45 ! convp [mm/h] convective precipitation rate * 46 ! grfraction [0-1] fraction of grid, for which precipitation occurs * 47 ! ix,jy indices of output grid cell for each particle * 48 ! itime [s] actual simulation time [s] * 49 ! jpart particle index * 50 ! ldeltat [s] interval since radioactive decay was computed * 51 ! lfr, cfr area fraction covered by precipitation for large scale * 52 ! and convective precipitation (dependent on prec. rate) * 53 ! loutnext [s] time for which gridded deposition is next output * 54 ! loutstep [s] interval at which gridded deposition is output * 55 ! lsp [mm/h] large scale precipitation rate * 56 ! ltsample [s] interval over which mass is deposited * 57 ! prec [mm/h] precipitation rate in subgrid, where precipitation occurs* 58 ! wetdeposit mass that is wet deposited * 59 ! wetgrid accumulated deposited mass on output grid * 60 ! wetscav scavenging coefficient * 61 ! * 62 ! Constants: * 63 ! * 64 !***************************************************************************** 65 65 66 66 use point_mod … … 71 71 72 72 integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy 73 integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v 73 integer :: ngrid,itage,nage,hz,il,interp_time, n 74 integer(kind=1) :: clouds_v 74 75 integer :: ks, kp 75 76 ! integer :: n1,n2, icbot,ictop, indcloud !TEST … … 79 80 real :: wetdeposit(maxspec),restmass 80 81 real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled 81 82 !save lfr,cfr 82 83 83 84 real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/) … … 91 92 integer :: blc_count, inc_count 92 93 real :: Si_dummy, wetscav_dummy 93 94 95 96 97 94 logical :: readclouds_this_nest 95 96 97 ! Compute interval since radioactive decay of deposited mass was computed 98 !************************************************************************ 98 99 99 100 if (itime.le.loutnext) then … … 103 104 endif 104 105 105 106 106 ! Loop over all particles 107 !************************ 107 108 108 109 blc_count=0 … … 118 119 endif 119 120 120 121 ! Determine age class of the particle 121 122 itage=abs(itra1(jpart)-itramem(jpart)) 122 123 do nage=1,nageclass … … 126 127 127 128 128 129 129 ! Determine which nesting level to be used 130 !***************************************** 130 131 131 132 ngrid=0 132 133 do j=numbnests,1,-1 133 134 if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. & 134 (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then135 (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then 135 136 ngrid=j 136 137 goto 23 … … 140 141 141 142 142 143 143 ! Determine nested grid coordinates 144 !********************************** 144 145 145 146 if (ngrid.gt.0) then … … 148 149 ix=int(xtn) 149 150 jy=int(ytn) 151 if (readclouds_nest(ngrid)) then 152 readclouds_this_nest=.true. 153 else 154 readclouds_this_nest=.false. 155 end if 150 156 else 151 157 ix=int(xtra1(jpart)) … … 154 160 155 161 156 157 158 159 162 ! Interpolate large scale precipitation, convective precipitation and 163 ! total cloud cover 164 ! Note that interpolated time refers to itime-0.5*ltsample [PS] 165 !******************************************************************** 160 166 interp_time=nint(itime-0.5*ltsample) 161 167 162 168 if (ngrid.eq.0) then 163 169 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, & 164 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &165 memtime(1),memtime(2),interp_time,lsp,convp,cc)170 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, & 171 memtime(1),memtime(2),interp_time,lsp,convp,cc) 166 172 else 167 173 call interpol_rain_nests(lsprecn,convprecn,tccn, & 168 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &169 memtime(1),memtime(2),interp_time,lsp,convp,cc)174 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 175 memtime(1),memtime(2),interp_time,lsp,convp,cc) 170 176 endif 171 177 … … 173 179 if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 174 180 175 181 ! get the level were the actual particle is in 176 182 do il=2,nz 177 183 if (height(il).gt.ztra1(jpart)) then 178 184 hz=il-1 179 goto 26 185 ! goto 26 186 exit 180 187 endif 181 188 end do 182 26 continue189 !26 continue 183 190 184 191 n=memind(2) … … 190 197 clouds_h=cloudsh(ix,jy,n) 191 198 else 192 ! new removal not implemented for nests yet193 199 clouds_v=cloudsn(ix,jy,hz,n,ngrid) 194 clouds_h=clouds nh(ix,jy,n,ngrid)195 endif 196 197 198 200 clouds_h=cloudshn(ix,jy,n,ngrid) 201 endif 202 203 ! if there is no precipitation or the particle is above the clouds no 204 ! scavenging is done 199 205 200 206 if (clouds_v.le.1) goto 20 201 202 203 204 205 206 207 207 208 ! 1) Parameterization of the the area fraction of the grid cell where the 209 ! precipitation occurs: the absolute limit is the total cloud cover, but 210 ! for low precipitation rates, an even smaller fraction of the grid cell 211 ! is used. Large scale precipitation occurs over larger areas than 212 ! convective precipitation. 213 !************************************************************************** 208 214 209 215 if (lsp.gt.20.) then … … 232 238 233 239 234 235 236 240 !ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp 241 ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms 242 ! for now they are treated the same 237 243 grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp)) 238 244 grfraction(2)=max(0.05,cc*(lfr(i))) … … 240 246 241 247 242 243 248 ! 2) Computation of precipitation rate in sub-grid cell 249 !****************************************************** 244 250 prec(1)=(lsp+convp)/grfraction(1) 245 251 prec(2)=(lsp)/grfraction(2) … … 247 253 248 254 249 250 251 255 ! 3) Computation of scavenging coefficients for all species 256 ! Computation of wet deposition 257 !********************************************************** 252 258 253 259 do ks=1,nspec ! loop over species … … 255 261 wetscav=0. 256 262 257 263 !ZHG test if it nested? 258 264 if (ngrid.gt.0) then 259 265 act_temp=ttn(ix,jy,hz,n,ngrid) … … 261 267 act_temp=tt(ix,jy,hz,n) 262 268 endif 263 264 265 266 267 269 270 271 !*********************** 272 ! BELOW CLOUD SCAVENGING 273 !*********************** 268 274 if (clouds_v.ge.4) then !below cloud 269 275 … … 274 280 if (dquer(ks) .le. 0.) then !is gas 275 281 wetscav=weta(ks)*prec(1)**wetb(ks) 276 282 277 283 !AEROSOL 278 284 else !is particle … … 281 287 ! for particles larger than 10 um use the largest size defined in the parameterizations (10um) 282 288 dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m 283 if (act_temp .ge. 273 .and. weta(ks).gt.0.) then !Rain284 285 286 289 if (act_temp .ge. 273. .and. weta(ks).gt.0.) then !Rain 290 ! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003, 291 ! the below-cloud scavenging (rain efficienty) 292 ! parameter A (=weta) from SPECIES file 287 293 wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* & 288 294 (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5)) 289 295 290 elseif (act_temp .lt. 273 .and. wetb(ks).gt.0.) then ! Snow291 292 293 296 elseif (act_temp .lt. 273. .and. wetb(ks).gt.0.) then ! Snow 297 ! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009, 298 ! the below-cloud scavenging (Snow efficiency) 299 ! parameter B (=wetb) from SPECIES file 294 300 wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* & 295 301 (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5)) … … 304 310 305 311 306 307 308 !******************************************************312 !******************** 313 ! IN CLOUD SCAVENGING 314 !******************** 309 315 if (clouds_v.lt.4) then ! In-cloud 310 316 ! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file … … 315 321 if (wetb_in(ks).lt.0.) wetb_in(ks)=0. 316 322 317 !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF 318 if (readclouds) then !icloud_stats(ix,jy,4,n) has units kg/m2 319 ! cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc) 320 ! ESO: stop using icloud_stats 321 cl =clw4(ix,jy,n)*(grfraction(1)/cc) 323 !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF 324 ! nested fields 325 if (ngrid.gt.0.and.readclouds_this_nest) then 326 cl=clw4n(ix,jy,n,ngrid)*(grfraction(1)/cc) 327 else if (ngrid.eq.0.and.readclouds) then 328 cl=clw4(ix,jy,n)*(grfraction(1)/cc) 322 329 else !parameterize cloudwater m2/m3 323 330 !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF 324 331 cl=1.6E-6*prec(1)**0.36 325 332 endif 326 333 327 328 if (act_temp .le. 253) then329 330 else if (act_temp .ge. 273) then331 332 333 liq_frac =((act_temp-273)/(273-253))**2334 endif335 336 337 338 339 340 341 334 !ZHG: Calculate the partition between liquid and water phase water. 335 if (act_temp .le. 253.) then 336 liq_frac=0 337 else if (act_temp .ge. 273.) then 338 liq_frac=1 339 else 340 liq_frac =((act_temp-273.)/(273.-253.))**2. 341 end if 342 ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi 343 frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks) 344 345 !ZHG Use the activated fraction and the liqid water to calculate the washout 346 347 ! AEROSOL 348 !************************************************** 342 349 if (dquer(ks).gt. 0.) then ! is particle 343 350 344 351 S_i= frac_act/cl 345 352 346 347 353 !********************* 354 ! GAS 348 355 else ! is gas 349 356 350 357 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 351 352 358 !REPLACE to switch old/ new scheme 359 ! S_i=frac_act/cle 353 360 S_i=1/cle 354 361 endif ! gas or particle 355 362 356 357 358 if ( readclouds) then363 ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol 364 !OLD 365 if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then 359 366 wetscav=S_i*(prec(1)/3.6E6) 360 367 else … … 362 369 endif 363 370 364 !ZHG 2015 TEST365 ! Si_dummy=frac_act/2E-7*prec(1)**0.36366 ! wetscav_dummy=Si_dummy*(prec(1)/3.6E6)/clouds_h367 ! if (clouds_v.lt.4) then368 ! talltest=talltest+1369 !if(talltest .eq. 1) OPEN(UNIT=199, FILE=utfil,FORM='FORMATTED',STATUS = 'UNKNOWN')370 !if(talltest .lt. 100001) write(199,*) prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90371 !if(talltest .lt. 100001) write(199,*) wetscav, wetscav_dummy,prec(1),ytra1(jpart)-90,clouds_v,cl372 !if(talltest .eq. 100001) CLOSE(199)373 !if(talltest .eq. 100001) STOP374 !375 !write(*,*) 'PREC kg/m2s CLOUD kg/m2', (prec(1)/3.6E6), cl !, '2E-7*prec(1)**0.36', 2E-7*prec(1)**0.36,'2E-7*prec(1)**0.36*clouds_h',2E-7*prec(1)**0.36*clouds_h376 !write(*,*) 'PREC kg/m2s LSP+convp kg/m2', prec(1), convp+lsp377 !write(*,*) wetscav, wetscav_dummy378 !write(*,*) cc, grfraction(1), cc/grfraction(1)379 380 !write(*,*) 'Lmbda_old', (prec(1)/3.6E6)/(clouds_h*2E-7*prec(1)**0.36)381 382 383 !write(*,*) '**************************************************'384 !write(*,*) 'clouds_h', clouds_h, 'clouds_v',clouds_v,'abs(ltsample)', abs(ltsample)385 !write(*,*) 'readclouds', readclouds, 'wetscav',wetscav, 'wetscav_dummy', wetscav_dummy386 !write(*,*) 'S_i', S_i , 'Si_dummy', Si_dummy, 'prec(1)', prec(1)387 388 389 ! write(*,*) 'PRECIPITATION ,cl ECMWF , cl PARAMETIZED, clouds_v, lat' &390 ! ,prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90391 392 !endif393 371 394 372 endif ! positive in-cloud scavenging parameters given in Species file 395 373 endif !incloud 396 374 !END ZHG TEST 397 398 399 400 401 402 375 376 !************************************************** 377 ! CALCULATE DEPOSITION 378 !************************************************** 379 ! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!' 380 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v 403 381 404 382 if (wetscav.gt.0.) then … … 421 399 if (restmass .gt. smallnum) then 422 400 xmass1(jpart,ks)=restmass 423 424 425 401 ! depostatistic 402 ! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 403 ! depostatistic 426 404 else 427 405 xmass1(jpart,ks)=0. 428 406 endif 429 430 407 ! Correct deposited mass to the last time step when radioactive decay of 408 ! gridded deposited mass was calculated 431 409 if (decay(ks).gt.0.) then 432 410 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks)) … … 436 414 end do !all species 437 415 438 439 440 416 ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs 417 ! Add the wet deposition to accumulated amount on output grid and nested output grid 418 !***************************************************************************** 441 419 442 420 if (ldirect.eq.1) then … … 450 428 end do ! all particles 451 429 452 430 ! count the total number of below-cloud and in-cloud occurences: 453 431 tot_blc_count=tot_blc_count+blc_count 454 432 tot_inc_count=tot_inc_count+inc_count -
src/wetdepokernel_nest.f90
re200b7a rdb712a8 20 20 !********************************************************************** 21 21 22 subroutine wetdepokernel_nest & 23 (nunc,deposit,x,y,nage,kp) 24 ! i i i i i i 22 subroutine wetdepokernel_nest(nunc,deposit,x,y,nage,kp) 23 ! i i i i i i 25 24 !***************************************************************************** 26 25 ! * … … 55 54 integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage 56 55 57 56 real :: dbg_dx, dbg_dy, dbg_xoutshiftn, dbg_youtshiftn, dbg_dxoutn, dbg_dyoutn,dbg_t 58 57 59 58 xl=(x*dx+xoutshiftn)/dxoutn 60 59 yl=(y*dy+youtshiftn)/dyoutn 61 ix=int(xl) 62 jy=int(yl) 60 61 ! old: 62 ! ix=int(xl) 63 ! jy=int(yl) 64 65 ! ESO: for xl,yl in range <-.5,-1> we get ix,jy=0 and thus negative 66 ! wx,wy as the int function rounds upwards for negative numbers. 67 ! Either use the floor function, or (perhaps more correctly?) use "(xl.gt.-0.5)" 68 ! in place of "(ix.ge.0)" and similar for the upper boundary. 69 70 ! new: 71 ix=floor(xl) 72 jy=floor(yl) 73 63 74 ddx=xl-real(ix) ! distance to left cell border 64 75 ddy=yl-real(jy) ! distance to lower cell border 76 65 77 66 78 if (ddx.gt.0.5) then … … 80 92 endif 81 93 82 83 ! Determine mass fractions for four grid points 84 !********************************************** 94 ! Determine mass fractions for four grid points 95 !********************************************** 85 96 86 97 do ks=1,nspec 87 88 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & 89 (jy.le.numygridn-1)) then 90 w=wx*wy 98 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & 99 (jy.le.numygridn-1)) then 100 w=wx*wy 91 101 wetgriduncn(ix,jy,ks,kp,nunc,nage)= & 92 102 wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w 93 endif103 endif 94 104 95 if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. &96 (jyp.le.numygridn-1)) then97 w=(1.-wx)*(1.-wy)105 if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. & 106 (jyp.le.numygridn-1)) then 107 w=(1.-wx)*(1.-wy) 98 108 wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= & 99 109 wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w 100 endif110 endif 101 111 102 if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. &103 (jy.le.numygridn-1)) then104 w=(1.-wx)*wy112 if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. & 113 (jy.le.numygridn-1)) then 114 w=(1.-wx)*wy 105 115 wetgriduncn(ixp,jy,ks,kp,nunc,nage)= & 106 116 wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w 107 endif117 endif 108 118 109 if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. &110 (jyp.le.numygridn-1)) then111 w=wx*(1.-wy)119 if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. & 120 (jyp.le.numygridn-1)) then 121 w=wx*(1.-wy) 112 122 wetgriduncn(ix,jyp,ks,kp,nunc,nage)= & 113 123 wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w 114 endif124 endif 115 125 116 126 end do
Note: See TracChangeset
for help on using the changeset viewer.