Context Navigation

surfini.F @ 2156

Last change on this file since 2156 was 1944, checked in by emillour, 6 years ago

Mars GCM:
A first step towards 1+1=2 (for now only works without tracers):

store and load "albedo" (surface albedo) and wstar (thermals' max vertical velocity) in physics (re)start file.
turn phyetat0 into a module in the process.

File size: 16.8 KB

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

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

Download in other formats:

Original Format