Context Navigation

surfini.F @ 1087

Last change on this file since 1087 was 1087, checked in by tnavarro, 11 years ago
bug when running the GCM without the water cycle : random ice patches could appear on the surface on the whole planet, thus affecting surface temperatures. Scary.
File size: 18.4 KB

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

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

Download in other formats:

Original Format