source: trunk/LMDZ.MARS/libf/phymars/physdem1.F @ 308

Last change on this file since 308 was 224, checked in by emillour, 13 years ago

Mars GCM:

Implemented using 'z0' roughness length map (important: 'z0' reference

field is in datafile surface.nc, which has also been updated).

  • made z0 a z0(ngridmx) array and moved 'z0' from 'planete.h' to 'surfdat.h'; added a 'z0_default' (common in surfdat.h) corresponding to the 'control' array value (contole(19) in startfi.nc).
  • adapted 'tabfi.F' to use 'z0_default'.
  • adapted 'phyetat0.F' to look for a 'z0' field in startfi.nc. If 'z0' is not found in the startfi.nc file, then the uniform default value (z0_default) is used.
  • modified 'physdem1.F' to write 'z0' field to restart.nc
  • adapted use of z0() in 'physiq.F' (diagnostic computation of surface stress), 'vdifc.F' and 'vdif_cd.F'.
  • adapted 'dustdevil.F' to use 'z0_default'.
  • 'testphys1d.F' now uses 'z0_default', and the value to use can be set in run.def (with "z0=TheValueYouWant?").
  • modified 'datareadnc.F' to load reference map of 'z0' from surface.nc, and added a 'z0' option in 'newstart.F' to force a uniform value of z0. Note that the use of the z0 map is automatic when using newstart, but only when it loads a start_archive.nc file.
File size: 21.1 KB
Line 
1      subroutine physdem1(filename,lonfi,latfi,nsoil,nq,
2     .                   phystep,day_ini,
3     .                   time,tsurf,tsoil,co2ice,emis,q2,qsurf,
4     .                   airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
5
6      implicit none
7c-------------------------------------------------------------
8c
9c create physics (re-)start data file "restartfi.nc"
10c
11c
12c
13#include "dimensions.h"
14#include "paramet.h"
15c-----------------------------------------------------------------------
16#include "comvert.h"
17#include "comgeom2.h"
18#include "control.h"
19#include "comdissnew.h"
20#include "logic.h"
21#include "ener.h"
22#include "netcdf.inc"
23#include "dimphys.h"
24#include"advtrac.h"
25#include"callkeys.h"
26c
27      INTEGER nid,iq
28      INTEGER, parameter :: ivap=1
29      REAL, parameter :: qsolmax= 150.0
30      character (len=*) :: filename
31      character (len=7) :: str7
32
33      REAL day_ini
34      INTEGER nsoil,nq
35      integer ierr,idim1,idim2,idim3,idim4,idim5,nvarid
36
37c
38      REAL phystep,time
39      REAL latfi(ngridmx), lonfi(ngridmx)
40!      REAL champhys(ngridmx)
41      REAL tsurf(ngridmx)
42      INTEGER length
43      PARAMETER (length=100)
44      REAL tab_cntrl(length)
45
46c
47
48!      EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression
49!      EXTERNAL exner_hyb , SSUM
50c
51#include "serre.h"
52#include "clesph0.h"
53#include "fxyprim.h"
54#include "comgeomfi.h"
55#include "surfdat.h"
56#include "comsoil.h"
57#include "planete.h"
58#include "dimradmars.h"
59#include "yomaer.h"
60#include "comcstfi.h"
61
62      real co2ice(ngridmx),tsoil(ngridmx,nsoil),emis(ngridmx)
63      real q2(ngridmx, llm+1),qsurf(ngridmx,nq)
64      real airefi(ngridmx)
65      real alb(ngridmx),ith(ngridmx,nsoil)
66      real pzmea(ngridmx),pzstd(ngridmx)
67      real pzsig(ngridmx),pzgam(ngridmx),pzthe(ngridmx)
68      integer ig
69
70! flag which identifies if we are using old tracer names (qsurf01,...)
71      logical :: oldtracernames=.false.
72      integer :: count
73      character(len=30) :: txt ! to store some text
74! indexes of water vapour & water ice tracers (if any):
75      integer :: i_h2o_vap=0
76      integer :: i_h2o_ice=0
77c-----------------------------------------------------------------------
78
79      ! copy airefi(:) to area(:)
80      CALL SCOPY(ngridmx,airefi,1,area,1)
81      ! note: area() is defined in comgeomfi.h
82
83      DO ig=1,ngridmx
84         albedodat(ig)=alb(ig) ! note: albedodat() is defined in surfdat.h
85         zmea(ig)=pzmea(ig) ! note: zmea() is defined in surfdat.h
86         zstd(ig)=pzstd(ig) ! note: zstd() is defined in surfdat.h
87         zsig(ig)=pzsig(ig) ! note: zsig() is defined in surfdat.h
88         zgam(ig)=pzgam(ig) ! note: zgam() is defined in surfdat.h
89         zthe(ig)=pzthe(ig) ! note: zthe() is defined in surfdat.h
90      ENDDO
91
92      inertiedat(:,:)=ith(:,:) ! note inertiedat() is defined in comsoil.h
93c
94c  things to store in the physics start file:
95c
96      ierr = NF_CREATE(adjustl(filename),NF_CLOBBER, nid)
97      IF (ierr.NE.NF_NOERR) THEN
98        WRITE(6,*)'physdem1: Problem creating file ',adjustl(filename)
99        write(6,*) NF_STRERROR(ierr)
100        CALL ABORT
101      ENDIF
102c
103      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 18,
104     .                       "Physics start file")
105c
106      ierr = NF_DEF_DIM (nid,"index",length,idim1)
107      if (ierr.ne.NF_NOERR) then
108        WRITE(6,*)'physdem1: Problem defining index dimension'
109        write(6,*) NF_STRERROR(ierr)
110        call abort
111      endif
112c
113      ierr = NF_DEF_DIM (nid,"physical_points",ngridmx,idim2)
114      if (ierr.ne.NF_NOERR) then
115        WRITE(6,*)'physdem1: Problem defining physical_points dimension'
116        write(6,*) NF_STRERROR(ierr)
117        call abort
118      endif
119c
120      ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoil,idim3)
121      if (ierr.ne.NF_NOERR) then
122      WRITE(6,*)'physdem1: Problem defining subsurface_layers dimension'
123        write(6,*) NF_STRERROR(ierr)
124        call abort
125      endif
126c
127!      ierr = NF_DEF_DIM (nid,"nlayer+1",llm+1,idim4)
128      ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
129      if (ierr.ne.NF_NOERR) then
130        WRITE(6,*)'physdem1: Problem defining nlayer+1 dimension'
131        write(6,*) NF_STRERROR(ierr)
132        call abort
133      endif
134c
135      ierr = NF_DEF_DIM (nid,"number_of_advected_fields",nq,idim5)
136      if (ierr.ne.NF_NOERR) then
137        WRITE(6,*)'physdem1: Problem defining advected fields dimension'
138        WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
139        write(6,*) NF_STRERROR(ierr)
140      endif
141
142      ierr = NF_ENDDEF(nid) ! exit NetCDF define mode
143
144c clear tab_cntrl(:) array
145      DO ierr = 1, length
146         tab_cntrl(ierr) = 0.0
147      ENDDO
148
149      write(*,*) "physdem1: ngridmx: ",ngridmx
150
151ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
152c Fill control array tab_cntrl(:) with paramleters for this run
153ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
154c Informations on the physics grid
155      tab_cntrl(1) = float(ngridmx)  ! number of nodes on physics grid
156      tab_cntrl(2) = float(nlayermx) ! number of atmospheric layers
157      tab_cntrl(3) = day_ini + int(time)         ! initial day
158      tab_cntrl(4) = time -int(time)            ! initiale time of day
159
160c Informations about Mars, used by dynamics and physics
161      tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
162      tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
163      tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
164      tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
165      tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
166      tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775
167
168      tab_cntrl(11) = phystep  ! time step in the physics
169      tab_cntrl(12) = 0.
170      tab_cntrl(13) = 0.
171
172c Informations about Mars, only for physics
173      tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
174      tab_cntrl(15) = periheli  ! min. Sun-Mars distance (Mkm) ~206.66
175      tab_cntrl(16) = aphelie   ! max. SUn-Mars distance (Mkm) ~249.22
176      tab_cntrl(17) = peri_day  ! date of perihelion (sols since N. spring)
177      tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98
178
179c Boundary layer and turbulence
180      tab_cntrl(19) = z0_default   ! default surface roughness (m) ~0.01
181      tab_cntrl(20) = lmixmin   ! mixing length ~100
182      tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
183
184c Optical properties of polar caps and ground emissivity
185      tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
186      tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
187      tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
188      tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
189      tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
190      tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
191      tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
192      tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
193      tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
194
195c dust aerosol properties
196      tab_cntrl(27) = tauvis      ! mean visible optical depth
197
198      tab_cntrl(28) = 0.
199      tab_cntrl(29) = 0.
200      tab_cntrl(30) = 0.
201
202! Soil properties:
203      tab_cntrl(35) = volcapa ! soil volumetric heat capacity
204     
205c
206!      write(*,*) "physdem1: tab_cntrl():",tab_cntrl
207     
208      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
209      IF (ierr.NE.NF_NOERR) THEN
210         PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
211         CALL abort
212      ENDIF
213      ! define variable
214#ifdef NC_DOUBLE
215      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
216#else
217      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
218#endif
219      IF (ierr.NE.NF_NOERR) THEN
220         PRINT*, 'physdem1: Failed to define controle'
221         CALL abort
222      ENDIF
223      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
224     .                        "Control parameters")
225      IF (ierr.NE.NF_NOERR) THEN
226         PRINT*, 'physdem1: Failed to define controle title attribute'
227         CALL abort
228      ENDIF
229      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
230      IF (ierr.NE.NF_NOERR) THEN
231         PRINT*, 'physdem1: Failed to swich out of NetCDF define mode'
232         CALL abort
233      ENDIF
234      ! write variable
235#ifdef NC_DOUBLE
236      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
237#else
238      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
239#endif
240      IF (ierr.NE.NF_NOERR) THEN
241         PRINT*, 'physdem1: Failed to write controle data'
242         CALL abort
243      ENDIF
244
245! write mid-layer depths mlayer() !known from comsoil.h
246
247      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
248      ! define variable
249#ifdef NC_DOUBLE
250      ierr = NF_DEF_VAR (nid,"soildepth",NF_DOUBLE,1,idim3,nvarid)
251#else
252      ierr = NF_DEF_VAR (nid,"soildepth",NF_FLOAT,1,idim3,nvarid)
253#endif
254      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
255     .                        "Soil mid-layer depth")
256      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
257      ! write variable
258#ifdef NC_DOUBLE
259      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
260#else
261      ierr = NF_PUT_VAR_REAL (nid,nvarid,mlayer)
262#endif
263
264c
265
266      ierr = NF_REDEF (nid)
267#ifdef NC_DOUBLE
268      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
269#else
270      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
271#endif
272      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 26,
273     .                        "Longitudes of physics grid")
274      ierr = NF_ENDDEF(nid)
275
276#ifdef NC_DOUBLE
277      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
278#else
279      ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
280#endif
281
282c
283
284      ierr = NF_REDEF (nid)
285#ifdef NC_DOUBLE
286      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
287#else
288      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
289#endif
290      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
291     .                        "Latitudes of physics grid")
292      ierr = NF_ENDDEF(nid)
293#ifdef NC_DOUBLE
294      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
295#else
296      ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
297#endif
298
299c
300
301      ierr = NF_REDEF (nid)
302#ifdef NC_DOUBLE
303      ierr = NF_DEF_VAR (nid, "area", NF_DOUBLE, 1, idim2,nvarid)
304#else
305      ierr = NF_DEF_VAR (nid, "area", NF_FLOAT, 1, idim2,nvarid)
306#endif
307      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
308     .                        "Mesh area")
309      ierr = NF_ENDDEF(nid)
310#ifdef NC_DOUBLE
311      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
312#else
313      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
314#endif
315
316c
317
318      ierr = NF_REDEF (nid)
319#ifdef NC_DOUBLE
320      ierr = NF_DEF_VAR (nid, "phisfi", NF_DOUBLE, 1, idim2,nvarid)
321#else
322      ierr = NF_DEF_VAR (nid, "phisfi", NF_FLOAT, 1, idim2,nvarid)
323#endif
324      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
325     .                        "Geopotential at the surface")
326      ierr = NF_ENDDEF(nid)
327#ifdef NC_DOUBLE
328      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phisfi)
329#else
330      ierr = NF_PUT_VAR_REAL (nid,nvarid,phisfi)
331#endif
332
333c
334
335      ierr = NF_REDEF (nid)
336#ifdef NC_DOUBLE
337      ierr = NF_DEF_VAR (nid, "albedodat", NF_DOUBLE, 1, idim2,nvarid)
338#else
339      ierr = NF_DEF_VAR (nid, "albedodat", NF_FLOAT, 1, idim2,nvarid)
340#endif
341      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
342     .                        "Albedo of bare ground")
343      ierr = NF_ENDDEF(nid)
344#ifdef NC_DOUBLE
345      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedodat)
346#else
347      ierr = NF_PUT_VAR_REAL (nid,nvarid,albedodat)
348#endif
349
350c
351c   some data for Francois Lott's programs
352c
353
354      ierr = NF_REDEF (nid)
355#ifdef NC_DOUBLE
356      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
357#else
358      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
359#endif
360      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
361     .                        "Relief: mean relief")
362      ierr = NF_ENDDEF(nid)
363#ifdef NC_DOUBLE
364      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
365#else
366      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
367#endif
368c
369      ierr = NF_REDEF (nid)
370#ifdef NC_DOUBLE
371      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
372#else
373      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
374#endif
375      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
376     .                        "Relief: standard deviation")
377      ierr = NF_ENDDEF(nid)
378#ifdef NC_DOUBLE
379      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
380#else
381      ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
382#endif
383c
384      ierr = NF_REDEF (nid)
385#ifdef NC_DOUBLE
386      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
387#else
388      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
389#endif
390      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
391     .                        "Relief: sigma parameter")
392      ierr = NF_ENDDEF(nid)
393#ifdef NC_DOUBLE
394      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
395#else
396      ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
397#endif
398c
399      ierr = NF_REDEF (nid)
400#ifdef NC_DOUBLE
401      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
402#else
403      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
404#endif
405      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
406     .                        "Relief: gamma parameter")
407      ierr = NF_ENDDEF(nid)
408#ifdef NC_DOUBLE
409      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
410#else
411      ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
412#endif
413c
414      ierr = NF_REDEF (nid)
415#ifdef NC_DOUBLE
416      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
417#else
418      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
419#endif
420      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
421     .                        "Relief: theta parameter")
422      ierr = NF_ENDDEF(nid)
423#ifdef NC_DOUBLE
424      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
425#else
426      ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
427#endif
428
429c Write the physical fields
430
431! CO2 Ice Cover
432
433      ierr = NF_REDEF (nid)
434#ifdef NC_DOUBLE
435      ierr = NF_DEF_VAR (nid, "co2ice", NF_DOUBLE, 1, idim2,nvarid)
436#else
437      ierr = NF_DEF_VAR (nid, "co2ice", NF_FLOAT, 1, idim2,nvarid)
438#endif
439      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 13,
440     .                        "CO2 ice cover")
441      ierr = NF_ENDDEF(nid)
442#ifdef NC_DOUBLE
443      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,co2ice)
444#else
445      ierr = NF_PUT_VAR_REAL (nid,nvarid,co2ice)
446#endif
447
448! Soil Thermal inertia
449
450      ierr = NF_REDEF (nid)
451#ifdef NC_DOUBLE
452      ierr = NF_DEF_VAR (nid,"inertiedat",NF_DOUBLE,
453     &                   2,(/idim2,idim3/),nvarid)
454#else
455      ierr = NF_DEF_VAR (nid,"inertiedat",NF_FLOAT,
456     &                   2,(/idim2,idim3/),nvarid)
457#endif
458      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 20,
459     .                        "Soil thermal inertia")
460      ierr = NF_ENDDEF(nid)
461#ifdef NC_DOUBLE
462      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,inertiedat)
463#else
464      ierr = NF_PUT_VAR_REAL (nid,nvarid,inertiedat)
465#endif
466
467!  Surface temperature
468
469      ierr = NF_REDEF (nid)
470#ifdef NC_DOUBLE
471      ierr = NF_DEF_VAR (nid, "tsurf", NF_DOUBLE, 1, idim2,nvarid)
472#else
473      ierr = NF_DEF_VAR (nid, "tsurf", NF_FLOAT, 1, idim2,nvarid)
474#endif
475      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
476     .                        "Surface temperature")
477      ierr = NF_ENDDEF(nid)
478#ifdef NC_DOUBLE
479      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsurf)
480#else
481      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsurf)
482#endif
483
484! Soil temperature
485
486      ierr = NF_REDEF (nid)
487#ifdef NC_DOUBLE
488      ierr = NF_DEF_VAR (nid,"tsoil",NF_DOUBLE,2,(/idim2,idim3/),nvarid)
489#else
490!      ierr = NF_DEF_VAR (nid, "tsoil", NF_FLOAT, 2, idim2,nvarid)
491      ierr = NF_DEF_VAR (nid,"tsoil",NF_FLOAT,2,(/idim2,idim3/),nvarid)
492#endif
493      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 16,
494     .                        "Soil temperature")
495      ierr = NF_ENDDEF(nid)
496#ifdef NC_DOUBLE
497      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsoil)
498#else
499      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsoil)
500#endif
501
502c emissivity
503
504      ierr = NF_REDEF (nid)
505#ifdef NC_DOUBLE
506      ierr = NF_DEF_VAR (nid, "emis", NF_DOUBLE, 1, idim2,nvarid)
507#else
508      ierr = NF_DEF_VAR (nid, "emis", NF_FLOAT, 1, idim2,nvarid)
509#endif
510      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 18,
511     .                        "Surface emissivity")
512      ierr = NF_ENDDEF(nid)
513#ifdef NC_DOUBLE
514      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,emis)
515#else
516      ierr = NF_PUT_VAR_REAL (nid,nvarid,emis)
517#endif
518
519! surface roughness length (z0 is a common in surfdat.h)
520
521      ierr = NF_REDEF (nid)
522#ifdef NC_DOUBLE
523      ierr = NF_DEF_VAR (nid, "z0", NF_DOUBLE, 1, idim2,nvarid)
524#else
525      ierr = NF_DEF_VAR (nid, "z0", NF_FLOAT, 1, idim2,nvarid)
526#endif
527      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 24,
528     .                        "Surface roughness length")
529      ierr = NF_ENDDEF(nid)
530#ifdef NC_DOUBLE
531      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,z0)
532#else
533      ierr = NF_PUT_VAR_REAL (nid,nvarid,z0)
534#endif
535
536
537c planetary boundary layer
538
539      ierr = NF_REDEF (nid)
540#ifdef NC_DOUBLE
541      ierr = NF_DEF_VAR (nid,"q2",NF_DOUBLE,2,(/idim2,idim4/),nvarid)
542#else
543      ierr = NF_DEF_VAR (nid,"q2",NF_FLOAT, 2,(/idim2,idim4/),nvarid)
544#endif
545      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
546     .                        "pbl wind variance")
547      ierr = NF_ENDDEF(nid)
548#ifdef NC_DOUBLE
549      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q2)
550#else
551      ierr = NF_PUT_VAR_REAL (nid,nvarid,q2)
552#endif
553      IF (ierr.NE.NF_NOERR) THEN
554        PRINT*, 'physdem1: Failed to write q2'
555        CALL abort
556      ENDIF
557
558c tracers
559
560! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
561!                    qsurf02, ...)
562      count=0
563      do iq=1,nqmx
564        txt= " "
565        write(txt,'(a1,i2.2)')'q',iq
566        if (txt.ne.tnom(iq)) then ! use tracer names stored in dynamics
567          ! did not find old tracer name
568          exit ! might as well stop here
569        else
570          ! found old tracer name
571          count=count+1
572        endif
573      enddo
574      if (count.eq.nqmx) then
575        write(*,*) "physdem1:tracers seem to follow old naming ",
576     &             "convention (qsurf01,qsurf02,...)"
577        write(*,*) "   => will work for now ... "
578        write(*,*) "      but you should run newstart to rename them"
579        oldtracernames=.true.
580        ! Moreover, if computing water cycle with ice, move surface ice
581        ! back to qsurf(nqmx)
582        IF (water) THEN
583          !"loop" to avoid potential out-of-bounds on arrays
584          write(*,*)'physdem1: moving surface water ice to index ',nqmx
585          do iq=nqmx,nqmx
586          qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq-1)
587          qsurf(1:ngridmx,iq-1)=0
588          enddo
589        ENDIF
590      endif ! of if (count.eq.nqmx)
591
592      IF(nq.GE.1) THEN
593! preliminary stuff: look for water vapour & water ice tracers (if any)
594        if (.not.oldtracernames) then
595         do iq=1,nq
596           if (tnom(iq).eq."h2o_vap") then
597             i_h2o_vap=iq
598           endif
599           if (tnom(iq).eq."h2o_ice") then
600             i_h2o_ice=iq
601           endif
602         enddo ! of iq=1,nq
603         ! handle special case of only water vapour tracer (no ice)
604         if ((i_h2o_vap.ne.0).and.(i_h2o_ice.eq.0)) then
605          ! then the index of (surface) ice is i_h2o_vap
606          i_h2o_ice=i_h2o_vap
607         endif
608        endif ! of if (.not.oldtracernames)
609
610         DO iq=1,nq
611           IF (oldtracernames) THEN
612             txt=" "
613             write(txt,'(a5,i2.2)')'qsurf',iq
614           ELSE
615             txt=tnom(iq)
616             ! Exception: there is no water vapour surface tracer
617             if (txt.eq."h2o_vap") then
618               write(*,*)"physdem1: skipping water vapour tracer"
619               if (i_h2o_ice.eq.i_h2o_vap) then
620               ! then there is no "water ice" tracer; but still
621               ! there is some water ice on the surface
622                 write(*,*)"          writting water ice instead"
623                 txt="h2o_ice"
624               else
625               ! there is a "water ice" tracer which has been / will be
626               ! delt with in due time
627                 cycle
628               endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
629             endif ! of if (txt.eq."h2o_vap")
630           ENDIF ! of IF (oldtracernames)
631
632           ierr=NF_REDEF(nid)
633           IF (ierr.NE.NF_NOERR) THEN
634             PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
635             CALL abort
636           ENDIF
637#ifdef NC_DOUBLE
638           ierr=NF_DEF_VAR(nid,txt,NF_DOUBLE,1,idim2,nvarid)
639#else
640           ierr=NF_DEF_VAR(nid,txt,NF_FLOAT,1,idim2,nvarid)
641#endif
642           IF (ierr.NE.NF_NOERR) THEN
643             PRINT*, 'physdem1: Failed to define ',trim(txt)
644             CALL abort
645           ENDIF
646           ierr=NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
647     &                        "tracer on surface")
648           IF (ierr.NE.NF_NOERR) THEN
649             PRINT*, 'physdem1: Failed to define ',trim(txt),
650     &               ' title attribute'
651             CALL abort
652           ENDIF
653           ierr=NF_ENDDEF(nid)
654           IF (ierr.NE.NF_NOERR) THEN
655             PRINT*, 'physdem1: Failed to swich out of define mode'
656             CALL abort
657           ENDIF
658           
659#ifdef NC_DOUBLE
660            ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,qsurf(1,iq))
661#else
662            ierr=NF_PUT_VAR_REAL (nid,nvarid,qsurf(1,iq))
663#endif
664           IF (ierr.NE.NF_NOERR) THEN
665             PRINT*, 'physdem1: Failed to write ',trim(txt)
666             CALL abort
667           ENDIF
668         ENDDO ! of DO iq=1,nq
669      ENDIF ! of IF(nq.GE.1)
670
671c close file
672      ierr = NF_CLOSE(nid)
673
674      RETURN
675
676      END
Note: See TracBrowser for help on using the repository browser.