Context Navigation

surfini.F @ 607

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

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

Download in other formats:

Original Format