Context Navigation

MISR_simulator.F @ 5086

Last change on this file since 5086 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

mise a jour des sources pour ISCCP, CALIPSO et PARASOL
prise en compte des changements de phases pour les nuages (Calipso)
rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

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

File size: 16.2 KB

Rev	Line
[1262]	1	!
	2	! Copyright (c) 2009, Roger Marchand, version 1.2
	3	! All rights reserved.
[2428]	4	! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $
	5	! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MISR_simulator/MISR_simulator.f $
[1262]	6	!
	7	! Redistribution and use in source and binary forms, with or without modification, are permitted
	8	! provided that the following conditions are met:
	9	!
	10	! * Redistributions of source code must retain the above copyright notice, this list of
	11	! conditions and the following disclaimer.
	12	! * Redistributions in binary form must reproduce the above copyright notice, this list
	13	! of conditions and the following disclaimer in the documentation and/or other materials
	14	! provided with the distribution.
	15	! * Neither the name of the University of Washington nor the names of its contributors may be used
	16	! to endorse or promote products derived from this software without specific prior written permission.
	17	!
	18	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
	19	! BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
	20	! SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
	21	! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
	22	! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
	23	! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	24	!
	25
	26	SUBROUTINE MISR_simulator(
	27	& npoints,
	28	& nlev,
	29	& ncol,
	30	& sunlit,
[2428]	31	& zfull,
	32	& at,
[1262]	33	& dtau_s,
[2428]	34	& dtau_c,
	35	& frac_out,
	36	& missing_value,
	37	& fq_MISR_TAU_v_CTH,
	38	& dist_model_layertops,
	39	& MISR_mean_ztop,
[1262]	40	& MISR_cldarea
	41	& )
[2428]	42
[1262]	43
	44	implicit none
	45	integer n_MISR_CTH
	46	parameter(n_MISR_CTH=16)
	47
	48	! -----
	49	! Input
	50	! -----
	51
	52	INTEGER npoints ! if ncol ==1, the number of model points in the horizontal grid
[2428]	53	! else the number of GCM grid points
	54
[1262]	55	INTEGER nlev ! number of model vertical levels
	56
	57	INTEGER ncol ! number of model sub columns
[2428]	58	! (must already be generated in via scops and passed to this
	59	! routine via the variable frac_out )
[1262]	60
	61	INTEGER sunlit(npoints) ! 1 for day points, 0 for night time
	62
[2428]	63	REAL zfull(npoints,nlev) ! height (in meters) of full model levels (i.e. midpoints)
[1262]	64	! zfull(npoints,1) is top level of model
	65	! zfull(npoints,nlev) is bottom level of model (closest point to surface)
	66
	67	REAL at(npoints,nlev) ! temperature in each model level (K)
	68
	69	REAL dtau_s(npoints,nlev) ! visible wavelength cloud optical depth ... for "stratiform" condensate
	70	! NOTE: this the cloud optical depth of only the
[2428]	71	! the model cell (i,j)
	72
[1262]	73	REAL dtau_c(npoints,nlev) ! visible wavelength cloud optical depth ... for "convective" condensate
	74	! NOTE: this the cloud optical depth of only the
[2428]	75	! the model cell (i,j)
[1262]	76
	77	REAL frac_out(npoints,ncol,nlev) ! NOTE: only need if columns>1 ... subgrid scheme in use.
[2428]	78
	79	REAL missing_value
[1262]	80
	81	! ------
	82	! Outputs
	83	! ------
[2428]	84
[1262]	85	REAL fq_MISR_TAU_v_CTH(npoints,7,n_MISR_CTH)
	86	REAL dist_model_layertops(npoints,n_MISR_CTH)
[2428]	87	REAL MISR_cldarea(npoints) ! fractional area coverged by clouds
	88	REAL MISR_mean_ztop(npoints) ! mean cloud top hieght(m) MISR would observe
	89	! NOTE: == 0 if area ==0
	90
[1262]	91
	92	! ------
	93	! Working variables
	94	! ------
	95
[2428]	96	REAL tau(npoints,ncol) ! total column optical depth ...
[1262]	97
[2428]	98	INTEGER j,ilev,ilev2,ibox,k
[1262]	99	INTEGER itau
	100
	101	LOGICAL box_cloudy(npoints,ncol)
	102
	103	real isccp_taumin
	104	real boxarea
	105	real tauchk
[2428]	106	REAL box_MISR_ztop(npoints,ncol) ! cloud top hieght(m) MISR would observe
[1262]	107
	108	integer thres_crossed_MISR
	109	integer loop,iMISR_ztop
	110
	111	real dtau, cloud_dtau, MISR_penetration_height,ztest
	112
	113	real MISR_CTH_boundaries(n_MISR_CTH+1)
	114
	115	DATA MISR_CTH_boundaries / -99, 0, 0.5, 1, 1.5, 2, 2.5, 3,
[2428]	116	c 4, 5, 7, 9, 11, 13, 15, 17, 99 /
[1262]	117
	118	DATA isccp_taumin / 0.3 /
	119
	120	tauchk = -1.*log(0.9999999)
[2428]	121
[1262]	122	!
[2428]	123	! For each GCM cell or horizontal model grid point ...
	124	!
	125	do j=1,npoints
[1262]	126
	127	!
[2428]	128	! estimate distribution of Model layer tops
	129	!
[1262]	130	dist_model_layertops(j,:)=0
	131
[2428]	132	do ilev=1,nlev
	133
	134	! define location of "layer top"
	135	if(ilev.eq.1 .or. ilev.eq.nlev) then
	136	ztest=zfull(j,ilev)
	137	else
	138	ztest=0.5*(zfull(j,ilev)+zfull(j,ilev-1))
	139	endif
[1262]	140
[2428]	141	! find MISR layer that contains this level
	142	! note, the first MISR level is "no height" level
	143	iMISR_ztop=2
	144	do loop=2,n_MISR_CTH
	145
	146	if ( ztest .gt.
	147	& 1000*MISR_CTH_boundaries(loop+1) ) then
	148
	149	iMISR_ztop=loop+1
	150	endif
	151	enddo
[1262]	152
[2428]	153	dist_model_layertops(j,iMISR_ztop)=
	154	& dist_model_layertops(j,iMISR_ztop)+1
	155	enddo
	156
	157
[1262]	158	!
	159	! compute total cloud optical depth for each column
	160	!
[2428]	161	do ibox=1,ncol
	162
	163	! Initialize tau to zero in each subcolum
	164	tau(j,ibox)=0.
	165	box_cloudy(j,ibox)=.false.
	166	box_MISR_ztop(j,ibox)=0
	167
	168	! initialize threshold detection for each sub column
	169	thres_crossed_MISR=0;
	170
	171	do ilev=1,nlev
[1262]	172
[2428]	173	dtau=0
	174
	175	if (frac_out(j,ibox,ilev).eq.1) then
[1262]	176	dtau = dtau_s(j,ilev)
	177	endif
	178
	179	if (frac_out(j,ibox,ilev).eq.2) then
	180	dtau = dtau_c(j,ilev)
[2428]	181	end if
[1262]	182
[2428]	183	tau(j,ibox)=tau(j,ibox)+ dtau
	184
	185
	186	! NOW for MISR ..
	187	! if there a cloud ... start the counter ... store this height
	188	if(thres_crossed_MISR .eq. 0 .and. dtau .gt. 0.) then
	189
	190	! first encountered a "cloud"
	191	thres_crossed_MISR=1
	192	cloud_dtau=0
	193	endif
	194
	195	if( thres_crossed_MISR .lt. 99 .and.
	196	& thres_crossed_MISR .gt. 0 ) then
[1262]	197
[2428]	198	if( dtau .eq. 0.) then
	199
	200	! we have come to the end of the current cloud
	201	! layer without yet selecting a CTH boundary.
	202	! ... restart cloud tau counter
	203	cloud_dtau=0
	204	else
	205	! add current optical depth to count for
	206	! the current cloud layer
	207	cloud_dtau=cloud_dtau+dtau
	208	endif
	209
	210	! if the cloud is continuous but optically thin (< 1)
	211	! from above the current layer cloud top to the current level
	212	! then MISR will like see a top below the top of the current
	213	! layer
	214	if( dtau.gt.0 .and. (cloud_dtau-dtau) .lt. 1) then
	215
	216	if(dtau .lt. 1 .or. ilev.eq.1 .or. ilev.eq.nlev) then
[1262]	217
[2428]	218	! MISR will likely penetrate to some point
	219	! within this layer ... the middle
	220	MISR_penetration_height=zfull(j,ilev)
[1262]	221
[2428]	222	else
	223	! take the OD = 1.0 level into this layer
	224	MISR_penetration_height=
	225	& 0.5*(zfull(j,ilev)+zfull(j,ilev-1)) -
	226	& 0.5*(zfull(j,ilev-1)-zfull(j,ilev+1))
	227	& /dtau
	228	endif
[1262]	229
[2428]	230	box_MISR_ztop(j,ibox)=MISR_penetration_height
	231
	232	endif
	233
	234	! check for a distinctive water layer
	235	if(dtau .gt. 1 .and. at(j,ilev).gt.273 ) then
[1262]	236
[2428]	237	! must be a water cloud ...
	238	! take this as CTH level
	239	thres_crossed_MISR=99
	240	endif
	241
	242	! if the total column optical depth is "large" than
	243	! MISR can't seen anything else ... set current point as CTH level
	244	if(tau(j,ibox) .gt. 5) then
[1262]	245
[2428]	246	thres_crossed_MISR=99
	247	endif
[1262]	248
[2428]	249	endif ! MISR CTH booundary not set
	250
	251	enddo !ilev - loop over vertical levesl
	252
	253	! written by roj 5/2006
	254	! check to see if there was a cloud for which we didn't
	255	! set a MISR cloud top boundary
	256	if( thres_crossed_MISR .eq. 1) then
	257
	258	! if the cloud has a total optical depth of greater
	259	! than ~ 0.5 MISR will still likely pick up this cloud
	260	! with a height near the true cloud top
	261	! otherwise there should be no CTH
	262	if( tau(j,ibox) .gt. 0.5) then
[1262]	263
[2428]	264	! keep MISR detected CTH
	265
	266	elseif(tau(j,ibox) .gt. 0.2) then
[1262]	267
[2428]	268	! MISR may detect but wont likley have a good height
	269	box_MISR_ztop(j,ibox)=-1
	270
	271	else
	272	! MISR not likely to even detect.
	273	! so set as not cloudy
	274	box_MISR_ztop(j,ibox)=0
[1262]	275
[2428]	276	endif
	277
	278	endif
	279
	280	enddo ! loop of subcolumns
[1262]	281	enddo ! loop of gridpoints
	282
	283
	284	!
[2428]	285	! Modify MISR CTH for satellite spatial / pattern matcher effects
	286	!
	287	! Code in this region added by roj 5/2006 to account
	288	! for spatial effect of the MISR pattern matcher.
	289	! Basically, if a column is found between two neighbors
	290	! at the same CTH, and that column has no hieght or
	291	! a lower CTH, THEN misr will tend to but place the
	292	! odd column at the same height as it neighbors.
	293	!
	294	! This setup assumes the columns represent a about a 1 to 4 km scale
	295	! it will need to be modified significantly, otherwise
	296	if(ncol.eq.1) then
	297
	298	! adjust based on neightboring points ... i.e. only 2D grid was input
[1262]	299	do j=2,npoints-1
[2428]	300
	301	if(box_MISR_ztop(j-1,1).gt.0 .and.
	302	& box_MISR_ztop(j+1,1).gt.0 ) then
[1262]	303
[2428]	304	if( abs( box_MISR_ztop(j-1,1) -
	305	& box_MISR_ztop(j+1,1) ) .lt. 500
	306	& .and.
	307	& box_MISR_ztop(j,1) .lt.
	308	& box_MISR_ztop(j+1,1) ) then
	309
	310	box_MISR_ztop(j,1) =
	311	& box_MISR_ztop(j+1,1)
	312	endif
[1262]	313
[2428]	314	endif
[1262]	315	enddo
[2428]	316	else
[1262]	317
	318	! adjust based on neighboring subcolumns ....
	319	do ibox=2,ncol-1
[2428]	320
	321	if(box_MISR_ztop(1,ibox-1).gt.0 .and.
	322	& box_MISR_ztop(1,ibox+1).gt.0 ) then
[1262]	323
[2428]	324	if( abs( box_MISR_ztop(1,ibox-1) -
	325	& box_MISR_ztop(1,ibox+1) ) .lt. 500
	326	& .and.
	327	& box_MISR_ztop(1,ibox) .lt.
	328	& box_MISR_ztop(1,ibox+1) ) then
	329
	330	box_MISR_ztop(1,ibox) =
	331	& box_MISR_ztop(1,ibox+1)
	332	endif
[1262]	333
[2428]	334	endif
[1262]	335	enddo
	336
[2428]	337	endif
[1262]	338
	339	!
[2428]	340	! DETERMINE CLOUD TYPE FREQUENCIES
	341	!
	342	! Now that ztop and tau have been determined,
	343	! determine amount of each cloud type
	344	boxarea=1./real(ncol)
	345	do j=1,npoints
[1262]	346
	347	! reset frequencies -- modified loop structure, roj 5/2006
[2428]	348	do ilev=1,7 ! "tau loop"
	349	do ilev2=1,n_MISR_CTH
	350	fq_MISR_TAU_v_CTH(j,ilev,ilev2)=0.
[1262]	351	enddo
[2428]	352	enddo
	353
	354	MISR_cldarea(j)=0.
	355	MISR_mean_ztop(j)=0.
[1262]	356
	357	do ibox=1,ncol
	358
	359	if (tau(j,ibox) .gt. (tauchk)) then
	360	box_cloudy(j,ibox)=.true.
	361	endif
	362
[2428]	363	itau = 0
	364
[1262]	365	if (box_cloudy(j,ibox)) then
[2428]	366
	367	!determine optical depth category
[1262]	368	if (tau(j,ibox) .lt. isccp_taumin) then
	369	itau=1
	370	else if (tau(j,ibox) .ge. isccp_taumin
	371	& .and. tau(j,ibox) .lt. 1.3) then
	372	itau=2
	373	else if (tau(j,ibox) .ge. 1.3
	374	& .and. tau(j,ibox) .lt. 3.6) then
	375	itau=3
	376	else if (tau(j,ibox) .ge. 3.6
	377	& .and. tau(j,ibox) .lt. 9.4) then
	378	itau=4
	379	else if (tau(j,ibox) .ge. 9.4
	380	& .and. tau(j,ibox) .lt. 23.) then
	381	itau=5
	382	else if (tau(j,ibox) .ge. 23.
	383	& .and. tau(j,ibox) .lt. 60.) then
	384	itau=6
	385	else if (tau(j,ibox) .ge. 60.) then
	386	itau=7
	387	endif
	388
[2428]	389	endif
[1262]	390
[2428]	391	! update MISR histograms and summary metrics - roj 5/2005
	392	if (sunlit(j).eq.1) then
	393
[1262]	394	!if cloudy added by roj 5/2005
[2428]	395	if( box_MISR_ztop(j,ibox).eq.0) then
	396
	397	! no cloud detected
	398	iMISR_ztop=0
[1262]	399
[2428]	400	elseif( box_MISR_ztop(j,ibox).eq.-1) then
[1262]	401
[2428]	402	! cloud can be detected but too thin to get CTH
	403	iMISR_ztop=1
[1262]	404
[2428]	405	fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop)=
	406	& fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop) + boxarea
[1262]	407
[2428]	408	else
	409
	410	!
	411	! determine index for MISR bin set
	412	!
[1262]	413
[2428]	414	iMISR_ztop=2
	415
	416	do loop=2,n_MISR_CTH
	417
	418	if ( box_MISR_ztop(j,ibox) .gt.
	419	& 1000*MISR_CTH_boundaries(loop+1) ) then
	420
	421	iMISR_ztop=loop+1
[1262]	422
[2428]	423	endif
	424	enddo
	425
	426	if(box_cloudy(j,ibox)) then
	427
	428	! there is an isccp clouds so itau(j) is defined
	429	fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop)=
	430	& fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop) + boxarea
[1262]	431
[2428]	432	else
	433	! MISR CTH resolution is trying to fill in a
	434	! broken cloud scene where there is no condensate.
	435	! The MISR CTH-1D-OD product will only put in a cloud
	436	! if the MISR cloud mask indicates cloud.
	437	! therefore we will not include this column in the histogram
	438	! in reality aerosoal and 3D effects or bright surfaces
	439	! could fool the MISR cloud mask
[1262]	440
[2428]	441	! the alternative is to count as very thin cloud ??
	442	! fq_MISR_TAU_v_CTH(1,iMISR_ztop)=
	443	! & fq_MISR_TAU_v_CTH(1,iMISR_ztop) + boxarea
	444	endif
[1262]	445
	446
[2428]	447	MISR_mean_ztop(j)=MISR_mean_ztop(j)+
	448	& box_MISR_ztop(j,ibox)*boxarea
[1262]	449
[2428]	450	MISR_cldarea(j)=MISR_cldarea(j) + boxarea
[1262]	451
[2428]	452	endif
	453	else
	454	! Set to issing data. A. Bodas - 14/05/2010
	455	do loop=1,n_MISR_CTH
	456	do k=1,7
	457	fq_MISR_TAU_v_CTH(j,k,loop) = missing_value
	458	enddo
	459	dist_model_layertops(j,loop) = missing_value
	460	enddo
	461	MISR_cldarea(j) = missing_value
	462	MISR_mean_ztop(npoints) = missing_value
	463
	464	endif ! is sunlight ?
	465
	466	enddo ! ibox - loop over subcolumns
[1262]	467
[2428]	468	if( MISR_cldarea(j) .gt. 0.) then
	469	MISR_mean_ztop(j)= MISR_mean_ztop(j) / MISR_cldarea(j) ! roj 5/2006
	470	endif
[1262]	471
[2428]	472	enddo ! loop over grid points
[1262]	473
	474	return
	475	end

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

Download in other formats:

Original Format