Context Navigation

isccp_cloud_types.F90

Last change on this file was 1999, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1920:1997 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 45.5 KB

Rev	Line
[1992]	1
[1344]	2	! $Id $
[524]	3
[1992]	4	SUBROUTINE isccp_cloud_types(debug, debugcol, npoints, sunlit, nlev, ncol, &
	5	seed, pfull, phalf, qv, cc, conv, dtau_s, dtau_c, top_height, overlap, &
	6	tautab, invtau, skt, emsfc_lw, at, dem_s, dem_c, fq_isccp, totalcldarea, &
	7	meanptop, meantaucld, boxtau, boxptop)
[524]	8
	9
[1992]	10	! Copyright Steve Klein and Mark Webb 2002 - all rights reserved.
[524]	11
[1992]	12	! This code is available without charge with the following conditions:
[524]	13
[1992]	14	! 1. The code is available for scientific purposes and is not for
	15	! commercial use.
	16	! 2. Any improvements you make to the code should be made available
	17	! to the to the authors for incorporation into a future release.
	18	! 3. The code should not be used in any way that brings the authors
	19	! or their employers into disrepute.
[524]	20
[1992]	21	IMPLICIT NONE
[524]	22
[1992]	23	! NOTE: the maximum number of levels and columns is set by
	24	! the following parameter statement
[524]	25
[1992]	26	INTEGER ncolprint
[524]	27
[1992]	28	! -----
	29	! Input
	30	! -----
[524]	31
[1992]	32	INTEGER npoints ! number of model points in the horizontal
	33	! PARAMETER(npoints=6722)
	34	INTEGER nlev ! number of model levels in column
	35	INTEGER ncol ! number of subcolumns
[524]	36
[1992]	37	INTEGER sunlit(npoints) ! 1 for day points, 0 for night time
[524]	38
[1992]	39	INTEGER seed(npoints) ! seed value for random number generator
	40	! ! ( see Numerical Recipes Chapter 7)
	41	! ! It is recommended that the seed is set
	42	! ! to a different value for each model
	43	! ! gridbox it is called on, as it is
	44	! ! possible that the choice of the samec
	45	! ! seed value every time may introduce some
	46	! ! statistical bias in the results, particularly
	47	! ! for low values of NCOL.
[524]	48
[1992]	49	REAL pfull(npoints, nlev) ! pressure of full model levels (Pascals)
	50	! ! pfull(npoints,1) is top level of model
	51	! ! pfull(npoints,nlev) is bottom level of model
[524]	52
[1992]	53	REAL phalf(npoints, nlev+1) ! pressure of half model levels (Pascals)
	54	! ! phalf(npoints,1) is top of model
	55	! ! phalf(npoints,nlev+1) is the surface pressure
[524]	56
[1992]	57	REAL qv(npoints, nlev) ! water vapor specific humidity (kg vapor/ kg air)
	58	! ! on full model levels
[524]	59
[1992]	60	REAL cc(npoints, nlev) ! input cloud cover in each model level (fraction)
	61	! ! NOTE: This is the HORIZONTAL area of each
	62	! ! grid box covered by clouds
[524]	63
[1992]	64	REAL conv(npoints, nlev) ! input convective cloud cover in each model level (fraction)
	65	! ! NOTE: This is the HORIZONTAL area of each
	66	! ! grid box covered by convective clouds
[524]	67
[1992]	68	REAL dtau_s(npoints, nlev) ! mean 0.67 micron optical depth of stratiform
	69	! ! clouds in each model level
	70	! ! NOTE: this the cloud optical depth of only the
	71	! ! cloudy part of the grid box, it is not weighted
	72	! ! with the 0 cloud optical depth of the clear
	73	! ! part of the grid box
[524]	74
[1992]	75	REAL dtau_c(npoints, nlev) ! mean 0.67 micron optical depth of convective
	76	! ! clouds in each
	77	! ! model level. Same note applies as in dtau_s.
[524]	78
[1992]	79	INTEGER overlap ! overlap type
[524]	80
[1992]	81	! 1=max
[524]	82
[1992]	83	! 2=rand
	84	! 3=max/rand
[524]	85
[1992]	86	INTEGER top_height ! 1 = adjust top height using both a computed
	87	! ! infrared brightness temperature and the visible
	88	! ! optical depth to adjust cloud top pressure. Note
	89	! ! that this calculation is most appropriate to compare
	90	! ! to ISCCP data during sunlit hours.
	91	! ! 2 = do not adjust top height, that is cloud top
	92	! ! pressure is the actual cloud top pressure
	93	! ! in the model
	94	! ! 3 = adjust top height using only the computed
	95	! ! infrared brightness temperature. Note that this
	96	! ! calculation is most appropriate to compare to ISCCP
	97	! ! IR only algortihm (i.e. you can compare to nighttime
	98	! ! ISCCP data with this option)
[524]	99
[1992]	100	REAL tautab(0:255) ! ISCCP table for converting count value to
	101	! ! optical thickness
[524]	102
[1992]	103	INTEGER invtau(-20:45000) ! ISCCP table for converting optical thickness
	104	! ! to count value
[524]	105
[1992]	106	! The following input variables are used only if top_height = 1 or
	107	! top_height = 3
[524]	108
[1992]	109	REAL skt(npoints) ! skin Temperature (K)
	110	REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction)
	111	REAL at(npoints, nlev) ! temperature in each model level (K)
	112	REAL dem_s(npoints, nlev) ! 10.5 micron longwave emissivity of stratiform
	113	! ! clouds in each
	114	! ! model level. Same note applies as in dtau_s.
	115	REAL dem_c(npoints, nlev) ! 10.5 micron longwave emissivity of convective
	116	! ! clouds in each
	117	! ! model level. Same note applies as in dtau_s.
	118	! IM reg.dyn BEG
	119	REAL t1, t2
	120	! REAL w(npoints) !vertical wind at 500 hPa
	121	! LOGICAL pct_ocean(npoints) !TRUE if oceanic point, FALSE otherway
	122	! INTEGER iw(npoints) , nw
	123	! REAL wmin, pas_w
	124	! INTEGER k, l, iwmx
	125	! PARAMETER(wmin=-100.,pas_w=10.,iwmx=30)
	126	! REAL fq_dynreg(7,7,iwmx)
	127	! REAL nfq_dynreg(7,7,iwmx)
	128	! LOGICAL pctj(7,7,iwmx)
	129	! IM reg.dyn END
	130	! ------
	131	! Output
	132	! ------
[524]	133
[1992]	134	REAL fq_isccp(npoints, 7, 7) ! the fraction of the model grid box covered by
	135	! ! each of the 49 ISCCP D level cloud types
[524]	136
[1992]	137	REAL totalcldarea(npoints) ! the fraction of model grid box columns
	138	! ! with cloud somewhere in them. This should
	139	! ! equal the sum over all entries of fq_isccp
[524]	140
	141
[1992]	142	! ! The following three means are averages over the cloudy areas only. If
	143	! no
	144	! ! clouds are in grid box all three quantities should equal zero.
[524]	145
[1992]	146	REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging
	147	! ! in cloud top pressure.
[524]	148
[1992]	149	REAL meantaucld(npoints) ! mean optical thickness
	150	! ! linear averaging in albedo performed.
[524]	151
[1992]	152	REAL boxtau(npoints, ncol) ! optical thickness in each column
[524]	153
[1992]	154	REAL boxptop(npoints, ncol) ! cloud top pressure (mb) in each column
[524]	155
	156
	157
[1992]	158	! ------
	159	! Working variables added when program updated to mimic Mark Webb's PV-Wave
	160	! code
	161	! ------
[524]	162
[1992]	163	REAL frac_out(npoints, ncol, nlev) ! boxes gridbox divided up into
	164	! ! Equivalent of BOX in original version, but
	165	! ! indexed by column then row, rather than
	166	! ! by row then column
[524]	167
[1992]	168	REAL tca(npoints, 0:nlev) ! total cloud cover in each model level (fraction)
	169	! ! with extra layer of zeroes on top
	170	! ! in this version this just contains the values input
	171	! ! from cc but with an extra level
	172	REAL cca(npoints, nlev) ! convective cloud cover in each model level (fraction)
	173	! ! from conv
[524]	174
[1992]	175	REAL threshold(npoints, ncol) ! pointer to position in gridbox
	176	REAL maxocc(npoints, ncol) ! Flag for max overlapped conv cld
	177	REAL maxosc(npoints, ncol) ! Flag for max overlapped strat cld
[524]	178
[1992]	179	REAL boxpos(npoints, ncol) ! ordered pointer to position in gridbox
[524]	180
[1992]	181	REAL threshold_min(npoints, ncol) ! minimum value to define range in with new threshold
	182	! ! is chosen
[524]	183
[1992]	184	REAL dem(npoints, ncol), bb(npoints) ! working variables for 10.5 micron longwave
	185	! ! emissivity in part of
	186	! ! gridbox under consideration
[524]	187
[1992]	188	REAL ran(npoints) ! vector of random numbers
	189	REAL ptrop(npoints)
	190	REAL attrop(npoints)
	191	REAL attropmin(npoints)
	192	REAL atmax(npoints)
	193	REAL atmin(npoints)
	194	REAL btcmin(npoints)
	195	REAL transmax(npoints)
[524]	196
[1992]	197	INTEGER i, j, ilev, ibox, itrop(npoints)
	198	INTEGER ipres(npoints)
	199	INTEGER itau(npoints), ilev2
	200	INTEGER acc(nlev, ncol)
	201	INTEGER match(npoints, nlev-1)
	202	INTEGER nmatch(npoints)
	203	INTEGER levmatch(npoints, ncol)
[524]	204
[1992]	205	! !variables needed for water vapor continuum absorption
	206	REAL fluxtop_clrsky(npoints), trans_layers_above_clrsky(npoints)
	207	REAL taumin(npoints)
	208	REAL dem_wv(npoints, nlev), wtmair, wtmh20, navo, grav, pstd, t0
	209	REAL press(npoints), dpress(npoints), atmden(npoints)
	210	REAL rvh20(npoints), wk(npoints), rhoave(npoints)
	211	REAL rh20s(npoints), rfrgn(npoints)
	212	REAL tmpexp(npoints), tauwv(npoints)
[524]	213
[1992]	214	CHARACTER *1 cchar(6), cchar_realtops(6)
	215	INTEGER icycle
	216	REAL tau(npoints, ncol)
	217	LOGICAL box_cloudy(npoints, ncol)
	218	REAL tb(npoints, ncol)
	219	REAL ptop(npoints, ncol)
	220	REAL emcld(npoints, ncol)
	221	REAL fluxtop(npoints, ncol)
	222	REAL trans_layers_above(npoints, ncol)
	223	REAL isccp_taumin, fluxtopinit(npoints), tauir(npoints)
	224	REAL meanalbedocld(npoints)
	225	REAL albedocld(npoints, ncol)
	226	REAL boxarea
	227	INTEGER debug ! set to non-zero value to print out inputs
	228	! ! with step debug
	229	INTEGER debugcol ! set to non-zero value to print out column
	230	! ! decomposition with step debugcol
[524]	231
[1992]	232	INTEGER index1(npoints), num1, jj
	233	REAL rec2p13, tauchk
[524]	234
[1992]	235	CHARACTER *10 ftn09
[524]	236
[1992]	237	DATA isccp_taumin/0.3/
	238	DATA cchar/' ', '-', '1', '+', 'I', '+'/
	239	DATA cchar_realtops/' ', ' ', '1', '1', 'I', 'I'/
[524]	240
[1992]	241	tauchk = -1.*log(0.9999999)
	242	rec2p13 = 1./2.13
[524]	243
[1992]	244	ncolprint = 0
[524]	245
[1992]	246	! IM
	247	! PRINT*,' isccp_cloud_types npoints=',npoints
[524]	248
[1992]	249	! if ( debug.ne.0 ) then
	250	! j=1
	251	! write(6,'(a10)') 'j='
	252	! write(6,'(8I10)') j
	253	! write(6,'(a10)') 'debug='
	254	! write(6,'(8I10)') debug
	255	! write(6,'(a10)') 'debugcol='
	256	! write(6,'(8I10)') debugcol
	257	! write(6,'(a10)') 'npoints='
	258	! write(6,'(8I10)') npoints
	259	! write(6,'(a10)') 'nlev='
	260	! write(6,'(8I10)') nlev
	261	! write(6,'(a10)') 'ncol='
	262	! write(6,'(8I10)') ncol
	263	! write(6,'(a10)') 'top_height='
	264	! write(6,'(8I10)') top_height
	265	! write(6,'(a10)') 'overlap='
	266	! write(6,'(8I10)') overlap
	267	! write(6,'(a10)') 'emsfc_lw='
	268	! write(6,'(8f10.2)') emsfc_lw
	269	! write(6,'(a10)') 'tautab='
	270	! write(6,'(8f10.2)') tautab
	271	! write(6,'(a10)') 'invtau(1:100)='
	272	! write(6,'(8i10)') (invtau(i),i=1,100)
	273	! do j=1,npoints,debug
	274	! write(6,'(a10)') 'j='
	275	! write(6,'(8I10)') j
	276	! write(6,'(a10)') 'sunlit='
	277	! write(6,'(8I10)') sunlit(j)
	278	! write(6,'(a10)') 'seed='
	279	! write(6,'(8I10)') seed(j)
	280	! write(6,'(a10)') 'pfull='
	281	! write(6,'(8f10.2)') (pfull(j,i),i=1,nlev)
	282	! write(6,'(a10)') 'phalf='
	283	! write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1)
	284	! write(6,'(a10)') 'qv='
	285	! write(6,'(8f10.3)') (qv(j,i),i=1,nlev)
	286	! write(6,'(a10)') 'cc='
	287	! write(6,'(8f10.3)') (cc(j,i),i=1,nlev)
	288	! write(6,'(a10)') 'conv='
	289	! write(6,'(8f10.2)') (conv(j,i),i=1,nlev)
	290	! write(6,'(a10)') 'dtau_s='
	291	! write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev)
	292	! write(6,'(a10)') 'dtau_c='
	293	! write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev)
	294	! write(6,'(a10)') 'skt='
	295	! write(6,'(8f10.2)') skt(j)
	296	! write(6,'(a10)') 'at='
	297	! write(6,'(8f10.2)') (at(j,i),i=1,nlev)
	298	! write(6,'(a10)') 'dem_s='
	299	! write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev)
	300	! write(6,'(a10)') 'dem_c='
	301	! write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev)
	302	! enddo
	303	! endif
[524]	304
[1992]	305	! ---------------------------------------------------!
[524]	306
[1992]	307	! assign 2d tca array using 1d input array cc
[524]	308
[1992]	309	DO j = 1, npoints
	310	tca(j, 0) = 0
	311	END DO
[524]	312
[1992]	313	DO ilev = 1, nlev
	314	DO j = 1, npoints
	315	tca(j, ilev) = cc(j, ilev)
	316	END DO
	317	END DO
[524]	318
[1992]	319	! assign 2d cca array using 1d input array conv
[524]	320
[1992]	321	DO ilev = 1, nlev
	322	! IM pas besoin do ibox=1,ncol
	323	DO j = 1, npoints
	324	cca(j, ilev) = conv(j, ilev)
	325	END DO
	326	! IM enddo
	327	END DO
[524]	328
[1992]	329	! IM
	330	! do j=1, iwmx
	331	! do l=1, 7
	332	! do k=1, 7
	333	! fq_dynreg(k,l,j) =0.
	334	! nfq_dynreg(k,l,j) =0.
	335	! enddo !k
	336	! enddo !l
	337	! enddo !j
	338	! IM
	339	! IM
	340	! if (ncolprint.ne.0) then
	341	! do j=1,npoints,1000
	342	! write(6,'(a10)') 'j='
	343	! write(6,'(8I10)') j
	344	! write (6,'(a)') 'seed:'
	345	! write (6,'(I3.2)') seed(j)
[524]	346
[1992]	347	! write (6,'(a)') 'tca_pp_rev:'
	348	! write (6,'(8f5.2)')
	349	! & ((tca(j,ilev)),
	350	! & ilev=1,nlev)
[524]	351
[1992]	352	! write (6,'(a)') 'cca_pp_rev:'
	353	! write (6,'(8f5.2)')
	354	! & ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev)
	355	! enddo
	356	! endif
[524]	357
[1992]	358	IF (top_height==1 .OR. top_height==3) THEN
[524]	359
[1992]	360	DO j = 1, npoints
	361	ptrop(j) = 5000.
	362	atmin(j) = 400.
	363	attropmin(j) = 400.
	364	atmax(j) = 0.
	365	attrop(j) = 120.
	366	itrop(j) = 1
	367	END DO
[524]	368
[1992]	369	DO ilev = 1, nlev
	370	DO j = 1, npoints
	371	IF (pfull(j,ilev)<40000. .AND. pfull(j,ilev)>5000. .AND. &
	372	at(j,ilev)<attropmin(j)) THEN
	373	ptrop(j) = pfull(j, ilev)
	374	attropmin(j) = at(j, ilev)
	375	attrop(j) = attropmin(j)
	376	itrop(j) = ilev
	377	END IF
	378	IF (at(j,ilev)>atmax(j)) atmax(j) = at(j, ilev)
	379	IF (at(j,ilev)<atmin(j)) atmin(j) = at(j, ilev)
	380	END DO
	381	END DO
[524]	382
[1992]	383	END IF
[524]	384
[1992]	385	! -----------------------------------------------------!
[524]	386
[1992]	387	! ---------------------------------------------------!
[524]	388
[1992]	389	! IM
	390	! do 13 ilev=1,nlev
	391	! num1=0
	392	! do j=1,npoints
	393	! if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then
	394	! num1=num1+1
	395	! index1(num1)=j
	396	! end if
	397	! enddo
	398	! do jj=1,num1
	399	! j=index1(jj)
	400	! write(6,*) ' error = cloud fraction less than zero'
	401	! write(6,*) ' or '
	402	! write(6,*) ' error = cloud fraction greater than 1'
	403	! write(6,*) 'value at point ',j,' is ',cc(j,ilev)
	404	! write(6,*) 'level ',ilev
	405	! STOP
	406	! enddo
	407	! num1=0
	408	! do j=1,npoints
	409	! if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then
	410	! num1=num1+1
	411	! index1(num1)=j
	412	! end if
	413	! enddo
	414	! do jj=1,num1
	415	! j=index1(jj)
	416	! write(6,*)
	417	! & ' error = convective cloud fraction less than zero'
	418	! write(6,*) ' or '
	419	! write(6,*)
	420	! & ' error = convective cloud fraction greater than 1'
	421	! write(6,*) 'value at point ',j,' is ',conv(j,ilev)
	422	! write(6,*) 'level ',ilev
	423	! STOP
	424	! enddo
[524]	425
[1992]	426	! num1=0
	427	! do j=1,npoints
	428	! if (dtau_s(j,ilev) .lt. 0.) then
	429	! num1=num1+1
	430	! index1(num1)=j
	431	! end if
	432	! enddo
	433	! do jj=1,num1
	434	! j=index1(jj)
	435	! write(6,*)
	436	! & ' error = stratiform cloud opt. depth less than zero'
	437	! write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev)
	438	! write(6,*) 'level ',ilev
	439	! STOP
	440	! enddo
	441	! num1=0
	442	! do j=1,npoints
	443	! if (dtau_c(j,ilev) .lt. 0.) then
	444	! num1=num1+1
	445	! index1(num1)=j
	446	! end if
	447	! enddo
	448	! do jj=1,num1
	449	! j=index1(jj)
	450	! write(6,*)
	451	! & ' error = convective cloud opt. depth less than zero'
	452	! write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev)
	453	! write(6,*) 'level ',ilev
	454	! STOP
	455	! enddo
[524]	456
[1992]	457	! num1=0
	458	! do j=1,npoints
	459	! if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then
	460	! num1=num1+1
	461	! index1(num1)=j
	462	! end if
	463	! enddo
	464	! do jj=1,num1
	465	! j=index1(jj)
	466	! write(6,*)
	467	! & ' error = stratiform cloud emissivity less than zero'
	468	! write(6,*)'or'
	469	! write(6,*)
	470	! & ' error = stratiform cloud emissivity greater than 1'
	471	! write(6,*) 'value at point ',j,' is ',dem_s(j,ilev)
	472	! write(6,*) 'level ',ilev
	473	! STOP
	474	! enddo
[524]	475
[1992]	476	! num1=0
	477	! do j=1,npoints
	478	! if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then
	479	! num1=num1+1
	480	! index1(num1)=j
	481	! end if
	482	! enddo
	483	! do jj=1,num1
	484	! j=index1(jj)
	485	! write(6,*)
	486	! & ' error = convective cloud emissivity less than zero'
	487	! write(6,*)'or'
	488	! write(6,*)
	489	! & ' error = convective cloud emissivity greater than 1'
	490	! write (6,*)
	491	! & 'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev)
	492	! STOP
	493	! enddo
	494	! 13 continue
[524]	495
	496
[1992]	497	DO ibox = 1, ncol
	498	DO j = 1, npoints
	499	boxpos(j, ibox) = (ibox-.5)/ncol
	500	END DO
	501	END DO
[524]	502
[1992]	503	! ---------------------------------------------------!
	504	! Initialise working variables
	505	! ---------------------------------------------------!
[524]	506
[1992]	507	! Initialised frac_out to zero
[524]	508
[1992]	509	DO ilev = 1, nlev
	510	DO ibox = 1, ncol
	511	DO j = 1, npoints
	512	frac_out(j, ibox, ilev) = 0.0
	513	END DO
	514	END DO
	515	END DO
[524]	516
[1992]	517	! IM
	518	! if (ncolprint.ne.0) then
	519	! write (6,'(a)') 'frac_out_pp_rev:'
	520	! do j=1,npoints,1000
	521	! write(6,'(a10)') 'j='
	522	! write(6,'(8I10)') j
	523	! write (6,'(8f5.2)')
	524	! & ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
[524]	525
[1992]	526	! enddo
	527	! write (6,'(a)') 'ncol:'
	528	! write (6,'(I3)') ncol
	529	! endif
	530	! if (ncolprint.ne.0) then
	531	! write (6,'(a)') 'last_frac_pp:'
	532	! do j=1,npoints,1000
	533	! write(6,'(a10)') 'j='
	534	! write(6,'(8I10)') j
	535	! write (6,'(8f5.2)') (tca(j,0))
	536	! enddo
	537	! endif
[524]	538
[1992]	539	! ---------------------------------------------------!
	540	! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
	541	! frac_out is the array that contains the information
	542	! where 0 is no cloud, 1 is a stratiform cloud and 2 is a
	543	! convective cloud
[524]	544
[1992]	545	!loop over vertical levels
	546	DO ilev = 1, nlev
[524]	547
[1992]	548	! Initialise threshold
[524]	549
[1992]	550	IF (ilev==1) THEN
	551	! If max overlap
	552	IF (overlap==1) THEN
	553	! select pixels spread evenly
	554	! across the gridbox
	555	DO ibox = 1, ncol
	556	DO j = 1, npoints
	557	threshold(j, ibox) = boxpos(j, ibox)
	558	END DO
	559	END DO
	560	ELSE
	561	DO ibox = 1, ncol
	562	CALL ran0_vec(npoints, seed, ran)
	563	! select random pixels from the non-convective
	564	! part the gridbox ( some will be converted into
	565	! convective pixels below )
	566	DO j = 1, npoints
	567	threshold(j, ibox) = cca(j, ilev) + (1-cca(j,ilev))*ran(j)
	568	END DO
	569	END DO
	570	END IF
	571	! IM
	572	! IF (ncolprint.ne.0) then
	573	! write (6,'(a)') 'threshold_nsf2:'
	574	! do j=1,npoints,1000
	575	! write(6,'(a10)') 'j='
	576	! write(6,'(8I10)') j
	577	! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
	578	! enddo
	579	! ENDIF
	580	END IF
[524]	581
[1992]	582	! IF (ncolprint.ne.0) then
	583	! write (6,'(a)') 'ilev:'
	584	! write (6,'(I2)') ilev
	585	! ENDIF
[524]	586
[1992]	587	DO ibox = 1, ncol
[524]	588
[1992]	589	! All versions
	590	DO j = 1, npoints
	591	IF (boxpos(j,ibox)<=cca(j,ilev)) THEN
	592	! IM REAL maxocc(j,ibox) = 1
	593	maxocc(j, ibox) = 1.0
	594	ELSE
	595	! IM REAL maxocc(j,ibox) = 0
	596	maxocc(j, ibox) = 0.0
	597	END IF
	598	END DO
[524]	599
[1992]	600	! Max overlap
	601	IF (overlap==1) THEN
	602	DO j = 1, npoints
	603	threshold_min(j, ibox) = cca(j, ilev)
	604	! IM REAL maxosc(j,ibox)=1
	605	maxosc(j, ibox) = 1.0
	606	END DO
	607	END IF
[524]	608
[1992]	609	! Random overlap
	610	IF (overlap==2) THEN
	611	DO j = 1, npoints
	612	threshold_min(j, ibox) = cca(j, ilev)
	613	! IM REAL maxosc(j,ibox)=0
	614	maxosc(j, ibox) = 0.0
	615	END DO
	616	END IF
[524]	617
[1992]	618	! Max/Random overlap
	619	IF (overlap==3) THEN
	620	DO j = 1, npoints
	621	threshold_min(j, ibox) = max(cca(j,ilev), min(tca(j,ilev-1),tca(j, &
	622	ilev)))
	623	IF (threshold(j,ibox)<min(tca(j,ilev-1),tca(j, &
	624	ilev)) .AND. (threshold(j,ibox)>cca(j,ilev))) THEN
	625	! IM REAL maxosc(j,ibox)= 1
	626	maxosc(j, ibox) = 1.0
	627	ELSE
	628	! IM REAL maxosc(j,ibox)= 0
	629	maxosc(j, ibox) = 0.0
	630	END IF
	631	END DO
	632	END IF
[524]	633
[1992]	634	! Reset threshold
	635	CALL ran0_vec(npoints, seed, ran)
	636
	637	DO j = 1, npoints
	638	threshold(j, ibox) = & !if max overlapped conv cloud
	639	maxocc(j, ibox)*(boxpos(j,ibox)) + & !else
	640	(1-maxocc(j,ibox))*( & !if max overlapped strat cloud
	641	(maxosc(j,ibox))*( & !threshold=boxpos
	642	threshold(j,ibox))+ & !else
	643	(1-maxosc(j,ibox))*( & !threshold_min=random[thrmin,1]
	644	threshold_min(j,ibox)+(1-threshold_min(j,ibox))*ran(j)))
	645	END DO
	646
	647	END DO ! ibox
	648
	649	! Fill frac_out with 1's where tca is greater than the threshold
	650
	651	DO ibox = 1, ncol
	652	DO j = 1, npoints
	653	IF (tca(j,ilev)>threshold(j,ibox)) THEN
	654	! IM REAL frac_out(j,ibox,ilev)=1
	655	frac_out(j, ibox, ilev) = 1.0
	656	ELSE
	657	! IM REAL frac_out(j,ibox,ilev)=0
	658	frac_out(j, ibox, ilev) = 0.0
	659	END IF
	660	END DO
	661	END DO
	662
	663	! Code to partition boxes into startiform and convective parts
	664	! goes here
	665
	666	DO ibox = 1, ncol
	667	DO j = 1, npoints
	668	IF (threshold(j,ibox)<=cca(j,ilev)) THEN
	669	! = 2 IF threshold le cca(j)
	670	! IM REAL frac_out(j,ibox,ilev) = 2
	671	frac_out(j, ibox, ilev) = 2.0
	672	ELSE
	673	! = the same IF NOT threshold le cca(j)
	674	frac_out(j, ibox, ilev) = frac_out(j, ibox, ilev)
	675	END IF
	676	END DO
	677	END DO
	678
	679	! Set last_frac to tca at this level, so as to be tca
	680	! from last level next time round
	681
	682	! IM
	683	! if (ncolprint.ne.0) then
	684
	685	! do j=1,npoints ,1000
	686	! write(6,'(a10)') 'j='
	687	! write(6,'(8I10)') j
	688	! write (6,'(a)') 'last_frac:'
	689	! write (6,'(8f5.2)') (tca(j,ilev-1))
	690
	691	! write (6,'(a)') 'cca:'
	692	! write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint)
	693
	694	! write (6,'(a)') 'max_overlap_cc:'
	695	! write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
	696
	697	! write (6,'(a)') 'max_overlap_sc:'
	698	! write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
	699
	700	! write (6,'(a)') 'threshold_min_nsf2:'
	701	! write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
	702
	703	! write (6,'(a)') 'threshold_nsf2:'
	704	! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
	705
	706	! write (6,'(a)') 'frac_out_pp_rev:'
	707	! write (6,'(8f5.2)')
	708	! & ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
	709	! enddo
	710	! endif
	711
	712
	713	END DO
	714
	715	! ---------------------------------------------------!
	716
	717
	718
	719	! ---------------------------------------------------!
	720	! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
	721	! put into vector tau
	722
	723	!initialize tau and albedocld to zero
	724	! loop over nlev
	725	DO ibox = 1, ncol
	726	DO j = 1, npoints
	727	tau(j, ibox) = 0.
	728	albedocld(j, ibox) = 0.
	729	boxtau(j, ibox) = 0.
	730	boxptop(j, ibox) = 0.
	731	box_cloudy(j, ibox) = .FALSE.
	732	END DO
	733	END DO
	734
	735	!compute total cloud optical depth for each column
	736	DO ilev = 1, nlev
	737	!increment tau for each of the boxes
	738	DO ibox = 1, ncol
	739	DO j = 1, npoints
	740	! IM REAL if (frac_out(j,ibox,ilev).eq.1) then
	741	IF (frac_out(j,ibox,ilev)==1.0) THEN
	742	tau(j, ibox) = tau(j, ibox) + dtau_s(j, ilev)
	743	END IF
	744	! IM REAL if (frac_out(j,ibox,ilev).eq.2) then
	745	IF (frac_out(j,ibox,ilev)==2.0) THEN
	746	tau(j, ibox) = tau(j, ibox) + dtau_c(j, ilev)
	747	END IF
	748	END DO
	749	END DO ! ibox
	750	END DO ! ilev
	751	! IM
	752	! if (ncolprint.ne.0) then
	753
	754	! do j=1,npoints ,1000
	755	! write(6,'(a10)') 'j='
	756	! write(6,'(8I10)') j
	757	! write(6,'(i2,1X,8(f7.2,1X))')
	758	! & ilev,
	759	! & (tau(j,ibox),ibox=1,ncolprint)
	760	! enddo
	761	! endif
	762
	763	! ---------------------------------------------------!
	764
	765
	766
	767
	768	! ---------------------------------------------------!
	769	! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
	770	! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
	771
	772	! again this is only done if top_height = 1 or 3
	773
	774	! fluxtop is the 10.5 micron radiance at the top of the
	775	! atmosphere
	776	! trans_layers_above is the total transmissivity in the layers
	777	! above the current layer
	778	! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
	779	! sky versions of these quantities.
	780
	781	IF (top_height==1 .OR. top_height==3) THEN
	782
	783
	784	!----------------------------------------------------------------------
	785	!
	786	! DO CLEAR SKY RADIANCE CALCULATION FIRST
	787	!
	788	!compute water vapor continuum emissivity
	789	!this treatment follows Schwarkzopf and Ramasamy
	790	!JGR 1999,vol 104, pages 9467-9499.
	791	!the emissivity is calculated at a wavenumber of 955 cm-1,
	792	!or 10.47 microns
	793	wtmair = 28.9644
	794	wtmh20 = 18.01534
	795	navo = 6.023E+23
	796	grav = 9.806650E+02
	797	pstd = 1.013250E+06
	798	t0 = 296.
	799	! IM
	800	! if (ncolprint .ne. 0)
	801	! & write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv'
	802	DO ilev = 1, nlev
	803	DO j = 1, npoints
	804	!press and dpress are dyne/cm2 = Pascals *10
	805	press(j) = pfull(j, ilev)*10.
	806	dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
	807	!atmden = g/cm2 = kg/m2 / 10
	808	atmden(j) = dpress(j)/grav
	809	rvh20(j) = qv(j, ilev)*wtmair/wtmh20
	810	wk(j) = rvh20(j)navoatmden(j)/wtmair
	811	rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
	812	rh20s(j) = rvh20(j)*rhoave(j)
	813	rfrgn(j) = rhoave(j) - rh20s(j)
	814	tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
	815	tauwv(j) = wk(j)1.E-20((0.0224697rh20s(j)tmpexp(j))+(3.41817E-7* &
	816	rfrgn(j)))*0.98
	817	dem_wv(j, ilev) = 1. - exp(-1.*tauwv(j))
	818	END DO
	819	! IM
	820	! if (ncolprint .ne. 0) then
	821	! do j=1,npoints ,1000
	822	! write(6,'(a10)') 'j='
	823	! write(6,'(8I10)') j
	824	! write(6,'(i2,1X,3(f8.3,3X))') ilev,
	825	! & qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
	826	! & tauwv(j),dem_wv(j,ilev)
	827	! enddo
	828	! endif
	829	END DO
	830
	831	!initialize variables
	832	DO j = 1, npoints
	833	fluxtop_clrsky(j) = 0.
	834	trans_layers_above_clrsky(j) = 1.
	835	END DO
	836
	837	DO ilev = 1, nlev
	838	DO j = 1, npoints
	839
	840	! Black body emission at temperature of the layer
	841
	842	bb(j) = 1/(exp(1307.27/at(j,ilev))-1.)
	843	!bb(j)= 5.67e-8at(j,ilev)*4
	844
	845	! increase TOA flux by flux emitted from layer
	846	! times total transmittance in layers above
	847
	848	fluxtop_clrsky(j) = fluxtop_clrsky(j) + dem_wv(j, ilev)bb(j) &
	849	trans_layers_above_clrsky(j)
	850
	851	! update trans_layers_above with transmissivity
	852	! from this layer for next time around loop
	853
	854	trans_layers_above_clrsky(j) = trans_layers_above_clrsky(j)* &
	855	(1.-dem_wv(j,ilev))
	856
	857
	858	END DO
	859	! IM
	860	! if (ncolprint.ne.0) then
	861	! do j=1,npoints ,1000
	862	! write(6,'(a10)') 'j='
	863	! write(6,'(8I10)') j
	864	! write (6,'(a)') 'ilev:'
	865	! write (6,'(I2)') ilev
	866
	867	! write (6,'(a)')
	868	! & 'emiss_layer,100.bb(j),100.f,total_trans:'
	869	! write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
	870	! & 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
	871	! enddo
	872	! endif
	873
	874	END DO !loop over level
	875
	876	DO j = 1, npoints
	877	!add in surface emission
	878	bb(j) = 1/(exp(1307.27/skt(j))-1.)
	879	!bb(j)=5.67e-8skt(j)*4
	880
	881	fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lwbb(j) &
	882	trans_layers_above_clrsky(j)
	883	END DO
	884
	885	! IM
	886	! if (ncolprint.ne.0) then
	887	! do j=1,npoints ,1000
	888	! write(6,'(a10)') 'j='
	889	! write(6,'(8I10)') j
	890	! write (6,'(a)') 'id:'
	891	! write (6,'(a)') 'surface'
	892
	893	! write (6,'(a)') 'emsfc,100.bb(j),100.f,total_trans:'
	894	! write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j),
	895	! & 100.*fluxtop_clrsky(j),
	896	! & trans_layers_above_clrsky(j)
	897	! enddo
	898	! endif
	899
	900
	901	!
	902	! END OF CLEAR SKY CALCULATION
	903	!
	904	!----------------------------------------------------------------
	905
	906
	907	! IM
	908	! if (ncolprint.ne.0) then
	909
	910	! do j=1,npoints ,1000
	911	! write(6,'(a10)') 'j='
	912	! write(6,'(8I10)') j
	913	! write (6,'(a)') 'ts:'
	914	! write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
	915
	916	! write (6,'(a)') 'ta_rev:'
	917	! write (6,'(8f7.2)')
	918	! & ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
	919
	920	! enddo
	921	! endif
	922	!loop over columns
	923	DO ibox = 1, ncol
	924	DO j = 1, npoints
	925	fluxtop(j, ibox) = 0.
	926	trans_layers_above(j, ibox) = 1.
	927	END DO
	928	END DO
	929
	930	DO ilev = 1, nlev
	931	DO j = 1, npoints
	932	! Black body emission at temperature of the layer
	933
	934	bb(j) = 1/(exp(1307.27/at(j,ilev))-1.)
	935	!bb(j)= 5.67e-8at(j,ilev)*4
	936	END DO
	937
	938	DO ibox = 1, ncol
	939	DO j = 1, npoints
	940
	941	! emissivity for point in this layer
	942	! IM REAL if (frac_out(j,ibox,ilev).eq.1) then
	943	IF (frac_out(j,ibox,ilev)==1.0) THEN
	944	dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_s(j,ilev)))
	945	! IM REAL else if (frac_out(j,ibox,ilev).eq.2) then
	946	ELSE IF (frac_out(j,ibox,ilev)==2.0) THEN
	947	dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_c(j,ilev)))
	948	ELSE
	949	dem(j, ibox) = dem_wv(j, ilev)
	950	END IF
	951
	952
	953	! increase TOA flux by flux emitted from layer
	954	! times total transmittance in layers above
	955
	956	fluxtop(j, ibox) = fluxtop(j, ibox) + dem(j, ibox)bb(j) &
	957	trans_layers_above(j, ibox)
	958
	959	! update trans_layers_above with transmissivity
	960	! from this layer for next time around loop
	961
	962	trans_layers_above(j, ibox) = trans_layers_above(j, ibox)* &
	963	(1.-dem(j,ibox))
	964
	965	END DO ! j
	966	END DO ! ibox
	967
	968	! IM
	969	! if (ncolprint.ne.0) then
	970	! do j=1,npoints,1000
	971	! write (6,'(a)') 'ilev:'
	972	! write (6,'(I2)') ilev
	973
	974	! write(6,'(a10)') 'j='
	975	! write(6,'(8I10)') j
	976	! write (6,'(a)') 'emiss_layer:'
	977	! write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
	978
	979	! write (6,'(a)') '100.*bb(j):'
	980	! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
	981
	982	! write (6,'(a)') '100.*f:'
	983	! write (6,'(8f7.2)')
	984	! & (100.*fluxtop(j,ibox),ibox=1,ncolprint)
	985
	986	! write (6,'(a)') 'total_trans:'
	987	! write (6,'(8f7.2)')
	988	! & (trans_layers_above(j,ibox),ibox=1,ncolprint)
	989	! enddo
	990	! endif
	991
	992	END DO ! ilev
	993
	994
	995	DO j = 1, npoints
	996	!add in surface emission
	997	bb(j) = 1/(exp(1307.27/skt(j))-1.)
	998	!bb(j)=5.67e-8skt(j)*4
	999	END DO
	1000
	1001	DO ibox = 1, ncol
	1002	DO j = 1, npoints
	1003
[524]	1004	!add in surface emission
	1005
[1992]	1006	fluxtop(j, ibox) = fluxtop(j, ibox) + emsfc_lwbb(j) &
	1007	trans_layers_above(j, ibox)
[524]	1008
[1992]	1009	END DO
	1010	END DO
[524]	1011
[1992]	1012	! IM
	1013	! if (ncolprint.ne.0) then
[524]	1014
[1992]	1015	! do j=1,npoints ,1000
	1016	! write(6,'(a10)') 'j='
	1017	! write(6,'(8I10)') j
	1018	! write (6,'(a)') 'id:'
	1019	! write (6,'(a)') 'surface'
[524]	1020
[1992]	1021	! write (6,'(a)') 'emiss_layer:'
	1022	! write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
[524]	1023
[1992]	1024	! write (6,'(a)') '100.*bb(j):'
	1025	! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
[524]	1026
[1992]	1027	! write (6,'(a)') '100.*f:'
	1028	! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
	1029	! end do
	1030	! endif
[524]	1031
[1992]	1032	!now that you have the top of atmosphere radiance account
	1033	!for ISCCP procedures to determine cloud top temperature
[524]	1034
[1992]	1035	!account for partially transmitting cloud recompute flux
	1036	!ISCCP would see assuming a single layer cloud
	1037	!note choice here of 2.13, as it is primarily ice
	1038	!clouds which have partial emissivity and need the
	1039	!adjustment performed in this section
	1040	!
	1041	!If it turns out that the cloud brightness temperature
	1042	!is greater than 260K, then the liquid cloud conversion
	1043	!factor of 2.56 is used.
	1044	!
	1045	!Note that this is discussed on pages 85-87 of
	1046	!the ISCCP D level documentation (Rossow et al. 1996)
[524]	1047
[1992]	1048	DO j = 1, npoints
	1049	!compute minimum brightness temperature and optical depth
	1050	btcmin(j) = 1./(exp(1307.27/(attrop(j)-5.))-1.)
	1051	END DO
	1052	DO ibox = 1, ncol
	1053	DO j = 1, npoints
	1054	transmax(j) = (fluxtop(j,ibox)-btcmin(j))/(fluxtop_clrsky(j)-btcmin(j &
	1055	))
	1056	!note that the initial setting of tauir(j) is needed so that
	1057	!tauir(j) has a realistic value should the next if block be
	1058	!bypassed
	1059	tauir(j) = tau(j, ibox)*rec2p13
	1060	taumin(j) = -1.*log(max(min(transmax(j),0.9999999),0.001))
[524]	1061
[1992]	1062	END DO
[524]	1063
[1992]	1064	IF (top_height==1) THEN
	1065	DO j = 1, npoints
	1066	IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN
	1067	fluxtopinit(j) = fluxtop(j, ibox)
	1068	tauir(j) = tau(j, ibox)*rec2p13
	1069	END IF
	1070	END DO
	1071	DO icycle = 1, 2
	1072	DO j = 1, npoints
	1073	IF (tau(j,ibox)>(tauchk)) THEN
	1074	IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN
	1075	emcld(j, ibox) = 1. - exp(-1.*tauir(j))
	1076	fluxtop(j, ibox) = fluxtopinit(j) - ((1.-emcld(j, &
	1077	ibox))*fluxtop_clrsky(j))
	1078	fluxtop(j, ibox) = max(1.E-06, (fluxtop(j,ibox)/emcld(j, &
	1079	ibox)))
	1080	tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox))))
	1081	IF (tb(j,ibox)>260.) THEN
	1082	tauir(j) = tau(j, ibox)/2.56
	1083	END IF
	1084	END IF
	1085	END IF
	1086	END DO
	1087	END DO
[524]	1088
[1992]	1089	END IF
[524]	1090
[1992]	1091	DO j = 1, npoints
	1092	IF (tau(j,ibox)>(tauchk)) THEN
	1093	!cloudy box
	1094	tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox))))
	1095	IF (top_height==1 .AND. tauir(j)<taumin(j)) THEN
	1096	tb(j, ibox) = attrop(j) - 5.
	1097	tau(j, ibox) = 2.13*taumin(j)
	1098	END IF
	1099	ELSE
	1100	!clear sky brightness temperature
	1101	tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
	1102	END IF
	1103	END DO ! j
	1104	END DO ! ibox
[524]	1105
[1992]	1106	! IM
	1107	! if (ncolprint.ne.0) then
[524]	1108
[1992]	1109	! do j=1,npoints,1000
	1110	! write(6,'(a10)') 'j='
	1111	! write(6,'(8I10)') j
[524]	1112
[1992]	1113	! write (6,'(a)') 'attrop:'
	1114	! write (6,'(8f7.2)') (attrop(j))
[524]	1115
[1992]	1116	! write (6,'(a)') 'btcmin:'
	1117	! write (6,'(8f7.2)') (btcmin(j))
[524]	1118
[1992]	1119	! write (6,'(a)') 'fluxtop_clrsky*100:'
	1120	! write (6,'(8f7.2)')
	1121	! & (100.*fluxtop_clrsky(j))
[524]	1122
[1992]	1123	! write (6,'(a)') '100.*f_adj:'
	1124	! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
[524]	1125
[1992]	1126	! write (6,'(a)') 'transmax:'
	1127	! write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
[524]	1128
[1992]	1129	! write (6,'(a)') 'tau:'
	1130	! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
[524]	1131
[1992]	1132	! write (6,'(a)') 'emcld:'
	1133	! write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
[524]	1134
[1992]	1135	! write (6,'(a)') 'total_trans:'
	1136	! write (6,'(8f7.2)')
	1137	! & (trans_layers_above(j,ibox),ibox=1,ncolprint)
[524]	1138
[1992]	1139	! write (6,'(a)') 'total_emiss:'
	1140	! write (6,'(8f7.2)')
	1141	! & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
[524]	1142
[1992]	1143	! write (6,'(a)') 'total_trans:'
	1144	! write (6,'(8f7.2)')
	1145	! & (trans_layers_above(j,ibox),ibox=1,ncolprint)
[524]	1146
[1992]	1147	! write (6,'(a)') 'ppout:'
	1148	! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
	1149	! enddo ! j
	1150	! endif
[524]	1151
[1992]	1152	END IF
[524]	1153
[1992]	1154	! ---------------------------------------------------!
[524]	1155
	1156
[1992]	1157	! ---------------------------------------------------!
	1158	! DETERMINE CLOUD TOP PRESSURE
[524]	1159
[1992]	1160	! again the 2 methods differ according to whether
	1161	! or not you use the physical cloud top pressure (top_height = 2)
	1162	! or the radiatively determined cloud top pressure (top_height = 1 or 3)
[524]	1163
	1164
[1992]	1165	!compute cloud top pressure
	1166	DO ibox = 1, ncol
	1167	!segregate according to optical thickness
	1168	IF (top_height==1 .OR. top_height==3) THEN
	1169	!find level whose temperature
	1170	!most closely matches brightness temperature
	1171	DO j = 1, npoints
	1172	nmatch(j) = 0
	1173	END DO
	1174	DO ilev = 1, nlev - 1
	1175	! cdir nodep
	1176	DO j = 1, npoints
	1177	IF ((at(j,ilev)>=tb(j,ibox) .AND. at(j,ilev+1)<tb(j, &
	1178	ibox)) .OR. (at(j,ilev)<=tb(j,ibox) .AND. at(j,ilev+1)>tb(j, &
	1179	ibox))) THEN
[524]	1180
[1992]	1181	nmatch(j) = nmatch(j) + 1
	1182	IF (abs(at(j,ilev)-tb(j,ibox))<abs(at(j,ilev+1)-tb(j,ibox))) THEN
	1183	match(j, nmatch(j)) = ilev
	1184	ELSE
	1185	match(j, nmatch(j)) = ilev + 1
	1186	END IF
	1187	END IF
	1188	END DO
	1189	END DO
[524]	1190
[1992]	1191	DO j = 1, npoints
	1192	IF (nmatch(j)>=1) THEN
	1193	ptop(j, ibox) = pfull(j, match(j,nmatch(j)))
	1194	levmatch(j, ibox) = match(j, nmatch(j))
	1195	ELSE
	1196	IF (tb(j,ibox)<atmin(j)) THEN
	1197	ptop(j, ibox) = ptrop(j)
	1198	levmatch(j, ibox) = itrop(j)
	1199	END IF
	1200	IF (tb(j,ibox)>atmax(j)) THEN
	1201	ptop(j, ibox) = pfull(j, nlev)
	1202	levmatch(j, ibox) = nlev
	1203	END IF
	1204	END IF
	1205	END DO ! j
[524]	1206
[1992]	1207	ELSE ! if (top_height .eq. 1 .or. top_height .eq. 3)
[524]	1208
[1992]	1209	DO j = 1, npoints
	1210	ptop(j, ibox) = 0.
	1211	END DO
	1212	DO ilev = 1, nlev
	1213	DO j = 1, npoints
	1214	IF ((ptop(j,ibox)==0.) & ! IM &
	1215	! .and.(frac_out(j,ibox,ilev) .ne. 0))
	1216	! then
	1217	.AND. (frac_out(j,ibox,ilev)/=0.0)) THEN
	1218	ptop(j, ibox) = pfull(j, ilev)
	1219	levmatch(j, ibox) = ilev
	1220	END IF
	1221	END DO
	1222	END DO
	1223	END IF
[524]	1224
[1992]	1225	DO j = 1, npoints
	1226	IF (tau(j,ibox)<=(tauchk)) THEN
	1227	ptop(j, ibox) = 0.
	1228	levmatch(j, ibox) = 0
	1229	END IF
	1230	END DO
[524]	1231
[1992]	1232	END DO
[524]	1233
	1234
	1235
[1992]	1236	! ---------------------------------------------------!
[524]	1237
	1238
	1239
[1992]	1240	! ---------------------------------------------------!
	1241	! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
[524]	1242
[1992]	1243	! Now that ptop and tau have been determined,
	1244	! determine amount of each of the 49 ISCCP cloud
	1245	! types
[524]	1246
[1992]	1247	! Also compute grid box mean cloud top pressure and
	1248	! optical thickness. The mean cloud top pressure and
	1249	! optical thickness are averages over the cloudy
	1250	! area only. The mean cloud top pressure is a linear
	1251	! average of the cloud top pressures. The mean cloud
	1252	! optical thickness is computed by converting optical
	1253	! thickness to an albedo, averaging in albedo units,
	1254	! then converting the average albedo back to a mean
	1255	! optical thickness.
[524]	1256
	1257
[1992]	1258	!compute isccp frequencies
[524]	1259
[1992]	1260	!reset frequencies
	1261	DO ilev = 1, 7
	1262	DO ilev2 = 1, 7
	1263	DO j = 1, npoints !
	1264	fq_isccp(j, ilev, ilev2) = 0.
	1265	END DO
	1266	END DO
	1267	END DO
[524]	1268
[1992]	1269	!reset variables need for averaging cloud properties
	1270	DO j = 1, npoints
	1271	totalcldarea(j) = 0.
	1272	meanalbedocld(j) = 0.
	1273	meanptop(j) = 0.
	1274	meantaucld(j) = 0.
	1275	END DO ! j
[524]	1276
[1992]	1277	boxarea = 1./real(ncol)
[524]	1278
[1992]	1279	!determine optical depth category
	1280	! IM do 39 j=1,npoints
	1281	! IM do ibox=1,ncol
	1282	DO ibox = 1, ncol
	1283	DO j = 1, npoints
[524]	1284
[1992]	1285	! IM
	1286	! CALL CPU_time(t1)
	1287	! IM
[524]	1288
[1992]	1289	IF (tau(j,ibox)>(tauchk) .AND. ptop(j,ibox)>0.) THEN
	1290	box_cloudy(j, ibox) = .TRUE.
	1291	END IF
[524]	1292
[1992]	1293	! IM
	1294	! CALL CPU_time(t2)
	1295	! print*,'IF tau t2 - t1',t2 - t1
[524]	1296
[1992]	1297	! CALL CPU_time(t1)
	1298	! IM
[524]	1299
[1992]	1300	IF (box_cloudy(j,ibox)) THEN
[524]	1301
[1992]	1302	! totalcldarea always diagnosed day or night
	1303	totalcldarea(j) = totalcldarea(j) + boxarea
[524]	1304
[1992]	1305	IF (sunlit(j)==1) THEN
[524]	1306
[1992]	1307	! tau diagnostics only with sunlight
[524]	1308
[1992]	1309	boxtau(j, ibox) = tau(j, ibox)
[524]	1310
[1992]	1311	!convert optical thickness to albedo
	1312	albedocld(j, ibox) = real(invtau(min(nint(100.*tau(j,ibox)), &
	1313	45000)))
[524]	1314
[1992]	1315	!contribute to averaging
	1316	meanalbedocld(j) = meanalbedocld(j) + albedocld(j, ibox)*boxarea
[524]	1317
[1992]	1318	END IF
[524]	1319
[1992]	1320	END IF
[524]	1321
[1992]	1322	! IM
	1323	! CALL CPU_time(t2)
	1324	! print*,'IF box_cloudy t2 - t1',t2 - t1
[524]	1325
[1992]	1326	! CALL CPU_time(t1)
	1327	! IM BEG
	1328	! IM !convert ptop to millibars
	1329	ptop(j, ibox) = ptop(j, ibox)/100.
[524]	1330
[1992]	1331	! IM !save for output cloud top pressure and optical
	1332	! thickness
	1333	boxptop(j, ibox) = ptop(j, ibox)
	1334	! IM END
[524]	1335
[1992]	1336	! IM BEG
	1337	!reset itau(j), ipres(j)
	1338	itau(j) = 0
	1339	ipres(j) = 0
[524]	1340
[1992]	1341	IF (tau(j,ibox)<isccp_taumin) THEN
	1342	itau(j) = 1
	1343	ELSE IF (tau(j,ibox)>=isccp_taumin .AND. tau(j,ibox)<1.3) THEN
	1344	itau(j) = 2
	1345	ELSE IF (tau(j,ibox)>=1.3 .AND. tau(j,ibox)<3.6) THEN
	1346	itau(j) = 3
	1347	ELSE IF (tau(j,ibox)>=3.6 .AND. tau(j,ibox)<9.4) THEN
	1348	itau(j) = 4
	1349	ELSE IF (tau(j,ibox)>=9.4 .AND. tau(j,ibox)<23.) THEN
	1350	itau(j) = 5
	1351	ELSE IF (tau(j,ibox)>=23. .AND. tau(j,ibox)<60.) THEN
	1352	itau(j) = 6
	1353	ELSE IF (tau(j,ibox)>=60.) THEN
	1354	itau(j) = 7
	1355	END IF
[524]	1356
[1992]	1357	!determine cloud top pressure category
	1358	IF (ptop(j,ibox)>0. .AND. ptop(j,ibox)<180.) THEN
	1359	ipres(j) = 1
	1360	ELSE IF (ptop(j,ibox)>=180. .AND. ptop(j,ibox)<310.) THEN
	1361	ipres(j) = 2
	1362	ELSE IF (ptop(j,ibox)>=310. .AND. ptop(j,ibox)<440.) THEN
	1363	ipres(j) = 3
	1364	ELSE IF (ptop(j,ibox)>=440. .AND. ptop(j,ibox)<560.) THEN
	1365	ipres(j) = 4
	1366	ELSE IF (ptop(j,ibox)>=560. .AND. ptop(j,ibox)<680.) THEN
	1367	ipres(j) = 5
	1368	ELSE IF (ptop(j,ibox)>=680. .AND. ptop(j,ibox)<800.) THEN
	1369	ipres(j) = 6
	1370	ELSE IF (ptop(j,ibox)>=800.) THEN
	1371	ipres(j) = 7
	1372	END IF
	1373	! IM END
[524]	1374
[1992]	1375	IF (sunlit(j)==1 .OR. top_height==3) THEN
[524]	1376
[1992]	1377	! IM !convert ptop to millibars
	1378	! IM ptop(j,ibox)=ptop(j,ibox) / 100.
[524]	1379
[1992]	1380	! IM !save for output cloud top pressure and optical
	1381	! thickness
	1382	! IM boxptop(j,ibox) = ptop(j,ibox)
[524]	1383
[1992]	1384	IF (box_cloudy(j,ibox)) THEN
[524]	1385
[1992]	1386	meanptop(j) = meanptop(j) + ptop(j, ibox)*boxarea
[524]	1387
[1992]	1388	! IM !reset itau(j), ipres(j)
	1389	! IM itau(j) = 0
	1390	! IM ipres(j) = 0
[524]	1391
[1992]	1392	! if (tau(j,ibox) .lt. isccp_taumin) then
	1393	! itau(j)=1
	1394	! else if (tau(j,ibox) .ge. isccp_taumin
	1395	! &
	1396	! & .and. tau(j,ibox) .lt. 1.3) then
	1397	! itau(j)=2
	1398	! else if (tau(j,ibox) .ge. 1.3
	1399	! & .and. tau(j,ibox) .lt. 3.6) then
	1400	! itau(j)=3
	1401	! else if (tau(j,ibox) .ge. 3.6
	1402	! & .and. tau(j,ibox) .lt. 9.4) then
	1403	! itau(j)=4
	1404	! else if (tau(j,ibox) .ge. 9.4
	1405	! & .and. tau(j,ibox) .lt. 23.) then
	1406	! itau(j)=5
	1407	! else if (tau(j,ibox) .ge. 23.
	1408	! & .and. tau(j,ibox) .lt. 60.) then
	1409	! itau(j)=6
	1410	! else if (tau(j,ibox) .ge. 60.) then
	1411	! itau(j)=7
	1412	! end if
[524]	1413
[1992]	1414	! !determine cloud top pressure category
	1415	! if ( ptop(j,ibox) .gt. 0.
	1416	! & .and.ptop(j,ibox) .lt. 180.) then
	1417	! ipres(j)=1
	1418	! else if(ptop(j,ibox) .ge. 180.
	1419	! & .and.ptop(j,ibox) .lt. 310.) then
	1420	! ipres(j)=2
	1421	! else if(ptop(j,ibox) .ge. 310.
	1422	! & .and.ptop(j,ibox) .lt. 440.) then
	1423	! ipres(j)=3
	1424	! else if(ptop(j,ibox) .ge. 440.
	1425	! & .and.ptop(j,ibox) .lt. 560.) then
	1426	! ipres(j)=4
	1427	! else if(ptop(j,ibox) .ge. 560.
	1428	! & .and.ptop(j,ibox) .lt. 680.) then
	1429	! ipres(j)=5
	1430	! else if(ptop(j,ibox) .ge. 680.
	1431	! & .and.ptop(j,ibox) .lt. 800.) then
	1432	! ipres(j)=6
	1433	! else if(ptop(j,ibox) .ge. 800.) then
	1434	! ipres(j)=7
	1435	! end if
	1436
	1437	!update frequencies
	1438	IF (ipres(j)>0 .AND. itau(j)>0) THEN
	1439	fq_isccp(j, itau(j), ipres(j)) = fq_isccp(j, itau(j), ipres(j)) + &
	1440	boxarea
	1441	END IF
	1442
	1443	! IM calcul stats regime dynamique BEG
	1444	! iw(j) = int((w(j)-wmin)/pas_w) +1
	1445	! pctj(itau(j),ipres(j),iw(j))=.FALSE.
	1446	! !update frequencies W500
	1447	! if (pct_ocean(j)) then
	1448	! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
	1449	! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
	1450	! print*,' ISCCP iw=',iw(j),j
	1451	! fq_dynreg(itau(j),ipres(j),iw(j))=
	1452	! & fq_dynreg(itau(j),ipres(j),iw(j))+
	1453	! & boxarea
	1454	! & fq_isccp(j,itau(j),ipres(j))
	1455	! pctj(itau(j),ipres(j),iw(j))=.TRUE.
	1456	! nfq_dynreg(itau(j),ipres(j),iw(j))=
	1457	! & nfq_dynreg(itau(j),ipres(j),iw(j))+1.
	1458	! end if
	1459	! end if
	1460	! end if
	1461	! IM calcul stats regime dynamique END
	1462	END IF !IM boxcloudy
	1463
	1464	END IF !IM sunlit
	1465
	1466	! IM
	1467	! CALL CPU_time(t2)
	1468	! print*,'IF sunlit boxcloudy t2 - t1',t2 - t1
	1469	! IM
	1470	END DO !IM ibox/j
	1471
	1472
	1473	! IM ajout stats s/ W500 BEG
	1474	! IM ajout stats s/ W500 END
	1475
	1476	! if(j.EQ.6722) then
	1477	! print*,' ISCCP',w(j),iw(j),ipres(j),itau(j)
	1478	! endif
	1479
	1480	! if (pct_ocean(j)) then
	1481	! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
	1482	! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
	1483	! if(pctj(itau(j),ipres(j),iw(j))) THEN
	1484	! nfq_dynreg(itau(j),ipres(j),iw(j))=
	1485	! & nfq_dynreg(itau(j),ipres(j),iw(j))+1.
	1486	! if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND.
	1487	! & iw(j).EQ.10) then
	1488	! PRINT*,' isccp AVANT',
	1489	! & nfq_dynreg(itau(j),ipres(j),iw(j)),
	1490	! & fq_dynreg(itau(j),ipres(j),iw(j))
	1491	! endif
	1492	! endif
	1493	! endif
	1494	! endif
	1495	! endif
	1496
	1497	END DO
	1498	!compute mean cloud properties
	1499	! IM j/ibox
	1500	DO j = 1, npoints
	1501	IF (totalcldarea(j)>0.) THEN
	1502	meanptop(j) = meanptop(j)/totalcldarea(j)
	1503	IF (sunlit(j)==1) THEN
	1504	meanalbedocld(j) = meanalbedocld(j)/totalcldarea(j)
	1505	meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j)))))
	1506	END IF
	1507	END IF
	1508	END DO ! j
	1509
	1510	! IM ajout stats s/ W500 BEG
	1511	! do nw = 1, iwmx
	1512	! do l = 1, 7
	1513	! do k = 1, 7
	1514	! if (nfq_dynreg(k,l,nw).GT.0.) then
	1515	! fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw)
	1516	! if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then
	1517	! print*,' isccp APRES',nfq_dynreg(k,l,nw),
	1518	! & fq_dynreg(k,l,nw)
	1519	! endif
	1520	! else
	1521	! if(fq_dynreg(k,l,nw).NE.0.) then
	1522	! print*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw)
	1523	! endif
	1524	! fq_dynreg(k,l,nw) = -1.E+20
	1525	! nfq_dynreg(k,l,nw) = 1.E+20
	1526	! end if
	1527	! enddo !k
	1528	! enddo !l
	1529	! enddo !nw
	1530	! IM ajout stats s/ W500 END
	1531	! ---------------------------------------------------!
	1532
	1533	! ---------------------------------------------------!
	1534	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
	1535
	1536	! cIM
	1537	! if (debugcol.ne.0) then
	1538
	1539	! do j=1,npoints,debugcol
	1540
	1541	! !produce character output
	1542	! do ilev=1,nlev
	1543	! do ibox=1,ncol
	1544	! acc(ilev,ibox)=0
	1545	! enddo
	1546	! enddo
	1547
	1548	! do ilev=1,nlev
	1549	! do ibox=1,ncol
	1550	! acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
	1551	! if (levmatch(j,ibox) .eq. ilev)
	1552	! & acc(ilev,ibox)=acc(ilev,ibox)+1
	1553	! enddo
	1554	! enddo
	1555
	1556	!print test
	1557
	1558	! write(ftn09,11) j
	1559	! 11 format('ftn09.',i4.4)
	1560	! open(9, FILE=ftn09, FORM='FORMATTED')
	1561
	1562	! write(9,'(a1)') ' '
	1563	! write(9,'(10i5)')
	1564	! & (ilev,ilev=5,nlev,5)
	1565	! write(9,'(a1)') ' '
	1566
	1567	! do ibox=1,ncol
	1568	! write(9,'(40(a1),1x,40(a1))')
	1569	! & (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
	1570	! & ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
	1571	! end do
	1572	! close(9)
	1573
	1574	! IM
	1575	! if (ncolprint.ne.0) then
	1576	! write(6,'(a1)') ' '
	1577	! write(6,'(a2,1X,5(a7,1X),a50)')
	1578	! & 'ilev',
	1579	! & 'pfull','at',
	1580	! & 'cc*100','dem_s','dtau_s',
	1581	! & 'cchar'
	1582
	1583	! do 4012 ilev=1,nlev
	1584	! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
	1585	! write(6,'(i2,1X,5(f7.2,1X),50(a1))')
	1586	! & ilev,
	1587	! & pfull(j,ilev)/100.,at(j,ilev),
	1588	! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
	1589	! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
	1590	! 4012 continue
	1591	! write (6,'(a)') 'skt(j):'
	1592	! write (6,'(8f7.2)') skt(j)
	1593
	1594	! write (6,'(8I7)') (ibox,ibox=1,ncolprint)
	1595
	1596	! write (6,'(a)') 'tau:'
	1597	! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
	1598
	1599	! write (6,'(a)') 'tb:'
	1600	! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
	1601
	1602	! write (6,'(a)') 'ptop:'
	1603	! write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
	1604	! endif
	1605
	1606	! enddo
	1607
	1608	! end if
	1609
	1610	RETURN
	1611	END SUBROUTINE isccp_cloud_types

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format