Context Navigation

surfini.F @ 652

Last change on this file since 652 was 643, checked in by tnavarro, 13 years ago
bug for mesocale in surfini corrected
File size: 14.6 KB

Rev	Line
[38]	1	SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
[520]	2	! to use 'getin'
	3	USE ioipsl_getincom
[38]	4	IMPLICIT NONE
	5	c=======================================================================
	6	c
	7	c creation des calottes pour l'etat initial
	8	c
	9	c=======================================================================
	10	c-----------------------------------------------------------------------
	11	c Declarations:
	12	c -------------
	13	#include "dimensions.h"
	14	#include "dimphys.h"
	15	#include "surfdat.h"
	16	#include "callkeys.h"
	17	#include "tracer.h"
[283]	18	#include "comgeomfi.h"
[285]	19	#include "comcstfi.h"
[38]	20	c
[283]	21	INTEGER ngrid,ig,icap,iq,alternate
[38]	22	REAL piceco2(ngrid),psolaralb(ngrid,2)
	23	REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
[520]	24	REAL icedryness ! ice dryness
[633]	25
	26	! longwatercaptag is watercaptag. Trick for some compilers
	27	LOGICAL, DIMENSION(100000) :: longwatercaptag
[38]	28
	29	EXTERNAL ISMIN,ISMAX
	30	INTEGER ISMIN,ISMAX
[520]	31
	32	! Choose false to have a somewhat non resolution dependant water ice caps distribution,
	33	! i.e based only on lat & lon values of each physical point.
	34	! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48
[633]	35	! For vizualisation : /u/tnalmd/bin/watercaps gcm_txt_output
[520]	36	LOGICAL,SAVE :: improvedicelocation = .true.
[38]	37	c
	38	c=======================================================================
	39
[283]	40	c water ice outliers
[38]	41	c ------------------------------------------
	42
[283]	43	IF ((water) .and. (caps)) THEN
[285]	44
[283]	45	c Perennial H20 north cap defined by watercaptag=true (allows surface to be
	46	c hollowed by sublimation in vdifc).
[38]	47
[283]	48	c We might not want albedodat to be modified because it is used to write
	49	c restart files. Instead, albedo is directly modified when needed (i.e.
	50	c if we have watercaptag and no co2 ice), below and in albedocaps.F90
	51
	52	c "Dryness coefficient" controlling the evaporation and
	53	c sublimation from the ground water ice (close to 1)
	54	c HERE, the goal is to correct for the fact
	55	c that the simulated permanent water ice polar caps
	56	c is larger than the actual cap and the atmospheric
	57	c opacity not always realistic.
	58
	59	alternate = 0
[520]	60
[633]	61	if (ngridmx .ne. 1) then
	62	watercaptag(:) = .false.
	63	longwatercaptag(:) = .false.
	64	endif
	65
[520]	66	write(,) "Ice dryness ?"
	67	icedryness=1. ! default value
	68	call getin("icedryness",icedryness)
	69	write(,) " icedryness = ",icedryness
[633]	70
[285]	71
	72	#ifdef MESOSCALE
[633]	73
[520]	74	do ig=1,ngridmx
[633]	75
[317]	76	!write(,) "all qsurf to zero. dirty."
[285]	77	do iq=1,nqmx
	78	qsurf(ig,iq)=0. !! on jette les inputs GCM
	79	!! on regle juste watercaptag
	80	!! il faudrait garder les inputs GCM
	81	!! si elles sont consequentes
	82	enddo
	83	if ( ( lati(ig)*180./pi .gt. 70. ) .and.
	84	. ( albedodat(ig) .ge. 0.26 ) ) then
	85	write(,)"outlier ",ig
	86	watercaptag(ig) = .true.
	87	dryness(ig) = 1.
	88	albedodat(ig) = albedo_h2o_ice !! pour output
	89	else
	90	watercaptag(ig) = .false.
	91	dryness(ig) = 1.
[633]	92	endif
	93
[643]	94	enddo
	95
[285]	96	#else
	97
[520]	98	if (improvedicelocation) then
	99
	100	if (ngridmx .eq. 738) then ! hopefully 32x24
[633]	101
[520]	102	print*,'water ice caps distribution for 32x24 resolution:'
[633]	103	longwatercaptag(1:9) = .true. ! central cap - core
	104	longwatercaptag(26:33) = .true. ! central cap
[520]	105
	106	else if (ngridmx .eq. 3010) then ! hopefully 64x48
[633]	107
	108	! For latitudes:
	109	! [ 67.5 71.25 75. 78.75 82.5 86.25]
	110	! The water ice caps should be (according to TES temperatures):
	111	! [ 8.63e-03 1.02e+00 5.99e+00 2.66e+01 5.65e+01]
	112
[520]	113	print*,'water ice caps distribution for 64x48 resolution:'
[633]	114	longwatercaptag(1:65) = .true. ! central cap - core
	115	longwatercaptag(75:85) = .true. ! central cap
	116	longwatercaptag(93:114) = .true. ! central cap
	117	!--------------------- OUTLIERS ----------------------------
	118	longwatercaptag(136) = .true. ! outlier, lat = 78.75
	119	longwatercaptag(138) = .true. ! outlier, lat = 78.75
	120	longwatercaptag(140) = .true. ! outlier, lat = 78.75
	121	longwatercaptag(142) = .true. ! outlier, lat = 78.75
	122	longwatercaptag(161) = .true. ! outlier, lat = 78.75
	123	longwatercaptag(163) = .true. ! outlier, lat = 78.75
	124	longwatercaptag(165) = .true. ! outlier, lat = 78.75
	125	longwatercaptag(183) = .true. ! outlier, lat = 78.75
	126	longwatercaptag(185) = .true. ! outlier, lat = 78.75
	127	longwatercaptag(187) = .true. ! outlier, lat = 78.75
	128	longwatercaptag(189) = .true. ! outlier, lat = 78.75
	129	longwatercaptag(191) = .true. ! outlier, lat = 78.75
	130	longwatercaptag(193) = .true. ! outlier, lat = 78.75
	131	longwatercaptag(194) = .true. ! outlier, lat = 75
	132	longwatercaptag(203) = .true. ! outlier, lat = 75
	133	longwatercaptag(207) = .true. ! outlier, lat = 75
	134	longwatercaptag(242) = .true. ! outlier, lat = 75
	135	longwatercaptag(244) = .true. ! outlier, lat = 75
	136	longwatercaptag(246) = .true. ! outlier, lat = 75
	137	longwatercaptag(248) = .true. ! outlier, lat = 75
	138	longwatercaptag(250) = .true. ! outlier, lat = 75
	139	longwatercaptag(252) = .true. ! outlier, lat = 75
	140	longwatercaptag(254) = .true. ! outlier, lat = 75
	141	longwatercaptag(256:257)= .true. ! outlier, lat = 75
	142	!!------------------------- OLD ----------------------------
	143	!! longwatercaptag(83:85) = .true.
	144	!! longwatercaptag(135:142) = .true.
	145	!! longwatercaptag(181:193) = .true.
	146	!! longwatercaptag(242:257) = .true.
	147
[520]	148
	149	else if (ngridmx .ne. 1) then
	150	print*,'No improved ice location for this resolution!'
	151	print*,'Please set improvedicelocation to false in surfini.'
	152	print*,'And check lat and lon conditions for ice caps in code.'
	153	call exit(1)
	154
	155	endif
	156
	157	! print caps locations
	158	print*,'latitude \| longitude \| ig'
	159	do ig=1,ngridmx
	160	dryness (ig) = icedryness
[633]	161
	162	if (longwatercaptag(ig)) then
	163	watercaptag(ig) = .True.
[520]	164	print,'ice water cap', lati(ig)180./pi,
	165	. long(ig)*180./pi, ig
	166	endif
	167	enddo
[283]	168
[633]	169
[520]	170	else
	171
	172	do ig=1,ngridmx
	173
	174	dryness (ig) = icedryness
	175
[285]	176	!!c Towards olympia planitia water caps ('relatively' low latitude ones)
	177	!!c---------------- proposition par AS --------------------
	178	!!c--------------------------------------------------------
	179	!!c if ( ( lati(ig)*180./pi .ge. 75 ) .and.
	180	!!c . ( lati(ig)*180./pi .le. 77 ) .and.
	181	!!c . ( ( ( long(ig)*180./pi .ge. 000. ) .and.
	182	!!c . ( long(ig)*180./pi .le. 120. ) )
	183	!!c . .or.
	184	!!c . ( ( long(ig)*180./pi .ge. -130. ) .and.
	185	!!c . ( long(ig)*180./pi .le. -115. ) ) ) ) then
	186	!!c---------------- proposition par TN --------------------
	187	!!c---------------- HIGHLY EXPERIMENTAL -------------------
	188	!!c--------------------------------------------------------
	189	!! if ( ( ( lati(ig)*180./pi .ge. 73. ) .and. ! cap #1
	190	!! . ( lati(ig)*180./pi .le. 75.1 ) .and.
	191	!! . ( long(ig)*180./pi .ge. 95. ) .and.
	192	!! . ( long(ig)*180./pi .le. 110. ) )
	193	!! . .or.
	194	!! . ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2
	195	!! . ( lati(ig)*180./pi .le. 80. ) .and.
	196	!! . ( long(ig)*180./pi .ge. 110. ) .and.
	197	!! . ( long(ig)*180./pi .le. 140. ) )
	198	!! . .or.
	199	!! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3
	200	!! . ( lati(ig)*180./pi .le. 78. ) .and.
	201	!! . ( long(ig)*180./pi .ge. 155. ) .and.
	202	!! . ( long(ig)*180./pi .le. 180. ) )
[520]	203	!! . .or.
	204	!! . ( ( lati(ig)*180./pi .ge. 71. ) .and. ! cap #4 (Korolev crater)
	205	!! . ( lati(ig)*180./pi .le. 72. ) .and.
	206	!! . ( long(ig)*180./pi .ge. 163. ) .and.
	207	!! . ( long(ig)*180./pi .le. 168. ) )
[285]	208	!! . .or.
	209	!! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5
	210	!! . ( lati(ig)*180./pi .le. 78. ) .and.
	211	!! . ( long(ig)*180./pi .ge. -160.) .and.
	212	!! . ( long(ig)*180./pi .le. -120.) ) )
	213	!! . then
[520]	214	if ( ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2
	215	. ( lati(ig)*180./pi .le. 80. ) .and.
	216	. ( long(ig)*180./pi .ge. 110. ) .and.
	217	. ( long(ig)*180./pi .le. 181. ) )
	218	. .or.
[283]	219
[520]	220	. ( ( lati(ig)*180./pi .ge. 75. ) .and. ! cap #4 (Korolev crater)
	221	. ( lati(ig)*180./pi .le. 76. ) .and.
	222	. ( long(ig)*180./pi .ge. 150. ) .and.
	223	. ( long(ig)*180./pi .le. 168. ) )
	224	. .or.
	225	. ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
	226	. ( lati(ig)*180./pi .le. 80. ) .and.
	227	. ( long(ig)*180./pi .ge. -150.) .and.
	228	. ( long(ig)*180./pi .le. -110.) ) )
	229	. then
	230
	231
	232	if ((alternate .eq. 0)) then ! 1/2 en 64x48 sinon trop large en lat
	233	if (ngridmx.ne.1) watercaptag(ig)=.true.
[633]	234	print,'ice water cap', lati(ig)180./pi,
	235	. long(ig)*180./pi, ig
	236	dryness(ig) = 1.
[520]	237	alternate = 1
	238	else
	239	alternate = 0
	240	endif !end if alternate = 0
[633]	241
[520]	242
	243	endif
	244
	245
	246	c Opposite olympia planitia water cap
	247	c---------------- proposition par AS --------------------
	248	c--------------------------------------------------------
	249	c if ( ( lati(ig)*180./pi .ge. 82 ) .and.
	250	c . ( lati(ig)*180./pi .le. 84 ) .and.
	251	c . ( ( long(ig)*180./pi .gt. -030. ) .and.
	252	c . ( long(ig)*180./pi .lt. 090. ) ) ) then
	253	c---------------- proposition par TN --------------------
	254	c--------------------------------------------------------
	255	if ( ( ( lati(ig)*180./pi .ge. 80 ) .and.
	256	. ( lati(ig)*180./pi .le. 84 ) )
	257	. .and.
	258	. ( ( long(ig)*180./pi .lt. -95. ) .or. !!! 32x24
	259	. ( long(ig)*180./pi .gt. 85. ) ) ) then !!! 32x24
	260	! . ( ( ( long(ig)*180./pi .ge. -29. ) .and. !!! 64x48
	261	! . ( long(ig)*180./pi .le. 90. ) ) .or. !!! 64x48
	262	! . ( ( long(ig)*180./pi .ge. -77. ) .and. !!! 64x48
	263	! . ( long(ig)*180./pi .le. -70. ) ) ) ) then !!! 64x48
	264	if (ngridmx.ne.1) then
	265	watercaptag(ig)=.true.
	266	print,'ice water cap', lati(ig)180./pi,
	267	. long(ig)*180./pi, ig
	268	endif
[633]	269	dryness(ig) = 1.
[520]	270	endif
	271
[285]	272	c Central cap
	273	c---------------- anciens reglages FF -------------------
[283]	274	c--------------------------------------------------------
	275
[633]	276	if (abs(lati(ig)*180./pi).gt.80) then
[520]	277	print,'ice water cap', lati(ig)180./pi,
	278	. long(ig)*180./pi, ig
[283]	279	if (ngridmx.ne.1) watercaptag(ig)=.true.
	280	!dryness(ig) = 1.
	281	c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
[285]	282	c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
	283	c if (ngridmx.ne.1) watercaptag(ig)=.true.
[283]	284	c dryness(ig) = 1.
[285]	285	c endif
[283]	286	c if (lati(ig)*180./pi.ge.85) then
	287	c if (ngridmx.ne.1) watercaptag(ig)=.true.
	288	c dryness(ig) = 1.
[285]	289	c endif
	290	endif ! (lati>80 deg)
[520]	291
	292	end do ! (ngridmx)
	293
	294	endif ! of if (improvedicelocation)
[285]	295	#endif
[283]	296
[520]	297	ENDIF ! (caps & water)
	298
[283]	299	c ===============================================================
	300	c INITIAL ALBEDO
	301	c ===============================================================
	302
	303	write(,)"surfini: water frost thickness",
	304	s frost_albedo_threshold
	305	write(,)"surfini: water ice albedo:", albedo_h2o_ice
	306	write(,)"surfini: water ice TI:", inert_h2o_ice
	307
	308	c To start with : Initial albedo = observed dataset
	309	c -------------------------------------------------
	310	DO ig=1,ngrid
	311	psolaralb(ig,1)=albedodat(ig)
	312	psolaralb(ig,2)=albedodat(ig)
	313	END DO
	314	PRINT*,'minimum albedo sans water caps',
[38]	315	s albedodat(ISMIN(ngrid,albedodat,1))
[283]	316	PRINT*,'maximum albedo sans water caps',
[38]	317	s albedodat(ISMAX(ngrid,albedodat,1))
	318
[283]	319	c initial albedo if permanent H2O ice is present
	320	c ------------------------------------------------
	321	IF ((water) .and. (caps)) THEN
	322	DO ig=1,ngrid
	323	IF (watercaptag(ig)) THEN
	324	psolaralb(ig,1) = albedo_h2o_ice
	325	psolaralb(ig,2) = albedo_h2o_ice
	326	ENDIF
	327	END DO
	328	PRINT*,'minimum albedo avec water caps',
	329	s psolaralb(ISMIN(ngrid,psolaralb,1),1)
	330	PRINT*,'maximum albedo avec water caps',
	331	s psolaralb(ISMAX(ngrid,psolaralb,1),1)
[520]	332
[283]	333	ENDIF
	334
	335	c changing initial albedo if CO2 ice is present
	336	c -------------------------------------------
	337
	338	DO ig=1,ngrid
[38]	339	IF (piceco2(ig) .GT. 0.) THEN
	340	IF(ig.GT.ngrid/2+1) THEN
	341	icap=2
	342	ELSE
	343	icap=1
	344	ENDIF
	345	psolaralb(ig,1) = albedice(icap)
[283]	346	psolaralb(ig,2) = albedice(icap)
[38]	347	END IF
[283]	348	END DO
[38]	349
[283]	350	c changing initial albedo if water ice frost is present
	351	c -------------------------------------------
	352	IF (water) THEN
	353	do iq=1,nqmx
	354	c if there is frost and surface albedo is set to albedo_h2o_ice
	355	if(noms(iq).eq."h2o_ice") then
	356	do ig=1,ngrid
	357	if ((piceco2(ig) .eq. 0.).and.
	358	& (qsurf(ig,iq).gt.frost_albedo_threshold)) then
	359	psolaralb(ig,1) = albedo_h2o_ice
	360	psolaralb(ig,2) = albedo_h2o_ice
	361	c PRINT*,'surfini.F frost',
	362	c & lati(ig)180./pi, long(ig)180./pi
	363	endif
	364	enddo
	365	endif
	366	end do
	367	PRINT*,'minimum albedo avec givre et co2',
[38]	368	s psolaralb(ISMIN(ngrid,psolaralb,1),1)
[283]	369	PRINT*,'maximum albedo avec givre et co2',
[38]	370	s psolaralb(ISMAX(ngrid,psolaralb,1),1)
[283]	371	END IF
	372
[38]	373
	374	RETURN
	375	END

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

Download in other formats:

Original Format