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

Line
1	subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq,
2	. phystep,day_ini,time,airefi,
3	. alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
4
5	use infotrac, only: nqtot, tnom
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
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"
40	!#include "dimphys.h"
41	!#include "advtrac.h"
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
51	INTEGER,INTENT(IN) :: nsoil,ngrid,nlay,nq
52	integer ierr,idim1,idim2,idim3,idim4,idim5,idim6,nvarid
53
54	c
55	REAL phystep,time
56	REAL latfi(ngrid), lonfi(ngrid)
57	! REAL champhys(ngrid)
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"
70	!#include "comgeomfi.h"
71	!#include "surfdat.h"
72	!#include "comsoil.h"
73	#include "planete.h"
74	!#include "dimradmars.h"
75	!#include "yomaer.h"
76	#include "comcstfi.h"
77
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)
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(:)
94	CALL SCOPY(ngrid,airefi,1,area,1)
95	! note: area() is defined in comgeomfi_h
96
97
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
105	ENDDO
106
107	inertiedat(:,:)=ith(:,:) ! note inertiedat() is defined in comsoil_h
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
129	ierr = NF_DEF_DIM (nid,"physical_points",ngrid,idim2)
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
172	write(,) "physdem0: ngrid: ",ngrid
173
174	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
175	c Fill control array tab_cntrl(:) with parameters for this run
176	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
177	c Informations on the physics grid
178	tab_cntrl(1) = float(ngrid) ! number of nodes on physics grid
179	tab_cntrl(2) = float(nlay) ! number of atmospheric layers
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
268	! write mid-layer depths mlayer() !known from comsoil_h
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
472	! surface roughness length (z0 is a common in surfdat_h)
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
594	do iq=1,nqtot
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
605	if (count.eq.nqtot) then
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.
611	endif ! of if (count.eq.nqtot)
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
699	subroutine physdem1(filename,nsoil,ngrid,nlay,nq,
700	. phystep,time,
701	. tsurf,tsoil,co2ice,emis,q2,qsurf)
702	use infotrac, only: nqtot, tnom
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"
729	!#include "dimphys.h"
730	!#include "advtrac.h"
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
740	INTEGER,INTENT(IN) :: nsoil,ngrid,nlay,nq
741	integer ierr,nvarid
742
743	c
744	REAL phystep,time
745	! REAL champhys(ngrid)
746	REAL tsurf(ngrid)
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"
761	!#include "comgeomfi.h"
762	!#include "surfdat.h"
763	!#include "comsoil.h"
764	#include "planete.h"
765	!#include "dimradmars.h"
766	!#include "yomaer.h"
767	#include "comcstfi.h"
768
769	real co2ice(ngrid),tsoil(ngrid,nsoil),emis(ngrid)
770	real q2(ngrid,nlay+1),qsurf(ngrid,nq)
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
823	edges(1)=ngrid
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
846	edges(1)=ngrid
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
869	edges(1)=ngrid
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
928	do iq=1,nqtot
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
939	if (count.eq.nqtot) then
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.
945	endif ! of if (count.eq.nqtot)
946	ENDIF ! of if(firstcall)
947
948	! If computing water cycle with ice, move surface ice
949	! back to qsurf(nqtot)
950	IF (oldtracernames .and. water) THEN
951	!"loop" to avoid potential out-of-bounds on arrays
952	write(,)'physdem1: moving surface water ice to index ',nqtot
953	do iq=nqtot,nqtot
954	qsurf(1:ngrid,iq)=qsurf(1:ngrid,iq-1)
955	qsurf(1:ngrid,iq-1)=0
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