Context Navigation

surfini.F @ 804

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

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

Download in other formats:

Original Format