source: trunk/LMDZ.GENERIC/libf/phystd/physdem1.F @ 764

Last change on this file since 764 was 726, checked in by jleconte, 13 years ago

17/07/2012 == JL for LK

  • Generalization of aerosol scheme:
    • any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
    • addition of a module with the id numbers for aerosols (aerosol_mod.F90).
    • initialization of aerosols id numbers in iniaerosol.F90
    • compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
  • For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

  • If starting from an old start file, recreate start file with the q=0 option in newstart.e.
  • update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

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