Changeset 5082 for LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90
- Timestamp:
- Jul 19, 2024, 5:41:58 PM (4 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90
r3491 r5082 134 134 ! ########################################################################## 135 135 136 if (debugcol .ne.0) then136 if (debugcol/=0) then 137 137 do j=1,npoints,debugcol 138 138 … … 140 140 do ilev=1,nlev 141 141 acc(ilev,1:ncol)=frac_out(j,1:ncol,ilev)*2 142 where(levmatch(j,1:ncol) .eq.ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1142 where(levmatch(j,1:ncol) == ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1 143 143 enddo 144 144 … … 224 224 225 225 ! Set tropopause values 226 if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then226 if (isccp_top_height == 1 .or. isccp_top_height == 3) then 227 227 ptrop(1:npoints) = 5000._wp 228 228 attropmin(1:npoints) = 400._wp … … 232 232 233 233 do ilev=1,nlev 234 where(pfull(1:npoints,ilev) .lt.40000. .and. &235 pfull(1:npoints,ilev) .gt.5000. .and. &236 at(1:npoints,ilev) .lt.attropmin(1:npoints))234 where(pfull(1:npoints,ilev) < 40000. .and. & 235 pfull(1:npoints,ilev) > 5000. .and. & 236 at(1:npoints,ilev) < attropmin(1:npoints)) 237 237 ptrop(1:npoints) = pfull(1:npoints,ilev) 238 238 attropmin(1:npoints) = at(1:npoints,ilev) … … 244 244 do ilev=1,nlev 245 245 atmax(1:npoints) = merge(at(1:npoints,ilev),atmax(1:npoints),& 246 at(1:npoints,ilev) .gt. atmax(1:npoints) .and. ilev .ge.itrop(1:npoints))246 at(1:npoints,ilev) > atmax(1:npoints) .and. ilev >= itrop(1:npoints)) 247 247 enddo 248 248 end if 249 249 250 if (isccp_top_height .eq. 1 .or. isccp_top_height .eq.3) then250 if (isccp_top_height == 1 .or. isccp_top_height == 3) then 251 251 ! ############################################################################ 252 252 ! Clear-sky radiance calculation … … 308 308 dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), & 309 309 1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), & 310 demIN(1:npoints,ibox,ilev) .eq.0)310 demIN(1:npoints,ibox,ilev) == 0) 311 311 312 312 ! Increase TOA flux emitted from layer … … 348 348 tauir(1:npoints) = tau(1:npoints,ibox)/2.13_wp 349 349 taumin(1:npoints) = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp)) 350 if (isccp_top_height .eq.1) then350 if (isccp_top_height == 1) then 351 351 do j=1,npoints 352 if (transmax(j) .gt. 0.001 .and. transmax(j) .le.0.9999999) then352 if (transmax(j) > 0.001 .and. transmax(j) <= 0.9999999) then 353 353 fluxtopinit(j) = fluxtop(j,ibox) 354 354 tauir(j) = tau(j,ibox)/2.13_wp … … 357 357 do icycle=1,2 358 358 do j=1,npoints 359 if (tau(j,ibox) .gt. (tauchk)) then360 if (transmax(j) .gt. 0.001 .and. transmax(j) .le.0.9999999) then359 if (tau(j,ibox) > (tauchk)) then 360 if (transmax(j) > 0.001 .and. transmax(j) <= 0.9999999) then 361 361 emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j) ) 362 362 fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j)) 363 363 fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox))) 364 364 tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox)))) 365 if (tb(j,ibox) .gt.260.) then365 if (tb(j,ibox) > 260.) then 366 366 tauir(j) = tau(j,ibox) / 2.56_wp 367 367 end if … … 373 373 374 374 ! Cloud-top temperature 375 where(tau(1:npoints,ibox) .gt.tauchk)375 where(tau(1:npoints,ibox) > tauchk) 376 376 tb(1:npoints,ibox)= 1307.27_wp/ (log(1. + (1._wp/fluxtop(1:npoints,ibox)))) 377 where (isccp_top_height .eq. 1 .and. tauir(1:npoints) .lt.taumin(1:npoints))377 where (isccp_top_height == 1 .and. tauir(1:npoints) < taumin(1:npoints)) 378 378 tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp 379 379 tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints) … … 382 382 383 383 ! Clear-sky brightness temperature 384 where(tau(1:npoints,ibox) .le. tauchk)384 where(tau(1:npoints,ibox) <= tauchk) 385 385 tb(1:npoints,ibox) = meantbclr(1:npoints) 386 386 endwhere … … 399 399 do ibox=1,ncol 400 400 !segregate according to optical thickness 401 if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then401 if (isccp_top_height == 1 .or. isccp_top_height == 3) then 402 402 403 403 ! Find level whose temperature most closely matches brightness temperature 404 404 nmatch(1:npoints)=0 405 405 do k1=1,nlev-1 406 ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2)406 ilev = merge(nlev-k1,k1,isccp_top_height_direction == 2) 407 407 do j=1,npoints 408 if (ilev .ge.itrop(j) .and. &409 ((at(j,ilev) .ge. tb(j,ibox) .and. &410 at(j,ilev+1) .le.tb(j,ibox)) .or. &411 (at(j,ilev) .le.tb(j,ibox) .and. &412 at(j,ilev+1) .ge. tb(j,ibox)))) then408 if (ilev >= itrop(j) .and. & 409 ((at(j,ilev) >= tb(j,ibox) .and. & 410 at(j,ilev+1) <= tb(j,ibox)) .or. & 411 (at(j,ilev) <= tb(j,ibox) .and. & 412 at(j,ilev+1) >= tb(j,ibox)))) then 413 413 nmatch(j)=nmatch(j)+1 414 414 match(j,nmatch(j))=ilev … … 418 418 419 419 do j=1,npoints 420 if (nmatch(j) .ge.1) then420 if (nmatch(j) >= 1) then 421 421 k1 = match(j,nmatch(j)) 422 422 k2 = k1 + 1 … … 426 426 logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd 427 427 ptop(j,ibox) = exp(logp) 428 levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) .lt.abs(pfull(j,k2)-ptop(j,ibox)))428 levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) < abs(pfull(j,k2)-ptop(j,ibox))) 429 429 else 430 if (tb(j,ibox) .le.attrop(j)) then430 if (tb(j,ibox) <= attrop(j)) then 431 431 ptop(j,ibox)=ptrop(j) 432 432 levmatch(j,ibox)=itrop(j) 433 433 end if 434 if (tb(j,ibox) .ge.atmax(j)) then434 if (tb(j,ibox) >= atmax(j)) then 435 435 ptop(j,ibox)=pfull(j,nlev) 436 436 levmatch(j,ibox)=nlev … … 441 441 ptop(1:npoints,ibox)=0. 442 442 do ilev=1,nlev 443 where((ptop(1:npoints,ibox) .eq. 0. ) .and.(frac_out(1:npoints,ibox,ilev) .ne.0))443 where((ptop(1:npoints,ibox) == 0. ) .and.(frac_out(1:npoints,ibox,ilev) /= 0)) 444 444 ptop(1:npoints,ibox)=phalf(1:npoints,ilev) 445 445 levmatch(1:npoints,ibox)=ilev … … 447 447 end do 448 448 end if 449 where(tau(1:npoints,ibox) .le.tauchk)449 where(tau(1:npoints,ibox) <= tauchk) 450 450 ptop(1:npoints,ibox)=0._wp 451 451 levmatch(1:npoints,ibox)=0._wp … … 460 460 do ibox=1,ncol 461 461 do j=1,npoints 462 if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt.0.) then463 if (sunlit(j) .eq.1 .or. isccp_top_height .eq.3) then462 if (tau(j,ibox) > (tauchk) .and. ptop(j,ibox) > 0.) then 463 if (sunlit(j)==1 .or. isccp_top_height == 3) then 464 464 boxtau(j,ibox) = tau(j,ibox) 465 465 boxptop(j,ibox) = ptop(j,ibox)!/100._wp … … 508 508 ! Brightness Temperature 509 509 ! #################################################################################### 510 if (isccp_top_height .eq. 1 .or. isccp_top_height .eq.3) then510 if (isccp_top_height == 1 .or. isccp_top_height == 3) then 511 511 meantb(1:npoints)=sum(boxttop,2)/ncol 512 512 else … … 535 535 do ilev2=1,7 536 536 do j=1,npoints ! 537 if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then537 if (sunlit(j)==1 .or. isccp_top_height == 3) then 538 538 fq_isccp(j,ilev,ilev2)= 0. 539 539 else … … 546 546 547 547 ! Reset variables need for averaging cloud properties 548 where(sunlit .eq. 1 .or. isccp_top_height .eq.3)548 where(sunlit == 1 .or. isccp_top_height == 3) 549 549 totalcldarea(1:npoints) = 0._wp 550 550 meanalbedocld(1:npoints) = 0._wp … … 561 561 do j=1,npoints 562 562 ! Subcolumns that are cloudy(true) and not(false) 563 box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt.0.)563 box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) > tauchk .and. boxptop(j,1:ncol) > 0.) 564 564 565 565 ! Compute joint histogram and column quantities for points that are sunlit and cloudy 566 if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then566 if (sunlit(j) ==1 .or. isccp_top_height == 3) then 567 567 ! Joint-histogram 568 568 call hist2D(boxtau(j,1:ncol),boxptop(j,1:ncol),ncol,isccp_histTau,numISCCPTauBins, & … … 572 572 573 573 ! Column cloud area 574 totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt.isccp_taumin))/ncol574 totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) > isccp_taumin))/ncol 575 575 576 576 ! Subcolumn cloud albedo 577 577 !albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp),& 578 578 ! 0._wp,box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin) 579 where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt.isccp_taumin)579 where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) > isccp_taumin) 580 580 albedocld(j,1:ncol) = (boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp) 581 581 elsewhere … … 587 587 588 588 ! Column cloud top pressure 589 meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt.isccp_taumin)/ncol589 meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) > isccp_taumin)/ncol 590 590 endif 591 591 enddo 592 592 593 593 ! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0 594 where(totalcldarea(1:npoints) .gt.0)594 where(totalcldarea(1:npoints) > 0) 595 595 meanptop(1:npoints) = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints) 596 596 meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints) … … 609 609 610 610 ! Represent in percent 611 where(totalcldarea .ne.output_missing_value) totalcldarea = totalcldarea*100._wp612 where(fq_isccp .ne.output_missing_value) fq_isccp = fq_isccp*100._wp611 where(totalcldarea /= output_missing_value) totalcldarea = totalcldarea*100._wp 612 where(fq_isccp /= output_missing_value) fq_isccp = fq_isccp*100._wp 613 613 614 614 … … 634 634 varOUT(1:dim1,1:dim2,1:dim3) = 0._wp 635 635 do j=1,dim2 636 where(flag(:,j,:) .eq.1)636 where(flag(:,j,:) == 1) 637 637 varOUT(:,j,:) = varIN2 638 638 endwhere 639 where(flag(:,j,:) .eq.2)639 where(flag(:,j,:) == 2) 640 640 varOUT(:,j,:) = varIN1 641 641 endwhere
Note: See TracChangeset
for help on using the changeset viewer.
