Context Navigation

surfini.F @ 1242

Last change on this file since 1242 was 1226, checked in by aslmd, 11 years ago
LMDZ.MARS : Replaced comcstfi and planete includes by modules.
File size: 20.1 KB

Rev	Line
[38]	1	SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
[520]	2	! to use 'getin'
	3	USE ioipsl_getincom
[697]	4	use netcdf
[1224]	5	use tracer_mod, only: nqmx, noms
[1047]	6	use comgeomfi_h, only: long, lati
	7	use surfdat_h, only: watercaptag, frost_albedo_threshold,
	8	& albedo_h2o_ice, inert_h2o_ice, albedodat,
[1224]	9	& albedice, dryness
[1212]	10	#ifndef MESOSCALE
[1130]	11	use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
	12	use mod_phys_lmdz_para, only : is_master, gather, scatter
[1212]	13	#endif
[1226]	14	USE comcstfi_h
[38]	15	IMPLICIT NONE
	16	c=======================================================================
	17	c
	18	c creation des calottes pour l'etat initial
	19	c
	20	c=======================================================================
	21	c-----------------------------------------------------------------------
	22	c Declarations:
	23	c -------------
	24	#include "dimensions.h"
[1047]	25	!#include "dimphys.h"
	26	!#include "surfdat.h"
[38]	27	#include "callkeys.h"
[1036]	28	!#include "tracer.h"
[1047]	29	!#include "comgeomfi.h"
[669]	30
	31	#include "datafile.h"
	32
[1130]	33	integer,intent(in) :: ngrid ! number of atmospheric columns
	34	real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
	35	real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
	36	real,intent(out) :: psolaralb(ngrid,2) ! albedo
	37
	38	INTEGER ig,icap,iq,alternate
[520]	39	REAL icedryness ! ice dryness
[633]	40
	41	! longwatercaptag is watercaptag. Trick for some compilers
	42	LOGICAL, DIMENSION(100000) :: longwatercaptag
[520]	43
[669]	44	! There are 3 different modes for ice distribution:
	45	! icelocationmode = 1 ---> based on data from surface.nc
	46	! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
	47	! icelocationmode = 3 ---> based on logical relations for latitude and longitude
	48	! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
	49	INTEGER,SAVE :: icelocationmode = 2
	50
	51
	52	!in case icelocationmode == 1
	53	INTEGER i,j
	54	INTEGER imd,jmd
	55	PARAMETER (imd=360,jmd=180)
[712]	56	REAL zdata(imd,jmd)
[669]	57	REAL zelat,zelon
	58
[1212]	59	#ifndef MESOSCALE
[1207]	60	INTEGER nb_ice(klon_glo,2) ! number of counts \| detected ice for GCM grid
[1212]	61	#endif
[669]	62	INTEGER latice(jjm,2),lonice (iim,2) ! number of counts \| detected ice along lat & lon axis
	63
	64	REAL step,count,ratiolat
	65
	66	INTEGER ierr,nid,nvarid
	67
	68	REAL,SAVE :: min_icevalue = 500.
[697]	69	character(len=50) :: string = 'thermal'
[669]	70
	71	character (len=100) :: zedatafile
[1130]	72
[1212]	73	#ifdef MESOSCALE
	74
	75	do ig=1,ngrid
	76
	77	!write(,) "all qsurf to zero. dirty."
	78	do iq=1,nqmx
	79	qsurf(ig,iq)=0. !! on jette les inputs GCM
	80	!! on regle juste watercaptag
	81	!! il faudrait garder les inputs GCM
	82	!! si elles sont consequentes
	83	enddo
	84	if ( ( lati(ig)*180./pi .gt. 70. ) .and.
	85	. ( albedodat(ig) .ge. 0.26 ) ) then
	86	write(,)"outlier ",ig
	87	watercaptag(ig) = .true.
	88	dryness(ig) = 1.
	89	albedodat(ig) = albedo_h2o_ice !! pour output
	90	else
	91	watercaptag(ig) = .false.
	92	dryness(ig) = 1.
	93	endif
	94
	95	enddo
	96	#endif
	97	! problem with nested precompiling flags
	98
	99	#ifndef MESOSCALE
[1130]	100	! to handle parallel cases
	101	#if CPP_PARA
	102	logical watercaptag_glo(klon_glo)
	103	real dryness_glo(klon_glo)
	104	real lati_glo(klon_glo)
	105	real long_glo(klon_glo)
	106	#else
	107	logical watercaptag_glo(ngrid)
	108	real dryness_glo(ngrid)
	109	real lati_glo(ngrid)
	110	real long_glo(ngrid)
	111	#endif
[1212]	112	#endif
[1130]	113
[1212]	114	#ifndef MESOSCALE
	115
[38]	116	c
	117	c=======================================================================
[1130]	118	! Initialize watercaptag (default is false)
	119	watercaptag_glo(:)=.false.
[38]	120
[283]	121	c water ice outliers
[38]	122	c ------------------------------------------
	123
[283]	124	IF ((water) .and. (caps)) THEN
[285]	125
[283]	126	c Perennial H20 north cap defined by watercaptag=true (allows surface to be
	127	c hollowed by sublimation in vdifc).
[38]	128
[283]	129	c We might not want albedodat to be modified because it is used to write
	130	c restart files. Instead, albedo is directly modified when needed (i.e.
	131	c if we have watercaptag and no co2 ice), below and in albedocaps.F90
	132
	133	c "Dryness coefficient" controlling the evaporation and
	134	c sublimation from the ground water ice (close to 1)
	135	c HERE, the goal is to correct for the fact
	136	c that the simulated permanent water ice polar caps
	137	c is larger than the actual cap and the atmospheric
	138	c opacity not always realistic.
	139
	140	alternate = 0
[520]	141
[1047]	142	if (ngrid .ne. 1) then
[633]	143	watercaptag(:) = .false.
	144	longwatercaptag(:) = .false.
	145	endif
	146
[1130]	147	write(,) "surfini: Ice dryness ?"
[520]	148	icedryness=1. ! default value
	149	call getin("icedryness",icedryness)
[1130]	150	write(,) "surfini: icedryness = ",icedryness
[669]	151	dryness (:) = icedryness
[633]	152
[1130]	153	! To be able to run in parallel, we work on the full grid
	154	! and dispatch results afterwards
[669]	155
[1130]	156	! start by geting latitudes and logitudes on full grid
	157	! (in serial mode, this is just a copy)
	158	call gather(lati,lati_glo)
	159	call gather(long,long_glo)
[669]	160
[1130]	161	if (is_master) then
	162
	163	IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
[669]	164
[1047]	165	print*, 'ngrid = 1, do no put ice caps in surfini.F'
[669]	166
[1130]	167	ELSE IF (icelocationmode .eq. 1) THEN
[669]	168
	169	print*,'Surfini: ice caps defined from surface.nc'
	170
	171	! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
	172	! Typically, it is for thermal inertia above 500 tiu.
	173	! Two conditions are verified:
	174	! 1. GCM ice caps are defined such as area is conserved for a given latitude
	175	! (the approximation is that all points within the GCM latitude resolution have the same area).
	176	! 2. caps are placed to fill the GCM points with the most detected ice first.
	177
	178
	179
[1130]	180	zedatafile = trim(datafile)
[669]	181
	182
[1130]	183	ierr=nf90_open(trim(zedatafile)//'/surface.nc',
	184	& NF90_NOWRITE,nid)
[669]	185
[1130]	186	IF (ierr.NE.nf90_noerr) THEN
[669]	187	write(,)'Error : cannot open file surface.nc '
	188	write(,)'(in phymars/surfini.F)'
	189	write(,)'It should be in :',trim(zedatafile),'/'
	190	write(,)'1) You can set this path in the callphys.def file:'
	191	write(,)' datadir=/path/to/the/datafiles'
	192	write(,)'2) If necessary, surface.nc (and other datafiles)'
	193	write(,)' can be obtained online on:'
	194	write(,)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
	195	CALL ABORT
[1130]	196	ENDIF
[669]	197
	198
[1130]	199	ierr=nf90_inq_varid(nid, string, nvarid)
	200	if (ierr.ne.nf90_noerr) then
	201	write(,) 'surfini error, cannot find ',trim(string)
	202	write(,) ' in file ',trim(zedatafile),'/surface.nc'
	203	write(,)trim(nf90_strerror(ierr))
	204	stop
	205	endif
[697]	206
[1130]	207	ierr=nf90_get_var(nid, nvarid, zdata)
[697]	208
[1130]	209	if (ierr.ne.nf90_noerr) then
	210	write(,) 'surfini: error failed loading ',trim(string)
	211	write(,)trim(nf90_strerror(ierr))
	212	stop
	213	endif
[669]	214
	215
[1130]	216	ierr=nf90_close(nid)
[669]	217
[285]	218
[1130]	219	nb_ice(:,1) = 1 ! default: there is no ice
	220	latice(:,1) = 1
	221	lonice(:,1) = 1
	222	nb_ice(:,2) = 0
	223	latice(:,2) = 0
	224	lonice(:,2) = 0
	225	!print*,'jjm,iim',jjm,iim ! jjm = nb lati , iim = nb longi
[669]	226
[1130]	227	! loop over the GCM grid - except for poles (ig=1 and ngrid)
	228	do ig=2,klon_glo-1
[520]	229
[1130]	230	! loop over the surface file grid
	231	do i=1,imd
	232	do j=1,jmd
	233	zelon = i - 180.
	234	zelat = 90. - j
	235	if ((abs(lati_glo(ig)*180./pi-zelat).le.90./real(jjm)) .and.
	236	& (abs(long_glo(ig)*180./pi-zelon).le.180./real(iim))) then
[669]	237	! count all points in that GCM grid point
	238	nb_ice(ig,1) = nb_ice(ig,1) + 1
[712]	239	if (zdata(i,j) > min_icevalue)
[669]	240	! count all detected points in that GCM grid point
	241	& nb_ice(ig,2) = nb_ice(ig,2) + 1
[712]	242	endif
[1130]	243	enddo
	244	enddo
[712]	245
[669]	246	! projection of nb_ice on GCM lat and lon axes
[1130]	247	latice(1+(ig-2)/iim,:) =
[669]	248	& latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
[1130]	249	lonice(1+mod(ig-2,iim),:) =
[669]	250	& lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
	251
[1130]	252	enddo ! of do ig=2,klon_glo-1
[669]	253
	254
	255
[1130]	256	! special case for poles
	257	nb_ice(1,2) = 1 ! ice prescribed on north pole
	258	latice(1,:) = nb_ice(1,:)
	259	lonice(1,:) = nb_ice(1,:)
	260	latice(jjm,:) = nb_ice(ngrid,:)
	261	lonice(iim,:) = nb_ice(ngrid,:)
[669]	262
	263
[712]	264	! print*, 'latice TOT', latice(:,1)
	265	! print*, 'latice FOUND', latice(:,2)
	266	! print*, 'lonice TOT', lonice(:,1)
	267	! print*, 'lonice FOUND', lonice(:,2)
[669]	268
[712]	269	! print, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))iim)
	270	! print, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))jjm)
[669]	271
[712]	272	! print*,''
	273	! print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
	274	! print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
[669]	275
	276
[1130]	277	! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
	278	do i=1,jjm/2
	279	step = 1. ! threshold to add ice cap
	280	count = 0. ! number of ice GCM caps at this latitude
	281	! ratiolat is the ratio of area covered by ice within this GCM latitude range
	282	ratiolat = real(latice(i,2))/real(latice(i,1))
	283	!print,'i',i,(i-1)iim+2,i*iim+1
[669]	284
[1130]	285	! put ice caps while there is not enough ice,
	286	! as long as the threshold is above 20%
	287	do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
	288	count = 0.
	289	! loop over GCM longitudes
	290	do j=1,iim
[669]	291	! if the detected ice ratio in the GCM grid point
	292	! is more than 'step', then add ice
	293	if (real(nb_ice((i-1)*iim+1+j,2))
	294	& / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
[1130]	295	watercaptag_glo((i-1)*iim+1+j) = .true.
[669]	296	count = count + 1
	297	endif
[1130]	298	enddo ! of do j=1,iim
	299	!print, 'step',step,count,ratiolatiim
	300	step = step - 0.01
	301	enddo ! of do while
[669]	302	!print, 'step',step,count,ratiolatiim
	303
[1130]	304	enddo ! of do i=1,jjm/2
[669]	305
	306
[1130]	307	ELSE IF (icelocationmode .eq. 2) THEN
[669]	308
[1130]	309	print*,'Surfini: predefined ice caps'
[669]	310
[1130]	311	if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
[633]	312
[669]	313	print*,'water ice caps distribution for 32x24 resolution'
[633]	314	longwatercaptag(1:9) = .true. ! central cap - core
	315	longwatercaptag(26:33) = .true. ! central cap
[669]	316	longwatercaptag(1:33) = .true. ! central cap
	317	longwatercaptag(56) = .true. ! central cap
	318	longwatercaptag(58) = .true. ! central cap
	319	longwatercaptag(60) = .true. ! central cap
	320	longwatercaptag(62) = .true. ! central cap
	321	longwatercaptag(64) = .true. ! central cap
	322	!--------------------- OUTLIERS ----------------------------
[633]	323
[1130]	324	else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
[633]	325
[669]	326	print*,'water ice caps distribution for 64x48 resolution'
[633]	327	longwatercaptag(1:65) = .true. ! central cap - core
	328	longwatercaptag(75:85) = .true. ! central cap
	329	longwatercaptag(93:114) = .true. ! central cap
[669]	330	!--------------------- OUTLIERS ----------------------------
	331	if (.true.) then
[633]	332	longwatercaptag(136) = .true. ! outlier, lat = 78.75
	333	longwatercaptag(138) = .true. ! outlier, lat = 78.75
	334	longwatercaptag(140) = .true. ! outlier, lat = 78.75
	335	longwatercaptag(142) = .true. ! outlier, lat = 78.75
	336	longwatercaptag(161) = .true. ! outlier, lat = 78.75
	337	longwatercaptag(163) = .true. ! outlier, lat = 78.75
	338	longwatercaptag(165) = .true. ! outlier, lat = 78.75
	339	longwatercaptag(183) = .true. ! outlier, lat = 78.75
	340	longwatercaptag(185) = .true. ! outlier, lat = 78.75
	341	longwatercaptag(187) = .true. ! outlier, lat = 78.75
	342	longwatercaptag(189) = .true. ! outlier, lat = 78.75
	343	longwatercaptag(191) = .true. ! outlier, lat = 78.75
	344	longwatercaptag(193) = .true. ! outlier, lat = 78.75
	345	longwatercaptag(194) = .true. ! outlier, lat = 75
	346	longwatercaptag(203) = .true. ! outlier, lat = 75
	347	longwatercaptag(207) = .true. ! outlier, lat = 75
	348	longwatercaptag(244) = .true. ! outlier, lat = 75
	349	longwatercaptag(246) = .true. ! outlier, lat = 75
	350	longwatercaptag(250) = .true. ! outlier, lat = 75
	351	longwatercaptag(252) = .true. ! outlier, lat = 75
	352	longwatercaptag(254) = .true. ! outlier, lat = 75
[740]	353	longwatercaptag(256) = .true. ! outlier, lat = 75
[669]	354	endif
	355	!--------------------------------------------------------------
[633]	356
[669]	357
[520]	358
[1130]	359	else if (klon_glo .ne. 1) then
[669]	360
[1130]	361	print*,'No predefined ice location for this resolution :',
	362	& iim,jjm
	363	print*,'Please change icelocationmode in surfini.F'
	364	print*,'Or add some new definitions ...'
	365	call abort
[520]	366
[1130]	367	endif
[669]	368
[1130]	369	do ig=1,klon_glo
	370	if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
	371	enddo
[283]	372
[633]	373
[1130]	374	ELSE IF (icelocationmode .eq. 3) THEN
[669]	375
[1130]	376	print*,'Surfini: ice caps defined by lat and lon values'
[520]	377
[1130]	378	do ig=1,klon_glo
[669]	379
	380	c-------- Towards olympia planitia water caps -----------
	381	c--------------------------------------------------------
[520]	382
[1130]	383	if ( ( ( lati_glo(ig)*180./pi .ge. 77. ) .and. ! cap #2
	384	. ( lati_glo(ig)*180./pi .le. 80. ) .and.
	385	. ( long_glo(ig)*180./pi .ge. 110. ) .and.
	386	. ( long_glo(ig)*180./pi .le. 181. ) )
[520]	387	. .or.
[283]	388
[1130]	389	. ( ( lati_glo(ig)*180./pi .ge. 75. ) .and. ! cap #4 (Korolev crater)
	390	. ( lati_glo(ig)*180./pi .le. 76. ) .and.
	391	. ( long_glo(ig)*180./pi .ge. 150. ) .and.
	392	. ( long_glo(ig)*180./pi .le. 168. ) )
[520]	393	. .or.
[1130]	394	. ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
	395	. ( lati_glo(ig)*180./pi .le. 80. ) .and.
	396	. ( long_glo(ig)*180./pi .ge. -150.) .and.
	397	. ( long_glo(ig)*180./pi .le. -110.) ) )
[520]	398	. then
	399
	400	if ((alternate .eq. 0)) then ! 1/2 en 64x48 sinon trop large en lat
[669]	401	! watercaptag(ig)=.true.
[520]	402	alternate = 1
[669]	403	else
[520]	404	alternate = 0
[669]	405	endif !end if alternate = 0
	406
[1130]	407	endif
[633]	408
[669]	409	c----------- Opposite olympia planitia water cap --------
	410	c--------------------------------------------------------
[520]	411
[1130]	412	if ( ( ( lati_glo(ig)*180./pi .ge. 80 ) .and.
	413	. ( lati_glo(ig)*180./pi .le. 84 ) )
[669]	414	. .and.
[1130]	415	. ( ( long_glo(ig)*180./pi .lt. -95. ) .or. !!! 32x24
	416	. ( long_glo(ig)*180./pi .gt. 85. ) ) ) then !!! 32x24
	417	! . ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and. !!! 64x48
	418	! . ( long_glo(ig)*180./pi .le. 90. ) ) .or. !!! 64x48
	419	! . ( ( long_glo(ig)*180./pi .ge. -77. ) .and. !!! 64x48
	420	! . ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then !!! 64x48
	421	! watercaptag_glo(ig)=.true.
	422	endif
[520]	423
	424
[669]	425	c -------------------- Central cap ----------------------
[283]	426	c--------------------------------------------------------
	427
[1130]	428	if (abs(lati_glo(ig)*180./pi).gt.80)
	429	. watercaptag_glo(ig)=.true.
[669]	430
	431	c--------------------------------------------------------
	432	c--------------------------------------------------------
[1130]	433	end do ! of (klon_glo)
[669]	434
	435
[1130]	436	ELSE
[669]	437
	438	print*, 'In surfini.F, icelocationmode is ', icelocationmode
	439	print*, 'It should be 1, 2 or 3.'
	440	call abort
	441
[1130]	442	ENDIF ! of if (icelocation)
[669]	443
	444
	445	! print caps locations - useful for plots too
[1130]	446	print*,'surfini: latitude \| longitude \| ig'
	447	do ig=1,klon_glo
	448	dryness_glo(ig) = icedryness
[669]	449
[1130]	450	if (watercaptag_glo(ig)) then
	451	print,'surfini: ice water cap', lati_glo(ig)180./pi,
	452	& long_glo(ig)*180./pi, ig
[669]	453	endif
	454	enddo
	455
[1130]	456	endif !of if (is_master)
	457
	458	! Now scatter fields watercaptag and dryness from master to all
	459	! (is just a plain copy in serial mode)
	460	call scatter(dryness_glo,dryness)
	461	call scatter(watercaptag_glo,watercaptag)
	462
	463	! end of #else of #ifdef MESOSCALE
[520]	464	ENDIF ! (caps & water)
[1212]	465	#endif
[520]	466
[283]	467	c ===============================================================
	468	c INITIAL ALBEDO
	469	c ===============================================================
	470
	471	write(,)"surfini: water frost thickness",
	472	s frost_albedo_threshold
	473	write(,)"surfini: water ice albedo:", albedo_h2o_ice
	474	write(,)"surfini: water ice TI:", inert_h2o_ice
	475
	476	c To start with : Initial albedo = observed dataset
	477	c -------------------------------------------------
	478	DO ig=1,ngrid
	479	psolaralb(ig,1)=albedodat(ig)
	480	psolaralb(ig,2)=albedodat(ig)
	481	END DO
[1130]	482	PRINT*,'surfini: minimum albedo without water caps',
	483	& minval(albedodat)
	484	PRINT*,'surfini: maximum albedo without water caps',
	485	& maxval(albedodat)
[38]	486
[283]	487	c initial albedo if permanent H2O ice is present
	488	c ------------------------------------------------
	489	IF ((water) .and. (caps)) THEN
	490	DO ig=1,ngrid
	491	IF (watercaptag(ig)) THEN
	492	psolaralb(ig,1) = albedo_h2o_ice
	493	psolaralb(ig,2) = albedo_h2o_ice
	494	ENDIF
	495	END DO
[1130]	496	PRINT*,'surfini: minimum albedo with water caps',
	497	& minval(albedodat)
	498	PRINT*,'surfini: maximum albedo with water caps',
	499	& maxval(albedodat)
[520]	500
[283]	501	ENDIF
	502
	503	c changing initial albedo if CO2 ice is present
	504	c -------------------------------------------
	505
	506	DO ig=1,ngrid
[38]	507	IF (piceco2(ig) .GT. 0.) THEN
[1130]	508	IF(lati(ig).LT. 0.) THEN
	509	icap=2 ! Southern hemisphere
[38]	510	ELSE
[1130]	511	icap=1 ! Northern hemisphere
[38]	512	ENDIF
	513	psolaralb(ig,1) = albedice(icap)
[283]	514	psolaralb(ig,2) = albedice(icap)
[38]	515	END IF
[283]	516	END DO
[38]	517
[283]	518	c changing initial albedo if water ice frost is present
	519	c -------------------------------------------
	520	IF (water) THEN
	521	do iq=1,nqmx
	522	c if there is frost and surface albedo is set to albedo_h2o_ice
	523	if(noms(iq).eq."h2o_ice") then
	524	do ig=1,ngrid
[712]	525
	526	if ((watercaptag(ig).eqv..false.)
	527	& .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
	528	print*, ''
	529	print*, '!!! PROBLEM in SURFINI !!!!'
	530	print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
	531	& WATERCAPTAG IS FALSE'
	532	print*, ''
	533	print*, 'ig,qsurf,threshold' ,
	534	& ig, qsurf(ig,iq), -frost_albedo_threshold
	535	print*, ''
	536	print*, '1) Check h2o_ice in startfi and ice
	537	& distribution in surfini'
	538	print*, '2) Use ini_h2osurf option in newstart'
	539	print*, ''
[1041]	540	#ifndef MESOINI
[712]	541	CALL ABORT
[1041]	542	#else
	543	qsurf(ig,iq) = 0.
	544	#endif
[712]	545	endif
	546
[283]	547	if ((piceco2(ig) .eq. 0.).and.
	548	& (qsurf(ig,iq).gt.frost_albedo_threshold)) then
	549	psolaralb(ig,1) = albedo_h2o_ice
	550	psolaralb(ig,2) = albedo_h2o_ice
	551	c PRINT*,'surfini.F frost',
	552	c & lati(ig)180./pi, long(ig)180./pi
	553	endif
	554	enddo
	555	endif
	556	end do
[1130]	557	PRINT*,'surfini: minimum albedo with frost and co2',
	558	& minval(albedodat)
	559	PRINT*,'surfini: maximum albedo with frost and co2',
	560	& maxval(albedodat)
[1087]	561	ELSE
	562	watercaptag(:) = .false.
[1130]	563	END IF ! OF IF(water)
[283]	564
[38]	565	RETURN
	566	END

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

Download in other formats:

Original Format