- Timestamp:
- Aug 2, 2024, 2:12:03 PM (7 weeks ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lidar_simulator.F90
r5099 r5158 175 175 ! PLANE PARRALLEL FIELDS 176 176 ! #################################################################################### 177 doicol=1,ncolumns177 DO icol=1,ncolumns 178 178 ! ################################################################################# 179 179 ! *) Total Backscatter signal … … 200 200 ! #################################################################################### 201 201 if (lphaseoptics) then 202 doicol=1,ncolumns202 DO icol=1,ncolumns 203 203 ! ################################################################################# 204 204 ! *) Ice/Liq Perpendicular Backscatter signal … … 206 206 ! Computation of ATBperp,ice/liq from ATBice/liq including the multiple scattering 207 207 ! contribution (Cesana and Chepfer 2013, JGR) 208 dok=1,nlev208 DO k=1,nlev 209 209 ! Ice particles 210 210 pnorm_perp_ice(1:npoints,icol,k) = Alpha * pnorm_ice(1:npoints,icol,k) … … 242 242 243 243 ! Other layers 244 dok=2,nlev244 DO k=2,nlev 245 245 ! Optical thickness of layer k 246 246 tautot_lay(1:npoints) = tautot(1:npoints,icol,k)-tautot(1:npoints,icol,k-1) … … 398 398 ! Compute LIDAR scattering ratio 399 399 if (use_vgrid) then 400 doic = 1, ncol400 DO ic = 1, ncol 401 401 pnorm_c = pnormFlip(:,ic,:) 402 402 where ((pnorm_c .lt. xmax) .and. (betamolFlip(:,1,:) .lt. xmax) .and. & … … 427 427 endif 428 428 else 429 doic = 1, ncol429 DO ic = 1, ncol 430 430 pnorm_c = pnorm(:,ic,:) 431 431 where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) … … 458 458 if (ok_lidar_cfad) then 459 459 ! CFADs of subgrid-scale lidar scattering ratios 460 doi=1,Npoints461 doj=1,llm460 DO i=1,Npoints 461 DO j=1,llm 462 462 cfad2(i,:,j) = hist1D(ncol,x3d(i,:,j),SR_BINS,histBsct) 463 463 enddo … … 499 499 500 500 ! Other layers 501 dok=2,nlev501 DO k=2,nlev 502 502 tautot_lay(:) = tau(:,k)-tau(:,k-1) 503 503 WHERE (tautot_lay(:) .gt. 0.) … … 527 527 epsrealwp = epsilon(1._wp) 528 528 beta(:,1) = pnorm(:,1) * (2._wp*tau(:,1))/(1._wp-exp(-2._wp*tau(:,1))) 529 dok=2,nlev529 DO k=2,nlev 530 530 tautot_lay(:) = tau(:,k)-tau(:,k-1) 531 531 WHERE ( EXP(-2._wp*tau(:,k-1)) .gt. epsrealwp ) … … 569 569 zeta50 = -9.4776e-07_wp ! 570 570 571 571 ! Inputs 572 572 integer,intent(in) :: & 573 573 Npoints, & ! Number of gridpoints … … 576 576 Ncat, & ! Number of cloud layer types 577 577 Nphase ! Number of cloud layer phase types 578 578 ! [ice,liquid,undefined,false ice,false liquid,Percent of ice] 579 579 real(wp),intent(in) :: & 580 580 S_att, & ! … … 590 590 pplay ! Pressure 591 591 592 592 ! Outputs 593 593 real(wp),intent(out),dimension(Npoints,Ntemp,5) :: & 594 594 lidarcldtemp ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq … … 625 625 626 626 ! #################################################################################### 627 ! 1) Initialize 627 ! 1) Initialize 628 628 ! #################################################################################### 629 629 lidarcld = 0._wp … … 648 648 ! 2) Cloud detection 649 649 ! #################################################################################### 650 dok=1,Nlevels650 DO k=1,Nlevels 651 651 ! Cloud detection at subgrid-scale: 652 652 where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) ) … … 671 671 cldlay = 0._wp 672 672 nsublay = 0._wp 673 dok=1,Nlevels674 doic = 1, Ncolumns675 doip = 1, Npoints673 DO k=1,Nlevels 674 DO ic = 1, Ncolumns 675 DO ip = 1, Npoints 676 676 677 677 ! Computation of the cloud fraction as a function of the temperature instead 678 678 ! of height, for ice,liquid and all clouds 679 679 if(srok(ip,ic,k).gt.0.)then 680 doitemp=1,Ntemp680 DO itemp=1,Ntemp 681 681 if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then 682 682 lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1._wp … … 686 686 687 687 if(cldy(ip,ic,k).eq.1.)then 688 do itemp=1,Ntemp688 DO itemp=1,Ntemp 689 689 if( (tmp(ip,k) .ge. tempmod(itemp)).and.(tmp(ip,k) .lt. tempmod(itemp+1)) )then 690 690 lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1._wp … … 723 723 cldlayer = 0._wp 724 724 nsublayer = 0._wp 725 doiz = 1, Ncat726 doic = 1, Ncolumns725 DO iz = 1, Ncat 726 DO ic = 1, Ncolumns 727 727 cldlayer(:,iz) = cldlayer(:,iz) + cldlay(:,ic,iz) 728 728 nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz) … … 742 742 ! 4.1) For Cloudy pixels with 8.16km < z < 19.2km 743 743 ! #################################################################################### 744 doncol=1,Ncolumns745 do i=1,Npoints746 donlev=1,23 ! from 19.2km until 8.16km744 DO ncol=1,Ncolumns 745 DO i=1,Npoints 746 DO nlev=1,23 ! from 19.2km until 8.16km 747 747 p1 = pplay(1,nlev) 748 748 … … 856 856 ! ############################################################################## 857 857 toplvlsat = 0 858 donlev=24,Nlevels! from 8.16km until 0km858 DO nlev=24,Nlevels! from 8.16km until 0km 859 859 p1 = pplay(i,nlev) 860 860 … … 979 979 ! ############################################################################## 980 980 if(toplvlsat.ne.0) then 981 donlev = toplvlsat,Nlevels981 DO nlev = toplvlsat,Nlevels 982 982 p1 = pplay(i,nlev) 983 983 if(cldy(i,ncol,nlev).eq.1.)then … … 1031 1031 1032 1032 ! Compute Phase low mid high cloud fractions 1033 doiz = 1, Ncat1034 doi=1,Nphase-31035 doic = 1, Ncolumns1033 DO iz = 1, Ncat 1034 DO i=1,Nphase-3 1035 DO ic = 1, Ncolumns 1036 1036 cldlayerphase(:,iz,i) = cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i) 1037 1037 cldlayerphasesum(:,iz) = cldlayerphasesum(:,iz) + cldlayphase(:,ic,iz,i) … … 1039 1039 enddo 1040 1040 enddo 1041 doiz = 1, Ncat1042 doi=4,51043 doic = 1, Ncolumns1041 DO iz = 1, Ncat 1042 DO i=4,5 1043 DO ic = 1, Ncolumns 1044 1044 cldlayerphase(:,iz,i) = cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i) 1045 1045 enddo … … 1055 1055 ENDWHERE 1056 1056 1057 doi=1,Nphase-11057 DO i=1,Nphase-1 1058 1058 WHERE ( cldlayerphasesum(:,:).gt.0.0 ) 1059 1059 cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:) … … 1061 1061 enddo 1062 1062 1063 doi=1,Npoints1064 doiz=1,Ncat1063 DO i=1,Npoints 1064 DO iz=1,Ncat 1065 1065 checkcldlayerphase=0. 1066 1066 checkcldlayerphase2=0. 1067 1067 if (cldlayerphasesum(i,iz) .gt. 0.0 )then 1068 doic=1,Nphase-31068 DO ic=1,Nphase-3 1069 1069 checkcldlayerphase = checkcldlayerphase+cldlayerphase(i,iz,ic) 1070 1070 enddo … … 1075 1075 enddo 1076 1076 1077 doi=1,Nphase-11077 DO i=1,Nphase-1 1078 1078 WHERE (nsublayer(:,:) .eq. 0.0) 1079 1079 cldlayerphase(:,:,i) = undef … … 1082 1082 1083 1083 ! Compute Phase 3D as a function of temperature 1084 donlev=1,Nlevels1085 doncol=1,Ncolumns1086 doi=1,Npoints1087 doitemp=1,Ntemp1084 DO nlev=1,Nlevels 1085 DO ncol=1,Ncolumns 1086 DO i=1,Npoints 1087 DO itemp=1,Ntemp 1088 1088 if(tmpi(i,ncol,nlev).gt.0.)then 1089 1089 if((tmpi(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpi(i,ncol,nlev) .lt. tempmod(itemp+1)) )then … … 1105 1105 1106 1106 ! Check temperature cloud fraction 1107 doi=1,Npoints1108 doitemp=1,Ntemp1107 DO i=1,Npoints 1108 DO itemp=1,Ntemp 1109 1109 checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4) 1110 1110 !if(checktemp .NE. lidarcldtemp(i,itemp,1))then … … 1124 1124 ENDWHERE 1125 1125 1126 doi=1,41126 DO i=1,4 1127 1127 WHERE(lidarcldtempind(:,:) .gt. 0.) 1128 1128 lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:) … … 1191 1191 ! 2) Cloud detection 1192 1192 ! #################################################################################### 1193 dok=1,Nlevels1193 DO k=1,Nlevels 1194 1194 ! Cloud detection at subgrid-scale: 1195 1195 where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) ) … … 1210 1210 ! 3) Grid-box 3D cloud fraction and layered cloud fractions(ISCCP pressure categories) 1211 1211 ! #################################################################################### 1212 dok=1,Nlevels1213 doic = 1, Ncolumns1214 doip = 1, Npoints1212 DO k=1,Nlevels 1213 DO ic = 1, Ncolumns 1214 DO ip = 1, Npoints 1215 1215 1216 1216 iz=1 … … 1244 1244 cldlayer = 0._wp 1245 1245 nsublayer = 0._wp 1246 doiz = 1, Ncat1247 doic = 1, Ncolumns1246 DO iz = 1, Ncat 1247 DO ic = 1, Ncolumns 1248 1248 cldlayer(:,iz) = cldlayer(:,iz) + cldlay(:,ic,iz) 1249 1249 nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz) … … 1272 1272 eta = 0.6_wp ! Multiple-scattering factor (Vaillant de Guelis et al. 2017a, AMT) 1273 1273 1274 1274 ! Inputs 1275 1275 integer,intent(in) :: & 1276 1276 Npoints, & ! Number of gridpoints … … 1291 1291 surfelev ! Surface Elevation (SE) 1292 1292 1293 1293 ! Outputs 1294 1294 real(wp),intent(out),dimension(Npoints,Nlevels,Ntype+1) :: & 1295 1295 lidarcldtype ! 3D OPAQ product fraction (opaque clouds, thin clouds, z_opaque, opacity) … … 1324 1324 1325 1325 ! #################################################################################### 1326 ! 1) Initialize 1326 ! 1) Initialize 1327 1327 ! #################################################################################### 1328 1328 cldtype(:,:) = 0._wp … … 1342 1342 ! 2) Cloud detection and Fully attenuated layer detection 1343 1343 ! #################################################################################### 1344 dok=1,Nlevels1344 DO k=1,Nlevels 1345 1345 ! Cloud detection at subgrid-scale: 1346 1346 where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) ) … … 1375 1375 ! #################################################################################### 1376 1376 1377 dok=1,Nlevels1378 doic = 1, Ncolumns1379 doip = 1, Npoints1377 DO k=1,Nlevels 1378 DO ic = 1, Ncolumns 1379 DO ip = 1, Npoints 1380 1380 1381 1381 cldlay(ip,ic,1) = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque cloud … … 1393 1393 1394 1394 ! OPAQ variables 1395 doic = 1, Ncolumns1396 doip = 1, Npoints1395 DO ic = 1, Ncolumns 1396 DO ip = 1, Npoints 1397 1397 1398 1398 ! Declaring non-opaque cloudy profiles as thin cloud profiles 1399 1400 1401 1399 if ( cldlay(ip,ic,4).gt. 0. .and. cldlay(ip,ic,1) .eq. 0. ) then 1400 cldlay(ip,ic,2) = 1._wp 1401 endif 1402 1402 1403 1403 ! Filling in 3D and 2D variables 1404 1404 1405 1405 ! Opaque cloud profiles 1406 1407 1408 1409 dok=1,Nlevels-11406 if ( cldlay(ip,ic,1) .eq. 1. ) then 1407 zopac = 0._wp 1408 z_top = 0._wp 1409 DO k=1,Nlevels-1 1410 1410 ! Declaring z_opaque altitude and opaque cloud fraction for 3D and 2D variables 1411 1411 ! From SFC-2-TOA ( actually from vgrid_z(SFC+1) = vgrid_z(Nlevels-1) ) 1412 1413 1414 1415 1416 1417 1418 1419 1420 1412 if ( cldy(ip,ic,Nlevels-k) .eq. 1. .and. zopac .eq. 0. ) then 1413 lidarcldtype(ip,Nlevels-k + 1,3) = lidarcldtype(ip,Nlevels-k + 1,3) + 1._wp 1414 cldlay(ip,ic,3) = vgrid_z(Nlevels-k+1) ! z_opaque altitude 1415 nsublay(ip,ic,3) = 1._wp 1416 zopac = Nlevels-k+1 ! z_opaque vertical index on vgrid_z 1417 endif 1418 if ( cldy(ip,ic,Nlevels-k) .eq. 1. ) then 1419 lidarcldtype(ip,Nlevels-k ,1) = lidarcldtype(ip,Nlevels-k ,1) + 1._wp 1420 z_top = Nlevels-k ! top cloud layer vertical index on vgrid_z 1421 1421 endif 1422 1422 enddo 1423 1423 ! Summing opaque cloud mean temperatures and altitudes 1424 1424 ! as defined in Vaillant de Guelis et al. 2017a, AMT … … 1432 1432 cldlay(ip,ic,1) = 0 1433 1433 endif 1434 1434 endif 1435 1435 1436 1436 ! Thin cloud profiles 1437 1438 1439 1440 1441 dok=1,Nlevels1437 if ( cldlay(ip,ic,2) .eq. 1. ) then 1438 topcloud = 0._wp 1439 z_top = 0._wp 1440 z_base = 0._wp 1441 DO k=1,Nlevels 1442 1442 ! Declaring thin cloud fraction for 3D variable 1443 1443 ! From TOA-2-SFC 1444 1444 if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 1. ) then 1445 1445 lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp 1446 1446 z_base = k ! bottom cloud layer 1447 1447 endif 1448 1448 if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 0. ) then 1449 1449 lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp 1450 1451 1450 z_top = k ! top cloud layer 1451 z_base = k ! bottom cloud layer 1452 1452 topcloud = 1._wp 1453 1454 1453 endif 1454 enddo 1455 1455 ! Computing mean emissivity using layers below the bottom cloud layer to the surface 1456 1457 1458 1459 dok=z_base+1,Nlevels1460 1461 1462 1456 srmean = 0._wp 1457 srcount = 0._wp 1458 cloudemis = 0._wp 1459 DO k=z_base+1,Nlevels 1460 if ( (x(ip,ic,k) .gt. S_att_opaq) .and. (x(ip,ic,k) .lt. 1.0) .and. (x(ip,ic,k) .ne. undef) ) then 1461 srmean = srmean + x(ip,ic,k) 1462 srcount = srcount + 1. 1463 1463 endif 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1464 enddo 1465 ! If clear sky layers exist below bottom cloud layer 1466 if ( srcount .gt. 0. ) then 1467 trans2 = srmean/srcount ! thin cloud transmittance**2 1468 tau_app = -(log(trans2))/2. ! apparent cloud optical depth 1469 tau_vis = tau_app/eta ! cloud visible optical depth (multiple scat.) 1470 tau_ir = tau_vis/2. ! approx. relation between visible and IR ODs 1471 cloudemis = 1. - exp(-tau_ir) ! no diffusion in IR considered : emis = 1-T 1472 count_emis(ip) = count_emis(ip) + 1. 1473 endif 1474 1474 ! Summing thin cloud mean temperatures and altitudes 1475 1475 ! as defined in Vaillant de Guelis et al. 2017a, AMT … … 1500 1500 ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=1) to SFC(k=Nlevels) 1501 1501 lidarcldtype(:,1,4) = lidarcldtype(:,1,3) !top layer equal to 3D z_opaque fraction 1502 doip = 1, Npoints1503 dok = 2, Nlevels1502 DO ip = 1, Npoints 1503 DO k = 2, Nlevels 1504 1504 if ( (lidarcldtype(ip,k,3) .ne. undef) .and. (lidarcldtype(ip,k-1,4) .ne. undef) ) then 1505 lidarcldtype(ip,k,4) = lidarcldtype(ip,k,3) + lidarcldtype(ip,k-1,4) 1506 else 1507 lidarcldtype(ip,k,4) = undef 1508 endif 1509 enddo 1505 lidarcldtype(ip,k,4) = lidarcldtype(ip,k,3) + lidarcldtype(ip,k-1,4) 1506 else 1507 lidarcldtype(ip,k,4) = undef 1508 endif 1510 1509 enddo 1510 enddo 1511 1511 1512 1512 ! Layered cloud types (opaque, thin and z_opaque 2D variables) 1513 1513 1514 doiz = 1, Ntype1515 doic = 1, Ncolumns1514 DO iz = 1, Ntype 1515 DO ic = 1, Ncolumns 1516 1516 cldtype(:,iz) = cldtype(:,iz) + cldlay(:,ic,iz) 1517 1517 nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
Note: See TracChangeset
for help on using the changeset viewer.
