Context Navigation

physdem.F @ 1085

Last change on this file since 1085 was 1085, checked in by emillour, 11 years ago
Removed wrongly referenced module variables in lwb and physdem. Added option "-full" to makegcm_* to force a full recompilation from scratch (ie: removing libo/... and makefile before launching makefile building and compilation). EM
File size: 30.5 KB

Rev	Line
[1047]	1	subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq,
[999]	2	. phystep,day_ini,time,airefi,
	3	. alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
	4
[1036]	5	use infotrac, only: nqtot, tnom
[1047]	6	use comsoil_h, only: inertiedat, volcapa, mlayer
	7	use comgeomfi_h, only: area
	8	use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe,
	9	& z0_default, albedice, emisice, emissiv,
	10	& iceradius, dtemisice, phisfi, z0
	11	use yomaer_h, only: tauvis
[999]	12	implicit none
	13	c
	14	c Variables qui NE possèdent PAS un axe des temps :
	15	c airefi
	16	c alb
	17	c ith
	18	c pzmea
	19	c pzstd
	20	c pzsig
	21	c pzgam
	22	c pzthe
	23	c
	24	c-------------------------------------------------------------
	25	c
	26	c create physics (re-)start data file "restartfi.nc"
	27	c
	28	c
	29	c
	30	#include "dimensions.h"
	31	#include "paramet.h"
	32	c-----------------------------------------------------------------------
	33	#include "comvert.h"
	34	#include "comgeom2.h"
	35	#include "control.h"
	36	#include "comdissnew.h"
	37	#include "logic.h"
	38	#include "ener.h"
	39	#include "netcdf.inc"
[1047]	40	!#include "dimphys.h"
[1036]	41	!#include "advtrac.h"
[999]	42	#include "callkeys.h"
	43	c
	44	INTEGER nid,iq
	45	INTEGER, parameter :: ivap=1
	46	REAL, parameter :: qsolmax= 150.0
	47	character (len=*) :: filename
	48	character (len=7) :: str7
	49
	50	REAL day_ini
[1047]	51	INTEGER,INTENT(IN) :: nsoil,ngrid,nlay,nq
[999]	52	integer ierr,idim1,idim2,idim3,idim4,idim5,idim6,nvarid
	53
	54	c
	55	REAL phystep,time
[1047]	56	REAL latfi(ngrid), lonfi(ngrid)
	57	! REAL champhys(ngrid)
[999]	58	INTEGER length
	59	PARAMETER (length=100)
	60	REAL tab_cntrl(length)
	61
	62	c
	63
	64	! EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression
	65	! EXTERNAL exner_hyb , SSUM
	66	c
	67	#include "serre.h"
	68	#include "clesph0.h"
	69	#include "fxyprim.h"
[1047]	70	!#include "comgeomfi.h"
	71	!#include "surfdat.h"
	72	!#include "comsoil.h"
[999]	73	#include "planete.h"
[1047]	74	!#include "dimradmars.h"
	75	!#include "yomaer.h"
[999]	76	#include "comcstfi.h"
	77
[1047]	78	real airefi(ngrid)
	79	real alb(ngrid),ith(ngrid,nsoil)
	80	real pzmea(ngrid),pzstd(ngrid)
	81	real pzsig(ngrid),pzgam(ngrid),pzthe(ngrid)
[999]	82	integer ig
	83
	84	! flag which identifies if we are using old tracer names (qsurf01,...)
	85	logical :: oldtracernames=.false.
	86	integer :: count
	87	character(len=30) :: txt ! to store some text
	88	! indexes of water vapour & water ice tracers (if any):
	89	integer :: i_h2o_vap=0
	90	integer :: i_h2o_ice=0
	91	c-----------------------------------------------------------------------
	92
	93	! copy airefi(:) to area(:)
[1047]	94	CALL SCOPY(ngrid,airefi,1,area,1)
	95	! note: area() is defined in comgeomfi_h
[999]	96
	97
[1047]	98	DO ig=1,ngrid
	99	albedodat(ig)=alb(ig) ! note: albedodat() is defined in surfdat_h
	100	zmea(ig)=pzmea(ig) ! note: zmea() is defined in surfdat_h
	101	zstd(ig)=pzstd(ig) ! note: zstd() is defined in surfdat_h
	102	zsig(ig)=pzsig(ig) ! note: zsig() is defined in surfdat_h
	103	zgam(ig)=pzgam(ig) ! note: zgam() is defined in surfdat_h
	104	zthe(ig)=pzthe(ig) ! note: zthe() is defined in surfdat_h
[999]	105	ENDDO
	106
[1047]	107	inertiedat(:,:)=ith(:,:) ! note inertiedat() is defined in comsoil_h
[999]	108	c
	109	c things to store in the physics start file:
	110	c
	111	ierr = NF_CREATE(adjustl(filename),
	112	& IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
	113	IF (ierr.NE.NF_NOERR) THEN
	114	WRITE(6,*)'physdem0: Problem creating file ',adjustl(filename)
	115	write(6,*) NF_STRERROR(ierr)
	116	CALL ABORT
	117	ENDIF
	118	c
	119	ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 18,
	120	. "Physics start file")
	121	c
	122	ierr = NF_DEF_DIM (nid,"index",length,idim1)
	123	if (ierr.ne.NF_NOERR) then
	124	WRITE(6,*)'physdem0: Problem defining index dimension'
	125	write(6,*) NF_STRERROR(ierr)
	126	call abort
	127	endif
	128	c
[1047]	129	ierr = NF_DEF_DIM (nid,"physical_points",ngrid,idim2)
[999]	130	if (ierr.ne.NF_NOERR) then
	131	WRITE(6,*)'physdem0: Problem defining physical_points dimension'
	132	write(6,*) NF_STRERROR(ierr)
	133	call abort
	134	endif
	135	c
	136	ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoil,idim3)
	137	if (ierr.ne.NF_NOERR) then
	138	WRITE(6,*)'physdem0: Problem defining subsurface_layers dimension'
	139	write(6,*) NF_STRERROR(ierr)
	140	call abort
	141	endif
	142	c
	143	! ierr = NF_DEF_DIM (nid,"nlayer+1",llm+1,idim4)
	144	ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
	145	if (ierr.ne.NF_NOERR) then
	146	WRITE(6,*)'physdem0: Problem defining nlayer+1 dimension'
	147	write(6,*) NF_STRERROR(ierr)
	148	call abort
	149	endif
	150	c
	151	ierr = NF_DEF_DIM (nid,"number_of_advected_fields",nq,idim5)
	152	if (ierr.ne.NF_NOERR) then
	153	WRITE(6,*)'physdem0: Problem defining advected fields dimension'
	154	WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
	155	write(6,*) NF_STRERROR(ierr)
	156	endif
	157	c
	158	ierr = NF_DEF_DIM (nid,"Time",NF_UNLIMITED,idim6)
	159	if (ierr.ne.NF_NOERR) then
	160	WRITE(6,*)'physdem0: Problem defining time dimension'
	161	WRITE(6,*)' ierr = ', ierr
	162	write(6,*) NF_STRERROR(ierr)
	163	endif
	164
	165	ierr = NF_ENDDEF(nid) ! exit NetCDF define mode
	166
	167	c clear tab_cntrl(:) array
	168	DO ierr = 1, length
	169	tab_cntrl(ierr) = 0.0
	170	ENDDO
	171
[1047]	172	write(,) "physdem0: ngrid: ",ngrid
[999]	173
	174	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	175	c Fill control array tab_cntrl(:) with parameters for this run
	176	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	177	c Informations on the physics grid
[1047]	178	tab_cntrl(1) = float(ngrid) ! number of nodes on physics grid
	179	tab_cntrl(2) = float(nlay) ! number of atmospheric layers
[999]	180	tab_cntrl(3) = day_ini + int(time) ! initial day
	181	tab_cntrl(4) = time -int(time) ! initial time of day
	182
	183	c Informations about Mars, used by dynamics and physics
	184	tab_cntrl(5) = rad ! radius of Mars (m) ~3397200
	185	tab_cntrl(6) = omeg ! rotation rate (rad.s-1)
	186	tab_cntrl(7) = g ! gravity (m.s-2) ~3.72
	187	tab_cntrl(8) = mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49
	188	tab_cntrl(9) = rcp ! = r/cp ~0.256793 (=kappa dans dynamique)
	189	tab_cntrl(10) = daysec ! length of a sol (s) ~88775
	190
	191	tab_cntrl(11) = phystep ! time step in the physics
	192	tab_cntrl(12) = 0.
	193	tab_cntrl(13) = 0.
	194
	195	c Informations about Mars, only for physics
	196	tab_cntrl(14) = year_day ! length of year (sols) ~668.6
	197	tab_cntrl(15) = periheli ! min. Sun-Mars distance (Mkm) ~206.66
	198	tab_cntrl(16) = aphelie ! max. SUn-Mars distance (Mkm) ~249.22
	199	tab_cntrl(17) = peri_day ! date of perihelion (sols since N. spring)
	200	tab_cntrl(18) = obliquit ! Obliquity of the planet (deg) ~23.98
	201
	202	c Boundary layer and turbulence
	203	tab_cntrl(19) = z0_default ! default surface roughness (m) ~0.01
	204	tab_cntrl(20) = lmixmin ! mixing length ~100
	205	tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
	206
	207	c Optical properties of polar caps and ground emissivity
	208	tab_cntrl(22) = albedice(1) ! Albedo of northern cap ~0.5
	209	tab_cntrl(23) = albedice(2) ! Albedo of southern cap ~0.5
	210	tab_cntrl(24) = emisice(1) ! Emissivity of northern cap ~0.95
	211	tab_cntrl(25) = emisice(2) ! Emissivity of southern cap ~0.95
	212	tab_cntrl(26) = emissiv ! Emissivity of martian soil ~.95
	213	tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
	214	tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
	215	tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
	216	tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
	217
	218	c dust aerosol properties
	219	tab_cntrl(27) = tauvis ! mean visible optical depth
	220
	221	tab_cntrl(28) = 0.
	222	tab_cntrl(29) = 0.
	223	tab_cntrl(30) = 0.
	224
	225	! Soil properties:
	226	tab_cntrl(35) = volcapa ! soil volumetric heat capacity
	227
	228	c
	229
	230	ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
	231	IF (ierr.NE.NF_NOERR) THEN
	232	PRINT*, 'physdem0: Failed to swich to NetCDF define mode'
	233	CALL abort
	234	ENDIF
	235	! define variable
	236	#ifdef NC_DOUBLE
	237	ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
	238	#else
	239	ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
	240	#endif
	241	IF (ierr.NE.NF_NOERR) THEN
	242	PRINT*, 'physdem0: Failed to define controle'
	243	CALL abort
	244	ENDIF
	245	ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
	246	. "Control parameters")
	247	IF (ierr.NE.NF_NOERR) THEN
	248	PRINT*, 'physdem0: Failed to define controle title attribute'
	249	CALL abort
	250	ENDIF
	251	ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
	252	IF (ierr.NE.NF_NOERR) THEN
	253	PRINT*, 'physdem0: Failed to swich out of NetCDF define mode'
	254	CALL abort
	255	ENDIF
	256	! write variable
	257	#ifdef NC_DOUBLE
	258	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
	259	#else
	260	ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
	261	#endif
	262	IF (ierr.NE.NF_NOERR) THEN
	263	PRINT*, 'physdem0: Failed to write controle data'
	264	CALL abort
	265	ENDIF
	266
	267
[1047]	268	! write mid-layer depths mlayer() !known from comsoil_h
[999]	269
	270	ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
	271	! define variable
	272	#ifdef NC_DOUBLE
	273	ierr = NF_DEF_VAR (nid,"soildepth",NF_DOUBLE,1,idim3,nvarid)
	274	#else
	275	ierr = NF_DEF_VAR (nid,"soildepth",NF_FLOAT,1,idim3,nvarid)
	276	#endif
	277	ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
	278	. "Soil mid-layer depth")
	279	ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
	280	! write variable
	281	#ifdef NC_DOUBLE
	282	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
	283	#else
	284	ierr = NF_PUT_VAR_REAL (nid,nvarid,mlayer)
	285	#endif
	286
	287	c
	288
	289	ierr = NF_REDEF (nid)
	290	#ifdef NC_DOUBLE
	291	ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
	292	#else
	293	ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
	294	#endif
	295	ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 26,
	296	. "Longitudes of physics grid")
	297	ierr = NF_ENDDEF(nid)
	298
	299	#ifdef NC_DOUBLE
	300	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
	301	#else
	302	ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
	303	#endif
	304
	305	c
	306
	307	ierr = NF_REDEF (nid)
	308	#ifdef NC_DOUBLE
	309	ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
	310	#else
	311	ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
	312	#endif
	313	ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
	314	. "Latitudes of physics grid")
	315	ierr = NF_ENDDEF(nid)
	316	#ifdef NC_DOUBLE
	317	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
	318	#else
	319	ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
	320	#endif
	321
	322	c
	323
	324	ierr = NF_REDEF (nid)
	325	#ifdef NC_DOUBLE
	326	ierr = NF_DEF_VAR (nid, "area", NF_DOUBLE, 1, idim2,nvarid)
	327	#else
	328	ierr = NF_DEF_VAR (nid, "area", NF_FLOAT, 1, idim2,nvarid)
	329	#endif
	330	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
	331	. "Mesh area")
	332	ierr = NF_ENDDEF(nid)
	333	#ifdef NC_DOUBLE
	334	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
	335	#else
	336	ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
	337	#endif
	338
	339	c
	340
	341	ierr = NF_REDEF (nid)
	342	#ifdef NC_DOUBLE
	343	ierr = NF_DEF_VAR (nid, "phisfi", NF_DOUBLE, 1, idim2,nvarid)
	344	#else
	345	ierr = NF_DEF_VAR (nid, "phisfi", NF_FLOAT, 1, idim2,nvarid)
	346	#endif
	347	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
	348	. "Geopotential at the surface")
	349	ierr = NF_ENDDEF(nid)
	350	#ifdef NC_DOUBLE
	351	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phisfi)
	352	#else
	353	ierr = NF_PUT_VAR_REAL (nid,nvarid,phisfi)
	354	#endif
	355
	356	c
	357
	358	ierr = NF_REDEF (nid)
	359	#ifdef NC_DOUBLE
	360	ierr = NF_DEF_VAR (nid, "albedodat", NF_DOUBLE, 1, idim2,nvarid)
	361	#else
	362	ierr = NF_DEF_VAR (nid, "albedodat", NF_FLOAT, 1, idim2,nvarid)
	363	#endif
	364	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
	365	. "Albedo of bare ground")
	366	ierr = NF_ENDDEF(nid)
	367	#ifdef NC_DOUBLE
	368	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedodat)
	369	#else
	370	ierr = NF_PUT_VAR_REAL (nid,nvarid,albedodat)
	371	#endif
	372
	373	c
	374	c some data for Francois Lott's programs
	375	c
	376
	377	ierr = NF_REDEF (nid)
	378	#ifdef NC_DOUBLE
	379	ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
	380	#else
	381	ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
	382	#endif
	383	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
	384	. "Relief: mean relief")
	385	ierr = NF_ENDDEF(nid)
	386	#ifdef NC_DOUBLE
	387	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
	388	#else
	389	ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
	390	#endif
	391	c
	392	ierr = NF_REDEF (nid)
	393	#ifdef NC_DOUBLE
	394	ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
	395	#else
	396	ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
	397	#endif
	398	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
	399	. "Relief: standard deviation")
	400	ierr = NF_ENDDEF(nid)
	401	#ifdef NC_DOUBLE
	402	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
	403	#else
	404	ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
	405	#endif
	406	c
	407	ierr = NF_REDEF (nid)
	408	#ifdef NC_DOUBLE
	409	ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
	410	#else
	411	ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
	412	#endif
	413	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
	414	. "Relief: sigma parameter")
	415	ierr = NF_ENDDEF(nid)
	416	#ifdef NC_DOUBLE
	417	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
	418	#else
	419	ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
	420	#endif
	421	c
	422	ierr = NF_REDEF (nid)
	423	#ifdef NC_DOUBLE
	424	ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
	425	#else
	426	ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
	427	#endif
	428	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
	429	. "Relief: gamma parameter")
	430	ierr = NF_ENDDEF(nid)
	431	#ifdef NC_DOUBLE
	432	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
	433	#else
	434	ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
	435	#endif
	436	c
	437	ierr = NF_REDEF (nid)
	438	#ifdef NC_DOUBLE
	439	ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
	440	#else
	441	ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
	442	#endif
	443	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
	444	. "Relief: theta parameter")
	445	ierr = NF_ENDDEF(nid)
	446	#ifdef NC_DOUBLE
	447	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
	448	#else
	449	ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
	450	#endif
	451
	452
	453	! Soil Thermal inertia
	454
	455	ierr = NF_REDEF (nid)
	456	#ifdef NC_DOUBLE
	457	ierr = NF_DEF_VAR (nid,"inertiedat",NF_DOUBLE,
	458	& 2,(/idim2,idim3/),nvarid)
	459	#else
	460	ierr = NF_DEF_VAR (nid,"inertiedat",NF_FLOAT,
	461	& 2,(/idim2,idim3/),nvarid)
	462	#endif
	463	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 20,
	464	. "Soil thermal inertia")
	465	ierr = NF_ENDDEF(nid)
	466	#ifdef NC_DOUBLE
	467	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,inertiedat)
	468	#else
	469	ierr = NF_PUT_VAR_REAL (nid,nvarid,inertiedat)
	470	#endif
	471
[1047]	472	! surface roughness length (z0 is a common in surfdat_h)
[999]	473
	474	ierr = NF_REDEF (nid)
	475	#ifdef NC_DOUBLE
	476	ierr = NF_DEF_VAR (nid, "z0", NF_DOUBLE, 1, idim2,nvarid)
	477	#else
	478	ierr = NF_DEF_VAR (nid, "z0", NF_FLOAT, 1, idim2,nvarid)
	479	#endif
	480	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 24,
	481	. "Surface roughness length")
	482	ierr = NF_ENDDEF(nid)
	483	#ifdef NC_DOUBLE
	484	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,z0)
	485	#else
	486	ierr = NF_PUT_VAR_REAL (nid,nvarid,z0)
	487	#endif
	488
	489
	490
	491
	492	c--------------------- Time dependent variables ---------------------
	493	c-------------------------------------------------------------------
	494	c-------------------------------------------------------------------
	495
	496
	497	! Time
	498
	499	ierr = NF_REDEF (nid)
	500	#ifdef NC_DOUBLE
	501	ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1,
	502	. (/idim6/), nvarid)
	503	#else
	504	ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1,
	505	. (/idim6/), nvarid)
	506	#endif
	507	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
	508	. "Temps de simulation")
	509	ierr = NF_ENDDEF(nid)
	510
	511
	512	! CO2 Ice Cover
	513
	514	ierr = NF_REDEF (nid)
	515	#ifdef NC_DOUBLE
	516	ierr = NF_DEF_VAR (nid, "co2ice", NF_DOUBLE, 2,
	517	. (/idim2,idim6/), nvarid)
	518	#else
	519	ierr = NF_DEF_VAR (nid, "co2ice", NF_FLOAT, 2,
	520	. (/idim2,idim6/), nvarid)
	521	#endif
	522	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 13,
	523	. "CO2 ice cover")
	524	ierr = NF_ENDDEF(nid)
	525
	526
	527	! Surface temperature
	528
	529	ierr = NF_REDEF (nid)
	530	#ifdef NC_DOUBLE
	531	ierr = NF_DEF_VAR (nid, "tsurf", NF_DOUBLE, 2,
	532	. (/idim2,idim6/), nvarid)
	533	#else
	534	ierr = NF_DEF_VAR (nid, "tsurf", NF_FLOAT, 2,
	535	. (/idim2,idim6/), nvarid)
	536	#endif
	537	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
	538	. "Surface temperature")
	539	ierr = NF_ENDDEF(nid)
	540
	541
	542	! Soil temperature
	543
	544	ierr = NF_REDEF (nid)
	545	#ifdef NC_DOUBLE
	546	ierr = NF_DEF_VAR (nid, "tsoil", NF_DOUBLE, 3,
	547	. (/idim2,idim3,idim6/), nvarid)
	548	#else
	549	ierr = NF_DEF_VAR (nid, "tsoil", NF_FLOAT, 3,
	550	. (/idim2,idim3,idim6/), nvarid)
	551	#endif
	552	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 16,
	553	. "Soil temperature")
	554	ierr = NF_ENDDEF(nid)
	555
	556
	557	c emissivity
	558
	559	ierr = NF_REDEF (nid)
	560	#ifdef NC_DOUBLE
	561	ierr = NF_DEF_VAR (nid,"emis", NF_DOUBLE, 2,
	562	. (/idim2,idim6/), nvarid)
	563	#else
	564	ierr = NF_DEF_VAR (nid, "emis", NF_FLOAT, 2,
	565	. (/idim2,idim6/), nvarid)
	566	#endif
	567	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 18,
	568	. "Surface emissivity")
	569	ierr = NF_ENDDEF(nid)
	570
	571
	572	c planetary boundary layer
	573
	574	ierr = NF_REDEF (nid)
	575	#ifdef NC_DOUBLE
	576	ierr = NF_DEF_VAR (nid, "q2", NF_DOUBLE, 3,
	577	. (/idim2,idim4,idim6/), nvarid)
	578	#else
	579	ierr = NF_DEF_VAR (nid, "q2", NF_FLOAT, 3,
	580	. (/idim2,idim4,idim6/), nvarid)
	581	#endif
	582	ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
	583	. "pbl wind variance")
	584	ierr = NF_ENDDEF(nid)
	585
	586
	587
	588
	589	c tracers
	590	!
	591	! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
	592	! qsurf02, ...)
	593	count=0
[1036]	594	do iq=1,nqtot
[999]	595	txt= " "
	596	write(txt,'(a1,i2.2)')'q',iq
	597	if (txt.ne.tnom(iq)) then ! use tracer names stored in dynamics
	598	! did not find old tracer name
	599	exit ! might as well stop here
	600	else
	601	! found old tracer name
	602	count=count+1
	603	endif
	604	enddo
[1036]	605	if (count.eq.nqtot) then
[999]	606	write(,) "physdem0:tracers seem to follow old naming ",
	607	& "convention (qsurf01,qsurf02,...)"
	608	write(,) " => will work for now ... "
	609	write(,) " but you should run newstart to rename them"
	610	oldtracernames=.true.
[1036]	611	endif ! of if (count.eq.nqtot)
[999]	612
	613	IF(nq.GE.1) THEN
	614	! preliminary stuff: look for water vapour & water ice tracers (if any)
	615	if (.not.oldtracernames) then
	616	do iq=1,nq
	617	if (tnom(iq).eq."h2o_vap") then
	618	i_h2o_vap=iq
	619	endif
	620	if (tnom(iq).eq."h2o_ice") then
	621	i_h2o_ice=iq
	622	endif
	623	enddo ! of iq=1,nq
	624	! handle special case of only water vapour tracer (no ice)
	625	if ((i_h2o_vap.ne.0).and.(i_h2o_ice.eq.0)) then
	626	! then the index of (surface) ice is i_h2o_vap
	627	i_h2o_ice=i_h2o_vap
	628	endif
	629	endif ! of if (.not.oldtracernames)
	630
	631	DO iq=1,nq
	632	IF (oldtracernames) THEN
	633	txt=" "
	634	write(txt,'(a5,i2.2)')'qsurf',iq
	635	ELSE
	636	txt=tnom(iq)
	637	! Exception: there is no water vapour surface tracer
	638	if (txt.eq."h2o_vap") then
	639	write(,)"physdem0: skipping water vapour tracer"
	640	if (i_h2o_ice.eq.i_h2o_vap) then
	641	! then there is no "water ice" tracer; but still
	642	! there is some water ice on the surface
	643	write(,)" writing water ice instead"
	644	txt="h2o_ice"
	645	else
	646	! there is a "water ice" tracer which has been / will be
	647	! delt with in due time
	648	cycle
	649	endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
	650	endif ! of if (txt.eq."h2o_vap")
	651	ENDIF ! of IF (oldtracernames)
	652
	653	ierr=NF_REDEF(nid)
	654	IF (ierr.NE.NF_NOERR) THEN
	655	PRINT*, 'physdem0: Failed to swich to NetCDF define mode'
	656	CALL abort
	657	ENDIF
	658	#ifdef NC_DOUBLE
	659	ierr=NF_DEF_VAR(nid,txt,NF_DOUBLE,2,(/idim2,idim6/),nvarid)
	660	#else
	661	ierr=NF_DEF_VAR(nid,txt,NF_FLOAT,2,(/idim2,idim6/),nvarid)
	662	#endif
	663	IF (ierr.NE.NF_NOERR) THEN
	664	PRINT*, 'physdem0: Failed to define ',trim(txt)
	665	CALL abort
	666	ENDIF
	667	ierr=NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
	668	& "tracer on surface")
	669	IF (ierr.NE.NF_NOERR) THEN
	670	PRINT*, 'physdem0: Failed to define ',trim(txt),
	671	& ' title attribute'
	672	CALL abort
	673	ENDIF
	674	ierr=NF_ENDDEF(nid)
	675	IF (ierr.NE.NF_NOERR) THEN
	676	PRINT*, 'physdem0: Failed to swich out of define mode'
	677	CALL abort
	678	ENDIF
	679	ENDDO
	680	ENDIF
	681
	682	c close file
	683	ierr = NF_CLOSE(nid)
	684
	685
	686
	687	END
	688
	689
	690
	691	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	692	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	693	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	694	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	695
	696
	697
	698
[1047]	699	subroutine physdem1(filename,nsoil,ngrid,nlay,nq,
[999]	700	. phystep,time,
	701	. tsurf,tsoil,co2ice,emis,q2,qsurf)
[1036]	702	use infotrac, only: nqtot, tnom
[999]	703	implicit none
	704	c
	705	c Variables qui possèdent un axe des temps :
	706	c tsurf
	707	c tsoil
	708	c co2ice
	709	c q2
	710	c qsurf
	711	c emis (modifié dans newcondens --> co2snow)
	712	c
	713	c-------------------------------------------------------------
	714	c
	715	c write physics (re-)start data file "restartfi.nc"
	716	c
	717	c
	718	c
	719	#include "dimensions.h"
	720	#include "paramet.h"
	721	c-----------------------------------------------------------------------
	722	#include "comvert.h"
	723	#include "comgeom2.h"
	724	#include "control.h"
	725	#include "comdissnew.h"
	726	#include "logic.h"
	727	#include "ener.h"
	728	#include "netcdf.inc"
[1047]	729	!#include "dimphys.h"
[1036]	730	!#include "advtrac.h"
[999]	731	#include "callkeys.h"
	732	c
	733	INTEGER nid,iq
	734	INTEGER, parameter :: ivap=1
	735	REAL, parameter :: qsolmax= 150.0
	736	character (len=*) :: filename
	737	character (len=7) :: str7
	738
	739	REAL day_ini
[1047]	740	INTEGER,INTENT(IN) :: nsoil,ngrid,nlay,nq
[999]	741	integer ierr,nvarid
	742
	743	c
	744	REAL phystep,time
[1047]	745	! REAL champhys(ngrid)
	746	REAL tsurf(ngrid)
[999]	747
	748	INTEGER nb
	749	INTEGER edges(3),corner(3)
	750
	751
	752
	753	c
	754
	755	! EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression
	756	! EXTERNAL exner_hyb , SSUM
	757	c
	758	#include "serre.h"
	759	#include "clesph0.h"
	760	#include "fxyprim.h"
[1047]	761	!#include "comgeomfi.h"
	762	!#include "surfdat.h"
	763	!#include "comsoil.h"
[999]	764	#include "planete.h"
[1047]	765	!#include "dimradmars.h"
	766	!#include "yomaer.h"
[999]	767	#include "comcstfi.h"
	768
[1047]	769	real co2ice(ngrid),tsoil(ngrid,nsoil),emis(ngrid)
	770	real q2(ngrid,nlay+1),qsurf(ngrid,nq)
[999]	771	integer ig
	772
	773	! flag which identifies if we are using old tracer names (qsurf01,...)
	774	logical,save :: oldtracernames=.false.
	775	logical,save :: firstcall = .true.
	776	integer :: count
	777
	778	character(len=30) :: txt ! to store some text
	779	! indexes of water vapour & water ice tracers (if any):
	780	integer,save :: i_h2o_vap=0
	781	integer,save :: i_h2o_ice=0
	782	c-----------------------------------------------------------------------
	783
	784
	785	ierr = NF_OPEN(filename, NF_WRITE, nid)
	786	IF (ierr .NE. NF_NOERR) THEN
	787	PRINT*, "Pb. d ouverture "//filename
	788	CALL abort
	789	ENDIF
	790
	791	c On a single run, different files can be written with physdem1.
	792	c Therefore, get the last time index from the file itself:
	793	ierr = NF_INQ_DIMID(nid,"Time",nvarid)
	794	ierr = NF_INQ_DIMLEN(nid,nvarid,nb)
	795
	796	c Ecriture/extension de la coordonnee temps
	797
	798	nb = nb + 1
	799	ierr = NF_INQ_VARID(nid, "Time", nvarid)
	800	IF (ierr .NE. NF_NOERR) THEN
	801	print *, NF_STRERROR(ierr)
	802	print*,'Variable Time n est pas definie'
	803	CALL abort
	804	ENDIF
	805	#ifdef NC_DOUBLE
	806	ierr = NF_PUT_VAR1_DOUBLE (nid,nvarid,nb,time)
	807	#else
	808	ierr = NF_PUT_VAR1_REAL (nid,nvarid,nb,time)
	809	#endif
	810	IF (ierr .NE. NF_NOERR) THEN
	811	print*, "Erreur ecriture temps!!"
	812	print*, NF_STRERROR(ierr)
	813	ENDIF
	814	!PRINT*, "Enregistrement pour ", nb, time
	815
	816
	817	c Write the physical fields
	818
	819	! Soil temperature
	820	corner(1)=1
	821	corner(2)=1
	822	corner(3)=nb
[1047]	823	edges(1)=ngrid
[999]	824	edges(2)=nsoil
	825	edges(3)=1
	826	ierr = NF_INQ_VARID(nid, "tsoil", nvarid)
	827	IF (ierr .NE. NF_NOERR) THEN
	828	PRINT*, "Variable tsoil n est pas definie"
	829	CALL abort
	830	ENDIF
	831	#ifdef NC_DOUBLE
	832	ierr = NF_PUT_VARA_DOUBLE (nid,nvarid,corner,edges,tsoil)
	833	#else
	834	ierr = NF_PUT_VARA_REAL (nid,nvarid,corner,edges,tsoil)
	835	#endif
	836	IF (ierr .NE. NF_NOERR) THEN
	837	print*, "Erreur ecriture tsoil!!"
	838	print*, NF_STRERROR(ierr)
	839	ENDIF
	840
	841
	842	c planetary boundary layer
	843	corner(1)=1
	844	corner(2)=1
	845	corner(3)=nb
[1047]	846	edges(1)=ngrid
[999]	847	edges(2)=llm+1
	848	edges(3)=1
	849	ierr = NF_INQ_VARID(nid, "q2", nvarid)
	850	IF (ierr .NE. NF_NOERR) THEN
	851	PRINT*, "Variable q2 n est pas definie"
	852	CALL abort
	853	ENDIF
	854	#ifdef NC_DOUBLE
	855	ierr = NF_PUT_VARA_DOUBLE (nid,nvarid,corner,edges,q2)
	856	#else
	857	ierr = NF_PUT_VARA_REAL (nid,nvarid,corner,edges,q2)
	858	#endif
	859	IF (ierr .NE. NF_NOERR) THEN
	860	print*, "Erreur ecriture q2!!"
	861	print*, NF_STRERROR(ierr)
	862	ENDIF
	863
	864
	865
	866	c Following corner and egdes are the same for co2ice,tsurf,qusrf,emis and tracers
	867	corner(1)=1
	868	corner(2)=nb
[1047]	869	edges(1)=ngrid
[999]	870	edges(2)=1
	871
	872	! CO2 Ice Cover
	873	ierr = NF_INQ_VARID(nid, "co2ice", nvarid)
	874	IF (ierr .NE. NF_NOERR) THEN
	875	PRINT*, "Variable co2ice n est pas definie"
	876	CALL abort
	877	ENDIF
	878	#ifdef NC_DOUBLE
	879	ierr = NF_PUT_VARA_DOUBLE (nid,nvarid,corner,edges,co2ice)
	880	#else
	881	ierr = NF_PUT_VARA_REAL (nid,nvarid,corner,edges,co2ice)
	882	#endif
	883	IF (ierr .NE. NF_NOERR) THEN
	884	print*, "Erreur ecriture co2ice!!"
	885	print*, NF_STRERROR(ierr)
	886	ENDIF
	887
	888	! Surface temperature
	889	ierr = NF_INQ_VARID(nid, "tsurf", nvarid)
	890	IF (ierr .NE. NF_NOERR) THEN
	891	PRINT*, "Variable tsurf n est pas definie"
	892	CALL abort
	893	ENDIF
	894	#ifdef NC_DOUBLE
	895	ierr = NF_PUT_VARA_DOUBLE (nid,nvarid,corner,edges,tsurf)
	896	#else
	897	ierr = NF_PUT_VARA_REAL (nid,nvarid,corner,edges,tsurf)
	898	#endif
	899	IF (ierr .NE. NF_NOERR) THEN
	900	print*, "Erreur ecriture tsurf!!"
	901	print*, NF_STRERROR(ierr)
	902	ENDIF
	903
	904
	905	c emissivity
	906	ierr = NF_INQ_VARID(nid, "emis", nvarid)
	907	IF (ierr .NE. NF_NOERR) THEN
	908	PRINT*, "Variable emis n est pas definie"
	909	CALL abort
	910	ENDIF
	911	#ifdef NC_DOUBLE
	912	ierr = NF_PUT_VARA_DOUBLE (nid,nvarid,corner,edges,emis)
	913	#else
	914	ierr = NF_PUT_VARA_REAL (nid,nvarid,corner,edges,emis)
	915	#endif
	916	IF (ierr .NE. NF_NOERR) THEN
	917	print*, "Erreur ecriture emis!!"
	918	print*, NF_STRERROR(ierr)
	919	ENDIF
	920
	921
	922	c tracers
	923
	924	! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
	925	! qsurf02, ...)
	926	IF (firstcall) THEN
	927	count=0
[1036]	928	do iq=1,nqtot
[999]	929	txt= " "
	930	write(txt,'(a1,i2.2)')'q',iq
	931	if (txt.ne.tnom(iq)) then ! use tracer names stored in dynamics
	932	! did not find old tracer name
	933	exit ! might as well stop here
	934	else
	935	! found old tracer name
	936	count=count+1
	937	endif
	938	enddo
[1036]	939	if (count.eq.nqtot) then
[999]	940	write(,) "physdem1:tracers seem to follow old naming ",
	941	& "convention (qsurf01,qsurf02,...)"
	942	write(,) " => will work for now ... "
	943	write(,) " but you should run newstart to rename them"
	944	oldtracernames=.true.
[1036]	945	endif ! of if (count.eq.nqtot)
[999]	946	ENDIF ! of if(firstcall)
	947
	948	! If computing water cycle with ice, move surface ice
[1036]	949	! back to qsurf(nqtot)
[999]	950	IF (oldtracernames .and. water) THEN
	951	!"loop" to avoid potential out-of-bounds on arrays
[1036]	952	write(,)'physdem1: moving surface water ice to index ',nqtot
	953	do iq=nqtot,nqtot
[1047]	954	qsurf(1:ngrid,iq)=qsurf(1:ngrid,iq-1)
	955	qsurf(1:ngrid,iq-1)=0
[999]	956	enddo
	957	ENDIF
	958
	959
	960
	961	IF(nq.GE.1) THEN
	962
	963	IF (firstcall) THEN
	964	! preliminary stuff: look for water vapour & water ice tracers (if any)
	965	if (.not.oldtracernames) then
	966	do iq=1,nq
	967	if (tnom(iq).eq."h2o_vap") then
	968	i_h2o_vap=iq
	969	endif
	970	if (tnom(iq).eq."h2o_ice") then
	971	i_h2o_ice=iq
	972	endif
	973	enddo ! of iq=1,nq
	974	! handle special case of only water vapour tracer (no ice)
	975	if ((i_h2o_vap.ne.0).and.(i_h2o_ice.eq.0)) then
	976	! then the index of (surface) ice is i_h2o_vap
	977	i_h2o_ice=i_h2o_vap
	978	endif
	979	endif ! of if (.not.oldtracernames)
	980	ENDIF ! of if(firstcall)
	981	firstcall = .false.
	982
	983
	984	DO iq=1,nq
	985	IF (oldtracernames) THEN
	986	txt=" "
	987	write(txt,'(a5,i2.2)')'qsurf',iq
	988	ELSE
	989	txt=tnom(iq)
	990	! Exception: there is no water vapour surface tracer
	991	if (txt.eq."h2o_vap") then
	992	if (i_h2o_ice.eq.i_h2o_vap) then
	993	! then there is no "water ice" tracer; but still
	994	! there is some water ice on the surface
	995	txt="h2o_ice"
	996	else
	997	! there is a "water ice" tracer which has been / will be
	998	! delt with in due time
	999	cycle
	1000	endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
	1001	endif ! of if (txt.eq."h2o_vap")
	1002	ENDIF ! of IF (oldtracernames)
	1003
	1004
	1005
	1006	ierr = NF_INQ_VARID(nid, trim(txt), nvarid)
	1007	IF (ierr .NE. NF_NOERR) THEN
	1008	PRINT*, "Variable "//trim(txt)//" n est pas definie"
	1009	CALL abort
	1010	ENDIF
	1011	#ifdef NC_DOUBLE
	1012	ierr = NF_PUT_VARA_DOUBLE (nid,nvarid,corner,edges,
	1013	. qsurf(:,iq))
	1014	#else
	1015	ierr = NF_PUT_VARA_REAL (nid,nvarid,corner,edges,
	1016	. qsurf(:,iq))
	1017	#endif
	1018	IF (ierr .NE. NF_NOERR) THEN
	1019	print*, "Erreur ecriture "//trim(txt)//"!!"
	1020	print*, NF_STRERROR(ierr)
	1021	ENDIF
	1022
	1023	ENDDO ! of DO iq=1,nq
	1024	ENDIF ! of IF(nq.GE.1)
	1025
	1026	c close file
	1027	ierr = NF_CLOSE(nid)
	1028
	1029	RETURN
	1030
	1031	END

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

Download in other formats:

Original Format