! $Id $ SUBROUTINE isccp_cloud_types(debug, debugcol, npoints, sunlit, nlev, ncol, & seed, pfull, phalf, qv, cc, conv, dtau_s, dtau_c, top_height, overlap, & tautab, invtau, skt, emsfc_lw, at, dem_s, dem_c, fq_isccp, totalcldarea, & meanptop, meantaucld, boxtau, boxptop) ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved. ! This code is available without charge with the following conditions: ! 1. The code is available for scientific purposes and is not for ! commercial use. ! 2. Any improvements you make to the code should be made available ! to the to the authors for incorporation into a future release. ! 3. The code should not be used in any way that brings the authors ! or their employers into disrepute. IMPLICIT NONE ! NOTE: the maximum number of levels and columns is set by ! the following parameter statement INTEGER ncolprint ! ----- ! Input ! ----- INTEGER npoints ! number of model points in the horizontal ! PARAMETER(npoints=6722) INTEGER nlev ! number of model levels in column INTEGER ncol ! number of subcolumns INTEGER sunlit(npoints) ! 1 for day points, 0 for night time INTEGER seed(npoints) ! seed value for random number generator ! ! ( see Numerical Recipes Chapter 7) ! ! It is recommended that the seed is set ! ! to a different value for each model ! ! gridbox it is called on, as it is ! ! possible that the choice of the samec ! ! seed value every time may introduce some ! ! statistical bias in the results, particularly ! ! for low values of NCOL. REAL pfull(npoints, nlev) ! pressure of full model levels (Pascals) ! ! pfull(npoints,1) is top level of model ! ! pfull(npoints,nlev) is bottom level of model REAL phalf(npoints, nlev+1) ! pressure of half model levels (Pascals) ! ! phalf(npoints,1) is top of model ! ! phalf(npoints,nlev+1) is the surface pressure REAL qv(npoints, nlev) ! water vapor specific humidity (kg vapor/ kg air) ! ! on full model levels REAL cc(npoints, nlev) ! input cloud cover in each model level (fraction) ! ! NOTE: This is the HORIZONTAL area of each ! ! grid box covered by clouds REAL conv(npoints, nlev) ! input convective cloud cover in each model level (fraction) ! ! NOTE: This is the HORIZONTAL area of each ! ! grid box covered by convective clouds REAL dtau_s(npoints, nlev) ! mean 0.67 micron optical depth of stratiform ! ! clouds in each model level ! ! NOTE: this the cloud optical depth of only the ! ! cloudy part of the grid box, it is not weighted ! ! with the 0 cloud optical depth of the clear ! ! part of the grid box REAL dtau_c(npoints, nlev) ! mean 0.67 micron optical depth of convective ! ! clouds in each ! ! model level. Same note applies as in dtau_s. INTEGER overlap ! overlap type ! 1=max ! 2=rand ! 3=max/rand INTEGER top_height ! 1 = adjust top height using both a computed ! ! infrared brightness temperature and the visible ! ! optical depth to adjust cloud top pressure. Note ! ! that this calculation is most appropriate to compare ! ! to ISCCP data during sunlit hours. ! ! 2 = do not adjust top height, that is cloud top ! ! pressure is the actual cloud top pressure ! ! in the model ! ! 3 = adjust top height using only the computed ! ! infrared brightness temperature. Note that this ! ! calculation is most appropriate to compare to ISCCP ! ! IR only algortihm (i.e. you can compare to nighttime ! ! ISCCP data with this option) REAL tautab(0:255) ! ISCCP table for converting count value to ! ! optical thickness INTEGER invtau(-20:45000) ! ISCCP table for converting optical thickness ! ! to count value ! The following input variables are used only if top_height = 1 or ! top_height = 3 REAL skt(npoints) ! skin Temperature (K) REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction) REAL at(npoints, nlev) ! temperature in each model level (K) REAL dem_s(npoints, nlev) ! 10.5 micron longwave emissivity of stratiform ! ! clouds in each ! ! model level. Same note applies as in dtau_s. REAL dem_c(npoints, nlev) ! 10.5 micron longwave emissivity of convective ! ! clouds in each ! ! model level. Same note applies as in dtau_s. ! IM reg.dyn BEG REAL t1, t2 ! REAL w(npoints) !vertical wind at 500 hPa ! LOGICAL pct_ocean(npoints) !TRUE if oceanic point, FALSE otherway ! INTEGER iw(npoints) , nw ! REAL wmin, pas_w ! INTEGER k, l, iwmx ! PARAMETER(wmin=-100.,pas_w=10.,iwmx=30) ! REAL fq_dynreg(7,7,iwmx) ! REAL nfq_dynreg(7,7,iwmx) ! LOGICAL pctj(7,7,iwmx) ! IM reg.dyn END ! ------ ! Output ! ------ REAL fq_isccp(npoints, 7, 7) ! the fraction of the model grid box covered by ! ! each of the 49 ISCCP D level cloud types REAL totalcldarea(npoints) ! the fraction of model grid box columns ! ! with cloud somewhere in them. This should ! ! equal the sum over all entries of fq_isccp ! ! The following three means are averages over the cloudy areas only. If ! no ! ! clouds are in grid box all three quantities should equal zero. REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging ! ! in cloud top pressure. REAL meantaucld(npoints) ! mean optical thickness ! ! linear averaging in albedo performed. REAL boxtau(npoints, ncol) ! optical thickness in each column REAL boxptop(npoints, ncol) ! cloud top pressure (mb) in each column ! ------ ! Working variables added when program updated to mimic Mark Webb's PV-Wave ! code ! ------ REAL frac_out(npoints, ncol, nlev) ! boxes gridbox divided up into ! ! Equivalent of BOX in original version, but ! ! indexed by column then row, rather than ! ! by row then column REAL tca(npoints, 0:nlev) ! total cloud cover in each model level (fraction) ! ! with extra layer of zeroes on top ! ! in this version this just contains the values input ! ! from cc but with an extra level REAL cca(npoints, nlev) ! convective cloud cover in each model level (fraction) ! ! from conv REAL threshold(npoints, ncol) ! pointer to position in gridbox REAL maxocc(npoints, ncol) ! Flag for max overlapped conv cld REAL maxosc(npoints, ncol) ! Flag for max overlapped strat cld REAL boxpos(npoints, ncol) ! ordered pointer to position in gridbox REAL threshold_min(npoints, ncol) ! minimum value to define range in with new threshold ! ! is chosen REAL dem(npoints, ncol), bb(npoints) ! working variables for 10.5 micron longwave ! ! emissivity in part of ! ! gridbox under consideration REAL ran(npoints) ! vector of random numbers REAL ptrop(npoints) REAL attrop(npoints) REAL attropmin(npoints) REAL atmax(npoints) REAL atmin(npoints) REAL btcmin(npoints) REAL transmax(npoints) INTEGER i, j, ilev, ibox, itrop(npoints) INTEGER ipres(npoints) INTEGER itau(npoints), ilev2 INTEGER acc(nlev, ncol) INTEGER match(npoints, nlev-1) INTEGER nmatch(npoints) INTEGER levmatch(npoints, ncol) ! !variables needed for water vapor continuum absorption REAL fluxtop_clrsky(npoints), trans_layers_above_clrsky(npoints) REAL taumin(npoints) REAL dem_wv(npoints, nlev), wtmair, wtmh20, navo, grav, pstd, t0 REAL press(npoints), dpress(npoints), atmden(npoints) REAL rvh20(npoints), wk(npoints), rhoave(npoints) REAL rh20s(npoints), rfrgn(npoints) REAL tmpexp(npoints), tauwv(npoints) CHARACTER *1 cchar(6), cchar_realtops(6) INTEGER icycle REAL tau(npoints, ncol) LOGICAL box_cloudy(npoints, ncol) REAL tb(npoints, ncol) REAL ptop(npoints, ncol) REAL emcld(npoints, ncol) REAL fluxtop(npoints, ncol) REAL trans_layers_above(npoints, ncol) REAL isccp_taumin, fluxtopinit(npoints), tauir(npoints) REAL meanalbedocld(npoints) REAL albedocld(npoints, ncol) REAL boxarea INTEGER debug ! set to non-zero value to print out inputs ! ! with step debug INTEGER debugcol ! set to non-zero value to print out column ! ! decomposition with step debugcol INTEGER index1(npoints), num1, jj REAL rec2p13, tauchk CHARACTER *10 ftn09 DATA isccp_taumin/0.3/ DATA cchar/' ', '-', '1', '+', 'I', '+'/ DATA cchar_realtops/' ', ' ', '1', '1', 'I', 'I'/ tauchk = -1.*log(0.9999999) rec2p13 = 1./2.13 ncolprint = 0 ! IM ! PRINT*,' isccp_cloud_types npoints=',npoints ! if ( debug.ne.0 ) then ! j=1 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write(6,'(a10)') 'debug=' ! write(6,'(8I10)') debug ! write(6,'(a10)') 'debugcol=' ! write(6,'(8I10)') debugcol ! write(6,'(a10)') 'npoints=' ! write(6,'(8I10)') npoints ! write(6,'(a10)') 'nlev=' ! write(6,'(8I10)') nlev ! write(6,'(a10)') 'ncol=' ! write(6,'(8I10)') ncol ! write(6,'(a10)') 'top_height=' ! write(6,'(8I10)') top_height ! write(6,'(a10)') 'overlap=' ! write(6,'(8I10)') overlap ! write(6,'(a10)') 'emsfc_lw=' ! write(6,'(8f10.2)') emsfc_lw ! write(6,'(a10)') 'tautab=' ! write(6,'(8f10.2)') tautab ! write(6,'(a10)') 'invtau(1:100)=' ! write(6,'(8i10)') (invtau(i),i=1,100) ! do j=1,npoints,debug ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write(6,'(a10)') 'sunlit=' ! write(6,'(8I10)') sunlit(j) ! write(6,'(a10)') 'seed=' ! write(6,'(8I10)') seed(j) ! write(6,'(a10)') 'pfull=' ! write(6,'(8f10.2)') (pfull(j,i),i=1,nlev) ! write(6,'(a10)') 'phalf=' ! write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1) ! write(6,'(a10)') 'qv=' ! write(6,'(8f10.3)') (qv(j,i),i=1,nlev) ! write(6,'(a10)') 'cc=' ! write(6,'(8f10.3)') (cc(j,i),i=1,nlev) ! write(6,'(a10)') 'conv=' ! write(6,'(8f10.2)') (conv(j,i),i=1,nlev) ! write(6,'(a10)') 'dtau_s=' ! write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev) ! write(6,'(a10)') 'dtau_c=' ! write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev) ! write(6,'(a10)') 'skt=' ! write(6,'(8f10.2)') skt(j) ! write(6,'(a10)') 'at=' ! write(6,'(8f10.2)') (at(j,i),i=1,nlev) ! write(6,'(a10)') 'dem_s=' ! write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev) ! write(6,'(a10)') 'dem_c=' ! write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev) ! enddo ! endif ! ---------------------------------------------------! ! assign 2d tca array using 1d input array cc DO j = 1, npoints tca(j, 0) = 0 END DO DO ilev = 1, nlev DO j = 1, npoints tca(j, ilev) = cc(j, ilev) END DO END DO ! assign 2d cca array using 1d input array conv DO ilev = 1, nlev ! IM pas besoin do ibox=1,ncol DO j = 1, npoints cca(j, ilev) = conv(j, ilev) END DO ! IM enddo END DO ! IM ! do j=1, iwmx ! do l=1, 7 ! do k=1, 7 ! fq_dynreg(k,l,j) =0. ! nfq_dynreg(k,l,j) =0. ! enddo !k ! enddo !l ! enddo !j ! IM ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'seed:' ! write (6,'(I3.2)') seed(j) ! write (6,'(a)') 'tca_pp_rev:' ! write (6,'(8f5.2)') ! & ((tca(j,ilev)), ! & ilev=1,nlev) ! write (6,'(a)') 'cca_pp_rev:' ! write (6,'(8f5.2)') ! & ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev) ! enddo ! endif IF (top_height==1 .OR. top_height==3) THEN DO j = 1, npoints ptrop(j) = 5000. atmin(j) = 400. attropmin(j) = 400. atmax(j) = 0. attrop(j) = 120. itrop(j) = 1 END DO DO ilev = 1, nlev DO j = 1, npoints IF (pfull(j,ilev)<40000. .AND. pfull(j,ilev)>5000. .AND. & at(j,ilev)<attropmin(j)) THEN ptrop(j) = pfull(j, ilev) attropmin(j) = at(j, ilev) attrop(j) = attropmin(j) itrop(j) = ilev END IF IF (at(j,ilev)>atmax(j)) atmax(j) = at(j, ilev) IF (at(j,ilev)<atmin(j)) atmin(j) = at(j, ilev) END DO END DO END IF ! -----------------------------------------------------! ! ---------------------------------------------------! ! IM ! do 13 ilev=1,nlev ! num1=0 ! do j=1,npoints ! if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then ! num1=num1+1 ! index1(num1)=j ! end if ! enddo ! do jj=1,num1 ! j=index1(jj) ! write(6,*) ' error = cloud fraction less than zero' ! write(6,*) ' or ' ! write(6,*) ' error = cloud fraction greater than 1' ! write(6,*) 'value at point ',j,' is ',cc(j,ilev) ! write(6,*) 'level ',ilev ! STOP ! enddo ! num1=0 ! do j=1,npoints ! if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then ! num1=num1+1 ! index1(num1)=j ! end if ! enddo ! do jj=1,num1 ! j=index1(jj) ! write(6,*) ! & ' error = convective cloud fraction less than zero' ! write(6,*) ' or ' ! write(6,*) ! & ' error = convective cloud fraction greater than 1' ! write(6,*) 'value at point ',j,' is ',conv(j,ilev) ! write(6,*) 'level ',ilev ! STOP ! enddo ! num1=0 ! do j=1,npoints ! if (dtau_s(j,ilev) .lt. 0.) then ! num1=num1+1 ! index1(num1)=j ! end if ! enddo ! do jj=1,num1 ! j=index1(jj) ! write(6,*) ! & ' error = stratiform cloud opt. depth less than zero' ! write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev) ! write(6,*) 'level ',ilev ! STOP ! enddo ! num1=0 ! do j=1,npoints ! if (dtau_c(j,ilev) .lt. 0.) then ! num1=num1+1 ! index1(num1)=j ! end if ! enddo ! do jj=1,num1 ! j=index1(jj) ! write(6,*) ! & ' error = convective cloud opt. depth less than zero' ! write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev) ! write(6,*) 'level ',ilev ! STOP ! enddo ! num1=0 ! do j=1,npoints ! if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then ! num1=num1+1 ! index1(num1)=j ! end if ! enddo ! do jj=1,num1 ! j=index1(jj) ! write(6,*) ! & ' error = stratiform cloud emissivity less than zero' ! write(6,*)'or' ! write(6,*) ! & ' error = stratiform cloud emissivity greater than 1' ! write(6,*) 'value at point ',j,' is ',dem_s(j,ilev) ! write(6,*) 'level ',ilev ! STOP ! enddo ! num1=0 ! do j=1,npoints ! if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then ! num1=num1+1 ! index1(num1)=j ! end if ! enddo ! do jj=1,num1 ! j=index1(jj) ! write(6,*) ! & ' error = convective cloud emissivity less than zero' ! write(6,*)'or' ! write(6,*) ! & ' error = convective cloud emissivity greater than 1' ! write (6,*) ! & 'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev) ! STOP ! enddo ! 13 continue DO ibox = 1, ncol DO j = 1, npoints boxpos(j, ibox) = (ibox-.5)/ncol END DO END DO ! ---------------------------------------------------! ! Initialise working variables ! ---------------------------------------------------! ! Initialised frac_out to zero DO ilev = 1, nlev DO ibox = 1, ncol DO j = 1, npoints frac_out(j, ibox, ilev) = 0.0 END DO END DO END DO ! IM ! if (ncolprint.ne.0) then ! write (6,'(a)') 'frac_out_pp_rev:' ! do j=1,npoints,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(8f5.2)') ! & ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev) ! enddo ! write (6,'(a)') 'ncol:' ! write (6,'(I3)') ncol ! endif ! if (ncolprint.ne.0) then ! write (6,'(a)') 'last_frac_pp:' ! do j=1,npoints,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(8f5.2)') (tca(j,0)) ! enddo ! endif ! ---------------------------------------------------! ! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS ! frac_out is the array that contains the information ! where 0 is no cloud, 1 is a stratiform cloud and 2 is a ! convective cloud !loop over vertical levels DO ilev = 1, nlev ! Initialise threshold IF (ilev==1) THEN ! If max overlap IF (overlap==1) THEN ! select pixels spread evenly ! across the gridbox DO ibox = 1, ncol DO j = 1, npoints threshold(j, ibox) = boxpos(j, ibox) END DO END DO ELSE DO ibox = 1, ncol CALL ran0_vec(npoints, seed, ran) ! select random pixels from the non-convective ! part the gridbox ( some will be converted into ! convective pixels below ) DO j = 1, npoints threshold(j, ibox) = cca(j, ilev) + (1-cca(j,ilev))*ran(j) END DO END DO END IF ! IM ! IF (ncolprint.ne.0) then ! write (6,'(a)') 'threshold_nsf2:' ! do j=1,npoints,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint) ! enddo ! ENDIF END IF ! IF (ncolprint.ne.0) then ! write (6,'(a)') 'ilev:' ! write (6,'(I2)') ilev ! ENDIF DO ibox = 1, ncol ! All versions DO j = 1, npoints IF (boxpos(j,ibox)<=cca(j,ilev)) THEN ! IM REAL maxocc(j,ibox) = 1 maxocc(j, ibox) = 1.0 ELSE ! IM REAL maxocc(j,ibox) = 0 maxocc(j, ibox) = 0.0 END IF END DO ! Max overlap IF (overlap==1) THEN DO j = 1, npoints threshold_min(j, ibox) = cca(j, ilev) ! IM REAL maxosc(j,ibox)=1 maxosc(j, ibox) = 1.0 END DO END IF ! Random overlap IF (overlap==2) THEN DO j = 1, npoints threshold_min(j, ibox) = cca(j, ilev) ! IM REAL maxosc(j,ibox)=0 maxosc(j, ibox) = 0.0 END DO END IF ! Max/Random overlap IF (overlap==3) THEN DO j = 1, npoints threshold_min(j, ibox) = max(cca(j,ilev), min(tca(j,ilev-1),tca(j, & ilev))) IF (threshold(j,ibox)<min(tca(j,ilev-1),tca(j, & ilev)) .AND. (threshold(j,ibox)>cca(j,ilev))) THEN ! IM REAL maxosc(j,ibox)= 1 maxosc(j, ibox) = 1.0 ELSE ! IM REAL maxosc(j,ibox)= 0 maxosc(j, ibox) = 0.0 END IF END DO END IF ! Reset threshold CALL ran0_vec(npoints, seed, ran) DO j = 1, npoints threshold(j, ibox) = & !if max overlapped conv cloud maxocc(j, ibox)*(boxpos(j,ibox)) + & !else (1-maxocc(j,ibox))*( & !if max overlapped strat cloud (maxosc(j,ibox))*( & !threshold=boxpos threshold(j,ibox))+ & !else (1-maxosc(j,ibox))*( & !threshold_min=random[thrmin,1] threshold_min(j,ibox)+(1-threshold_min(j,ibox))*ran(j))) END DO END DO ! ibox ! Fill frac_out with 1's where tca is greater than the threshold DO ibox = 1, ncol DO j = 1, npoints IF (tca(j,ilev)>threshold(j,ibox)) THEN ! IM REAL frac_out(j,ibox,ilev)=1 frac_out(j, ibox, ilev) = 1.0 ELSE ! IM REAL frac_out(j,ibox,ilev)=0 frac_out(j, ibox, ilev) = 0.0 END IF END DO END DO ! Code to partition boxes into startiform and convective parts ! goes here DO ibox = 1, ncol DO j = 1, npoints IF (threshold(j,ibox)<=cca(j,ilev)) THEN ! = 2 IF threshold le cca(j) ! IM REAL frac_out(j,ibox,ilev) = 2 frac_out(j, ibox, ilev) = 2.0 ELSE ! = the same IF NOT threshold le cca(j) frac_out(j, ibox, ilev) = frac_out(j, ibox, ilev) END IF END DO END DO ! Set last_frac to tca at this level, so as to be tca ! from last level next time round ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints ,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'last_frac:' ! write (6,'(8f5.2)') (tca(j,ilev-1)) ! write (6,'(a)') 'cca:' ! write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint) ! write (6,'(a)') 'max_overlap_cc:' ! write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'max_overlap_sc:' ! write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'threshold_min_nsf2:' ! write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'threshold_nsf2:' ! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'frac_out_pp_rev:' ! write (6,'(8f5.2)') ! & ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev) ! enddo ! endif END DO ! ---------------------------------------------------! ! ---------------------------------------------------! ! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and ! put into vector tau !initialize tau and albedocld to zero ! loop over nlev DO ibox = 1, ncol DO j = 1, npoints tau(j, ibox) = 0. albedocld(j, ibox) = 0. boxtau(j, ibox) = 0. boxptop(j, ibox) = 0. box_cloudy(j, ibox) = .FALSE. END DO END DO !compute total cloud optical depth for each column DO ilev = 1, nlev !increment tau for each of the boxes DO ibox = 1, ncol DO j = 1, npoints ! IM REAL if (frac_out(j,ibox,ilev).eq.1) then IF (frac_out(j,ibox,ilev)==1.0) THEN tau(j, ibox) = tau(j, ibox) + dtau_s(j, ilev) END IF ! IM REAL if (frac_out(j,ibox,ilev).eq.2) then IF (frac_out(j,ibox,ilev)==2.0) THEN tau(j, ibox) = tau(j, ibox) + dtau_c(j, ilev) END IF END DO END DO ! ibox END DO ! ilev ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints ,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write(6,'(i2,1X,8(f7.2,1X))') ! & ilev, ! & (tau(j,ibox),ibox=1,ncolprint) ! enddo ! endif ! ---------------------------------------------------! ! ---------------------------------------------------! ! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES ! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE ! again this is only done if top_height = 1 or 3 ! fluxtop is the 10.5 micron radiance at the top of the ! atmosphere ! trans_layers_above is the total transmissivity in the layers ! above the current layer ! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear ! sky versions of these quantities. IF (top_height==1 .OR. top_height==3) THEN !---------------------------------------------------------------------- ! ! DO CLEAR SKY RADIANCE CALCULATION FIRST ! !compute water vapor continuum emissivity !this treatment follows Schwarkzopf and Ramasamy !JGR 1999,vol 104, pages 9467-9499. !the emissivity is calculated at a wavenumber of 955 cm-1, !or 10.47 microns wtmair = 28.9644 wtmh20 = 18.01534 navo = 6.023E+23 grav = 9.806650E+02 pstd = 1.013250E+06 t0 = 296. ! IM ! if (ncolprint .ne. 0) ! & write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv' DO ilev = 1, nlev DO j = 1, npoints !press and dpress are dyne/cm2 = Pascals *10 press(j) = pfull(j, ilev)*10. dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10 !atmden = g/cm2 = kg/m2 / 10 atmden(j) = dpress(j)/grav rvh20(j) = qv(j, ilev)*wtmair/wtmh20 wk(j) = rvh20(j)*navo*atmden(j)/wtmair rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev)) rh20s(j) = rvh20(j)*rhoave(j) rfrgn(j) = rhoave(j) - rh20s(j) tmpexp(j) = exp(-0.02*(at(j,ilev)-t0)) tauwv(j) = wk(j)*1.E-20*((0.0224697*rh20s(j)*tmpexp(j))+(3.41817E-7* & rfrgn(j)))*0.98 dem_wv(j, ilev) = 1. - exp(-1.*tauwv(j)) END DO ! IM ! if (ncolprint .ne. 0) then ! do j=1,npoints ,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write(6,'(i2,1X,3(f8.3,3X))') ilev, ! & qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.), ! & tauwv(j),dem_wv(j,ilev) ! enddo ! endif END DO !initialize variables DO j = 1, npoints fluxtop_clrsky(j) = 0. trans_layers_above_clrsky(j) = 1. END DO DO ilev = 1, nlev DO j = 1, npoints ! Black body emission at temperature of the layer bb(j) = 1/(exp(1307.27/at(j,ilev))-1.) !bb(j)= 5.67e-8*at(j,ilev)**4 ! increase TOA flux by flux emitted from layer ! times total transmittance in layers above fluxtop_clrsky(j) = fluxtop_clrsky(j) + dem_wv(j, ilev)*bb(j)* & trans_layers_above_clrsky(j) ! update trans_layers_above with transmissivity ! from this layer for next time around loop trans_layers_above_clrsky(j) = trans_layers_above_clrsky(j)* & (1.-dem_wv(j,ilev)) END DO ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints ,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'ilev:' ! write (6,'(I2)') ilev ! write (6,'(a)') ! & 'emiss_layer,100.*bb(j),100.*f,total_trans:' ! write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j), ! & 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j) ! enddo ! endif END DO !loop over level DO j = 1, npoints !add in surface emission bb(j) = 1/(exp(1307.27/skt(j))-1.) !bb(j)=5.67e-8*skt(j)**4 fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw*bb(j)* & trans_layers_above_clrsky(j) END DO ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints ,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'id:' ! write (6,'(a)') 'surface' ! write (6,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:' ! write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j), ! & 100.*fluxtop_clrsky(j), ! & trans_layers_above_clrsky(j) ! enddo ! endif ! ! END OF CLEAR SKY CALCULATION ! !---------------------------------------------------------------- ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints ,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'ts:' ! write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint) ! write (6,'(a)') 'ta_rev:' ! write (6,'(8f7.2)') ! & ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev) ! enddo ! endif !loop over columns DO ibox = 1, ncol DO j = 1, npoints fluxtop(j, ibox) = 0. trans_layers_above(j, ibox) = 1. END DO END DO DO ilev = 1, nlev DO j = 1, npoints ! Black body emission at temperature of the layer bb(j) = 1/(exp(1307.27/at(j,ilev))-1.) !bb(j)= 5.67e-8*at(j,ilev)**4 END DO DO ibox = 1, ncol DO j = 1, npoints ! emissivity for point in this layer ! IM REAL if (frac_out(j,ibox,ilev).eq.1) then IF (frac_out(j,ibox,ilev)==1.0) THEN dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_s(j,ilev))) ! IM REAL else if (frac_out(j,ibox,ilev).eq.2) then ELSE IF (frac_out(j,ibox,ilev)==2.0) THEN dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_c(j,ilev))) ELSE dem(j, ibox) = dem_wv(j, ilev) END IF ! increase TOA flux by flux emitted from layer ! times total transmittance in layers above fluxtop(j, ibox) = fluxtop(j, ibox) + dem(j, ibox)*bb(j)* & trans_layers_above(j, ibox) ! update trans_layers_above with transmissivity ! from this layer for next time around loop trans_layers_above(j, ibox) = trans_layers_above(j, ibox)* & (1.-dem(j,ibox)) END DO ! j END DO ! ibox ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints,1000 ! write (6,'(a)') 'ilev:' ! write (6,'(I2)') ilev ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'emiss_layer:' ! write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') '100.*bb(j):' ! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint) ! write (6,'(a)') '100.*f:' ! write (6,'(8f7.2)') ! & (100.*fluxtop(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'total_trans:' ! write (6,'(8f7.2)') ! & (trans_layers_above(j,ibox),ibox=1,ncolprint) ! enddo ! endif END DO ! ilev DO j = 1, npoints !add in surface emission bb(j) = 1/(exp(1307.27/skt(j))-1.) !bb(j)=5.67e-8*skt(j)**4 END DO DO ibox = 1, ncol DO j = 1, npoints !add in surface emission fluxtop(j, ibox) = fluxtop(j, ibox) + emsfc_lw*bb(j)* & trans_layers_above(j, ibox) END DO END DO ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints ,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'id:' ! write (6,'(a)') 'surface' ! write (6,'(a)') 'emiss_layer:' ! write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint) ! write (6,'(a)') '100.*bb(j):' ! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint) ! write (6,'(a)') '100.*f:' ! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) ! end do ! endif !now that you have the top of atmosphere radiance account !for ISCCP procedures to determine cloud top temperature !account for partially transmitting cloud recompute flux !ISCCP would see assuming a single layer cloud !note choice here of 2.13, as it is primarily ice !clouds which have partial emissivity and need the !adjustment performed in this section ! !If it turns out that the cloud brightness temperature !is greater than 260K, then the liquid cloud conversion !factor of 2.56 is used. ! !Note that this is discussed on pages 85-87 of !the ISCCP D level documentation (Rossow et al. 1996) DO j = 1, npoints !compute minimum brightness temperature and optical depth btcmin(j) = 1./(exp(1307.27/(attrop(j)-5.))-1.) END DO DO ibox = 1, ncol DO j = 1, npoints transmax(j) = (fluxtop(j,ibox)-btcmin(j))/(fluxtop_clrsky(j)-btcmin(j & )) !note that the initial setting of tauir(j) is needed so that !tauir(j) has a realistic value should the next if block be !bypassed tauir(j) = tau(j, ibox)*rec2p13 taumin(j) = -1.*log(max(min(transmax(j),0.9999999),0.001)) END DO IF (top_height==1) THEN DO j = 1, npoints IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN fluxtopinit(j) = fluxtop(j, ibox) tauir(j) = tau(j, ibox)*rec2p13 END IF END DO DO icycle = 1, 2 DO j = 1, npoints IF (tau(j,ibox)>(tauchk)) THEN IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN emcld(j, ibox) = 1. - exp(-1.*tauir(j)) fluxtop(j, ibox) = fluxtopinit(j) - ((1.-emcld(j, & ibox))*fluxtop_clrsky(j)) fluxtop(j, ibox) = max(1.E-06, (fluxtop(j,ibox)/emcld(j, & ibox))) tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox)))) IF (tb(j,ibox)>260.) THEN tauir(j) = tau(j, ibox)/2.56 END IF END IF END IF END DO END DO END IF DO j = 1, npoints IF (tau(j,ibox)>(tauchk)) THEN !cloudy box tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox)))) IF (top_height==1 .AND. tauir(j)<taumin(j)) THEN tb(j, ibox) = attrop(j) - 5. tau(j, ibox) = 2.13*taumin(j) END IF ELSE !clear sky brightness temperature tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j)))) END IF END DO ! j END DO ! ibox ! IM ! if (ncolprint.ne.0) then ! do j=1,npoints,1000 ! write(6,'(a10)') 'j=' ! write(6,'(8I10)') j ! write (6,'(a)') 'attrop:' ! write (6,'(8f7.2)') (attrop(j)) ! write (6,'(a)') 'btcmin:' ! write (6,'(8f7.2)') (btcmin(j)) ! write (6,'(a)') 'fluxtop_clrsky*100:' ! write (6,'(8f7.2)') ! & (100.*fluxtop_clrsky(j)) ! write (6,'(a)') '100.*f_adj:' ! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'transmax:' ! write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint) ! write (6,'(a)') 'tau:' ! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'emcld:' ! write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'total_trans:' ! write (6,'(8f7.2)') ! & (trans_layers_above(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'total_emiss:' ! write (6,'(8f7.2)') ! & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'total_trans:' ! write (6,'(8f7.2)') ! & (trans_layers_above(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'ppout:' ! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) ! enddo ! j ! endif END IF ! ---------------------------------------------------! ! ---------------------------------------------------! ! DETERMINE CLOUD TOP PRESSURE ! again the 2 methods differ according to whether ! or not you use the physical cloud top pressure (top_height = 2) ! or the radiatively determined cloud top pressure (top_height = 1 or 3) !compute cloud top pressure DO ibox = 1, ncol !segregate according to optical thickness IF (top_height==1 .OR. top_height==3) THEN !find level whose temperature !most closely matches brightness temperature DO j = 1, npoints nmatch(j) = 0 END DO DO ilev = 1, nlev - 1 ! cdir nodep DO j = 1, npoints IF ((at(j,ilev)>=tb(j,ibox) .AND. at(j,ilev+1)<tb(j, & ibox)) .OR. (at(j,ilev)<=tb(j,ibox) .AND. at(j,ilev+1)>tb(j, & ibox))) THEN nmatch(j) = nmatch(j) + 1 IF (abs(at(j,ilev)-tb(j,ibox))<abs(at(j,ilev+1)-tb(j,ibox))) THEN match(j, nmatch(j)) = ilev ELSE match(j, nmatch(j)) = ilev + 1 END IF END IF END DO END DO DO j = 1, npoints IF (nmatch(j)>=1) THEN ptop(j, ibox) = pfull(j, match(j,nmatch(j))) levmatch(j, ibox) = match(j, nmatch(j)) ELSE IF (tb(j,ibox)<atmin(j)) THEN ptop(j, ibox) = ptrop(j) levmatch(j, ibox) = itrop(j) END IF IF (tb(j,ibox)>atmax(j)) THEN ptop(j, ibox) = pfull(j, nlev) levmatch(j, ibox) = nlev END IF END IF END DO ! j ELSE ! if (top_height .eq. 1 .or. top_height .eq. 3) DO j = 1, npoints ptop(j, ibox) = 0. END DO DO ilev = 1, nlev DO j = 1, npoints IF ((ptop(j,ibox)==0.) & ! IM & ! .and.(frac_out(j,ibox,ilev) .ne. 0)) ! then .AND. (frac_out(j,ibox,ilev)/=0.0)) THEN ptop(j, ibox) = pfull(j, ilev) levmatch(j, ibox) = ilev END IF END DO END DO END IF DO j = 1, npoints IF (tau(j,ibox)<=(tauchk)) THEN ptop(j, ibox) = 0. levmatch(j, ibox) = 0 END IF END DO END DO ! ---------------------------------------------------! ! ---------------------------------------------------! ! DETERMINE ISCCP CLOUD TYPE FREQUENCIES ! Now that ptop and tau have been determined, ! determine amount of each of the 49 ISCCP cloud ! types ! Also compute grid box mean cloud top pressure and ! optical thickness. The mean cloud top pressure and ! optical thickness are averages over the cloudy ! area only. The mean cloud top pressure is a linear ! average of the cloud top pressures. The mean cloud ! optical thickness is computed by converting optical ! thickness to an albedo, averaging in albedo units, ! then converting the average albedo back to a mean ! optical thickness. !compute isccp frequencies !reset frequencies DO ilev = 1, 7 DO ilev2 = 1, 7 DO j = 1, npoints ! fq_isccp(j, ilev, ilev2) = 0. END DO END DO END DO !reset variables need for averaging cloud properties DO j = 1, npoints totalcldarea(j) = 0. meanalbedocld(j) = 0. meanptop(j) = 0. meantaucld(j) = 0. END DO ! j boxarea = 1./real(ncol) !determine optical depth category ! IM do 39 j=1,npoints ! IM do ibox=1,ncol DO ibox = 1, ncol DO j = 1, npoints ! IM ! CALL CPU_time(t1) ! IM IF (tau(j,ibox)>(tauchk) .AND. ptop(j,ibox)>0.) THEN box_cloudy(j, ibox) = .TRUE. END IF ! IM ! CALL CPU_time(t2) ! print*,'IF tau t2 - t1',t2 - t1 ! CALL CPU_time(t1) ! IM IF (box_cloudy(j,ibox)) THEN ! totalcldarea always diagnosed day or night totalcldarea(j) = totalcldarea(j) + boxarea IF (sunlit(j)==1) THEN ! tau diagnostics only with sunlight boxtau(j, ibox) = tau(j, ibox) !convert optical thickness to albedo albedocld(j, ibox) = real(invtau(min(nint(100.*tau(j,ibox)), & 45000))) !contribute to averaging meanalbedocld(j) = meanalbedocld(j) + albedocld(j, ibox)*boxarea END IF END IF ! IM ! CALL CPU_time(t2) ! print*,'IF box_cloudy t2 - t1',t2 - t1 ! CALL CPU_time(t1) ! IM BEG ! IM !convert ptop to millibars ptop(j, ibox) = ptop(j, ibox)/100. ! IM !save for output cloud top pressure and optical ! thickness boxptop(j, ibox) = ptop(j, ibox) ! IM END ! IM BEG !reset itau(j), ipres(j) itau(j) = 0 ipres(j) = 0 IF (tau(j,ibox)<isccp_taumin) THEN itau(j) = 1 ELSE IF (tau(j,ibox)>=isccp_taumin .AND. tau(j,ibox)<1.3) THEN itau(j) = 2 ELSE IF (tau(j,ibox)>=1.3 .AND. tau(j,ibox)<3.6) THEN itau(j) = 3 ELSE IF (tau(j,ibox)>=3.6 .AND. tau(j,ibox)<9.4) THEN itau(j) = 4 ELSE IF (tau(j,ibox)>=9.4 .AND. tau(j,ibox)<23.) THEN itau(j) = 5 ELSE IF (tau(j,ibox)>=23. .AND. tau(j,ibox)<60.) THEN itau(j) = 6 ELSE IF (tau(j,ibox)>=60.) THEN itau(j) = 7 END IF !determine cloud top pressure category IF (ptop(j,ibox)>0. .AND. ptop(j,ibox)<180.) THEN ipres(j) = 1 ELSE IF (ptop(j,ibox)>=180. .AND. ptop(j,ibox)<310.) THEN ipres(j) = 2 ELSE IF (ptop(j,ibox)>=310. .AND. ptop(j,ibox)<440.) THEN ipres(j) = 3 ELSE IF (ptop(j,ibox)>=440. .AND. ptop(j,ibox)<560.) THEN ipres(j) = 4 ELSE IF (ptop(j,ibox)>=560. .AND. ptop(j,ibox)<680.) THEN ipres(j) = 5 ELSE IF (ptop(j,ibox)>=680. .AND. ptop(j,ibox)<800.) THEN ipres(j) = 6 ELSE IF (ptop(j,ibox)>=800.) THEN ipres(j) = 7 END IF ! IM END IF (sunlit(j)==1 .OR. top_height==3) THEN ! IM !convert ptop to millibars ! IM ptop(j,ibox)=ptop(j,ibox) / 100. ! IM !save for output cloud top pressure and optical ! thickness ! IM boxptop(j,ibox) = ptop(j,ibox) IF (box_cloudy(j,ibox)) THEN meanptop(j) = meanptop(j) + ptop(j, ibox)*boxarea ! IM !reset itau(j), ipres(j) ! IM itau(j) = 0 ! IM ipres(j) = 0 ! if (tau(j,ibox) .lt. isccp_taumin) then ! itau(j)=1 ! else if (tau(j,ibox) .ge. isccp_taumin ! & ! & .and. tau(j,ibox) .lt. 1.3) then ! itau(j)=2 ! else if (tau(j,ibox) .ge. 1.3 ! & .and. tau(j,ibox) .lt. 3.6) then ! itau(j)=3 ! else if (tau(j,ibox) .ge. 3.6 ! & .and. tau(j,ibox) .lt. 9.4) then ! itau(j)=4 ! else if (tau(j,ibox) .ge. 9.4 ! & .and. tau(j,ibox) .lt. 23.) then ! itau(j)=5 ! else if (tau(j,ibox) .ge. 23. ! & .and. tau(j,ibox) .lt. 60.) then ! itau(j)=6 ! else if (tau(j,ibox) .ge. 60.) then ! itau(j)=7 ! end if ! !determine cloud top pressure category ! if ( ptop(j,ibox) .gt. 0. ! & .and.ptop(j,ibox) .lt. 180.) then ! ipres(j)=1 ! else if(ptop(j,ibox) .ge. 180. ! & .and.ptop(j,ibox) .lt. 310.) then ! ipres(j)=2 ! else if(ptop(j,ibox) .ge. 310. ! & .and.ptop(j,ibox) .lt. 440.) then ! ipres(j)=3 ! else if(ptop(j,ibox) .ge. 440. ! & .and.ptop(j,ibox) .lt. 560.) then ! ipres(j)=4 ! else if(ptop(j,ibox) .ge. 560. ! & .and.ptop(j,ibox) .lt. 680.) then ! ipres(j)=5 ! else if(ptop(j,ibox) .ge. 680. ! & .and.ptop(j,ibox) .lt. 800.) then ! ipres(j)=6 ! else if(ptop(j,ibox) .ge. 800.) then ! ipres(j)=7 ! end if !update frequencies IF (ipres(j)>0 .AND. itau(j)>0) THEN fq_isccp(j, itau(j), ipres(j)) = fq_isccp(j, itau(j), ipres(j)) + & boxarea END IF ! IM calcul stats regime dynamique BEG ! iw(j) = int((w(j)-wmin)/pas_w) +1 ! pctj(itau(j),ipres(j),iw(j))=.FALSE. ! !update frequencies W500 ! if (pct_ocean(j)) then ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then ! print*,' ISCCP iw=',iw(j),j ! fq_dynreg(itau(j),ipres(j),iw(j))= ! & fq_dynreg(itau(j),ipres(j),iw(j))+ ! & boxarea ! & fq_isccp(j,itau(j),ipres(j)) ! pctj(itau(j),ipres(j),iw(j))=.TRUE. ! nfq_dynreg(itau(j),ipres(j),iw(j))= ! & nfq_dynreg(itau(j),ipres(j),iw(j))+1. ! end if ! end if ! end if ! IM calcul stats regime dynamique END END IF !IM boxcloudy END IF !IM sunlit ! IM ! CALL CPU_time(t2) ! print*,'IF sunlit boxcloudy t2 - t1',t2 - t1 ! IM END DO !IM ibox/j ! IM ajout stats s/ W500 BEG ! IM ajout stats s/ W500 END ! if(j.EQ.6722) then ! print*,' ISCCP',w(j),iw(j),ipres(j),itau(j) ! endif ! if (pct_ocean(j)) then ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then ! if(pctj(itau(j),ipres(j),iw(j))) THEN ! nfq_dynreg(itau(j),ipres(j),iw(j))= ! & nfq_dynreg(itau(j),ipres(j),iw(j))+1. ! if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND. ! & iw(j).EQ.10) then ! PRINT*,' isccp AVANT', ! & nfq_dynreg(itau(j),ipres(j),iw(j)), ! & fq_dynreg(itau(j),ipres(j),iw(j)) ! endif ! endif ! endif ! endif ! endif END DO !compute mean cloud properties ! IM j/ibox DO j = 1, npoints IF (totalcldarea(j)>0.) THEN meanptop(j) = meanptop(j)/totalcldarea(j) IF (sunlit(j)==1) THEN meanalbedocld(j) = meanalbedocld(j)/totalcldarea(j) meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j))))) END IF END IF END DO ! j ! IM ajout stats s/ W500 BEG ! do nw = 1, iwmx ! do l = 1, 7 ! do k = 1, 7 ! if (nfq_dynreg(k,l,nw).GT.0.) then ! fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw) ! if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then ! print*,' isccp APRES',nfq_dynreg(k,l,nw), ! & fq_dynreg(k,l,nw) ! endif ! else ! if(fq_dynreg(k,l,nw).NE.0.) then ! print*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw) ! endif ! fq_dynreg(k,l,nw) = -1.E+20 ! nfq_dynreg(k,l,nw) = 1.E+20 ! end if ! enddo !k ! enddo !l ! enddo !nw ! IM ajout stats s/ W500 END ! ---------------------------------------------------! ! ---------------------------------------------------! ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM ! cIM ! if (debugcol.ne.0) then ! do j=1,npoints,debugcol ! !produce character output ! do ilev=1,nlev ! do ibox=1,ncol ! acc(ilev,ibox)=0 ! enddo ! enddo ! do ilev=1,nlev ! do ibox=1,ncol ! acc(ilev,ibox)=frac_out(j,ibox,ilev)*2 ! if (levmatch(j,ibox) .eq. ilev) ! & acc(ilev,ibox)=acc(ilev,ibox)+1 ! enddo ! enddo !print test ! write(ftn09,11) j ! 11 format('ftn09.',i4.4) ! open(9, FILE=ftn09, FORM='FORMATTED') ! write(9,'(a1)') ' ' ! write(9,'(10i5)') ! & (ilev,ilev=5,nlev,5) ! write(9,'(a1)') ' ' ! do ibox=1,ncol ! write(9,'(40(a1),1x,40(a1))') ! & (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev) ! & ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev) ! end do ! close(9) ! IM ! if (ncolprint.ne.0) then ! write(6,'(a1)') ' ' ! write(6,'(a2,1X,5(a7,1X),a50)') ! & 'ilev', ! & 'pfull','at', ! & 'cc*100','dem_s','dtau_s', ! & 'cchar' ! do 4012 ilev=1,nlev ! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint) ! write(6,'(i2,1X,5(f7.2,1X),50(a1))') ! & ilev, ! & pfull(j,ilev)/100.,at(j,ilev), ! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev) ! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint) ! 4012 continue ! write (6,'(a)') 'skt(j):' ! write (6,'(8f7.2)') skt(j) ! write (6,'(8I7)') (ibox,ibox=1,ncolprint) ! write (6,'(a)') 'tau:' ! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'tb:' ! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) ! write (6,'(a)') 'ptop:' ! write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint) ! endif ! enddo ! end if RETURN END SUBROUTINE isccp_cloud_types