Changeset 5081 for LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2
- Timestamp:
- Jul 19, 2024, 4:15:44 PM (4 months ago)
- Location:
- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2
- Files:
-
- 6 edited
-
cosp_optics.F90 (modified) (6 diffs)
-
mrgrnk.F90 (modified) (3 diffs)
-
optics_lib.F90 (modified) (7 diffs)
-
quickbeam.F90 (modified) (14 diffs)
-
quickbeam_optics.F90 (modified) (7 diffs)
-
scops.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90
r3491 r5081 72 72 varOUT(1:dim1,1:dim2,1:dim3) = 0._wp 73 73 do j=1,dim2 74 where(flag(:,j,:) .eq.1)74 where(flag(:,j,:) == 1) 75 75 varOUT(:,j,:) = varIN2 76 76 endwhere 77 where(flag(:,j,:) .eq.2)77 where(flag(:,j,:) == 2) 78 78 varOUT(:,j,:) = varIN1 79 79 endwhere … … 96 96 97 97 varOUT(1:dim1,1:dim2,1:dim3) = 0._wp 98 where(flag(:,:,:) .eq.1)98 where(flag(:,:,:) == 1) 99 99 varOUT(:,:,:) = varIN2 100 100 endwhere 101 where(flag(:,:,:) .eq.2)101 where(flag(:,:,:) == 2) 102 102 varOUT(:,:,:) = varIN1 103 103 endwhere … … 295 295 296 296 ! Which LIDAR frequency are we using? 297 if (lidar_freq .eq.355) then297 if (lidar_freq == 355) then 298 298 Cmol = Cmol_355nm 299 299 rdiffm = rdiffm_355nm 300 300 endif 301 if (lidar_freq .eq.532) then301 if (lidar_freq == 532) then 302 302 Cmol = Cmol_532nm 303 303 rdiffm = rdiffm_532nm … … 336 336 337 337 ! LS and CONV Ice water coefficients 338 if (ice_type .eq.0) then338 if (ice_type == 0) then 339 339 polpart(INDX_LSICE,1:5) = polpartLSICE0 340 340 polpart(INDX_CVICE,1:5) = polpartCVICE0 341 341 endif 342 if (ice_type .eq.1) then342 if (ice_type == 1) then 343 343 polpart(INDX_LSICE,1:5) = polpartLSICE1 344 344 polpart(INDX_CVICE,1:5) = polpartCVICE1 … … 393 393 ! Polynomials kp_lidar derived from Mie theory 394 394 do i = 1, npart 395 where (rad_part(1:npoints,1:nlev,i) .gt.0.0)395 where (rad_part(1:npoints,1:nlev,i) > 0.0) 396 396 kp_part(1:npoints,1:nlev,i) = & 397 397 polpart(i,1)*(rad_part(1:npoints,1:nlev,i)*1e6)**4 & … … 426 426 ! Alpha of particles in each subcolumn: 427 427 do i = 1, npart 428 where (rad_part(1:npoints,1:nlev,i) .gt.0.0)428 where (rad_part(1:npoints,1:nlev,i) > 0.0) 429 429 alpha_part(1:npoints,1:nlev,i) = 3._wp/4._wp * Qscat & 430 430 * rhoair(1:npoints,1:nlev) * qpart(1:npoints,1:nlev,i) & -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mrgrnk.F90
r3491 r5081 68 68 IRNGT (1) = 1 69 69 Return 70 Case Default71 Continue72 70 End Select 73 71 ! … … 268 266 IRNGT (1) = 1 269 267 Return 270 Case Default271 Continue272 268 End Select 273 269 ! … … 467 463 IRNGT (1) = 1 468 464 Return 469 Case Default470 Continue471 465 End Select 472 466 ! -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90
r3491 r5081 558 558 if(tk < temref(4)) tk=temref(4) 559 559 do i=2,4 560 if(tk .ge.temref(i)) go to 12560 if(tk>=temref(i)) go to 12 561 561 enddo 562 562 12 lt1 = i 563 563 lt2 = i-1 564 564 do i=2,nwlt 565 if(alam .le.wlt(i)) go to 14565 if(alam<=wlt(i)) go to 14 566 566 enddo 567 567 14 x1 = log(wlt(i-1)) … … 652 652 Complex(wp) :: A1 653 653 654 If ((Dx .Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then654 If ((Dx>Imaxx) .Or. (InP>ImaxNP)) Then 655 655 Error = 1 656 656 Return … … 659 659 Ir = 1 / Cm 660 660 Y = Dx * Cm 661 If (Dx .Lt.0.02) Then661 If (Dx<0.02) Then 662 662 NStop = 2 663 663 Else 664 If (Dx .Le.8.0) Then664 If (Dx<=8.0) Then 665 665 NStop = Dx + 4.00*Dx**(1./3.) + 2.0 666 666 Else 667 If (Dx .Lt.4200.0) Then667 If (Dx< 4200.0) Then 668 668 NStop = Dx + 4.05*Dx**(1./3.) + 2.0 669 669 Else … … 673 673 End If 674 674 NmX = Max(Real(NStop),Real(Abs(Y))) + 15. 675 If (Nmx .gt.Itermax) then675 If (Nmx > Itermax) then 676 676 Error = 1 677 677 Return … … 726 726 !ds Dqxt = Tnp1 * Dble(A + B) + Dqxt 727 727 Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc 728 If (N .Gt.1) then728 If (N>1) then 729 729 Dg = Dg + (dN*dN - 1) * (ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 *(ANM1*Conjg(BNM1)) / (dN*dN - dN) 730 730 !ds Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN) … … 735 735 AMB = A2 * (A - B) 736 736 Do I = 1,Inp2 737 If (I .GT.Inp) Then737 If (I>Inp) Then 738 738 S(I) = -Pi1(I) 739 739 Else … … 756 756 End Do 757 757 758 If (Dg .GT.0) Dg = 2 * Dg / Dqsc758 If (Dg >0) Dg = 2 * Dg / Dqsc 759 759 Dqsc = 2 * Dqsc / Dx**2 760 760 Dqxt = 2 * Dqxt / Dx**2 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90
r3491 r5081 179 179 180 180 ! Attenuation due to gaseous absorption between radar and volume 181 if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr .eq.1)) then181 if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr == 1)) then 182 182 if (d_gate==1) then 183 183 if (k>1) then … … 272 272 273 273 ! Which platforms to create diagnostics for? 274 if (platform .eq.'cloudsat') lcloudsat=.true.274 if (platform == 'cloudsat') lcloudsat=.true. 275 275 276 276 ! Create Cloudsat diagnostics. … … 289 289 enddo 290 290 enddo 291 where(cfad_ze .ne.R_UNDEF) cfad_ze = cfad_ze/Ncolumns291 where(cfad_ze /= R_UNDEF) cfad_ze = cfad_ze/Ncolumns 292 292 293 293 ! Compute cloudsat near-surface precipitation diagnostics … … 306 306 enddo 307 307 enddo 308 where(cfad_ze .ne.R_UNDEF) cfad_ze = cfad_ze/Ncolumns308 where(cfad_ze /= R_UNDEF) cfad_ze = cfad_ze/Ncolumns 309 309 endif 310 310 endif … … 402 402 do pr=1,Ncolumns 403 403 ! 1) Compute the PIA in all profiles containing hydrometeors 404 if ( (Ze_non_out(i,pr,cloudsat_preclvl) .gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl).gt.-100) ) then405 if ( (Ze_non_out(i,pr,cloudsat_preclvl) .lt.100) .and. (Ze_out(i,pr,cloudsat_preclvl).lt.100) ) then404 if ( (Ze_non_out(i,pr,cloudsat_preclvl)>-100) .and. (Ze_out(i,pr,cloudsat_preclvl)>-100) ) then 405 if ( (Ze_non_out(i,pr,cloudsat_preclvl)<100) .and. (Ze_out(i,pr,cloudsat_preclvl)<100) ) then 406 406 cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl) - Ze_out(i,pr,cloudsat_preclvl) 407 407 endif … … 412 412 ! 2a) Oceanic points. 413 413 ! ################################################################################ 414 if (land(i) .eq.0) then414 if (land(i) == 0) then 415 415 ! print*, 'aaa i, pr, fracPrecipIce(i,pr) : ', i, pr, fracPrecipIce(i,pr) !Artem 416 416 ! Snow 417 if(fracPrecipIce(i,pr) .gt.0.9) then418 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(2)) then417 if(fracPrecipIce(i,pr)>0.9) then 418 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(2)) then 419 419 cloudsat_pflag(i,pr) = pClass_Snow2 ! TSL: Snow certain 420 420 endif 421 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(4).and. &422 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(2)) then421 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. & 422 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(2)) then 423 423 cloudsat_pflag(i,pr) = pClass_Snow1 ! TSL: Snow possible 424 424 endif … … 426 426 427 427 ! Mixed 428 if(fracPrecipIce(i,pr) .gt.0.1.and.fracPrecipIce(i,pr).le.0.9) then429 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(2)) then428 if(fracPrecipIce(i,pr)>0.1.and.fracPrecipIce(i,pr)<=0.9) then 429 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(2)) then 430 430 cloudsat_pflag(i,pr) = pClass_Mixed2 ! TSL: Mixed certain 431 431 endif 432 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(4).and. &433 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(2)) then432 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. & 433 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(2)) then 434 434 cloudsat_pflag(i,pr) = pClass_Mixed1 ! TSL: Mixed possible 435 435 endif … … 437 437 438 438 ! Rain 439 if(fracPrecipIce(i,pr) .le.0.1) then440 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(1)) then439 if(fracPrecipIce(i,pr)<=0.1) then 440 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(1)) then 441 441 cloudsat_pflag(i,pr) = pClass_Rain3 ! TSL: Rain certain 442 442 endif 443 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(3).and. &444 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(1)) then443 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(3).and. & 444 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(1)) then 445 445 cloudsat_pflag(i,pr) = pClass_Rain2 ! TSL: Rain probable 446 446 endif 447 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(4).and. &448 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(3)) then447 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. & 448 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(3)) then 449 449 cloudsat_pflag(i,pr) = pClass_Rain1 ! TSL: Rain possible 450 450 endif 451 if(cloudsat_precip_pia(i,pr) .gt.40) then451 if(cloudsat_precip_pia(i,pr)>40) then 452 452 cloudsat_pflag(i,pr) = pClass_Rain4 ! TSL: Heavy Rain 453 453 endif … … 455 455 456 456 ! No precipitation 457 if(Ze_non_out(i,pr,cloudsat_preclvl) .le.-15) then457 if(Ze_non_out(i,pr,cloudsat_preclvl)<=-15) then 458 458 cloudsat_pflag(i,pr) = pClass_noPrecip ! TSL: Not Raining 459 459 endif … … 463 463 ! 2b) Land points. 464 464 ! ################################################################################ 465 if (land(i) .eq.1) then465 if (land(i) == 1) then 466 466 ! Find Zmax, the maximum reflectivity value in the attenuated profile (Ze_out); 467 467 Zmax=maxval(Ze_out(i,pr,:)) 468 468 469 469 ! Snow (T<273) 470 if(t2m(i) .lt.273._wp) then471 if(Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(5)) then470 if(t2m(i) < 273._wp) then 471 if(Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(5)) then 472 472 cloudsat_pflag(i,pr) = pClass_Snow2 ! JEK: Snow certain 473 473 endif 474 if(Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6) .and. &475 Ze_out(i,pr,cloudsat_preclvl) .le.Zbinvallnd(5)) then474 if(Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6) .and. & 475 Ze_out(i,pr,cloudsat_preclvl)<=Zbinvallnd(5)) then 476 476 cloudsat_pflag(i,pr) = pClass_Snow1 ! JEK: Snow possible 477 477 endif … … 479 479 480 480 ! Mized phase (273<T<275) 481 if(t2m(i) .ge. 273._wp .and. t2m(i) .le.275._wp) then482 if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &483 (Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(4))) then481 if(t2m(i) >= 273._wp .and. t2m(i) <= 275._wp) then 482 if ((Zmax > Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr)>30) .or. & 483 (Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(4))) then 484 484 cloudsat_pflag(i,pr) = pClass_Mixed2 ! JEK: Mixed certain 485 485 endif 486 if ((Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6) .and. &487 Ze_out(i,pr,cloudsat_preclvl) .le.Zbinvallnd(4)) .and. &488 (Zmax .gt.Zbinvallnd(5)) ) then486 if ((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6) .and. & 487 Ze_out(i,pr,cloudsat_preclvl) <= Zbinvallnd(4)) .and. & 488 (Zmax > Zbinvallnd(5)) ) then 489 489 cloudsat_pflag(i,pr) = pClass_Mixed1 ! JEK: Mixed possible 490 490 endif … … 492 492 493 493 ! Rain (T>275) 494 if(t2m(i) .gt.275) then495 if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &496 (Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(2))) then494 if(t2m(i) > 275) then 495 if ((Zmax > Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr)>30) .or. & 496 (Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(2))) then 497 497 cloudsat_pflag(i,pr) = pClass_Rain3 ! JEK: Rain certain 498 498 endif 499 if((Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6)) .and. &500 (Zmax .gt.Zbinvallnd(3))) then499 if((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6)) .and. & 500 (Zmax > Zbinvallnd(3))) then 501 501 cloudsat_pflag(i,pr) = pClass_Rain2 ! JEK: Rain probable 502 502 endif 503 if((Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6)) .and. &504 (Zmax .lt.Zbinvallnd(3))) then503 if((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6)) .and. & 504 (Zmax<Zbinvallnd(3))) then 505 505 cloudsat_pflag(i,pr) = pClass_Rain1 ! JEK: Rain possible 506 506 endif 507 if(cloudsat_precip_pia(i,pr) .gt.40) then507 if(cloudsat_precip_pia(i,pr)>40) then 508 508 cloudsat_pflag(i,pr) = pClass_Rain4 ! JEK: Heavy Rain 509 509 endif … … 511 511 512 512 ! No precipitation 513 if(Ze_out(i,pr,cloudsat_preclvl) .le.-15) then513 if(Ze_out(i,pr,cloudsat_preclvl)<=-15) then 514 514 cloudsat_pflag(i,pr) = pClass_noPrecip ! JEK: Not Precipitating 515 515 endif … … 526 526 ! Gridmean precipitation fraction for each precipitation type 527 527 do k=1,nCloudsatPrecipClass 528 if (any(cloudsat_pflag(i,:) .eq.k-1)) then529 cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) .eq.k-1)528 if (any(cloudsat_pflag(i,:) == k-1)) then 529 cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) == k-1) 530 530 endif 531 531 enddo -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90
r3491 r5081 172 172 173 173 ! Compute effective radius from number concentration and distribution parameters 174 if (Re_internal .eq.0) then174 if (Re_internal == 0) then 175 175 call calc_Re(hm_matrix(pr,k,tp),Np_matrix(pr,k,tp),rho_a, & 176 176 sd%dtype(tp),sd%apm(tp),sd%bpm(tp),sd%rho(tp),sd%p1(tp),sd%p2(tp),sd%p3(tp),Re) … … 187 187 ! Index into particle size dimension of scaling tables 188 188 iRe_type=1 189 if(Re .gt.0) then189 if(Re>0) then 190 190 ! Determine index in to scale LUT 191 191 ! Distance between Re points (defined by "base" and "step") for … … 197 197 base = rcfg%base_list(n+1) 198 198 iRe_type=Re/step 199 if (iRe_type .lt.1) iRe_type=1199 if (iRe_type<1) iRe_type=1 200 200 Re=step*(iRe_type+0.5_wp) ! set value of Re to closest value allowed in LUT. 201 201 iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step) 202 202 203 203 ! Make sure iRe_type is within bounds 204 if (iRe_type .ge.nRe_types) then204 if (iRe_type>=nRe_types) then 205 205 !write(*,*) 'Warning: size of Re exceed value permitted ', & 206 206 ! 'in Look-Up Table (LUT). Will calculate. ' … … 405 405 ! Exponential is same as modified gamma with vu =1 406 406 ! if Np is specified then we will just treat as modified gamma 407 if(dtype .eq. 2 .and. Np .gt.0) then407 if(dtype == 2 .and. Np > 0) then 408 408 local_dtype = 1 409 409 local_p3 = 1 … … 441 441 endif 442 442 443 if( Np .eq.0 .and. p2+1 > 1E-8) then ! use default value for MEAN diameter as first default443 if( Np==0 .and. p2+1 > 1E-8) then ! use default value for MEAN diameter as first default 444 444 dm = p2 ! by definition, should have units of microns 445 445 D0 = gamma(vu)/gamma(vu+1)*dm 446 446 else ! use value of Np 447 if(Np .eq.0) then447 if(Np==0) then 448 448 if( abs(p1+1) > 1E-8 ) then ! use default number concentration 449 449 local_Np = p1 ! total number concentration / pa --- units kg^-1 … … 525 525 526 526 ! get rg ... 527 if( Np .eq.0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg527 if( Np==0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg 528 528 rg = p2 529 529 else … … 826 826 log_sigma_g = p3 827 827 tmp2 = (bpm*log_sigma_g)*(bpm*log_sigma_g) 828 if(Re .le.0) then828 if(Re<=0) then 829 829 rg = p2 830 830 else -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90
r3491 r5081 75 75 76 76 ! Test for valid input overlap assumption 77 if (overlap .ne. 1 .and. overlap .ne. 2 .and. overlap .ne.3) then77 if (overlap /= 1 .and. overlap /= 2 .and. overlap /= 3) then 78 78 overlap=default_overlap 79 79 call errorMessage('ERROR(scops): Invalid overlap assumption provided. Using default overlap assumption (max/ran)') … … 92 92 tca(1:npoints,1:nlev) = cc(1:npoints,1:nlev) 93 93 94 if (ncolprint .ne.0) then94 if (ncolprint/=0) then 95 95 write (6,'(a)') 'frac_out_pp_rev:' 96 96 do j=1,npoints,1000 … … 102 102 write (6,'(I3)') ncol 103 103 endif 104 if (ncolprint .ne.0) then104 if (ncolprint/=0) then 105 105 write (6,'(a)') 'last_frac_pp:' 106 106 do j=1,npoints,1000 … … 122 122 123 123 ! Initialise threshold 124 IF (ilev .eq.1) then124 IF (ilev==1) then 125 125 ! If max overlap 126 IF (overlap .eq.1) then126 IF (overlap==1) then 127 127 ! Select pixels spread evenly across the gridbox 128 128 threshold(1:npoints,1:ncol)=boxpos(1:npoints,1:ncol) … … 137 137 enddo 138 138 ENDIF 139 IF (ncolprint .ne.0) then139 IF (ncolprint/=0) then 140 140 write (6,'(a)') 'threshold_nsf2:' 141 141 do j=1,npoints,1000 … … 147 147 ENDIF 148 148 149 IF (ncolprint .ne.0) then149 IF (ncolprint/=0) then 150 150 write (6,'(a)') 'ilev:' 151 151 write (6,'(I2)') ilev … … 157 157 !maxocc(1:npoints,ibox) = merge(1,0, conv(1:npoints,ilev) .gt. boxpos(1:npoints,ibox)) 158 158 do j=1,npoints 159 if (boxpos(j,ibox) .le.conv(j,ilev)) then159 if (boxpos(j,ibox)<=conv(j,ilev)) then 160 160 maxocc(j,ibox) = 1 161 161 else … … 165 165 166 166 ! Max overlap 167 if (overlap .eq.1) then167 if (overlap==1) then 168 168 threshold_min(1:npoints,ibox) = conv(1:npoints,ilev) 169 169 maxosc(1:npoints,ibox) = 1 … … 171 171 172 172 ! Random overlap 173 if (overlap .eq.2) then173 if (overlap==2) then 174 174 threshold_min(1:npoints,ibox) = conv(1:npoints,ilev) 175 175 maxosc(1:npoints,ibox) = 0 176 176 endif 177 177 ! Max/Random overlap 178 if (overlap .eq.3) then178 if (overlap==3) then 179 179 ! DS2014 START: The bounds on tca are not valid when ilev=1. 180 180 !threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev))) … … 182 182 ! min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. & 183 183 ! (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev))) 184 if (ilev .ne.1) then184 if (ilev /= 1) then 185 185 threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev))) 186 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt.&186 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) < & 187 187 min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. & 188 (threshold(1:npoints,ibox) .gt.conv(1:npoints,ilev)))188 (threshold(1:npoints,ibox)>conv(1:npoints,ilev))) 189 189 else 190 190 threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(0._wp,tca(1:npoints,ilev))) 191 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt.&191 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) < & 192 192 min(0._wp,tca(1:npoints,ilev)) .and. & 193 (threshold(1:npoints,ibox) .gt.conv(1:npoints,ilev)))193 (threshold(1:npoints,ibox)>conv(1:npoints,ilev))) 194 194 endif 195 195 endif … … 205 205 206 206 ! Fill frac_out with 1's where tca is greater than the threshold 207 frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev) .gt.threshold(1:npoints,ibox))207 frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev)>threshold(1:npoints,ibox)) 208 208 209 209 ! Code to partition boxes into startiform and convective parts goes here 210 where(threshold(1:npoints,ibox) .le.conv(1:npoints,ilev) .and. conv(1:npoints,ilev).gt.0.) frac_out(1:npoints,ibox,ilev)=2210 where(threshold(1:npoints,ibox)<=conv(1:npoints,ilev) .and. conv(1:npoints,ilev)>0.) frac_out(1:npoints,ibox,ilev)=2 211 211 ENDDO ! ibox 212 212 213 213 214 214 ! Set last_frac to tca at this level, so as to be tca from last level next time round 215 if (ncolprint .ne.0) then215 if (ncolprint/=0) then 216 216 do j=1,npoints ,1000 217 217 write(6,'(a10)') 'j='
Note: See TracChangeset
for help on using the changeset viewer.
