Context Navigation

surfini.F @ 700

Last change on this file since 700 was 697, checked in by emillour, 13 years ago

Mars GCM:

Minor cleanup in surfini.F
Corrected the polar mesh surface area which was wrong in the physics (changes in phyetat0.F, calfis.F and newstart.F)

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

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

Download in other formats:

Original Format