source: trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/lect_start_archive.F @ 3341

Last change on this file since 3341 was 3100, checked in by bhatnags, 13 months ago

Generic-PCM

This commit updates the slab ocean module to a parallelisable dynamic slab ocean module. This is particularly relevant if you want to be able to use oceanic heat transport in parallel mode.

It has the following features:

(a) Computes sea ice creation and evolution.
(b) Snow has thermodynamic properties.
(c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
(d) Can be used in parallel mode.

Required callphys.def flags:
The slab ocean and its dependencies can be activated with the following flags (already added to deftank):
## Ocean options
## ~
# Model slab-ocean (Main flag for slab ocean)
ok_slab_ocean = .true.
# The following flags can only be set to true if ok_slab_ocean is true
# Ekman transport
slab_ekman = .true.
# Ekman zonal advection
slab_ekman_zonadv = .true.
# Horizontal diffusion (default coef_hdiff=25000., can be changed)
slab_hdiff = .true.
# Slab-ocean timestep (in physics timesteps)
cpl_pas = 1
# Gent-McWilliams? Scheme (can only be true if slab_ekman is true)
slab_gm = .true.

Notes:
In the current state, the model crashes if moistadjustment = .true. Unsure whether this is due to the slab or is an inherent issue with moistadj (under investigation).

SB and EM
File size: 51.3 KB
Line 
1      SUBROUTINE lect_start_archive(ngrid,nlayer,
2     &     date,tsurf,tsoil,emis,q2,
3     &     t,ucov,vcov,ps,h,phisold_newgrid,
4     &     q,qsurf,surfith,nid,
5     &     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice,
6     &     du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
7
8      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
9      USE tracer_h, ONLY: igcm_co2_ice
10      USE infotrac, ONLY: tname, nqtot
11!      USE slab_ice_h, ONLY: noceanmx
12      USE ocean_slab_mod, ONLY: nslay
13      USE callkeys_mod, ONLY: ok_slab_ocean
14      USE comvert_mod, ONLY: ap,bp,aps,bps,preff
15      USE comconst_mod, ONLY: kappa,g,pi
16
17c=======================================================================
18c
19c    Routine to load variables from the "start_archive.nc" file
20c
21c=======================================================================
22
23      implicit none
24
25      include "dimensions.h"
26      include "paramet.h"
27      include "comgeom2.h"
28      include "netcdf.inc"
29
30c=======================================================================
31c   Declarations
32c=======================================================================
33
34      INTEGER,INTENT(IN) :: ngrid, nlayer
35
36c Old variables dimensions (from file)
37c------------------------------------
38      INTEGER   imold,jmold,lmold,nsoilold,nqold
39
40
41c Variables pour les lectures des fichiers "ini"
42c--------------------------------------------------
43!      INTEGER sizei,
44      integer timelen,dimid
45!      INTEGER length
46!      parameter (length = 100)
47      INTEGER tab0
48      INTEGER isoil,iq,iqmax
49      CHARACTER*2   str2
50
51      REAL,INTENT(OUT) :: date
52      INTEGER   memo
53!      character (len=50) :: tmpname
54
55c Variable histoire
56c------------------
57      REAL,INTENT(OUT) :: vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
58      REAL,INTENT(OUT) :: h(iip1,jjp1,llm),ps(iip1,jjp1)
59      REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot)
60
61c Physique sur grille scalaire
62c----------------------------
63
64c variable physique
65c------------------
66      REAL,INTENT(OUT) :: tsurf(ngrid) ! surface temperature
67      REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx) ! soil temperature
68      REAL co2ice(ngrid) ! CO2 ice layer
69      REAL,INTENT(OUT) :: emis(ngrid)
70      REAL,INTENT(OUT) :: q2(ngrid,llm+1),qsurf(ngrid,nqtot)
71      REAL,INTENT(OUT) :: tslab(ngrid,nslay)
72      REAL ,INTENT(OUT) ::rnat(ngrid),pctsrf_sic(ngrid)
73      REAL,INTENT(OUT) :: tsea_ice(ngrid),sea_ice(ngrid)
74c     REAL phisfi(ngrid)
75      REAL,INTENT(OUT):: du_nonoro_gwd(ngrid,llm)
76      REAL,INTENT(OUT):: dv_nonoro_gwd(ngrid,llm)
77      REAL,INTENT(OUT):: east_gwstress(ngrid,llm)
78      REAL,INTENT(OUT):: west_gwstress(ngrid,llm)
79
80      INTEGER i,j,l
81      INTEGER,INTENT(IN) :: nid
82      INTEGER :: nvarid
83c     REAL year_day,periheli,aphelie,peri_day
84c     REAL obliquit,z0,emin_turb,lmixmin
85c     REAL emissiv,emisice(2),albedice(2)
86c     REAL iceradius(2) , dtemisice(2)
87
88      integer ierr
89      integer, dimension(4) :: start,count
90
91c Variable nouvelle grille naturelle au point scalaire
92c------------------------------------------------------
93      real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
94      REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1)
95      REAL,INTENT(OUT) :: t(iip1,jjp1,llm)
96      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
97      real inertiedatS(iip1,jjp1,nsoilmx)
98      real co2iceS(iip1,jjp1)
99      real emisS(iip1,jjp1)
100      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
101      real tslabS(iip1,jjp1,nslay)
102      real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1)
103      real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1)
104      real du_nonoro_gwdS(iip1,jjp1,llm),dv_nonoro_gwdS(iip1,jjp1,llm)
105      real east_gwstressS(iip1,jjp1,llm),west_gwstressS(iip1,jjp1,llm)
106
107      real ptotal, co2icetotal
108
109c Var intermediaires : vent naturel, mais pas coord scalaire
110c-----------------------------------------------------------
111      real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
112
113
114c Variable de l'ancienne grille
115c---------------------------------------------------------
116
117      real, dimension(:), allocatable :: timelist
118      real, dimension(:), allocatable :: rlonuold, rlatvold
119      real, dimension(:), allocatable :: rlonvold, rlatuold
120      real, dimension(:), allocatable :: apsold,bpsold
121      real, dimension(:), allocatable :: mlayerold
122      real, dimension(:,:,:), allocatable :: uold,vold,told,q2old
123      real, dimension(:,:,:), allocatable :: tsoilold,qsurfold
124      real, dimension(:,:,:),allocatable :: tsoiloldnew
125! tsoiloldnew: old soil values, but along new subterranean grid
126      real, dimension(:,:,:), allocatable :: inertiedatold
127! inertiedatoldnew: old inertia values, but along new subterranean grid
128      real, dimension(:,:,:), allocatable :: inertiedatoldnew
129      real, dimension(:,:), allocatable :: psold,phisold
130      real, dimension(:,:), allocatable :: co2iceold
131      real, dimension(:,:), allocatable :: tsurfold
132      real, dimension(:,:), allocatable :: emisold
133      real, dimension(:,:,:,:), allocatable :: qold
134      real, dimension(:,:,:), allocatable :: tslabold
135      real, dimension(:,:), allocatable :: rnatold,pctsrf_sicold
136      real, dimension(:,:), allocatable :: tsea_iceold,sea_iceold
137      real,allocatable :: du_nonoro_gwdold(:,:,:)
138      real,allocatable :: dv_nonoro_gwdold(:,:,:)
139      real,allocatable :: east_gwstressold(:,:,:)
140      real,allocatable :: west_gwstressold(:,:,:)
141
142      real tab_cntrl(100)
143
144      real ptotalold, co2icetotalold
145
146      logical :: olddepthdef=.false. ! flag
147! olddepthdef=.true. if soil depths are in 'old' (unspecified) format
148      logical :: depthinterpol=.false. ! flag
149! depthinterpol=.true. if interpolation will be requiered
150      logical :: therminertia_3D=.true. ! flag
151! therminertia_3D=.true. if thermal inertia is 3D and read from datafile
152c Variable intermediaires iutilise pour l'extrapolation verticale
153c----------------------------------------------------------------
154      real, dimension(:,:,:), allocatable :: var,varp1
155      real, dimension(:), allocatable :: oldgrid, oldval
156      real, dimension(:), allocatable :: newval
157
158      real,intent(out) :: surfith(iip1,jjp1) ! surface thermal inertia
159      ! surface thermal inertia at old horizontal grid resolution
160      real, dimension(:,:), allocatable :: surfithold
161
162      character(len=30) :: txt ! to store some text
163
164c=======================================================================
165
166! 0. Preliminary stuff
167
168
169
170!-----------------------------------------------------------------------
171! 1. Read data dimensions (i.e. size and length)
172!-----------------------------------------------------------------------
173
174! 1.2 Read the various dimension lengths of data in file
175
176      ierr= NF_INQ_DIMID(nid,"Time",dimid)
177      if (ierr.ne.NF_NOERR) then
178         ierr= NF_INQ_DIMID(nid,"temps",dimid)
179      endif
180      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
181
182      ierr= NF_INQ_DIMID(nid,"latitude",dimid)
183      if (ierr.ne.NF_NOERR) then
184         ierr= NF_INQ_DIMID(nid,"rlatu",dimid)
185      endif
186      ierr= NF_INQ_DIMLEN(nid,dimid,jmold)
187      jmold=jmold-1
188
189      ierr= NF_INQ_DIMID(nid,"longitude",dimid)
190      if (ierr.ne.NF_NOERR) then
191         ierr= NF_INQ_DIMID(nid,"rlonv",dimid)
192      endif
193      ierr= NF_INQ_DIMLEN(nid,dimid,imold)
194      imold=imold-1
195
196      ierr= NF_INQ_DIMID(nid,"altitude",dimid)
197      if (ierr.ne.NF_NOERR) then
198         ierr= NF_INQ_DIMID(nid,"sig_s",dimid)
199      endif
200      ierr= NF_INQ_DIMLEN(nid,dimid,lmold)
201
202      nqold=0
203      do
204         write(str2,'(i2.2)') nqold+1
205         ierr= NF_INQ_VARID(nid,'q'//str2,dimid)
206!        write(*,*) 'q'//str2
207         if (ierr.eq.NF_NOERR) then
208            nqold=nqold+1
209         else
210            exit
211         endif
212      enddo
213
214! 1.2.1 find out the # of subsurface_layers
215      nsoilold=0 !dummy initialisation
216      ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid)
217      if (ierr.eq.NF_NOERR) then
218        ierr=NF_INQ_DIMLEN(nid,dimid,nsoilold)
219        if (ierr.ne.NF_NOERR) then
220         write(*,*)'lec_start_archive: ',
221     &              'Failed reading subsurface_layers length'
222        endif
223      else
224        write(*,*)"lec_start_archive: did not find subsurface_layers"
225      endif
226
227      if (nsoilold.eq.0) then ! 'old' archive format;
228      ! must use Tg//str2 fields to compute nsoilold
229      write(*,*)"lec_start_archive: building nsoilold from Tg fields"
230        do
231         write(str2,'(i2.2)') nsoilold+1
232         ierr=NF_INQ_VARID(nid,'Tg'//str2,dimid)
233         if (ierr.eq.NF_NOERR) then
234          nsoilold=nsoilold+1
235         else
236          exit
237         endif
238        enddo
239      endif
240
241
242      if (nsoilold.ne.nsoilmx) then ! interpolation will be required
243        depthinterpol=.true.
244      endif
245
246! 1.3 Report dimensions
247     
248      write(*,*) "Start_archive dimensions:"
249      write(*,*) "longitude: ",imold
250      write(*,*) "latitude: ",jmold
251      write(*,*) "altitude: ",lmold
252      write(*,*) "tracers: ",nqold
253      write(*,*) "subsurface_layers: ",nsoilold
254      if (depthinterpol) then
255      write(*,*) " => Warning, nsoilmx= ",nsoilmx
256      write(*,*) '    which implies that you want subterranean interpola
257     &tion.'
258      write(*,*) '  Otherwise, set nsoilmx -in dimphys.h- to: ',nsoilold
259      endif
260      write(*,*) "time lenght: ",timelen
261      write(*,*)
262
263!-----------------------------------------------------------------------
264! 2. Allocate arrays to store datasets
265!-----------------------------------------------------------------------
266
267      allocate(timelist(timelen))
268      allocate(rlonuold(imold+1), rlatvold(jmold))
269      allocate(rlonvold(imold+1), rlatuold(jmold+1))
270      allocate (apsold(lmold),bpsold(lmold))
271      allocate(uold(imold+1,jmold+1,lmold))
272      allocate(vold(imold+1,jmold+1,lmold))
273      allocate(told(imold+1,jmold+1,lmold))
274      allocate(psold(imold+1,jmold+1))
275      allocate(phisold(imold+1,jmold+1))
276      allocate(qold(imold+1,jmold+1,lmold,nqtot))
277      allocate(co2iceold(imold+1,jmold+1))
278      allocate(tsurfold(imold+1,jmold+1))
279      allocate(emisold(imold+1,jmold+1))
280      allocate(q2old(imold+1,jmold+1,lmold+1))
281!      allocate(tsoilold(imold+1,jmold+1,nsoilmx))
282      allocate(tsoilold(imold+1,jmold+1,nsoilold))
283      allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx))
284      allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia
285      allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx))
286      ! surface thermal inertia at old horizontal grid resolution
287      allocate(surfithold(imold+1,jmold+1))
288      allocate(mlayerold(nsoilold))
289      allocate(qsurfold(imold+1,jmold+1,nqtot))
290      allocate(tslabold(imold+1,jmold+1,nslay))
291      allocate(rnatold(imold+1,jmold+1))
292      allocate(pctsrf_sicold(imold+1,jmold+1))
293      allocate(tsea_iceold(imold+1,jmold+1))
294      allocate(sea_iceold(imold+1,jmold+1))
295
296      allocate(du_nonoro_gwdold(imold+1,jmold+1,lmold))
297      allocate(dv_nonoro_gwdold(imold+1,jmold+1,lmold))
298      allocate(east_gwstressold(imold+1,jmold+1,lmold))
299      allocate(west_gwstressold(imold+1,jmold+1,lmold))
300
301      allocate(var (imold+1,jmold+1,llm))
302      allocate(varp1 (imold+1,jmold+1,llm+1))
303
304      write(*,*) 'lect_start_archive: q2',ngrid,llm+1
305      write(*,*) 'lect_start_archive: q2S',iip1,jjp1,llm+1
306      write(*,*) 'lect_start_archive: q2old',imold+1,jmold+1,lmold+1
307
308!-----------------------------------------------------------------------
309! 3. Read time-independent data
310!-----------------------------------------------------------------------
311
312C-----------------------------------------------------------------------
313c 3.1. Lecture du tableau des parametres du run
314c     (pour  la lecture ulterieure de "ptotalold" et "co2icetotalold")
315c-----------------------------------------------------------------------
316c
317      ierr = NF_INQ_VARID (nid, "controle", nvarid)
318      IF (ierr .NE. NF_NOERR) THEN
319         PRINT*, "Lect_start_archive: champ <controle> not found"
320         CALL abort
321      ENDIF
322#ifdef NC_DOUBLE
323      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
324#else
325      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
326#endif
327      IF (ierr .NE. NF_NOERR) THEN
328         PRINT*, "lect_start_archive: Lecture echoue pour <controle>"
329         CALL abort
330      ENDIF
331c
332      tab0 = 50
333
334c-----------------------------------------------------------------------
335c 3.2 Lecture des longitudes et latitudes
336c-----------------------------------------------------------------------
337c
338      ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
339      IF (ierr .NE. NF_NOERR) THEN
340         PRINT*, "lect_start_archive: Field <rlonv> not found"
341         CALL abort
342      ENDIF
343#ifdef NC_DOUBLE
344      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
345#else
346      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
347#endif
348      IF (ierr .NE. NF_NOERR) THEN
349         PRINT*, "lect_start_archive: Failed loading <rlonv>"
350         CALL abort
351      ENDIF
352c
353      ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
354      IF (ierr .NE. NF_NOERR) THEN
355         PRINT*, "lect_start_archive: Field <rlatu> not found"
356         CALL abort
357      ENDIF
358#ifdef NC_DOUBLE
359      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
360#else
361      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
362#endif
363      IF (ierr .NE. NF_NOERR) THEN
364         PRINT*, "lect_start_archive: Failed loading <rlatu>"
365         CALL abort
366      ENDIF
367c
368      ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
369      IF (ierr .NE. NF_NOERR) THEN
370         PRINT*, "lect_start_archive: Field <rlonu> not found"
371         CALL abort
372      ENDIF
373#ifdef NC_DOUBLE
374      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
375#else
376      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
377#endif
378      IF (ierr .NE. NF_NOERR) THEN
379         PRINT*, "lect_start_archive: Failed loading <rlonu>"
380         CALL abort
381      ENDIF
382c
383      ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
384      IF (ierr .NE. NF_NOERR) THEN
385         PRINT*, "lect_start_archive: Field <rlatv> not found"
386         CALL abort
387      ENDIF
388#ifdef NC_DOUBLE
389      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
390#else
391      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
392#endif
393      IF (ierr .NE. NF_NOERR) THEN
394         PRINT*, "lect_start_archive: Failed loading <rlatv>"
395         CALL abort
396      ENDIF
397c
398
399c-----------------------------------------------------------------------
400c 3.3. Lecture des niveaux verticaux
401c-----------------------------------------------------------------------
402c
403      ierr = NF_INQ_VARID (nid, "aps", nvarid)
404      IF (ierr .NE. NF_NOERR) THEN
405         PRINT*, "lect_start_archive: Field <aps> not found"
406         apsold=0
407         PRINT*, "<aps> set to 0"
408      ELSE
409#ifdef NC_DOUBLE
410         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold)
411#else
412         ierr = NF_GET_VAR_REAL(nid, nvarid, apsold)
413#endif
414         IF (ierr .NE. NF_NOERR) THEN
415            PRINT*, "lect_start_archive: Failed loading <aps>"
416         ENDIF
417      ENDIF
418c
419      ierr = NF_INQ_VARID (nid, "bps", nvarid)
420      IF (ierr .NE. NF_NOERR) THEN
421         PRINT*, "lect_start_archive: Field <bps> not found"
422         PRINT*, "It must be an old start_archive, lets look for sig_s"
423         ierr = NF_INQ_VARID (nid, "sig_s", nvarid)
424         IF (ierr .NE. NF_NOERR) THEN
425            PRINT*, "Nothing to do..."
426            CALL abort
427         ENDIF
428      ENDIF
429#ifdef NC_DOUBLE
430      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold)
431#else
432      ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold)
433#endif
434      IF (ierr .NE. NF_NOERR) THEN
435         PRINT*, "lect_start_archive: Failed loading <bps>"
436         CALL abort
437      END IF
438
439c-----------------------------------------------------------------------
440c 3.4 Read Soil layers depths
441c-----------------------------------------------------------------------
442     
443      ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
444      if (ierr.ne.NF_NOERR) then
445       write(*,*)'lect_start_archive: Could not find <soildepth>'
446       write(*,*)' => Assuming this is an archive in old format'
447       olddepthdef=.true.
448       depthinterpol=.true.
449       ! this is how soil depth was defined in ye old days
450        do isoil=1,nsoilold
451          mlayerold(isoil)=sqrt(887.75/3.14)*((2.**(isoil-0.5))-1.)
452        enddo
453      else
454#ifdef NC_DOUBLE
455        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold)
456#else
457        ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold)
458#endif
459       if (ierr .NE. NF_NOERR) then
460         PRINT*, "lect_start_archive: Failed reading <soildepth>"
461         CALL abort
462       endif
463
464      endif !of if(ierr.ne.NF_NOERR)
465
466      ! Read (or build) mlayer()
467      if (depthinterpol) then
468       ! Build (default) new soil depths (mlayer(:) is in comsoil.h),
469       ! as in soil_settings.F
470       write(*,*)' => Building default soil depths'
471       do isoil=0,nsoilmx-1
472         mlayer(isoil)=2.e-4*(2.**(isoil-0.5))
473       enddo
474       write(*,*)' => mlayer: ',mlayer
475       ! Also build (default) new soil interlayer depth layer(:)
476       do isoil=1,nsoilmx
477         layer(isoil)=sqrt(mlayer(0)*mlayer(1))*
478     &                      ((mlayer(1)/mlayer(0))**(isoil-1))
479       enddo
480       write(*,*)' =>  layer: ',layer
481      else ! read mlayer() from file
482#ifdef NC_DOUBLE
483        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
484#else
485        ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
486#endif
487       if (ierr .NE. NF_NOERR) then
488         PRINT*, "lect_start_archive: Failed reading <soildepth>"
489         CALL abort
490       endif
491      endif ! of if (depthinterpol)
492
493c-----------------------------------------------------------------------
494c 3.5 Read Soil thermal inertia
495c-----------------------------------------------------------------------
496
497      ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
498      if (ierr.ne.NF_NOERR) then
499       write(*,*)'lect_start_archive: Could not find <inertiedat>'
500       write(*,*)' => Assuming this is an archive in old format'
501       therminertia_3D=.false.
502       write(*,*)' => Thermal inertia will be read from reference file'
503       volcapa=1.e6
504       write(*,*)'    and soil volumetric heat capacity is set to ',
505     &           volcapa
506      else
507#ifdef NC_DOUBLE
508        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold)
509#else
510        ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold)
511#endif
512       if (ierr .NE. NF_NOERR) then
513         PRINT*, "lect_start_archive: Failed reading <inertiedat>"
514         CALL abort
515       endif
516      endif
517
518c-----------------------------------------------------------------------
519c 3.6 Lecture geopotentiel au sol
520c-----------------------------------------------------------------------
521c
522      ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
523      IF (ierr .NE. NF_NOERR) THEN
524         PRINT*, "lect_start_archive: Field <phisinit> not found"
525         CALL abort
526      ENDIF
527#ifdef NC_DOUBLE
528      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
529#else
530      ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
531#endif
532      IF (ierr .NE. NF_NOERR) THEN
533         PRINT*, "lect_start_archive: Failed loading <phisinit>"
534         CALL abort
535      ENDIF
536
537C-----------------------------------------------------------------------
538c   lecture de "ptotalold" et "co2icetotalold"
539c-----------------------------------------------------------------------
540      ptotalold = tab_cntrl(tab0+49)
541      co2icetotalold = tab_cntrl(tab0+50)
542 
543c-----------------------------------------------------------------------
544c 4. Lecture du temps et choix
545c-----------------------------------------------------------------------
546 
547c  lecture du temps
548c
549      ierr = NF_INQ_DIMID (nid, "Time", nvarid)
550      IF (ierr .NE. NF_NOERR) THEN
551         ierr = NF_INQ_DIMID (nid, "temps", nvarid)
552         IF (ierr .NE. NF_NOERR) THEN
553            PRINT*, "lect_start_archive: Field <Time> not found"
554            CALL abort
555         endif
556      ENDIF
557
558      ierr = NF_INQ_VARID (nid, "Time", nvarid)
559      IF (ierr .NE. NF_NOERR) THEN
560         ierr = NF_INQ_VARID (nid, "temps", nvarid)
561      endif
562#ifdef NC_DOUBLE
563      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
564#else
565      ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
566#endif
567      IF (ierr .NE. NF_NOERR) THEN
568         PRINT*, "lect_start_archive: Failed loading <Time>"
569         CALL abort
570      ENDIF
571c
572      write(*,*)
573      write(*,*)
574      write(*,*) 'Available dates for the stored initial conditions:'
575      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
576      pi=2.*ASIN(1.)
577      do i=1,timelen
578c       call solarlong(timelist(i),sollong(i))
579c       sollong(i) = sollong(i)*180./pi
580        write(*,*) 'initial state for day ' ,int(timelist(i))
581c       write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
582c    .    sollong(i)
583      end do
584
585   6  FORMAT(i7,i7,f9.3)
586 
587      write(*,*)
588      write(*,*) 'Choice for the date'
589 123  read(*,*,iostat=ierr) date
590      if(ierr.ne.0) goto 123
591      memo = 0
592      do i=1,timelen
593        if (date.eq.int(timelist(i))) then
594            memo = i
595        endif
596      end do
597 
598      if (memo.eq.0) then
599        write(*,*)
600        write(*,*)
601        write(*,*) "Wrong value... can't you read !?!"
602        write(*,*)
603        write(*,*) 'Available dates for the stored initial conditions:'
604        write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
605        do i=1,timelen
606          write(*,*) 'initial state for day ' ,nint(timelist(i))
607c         write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
608        end do
609        goto 123
610      endif
611
612!-----------------------------------------------------------------------
613! 5. Read (time-dependent) data from datafile
614!-----------------------------------------------------------------------
615
616
617c-----------------------------------------------------------------------
618c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], qsurf)
619c-----------------------------------------------------------------------
620 
621      start=(/1,1,memo,0/)
622      count=(/imold+1,jmold+1,1,0/)
623       
624      ierr = NF_INQ_VARID (nid, "emis", nvarid)
625      IF (ierr .NE. NF_NOERR) THEN
626         PRINT*, "lect_start_archive: Field <emis> not found"
627         CALL abort
628      ENDIF
629#ifdef NC_DOUBLE
630      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
631#else
632      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
633#endif
634      IF (ierr .NE. NF_NOERR) THEN
635         PRINT*, "lect_start_archive: Failed loading <emis>"
636         CALL abort
637      ENDIF
638c
639      ierr = NF_INQ_VARID (nid, "ps", nvarid)
640      IF (ierr .NE. NF_NOERR) THEN
641         PRINT*, "lect_start_archive: Field <ps> not found"
642         CALL abort
643      ENDIF
644#ifdef NC_DOUBLE
645      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
646#else
647      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
648#endif
649      IF (ierr .NE. NF_NOERR) THEN
650         PRINT*, "lect_start_archive: Failed loading <ps>"
651         CALL abort
652      ENDIF
653c
654      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
655      IF (ierr .NE. NF_NOERR) THEN
656         PRINT*, "lect_start_archive: Field <tsurf> not found"
657         CALL abort
658      ENDIF
659#ifdef NC_DOUBLE
660      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
661#else
662      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
663#endif
664      IF (ierr .NE. NF_NOERR) THEN
665         PRINT*, "lect_start_archive: Failed loading <tsurf>"
666         CALL abort
667      ENDIF
668c
669      ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
670      IF (ierr .NE. NF_NOERR) THEN
671         PRINT*, "lect_start_archive: Field <q2surf> not found"
672         CALL abort
673      ENDIF
674#ifdef NC_DOUBLE
675      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
676#else
677      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
678#endif
679      IF (ierr .NE. NF_NOERR) THEN
680         PRINT*, "lect_start_archive: Failed loading <q2surf>"
681         CALL abort
682      ENDIF
683c
684cc Slab ocean
685      if(ok_slab_ocean) then
686      start=(/1,1,1,memo/)
687      count=(/imold+1,jmold+1,nslay,1/)
688
689       ierr=NF_INQ_VARID(nid,"tslab",nvarid)
690       IF (ierr.ne.NF_NOERR) then
691          PRINT*,"lect_start_archive: Cannot find <tslab>"
692       ENDIF
693#ifdef NC_DOUBLE
694      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tslabold)
695#else
696      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tslabold)
697#endif
698      IF (ierr .NE. NF_NOERR) THEN
699         PRINT*, "lect_start_archive: Failed loading <tslab>"
700      ENDIF
701
702
703c
704      start=(/1,1,memo,0/)
705      count=(/imold+1,jmold+1,1,0/)
706
707      ierr = NF_INQ_VARID (nid, "rnat", nvarid)
708      IF (ierr .NE. NF_NOERR) THEN
709         PRINT*, "lect_start_archive: Field <rnat> not found"
710      ENDIF
711#ifdef NC_DOUBLE
712      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,rnatold)
713#else
714      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,rnatold)
715#endif
716      IF (ierr .NE. NF_NOERR) THEN
717         PRINT*, "lect_start_archive: Failed loading <rnat>"
718      ENDIF
719c
720      ierr = NF_INQ_VARID (nid, "pctsrf_sic", nvarid)
721      IF (ierr .NE. NF_NOERR) THEN
722         PRINT*, "lect_start_archive: Field <pctsrf_sic> not found"
723      ENDIF
724#ifdef NC_DOUBLE
725      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,pctsrf_sicold)
726#else
727      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,pctsrf_sicold)
728#endif
729      IF (ierr .NE. NF_NOERR) THEN
730         PRINT*, "lect_start_archive: Failed loading <pctsrf_sic>"
731      ENDIF
732c
733      ierr = NF_INQ_VARID (nid, "tsea_ice", nvarid)
734      IF (ierr .NE. NF_NOERR) THEN
735         PRINT*, "lect_start_archive: Field <tsea_ice> not found"
736      ENDIF
737#ifdef NC_DOUBLE
738      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsea_iceold)
739#else
740      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsea_iceold)
741#endif
742      IF (ierr .NE. NF_NOERR) THEN
743         PRINT*, "lect_start_archive: Failed loading <tsea_ice>"
744      ENDIF
745c
746      ierr = NF_INQ_VARID (nid, "sea_ice", nvarid)
747      IF (ierr .NE. NF_NOERR) THEN
748         PRINT*, "lect_start_archive: Field <sea_ice> not found"
749      ENDIF
750#ifdef NC_DOUBLE
751      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,sea_iceold)
752#else
753      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,sea_iceold)
754#endif
755      IF (ierr .NE. NF_NOERR) THEN
756         PRINT*, "lect_start_archive: Failed loading <sea_ice>"
757      ENDIF
758 
759      ENDIF! ok_slab_ocean
760c
761      write(*,*)"lect_start_archive: rlonuold:"
762     &           ,rlonuold," rlatvold:",rlatvold
763      write(*,*)
764
765! Surface tracers:     
766      ! initialize all surface tracers to zero
767      qsurfold(1:imold+1,1:jmold+1,1:nqtot)=0
768
769      DO iq=1,nqtot
770          txt=trim(tname(iq))//"_surf"
771          if (txt.eq."h2o_vap") then
772            ! There is no surface tracer for h2o_vap;
773            ! "h2o_ice" should be loaded instead
774            txt="h2o_ice_surf"
775            write(*,*) 'lect_start_archive: loading surface tracer',
776     &                     ' h2o_ice instead of h2o_vap'
777          endif
778
779       
780        write(*,*) "lect_start_archive: loading tracer ",trim(txt)
781        ierr = NF_INQ_VARID (nid,txt,nvarid)
782        IF (ierr .NE. NF_NOERR) THEN
783          PRINT*, "lect_start_archive: ",
784     &              " Tracer <",trim(txt),"> not found"
785
786!          print*,'RDW has added hack to let me continue...'
787!          CALL abort
788        ENDIF
789#ifdef NC_DOUBLE
790        ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
791     &          qsurfold(1,1,iq))
792#else
793        ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
794     &          qsurfold(1,1,iq))
795#endif
796        IF (ierr .NE. NF_NOERR) THEN
797          PRINT*, "lect_start_archive: ",
798     &             " Failed loading <",trim(txt),">"
799          write (*,*) trim(txt),'    is set to 0'
800        ENDIF
801
802      ENDDO ! of DO iq=1,nqtot
803
804
805!-----------------------------------------------------------------------
806! 5.2 Read 3D subterranean fields
807!-----------------------------------------------------------------------
808
809      start=(/1,1,1,memo/)
810      count=(/imold+1,jmold+1,nsoilold,1/)
811!
812! Read soil temperatures
813!
814      if (olddepthdef) then ! tsoil stored using the 'old format'
815         start=(/1,1,memo,0/)
816         count=(/imold+1,jmold+1,1,0/) ! because the "Tg" are 2D datasets
817       do isoil=1,nsoilold
818         write(str2,'(i2.2)') isoil
819c
820         ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
821         IF (ierr .NE. NF_NOERR) THEN
822            PRINT*, "lect_start_archive: ",
823     &              "Field <","Tg"//str2,"> not found"
824            CALL abort
825         ENDIF
826#ifdef NC_DOUBLE
827         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
828     &          tsoilold(1,1,isoil))
829#else
830         ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
831     &          tsoilold(1,1,isoil))
832#endif
833         IF (ierr .NE. NF_NOERR) THEN
834            PRINT*, "lect_start_archive: ",
835     &            "Failed reading <","Tg"//str2,">"
836            CALL abort
837         ENDIF
838c
839       enddo ! of do isoil=1,nsoilold
840     
841      ! reset 'start' and 'count' to "3D" behaviour
842      start=(/1,1,1,memo/)
843      count=(/imold+1,jmold+1,nsoilold,1/)
844     
845      else
846       write(*,*) "lect_start_archive: loading tsoil "
847       ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
848       if (ierr.ne.NF_NOERR) then
849        write(*,*)"lect_start_archive: Cannot find <tsoil>"
850        call abort
851       else
852#ifdef NC_DOUBLE
853      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
854#else
855      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
856#endif
857       endif ! of if (ierr.ne.NF_NOERR)
858       
859      endif ! of if (olddepthdef)
860
861c-----------------------------------------------------------------------
862c 5.3   Lecture des champs 3D (t,u,v, q2atm,q)
863c-----------------------------------------------------------------------
864
865      start=(/1,1,1,memo/)
866      count=(/imold+1,jmold+1,lmold,1/)
867
868c
869      ierr = NF_INQ_VARID (nid,"temp", nvarid)
870      IF (ierr .NE. NF_NOERR) THEN
871         PRINT*, "lect_start_archive: Field <temp> not found"
872         CALL abort
873      ENDIF
874#ifdef NC_DOUBLE
875      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
876#else
877      ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
878#endif
879      IF (ierr .NE. NF_NOERR) THEN
880         PRINT*, "lect_start_archive: Failed loading <temp>"
881         CALL abort
882      ENDIF
883c
884      ierr = NF_INQ_VARID (nid,"u", nvarid)
885      IF (ierr .NE. NF_NOERR) THEN
886         PRINT*, "lect_start_archive: Field <u> not found"
887         CALL abort
888      ENDIF
889#ifdef NC_DOUBLE
890      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
891#else
892      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
893#endif
894      IF (ierr .NE. NF_NOERR) THEN
895         PRINT*, "lect_start_archive: Failed loading <u>"
896         CALL abort
897      ENDIF
898c
899      ierr = NF_INQ_VARID (nid,"v", nvarid)
900      IF (ierr .NE. NF_NOERR) THEN
901         PRINT*, "lect_start_archive: Field <v> not found"
902         CALL abort
903      ENDIF
904#ifdef NC_DOUBLE
905      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
906#else
907      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
908#endif
909      IF (ierr .NE. NF_NOERR) THEN
910         PRINT*, "lect_start_archive: Failed loading <v>"
911         CALL abort
912      ENDIF
913c
914      ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
915      IF (ierr .NE. NF_NOERR) THEN
916         PRINT*, "lect_start_archive: Field <q2atm> not found"
917         CALL abort
918      ENDIF
919#ifdef NC_DOUBLE
920      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
921#else
922      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
923#endif
924      IF (ierr .NE. NF_NOERR) THEN
925         PRINT*, "lect_start_archive: Failed loading <q2atm>"
926         CALL abort
927      ENDIF
928c
929
930! Tracers:     
931      qold(1:imold+1,1:jmold+1,1:lmold,1:nqtot)=0.
932
933      DO iq=1,nqtot
934        txt=tname(iq)
935        write(*,*)"lect_start_archive: loading tracer ",trim(txt)
936        ierr = NF_INQ_VARID (nid,txt,nvarid)
937        IF (ierr .NE. NF_NOERR) THEN
938            PRINT*, "lect_start_archive: ",
939     &              " Tracer <",trim(txt),"> not found"
940!            CALL abort
941        ENDIF
942#ifdef NC_DOUBLE
943        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
944#else
945        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
946#endif
947        IF (ierr .NE. NF_NOERR) THEN
948          PRINT*, "lect_start_archive: ",
949     &             "  Failed loading <",trim(txt),">"
950          write (*,*) trim(txt),'      set to 1.E-30'
951          do l=1,lmold
952            do j=1,jmold+1
953              do i=1,imold+1
954                 qold(i,j,l,iq)=1.e-30
955              end do
956            end do
957          end do
958        ENDIF
959
960      ENDDO ! of DO iq=1,nqtot
961
962! Non-orographic GWs:
963      write(*,*)"lect_start_archive: loading du_nonoro_gwd"
964      ierr = NF_INQ_VARID (nid,"du_nonoro_gwd", nvarid)
965      IF (ierr .NE. NF_NOERR) THEN
966         PRINT*, "lect_start_archive: Field <du_nonoro_gwd> not found"
967         PRINT*, "Setting it to zero"
968         du_nonoro_gwdold(:,:,:)=0
969      ELSE
970#ifdef NC_DOUBLE
971        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,du_nonoro_gwdold)
972#else
973        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,du_nonoro_gwdold)
974#endif
975        IF (ierr .NE. NF_NOERR) THEN
976          PRINT*, "lect_start_archive: Failed loading <du_nonoro_gwd>"
977          CALL abort
978        ENDIF
979      ENDIF
980
981      write(*,*)"lect_start_archive: loading dv_nonoro_gwd"
982      ierr = NF_INQ_VARID (nid,"dv_nonoro_gwd", nvarid)
983      IF (ierr .NE. NF_NOERR) THEN
984         PRINT*, "lect_start_archive: Field <dv_nonoro_gwd> not found"
985         PRINT*, "Setting it to zero"
986         dv_nonoro_gwdold(:,:,:)=0
987      ELSE
988#ifdef NC_DOUBLE
989        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,dv_nonoro_gwdold)
990#else
991        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,dv_nonoro_gwdold)
992#endif
993        IF (ierr .NE. NF_NOERR) THEN
994          PRINT*, "lect_start_archive: Failed loading <dv_nonoro_gwd>"
995          CALL abort
996        ENDIF
997      ENDIF
998
999      write(*,*)"lect_start_archive: loading east_gwstress"
1000      ierr = NF_INQ_VARID (nid,"east_gwstress", nvarid)
1001      IF (ierr .NE. NF_NOERR) THEN
1002         PRINT*, "lect_start_archive: Field <east_gwstress> not found"
1003         PRINT*, "Setting it to zero"
1004         east_gwstressold(:,:,:)=0
1005      ELSE
1006#ifdef NC_DOUBLE
1007        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,east_gwstressold)
1008#else
1009        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,east_gwstressold)
1010#endif
1011        IF (ierr .NE. NF_NOERR) THEN
1012          PRINT*, "lect_start_archive: Failed loading <east_gwstress>"
1013          CALL abort
1014        ENDIF
1015      ENDIF
1016
1017      write(*,*)"lect_start_archive: loading west_gwstress"
1018      ierr = NF_INQ_VARID (nid,"west_gwstress", nvarid)
1019      IF (ierr .NE. NF_NOERR) THEN
1020         PRINT*, "lect_start_archive: Field <west_gwstress> not found"
1021         PRINT*, "Setting it to zero"
1022         west_gwstressold(:,:,:)=0
1023      ELSE
1024#ifdef NC_DOUBLE
1025        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,west_gwstressold)
1026#else
1027        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,west_gwstressold)
1028#endif
1029        IF (ierr .NE. NF_NOERR) THEN
1030          PRINT*, "lect_start_archive: Failed loading <west_gwstress>"
1031          CALL abort
1032        ENDIF
1033      ENDIF
1034
1035!=======================================================================
1036! 6. Interpolation from old grid to new grid
1037!=======================================================================
1038
1039c=======================================================================
1040c   INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
1041c=======================================================================
1042c  Interpolation horizontale puis passage dans la grille physique pour
1043c  les variables physique
1044c  Interpolation verticale puis horizontale pour chaque variable 3D
1045c=======================================================================
1046
1047c-----------------------------------------------------------------------
1048c 6.1   Variable 2d :
1049c-----------------------------------------------------------------------
1050c Relief
1051      call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
1052     &                   rlonuold,rlatvold,rlonu,rlatv)
1053
1054! CO2 ice is now in qsurf(igcm_co2_ice)
1055!      call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
1056!     &                   rlonuold,rlatvold,rlonu,rlatv)
1057
1058c Temperature de surface
1059      call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
1060     &                   rlonuold,rlatvold,rlonu,rlatv)
1061      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs,tsurf)
1062c     write(44,*) 'tsurf', tsurf
1063
1064c Temperature du sous-sol
1065!      call interp_horiz(tsoilold,tsoils,
1066!     &                  imold,jmold,iim,jjm,nsoilmx,
1067!     &                   rlonuold,rlatvold,rlonu,rlatv)
1068!      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoils,tsoil)
1069c     write(45,*) 'tsoil',tsoil
1070
1071c Emissivite de la surface
1072      call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
1073     &                   rlonuold,rlatvold,rlonu,rlatv)
1074      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss,emis)
1075c     write(46,*) 'emis',emis
1076
1077
1078
1079c-----------------------------------------------------------------------
1080c 6.1.2 Traitement special de la pression au sol :
1081c-----------------------------------------------------------------------
1082
1083c  Extrapolation la pression dans la nouvelle grille
1084      call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
1085     &                   rlonuold,rlatvold,rlonu,rlatv)
1086
1087c-----------------------------------------------------------------------
1088c       On assure la conservation de la masse de l'atmosphere + calottes
1089c-----------------------------------------------------------------------
1090
1091      ptotal =  0.
1092      DO j=1,jjp1
1093         DO i=1,iim
1094            ptotal=ptotal+ps(i,j)*aire(i,j)/g
1095         END DO
1096      END DO
1097      co2icetotal = 0.
1098      if (igcm_co2_ice.ne.0) then
1099        ! recast surface CO2 ice on new grid
1100        call interp_horiz(qsurfold(1,1,igcm_co2_ice),
1101     &                  qsurfs(1,1,igcm_co2_ice),
1102     &                  imold,jmold,iim,jjm,1,
1103     &                  rlonuold,rlatvold,rlonu,rlatv)
1104       DO j=1,jjp1
1105         DO i=1,iim
1106           !co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
1107           co2icetotal=co2icetotal+qsurfS(i,j,igcm_co2_ice)*aire(i,j)
1108         END DO
1109       END DO
1110      else
1111        write(*,*) "Warning: No co2_ice surface tracer"
1112      endif
1113
1114      write(*,*)
1115      write(*,*)'Old grid: atmospheric mass :',ptotalold
1116      write(*,*)'New grid: atmospheric mass :',ptotal
1117      write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold
1118      write(*,*)
1119      write(*,*)'Old grid: mass of CO2 ice:',co2icetotalold
1120      write(*,*)'New grid: mass of CO2 ice:',co2icetotal
1121      if (co2icetotalold.ne.0.) then
1122      write(*,*)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
1123      endif
1124      write(*,*)
1125
1126
1127      DO j=1,jjp1
1128         DO i=1,iip1
1129            ps(i,j)=ps(i,j) * ptotalold/ptotal
1130         END DO
1131      END DO
1132
1133      if ( co2icetotalold.gt.0.) then
1134!         DO j=1,jjp1
1135!            DO i=1,iip1
1136!               co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
1137!            END DO
1138!         END DO
1139        write(*,*) "Not enforcing conservation of surface CO2 ice"
1140        write(*,*) " (should be OK when CO2 ice is a tracer)"
1141      end if
1142
1143c-----------------------------------------------------------------------
1144c 6.2 Subterranean 3d variables:
1145c-----------------------------------------------------------------------
1146
1147c-----------------------------------------------------------------------
1148c 6.2.1 Thermal Inertia
1149c       Note: recall that inertiedat is a common in "comsoil.h"
1150c-----------------------------------------------------------------------
1151
1152      ! depth-wise interpolation, if required
1153      if (depthinterpol.and.(.not.olddepthdef)) then
1154        allocate(oldval(nsoilold))
1155        allocate(newval(nsoilmx))
1156        write(*,*)'lect_start_archive: WARNING: vertical interpolation o
1157     &f soil thermal inertia; might be wiser to reset it.'
1158        write(*,*)
1159       
1160        do i=1,imold+1
1161         do j=1,jmold+1
1162           !copy old values
1163           oldval(1:nsoilold)=inertiedatold(i,j,1:nsoilold)
1164           !interpolate
1165           call interp_line(mlayerold,oldval,nsoilold,
1166     &                     mlayer,newval,nsoilmx)
1167           !copy interpolated values
1168           inertiedatoldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1169         enddo
1170        enddo
1171        ! cleanup
1172        deallocate(oldval)
1173        deallocate(newval)
1174      endif !of if (depthinterpol)
1175
1176      if (therminertia_3D) then
1177        ! We have inertiedatold
1178       if((imold.ne.iim).or.(jmold.ne.jjm)) then
1179       write(*,*)'lect_start_archive: WARNING: horizontal interpolation
1180     &of thermal inertia; might be better to reset it.'
1181       write(*,*)
1182       endif
1183       
1184        ! Do horizontal interpolation
1185        if (depthinterpol) then
1186          call interp_horiz(inertiedatoldnew,inertiedatS,
1187     &                  imold,jmold,iim,jjm,nsoilmx,
1188     &                   rlonuold,rlatvold,rlonu,rlatv)
1189        else
1190          call interp_horiz(inertiedatold,inertiedatS,
1191     &                  imold,jmold,iim,jjm,nsoilold,
1192     &                   rlonuold,rlatvold,rlonu,rlatv)
1193        endif ! of if (depthinterpol)
1194
1195      else ! no 3D thermal inertia data
1196       write(*,*)'lect_start_archive: using reference surface inertia'
1197        ! Use surface inertia (and extend it to all depths)
1198        do i=1,nsoilmx
1199         inertiedatS(1:iip1,1:jjp1,i)=surfith(1:iip1,1:jjp1)
1200        enddo
1201        ! Build an old resolution surface thermal inertia
1202        ! (will be needed for tsoil interpolation)
1203        call interp_horiz(surfith,surfithold,
1204     &                    iim,jjm,imold,jmold,1,
1205     &                    rlonu,rlatv,rlonuold,rlatvold)
1206      endif
1207
1208
1209      ! Reshape inertiedatS to scalar grid as inertiedat
1210      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,
1211     &                  inertiedatS,inertiedat)
1212     
1213c-----------------------------------------------------------------------
1214c 6.2.2 Soil temperature
1215c-----------------------------------------------------------------------
1216!      write(*,*) 'Soil'
1217
1218      !print*,'Problem in lect_start_archive interpolating'
1219      !print*,'to new resolution!!'
1220
1221      ! Recast temperatures along soil depth, if necessary
1222      if (olddepthdef) then
1223        allocate(oldgrid(nsoilold+1))
1224        allocate(oldval(nsoilold+1))
1225        allocate(newval(nsoilmx))
1226        do i=1,imold+1
1227         do j=1,jmold+1
1228
1229            !if(i.gt.iip1 .or. j.gt.jjp1)then
1230               !print*,'Problem in lect_start_archive interpolating'
1231               !print*,'to new resolution!!'
1232               !call abort
1233            !endif
1234
1235           ! copy values
1236           oldval(1)=tsurfold(i,j)
1237!          oldval(1)=tsurfS(i,j)
1238           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1239           ! build vertical coordinate
1240           oldgrid(1)=0. ! ground
1241           oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)*
1242     &                (surfithold(i,j)/1.e6)
1243          ! Note; at this stage, we impose volcapa=1.e6 above
1244          ! since volcapa isn't set in old soil definitions
1245
1246          ! interpolate
1247          call interp_line(oldgrid,oldval,nsoilold+1,
1248     &                     mlayer,newval,nsoilmx)
1249         ! copy result in tsoilold
1250         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1251         enddo
1252        enddo
1253        ! cleanup
1254        deallocate(oldgrid)
1255        deallocate(oldval)
1256        deallocate(newval)
1257
1258      else
1259       if (depthinterpol) then ! if vertical interpolation is required
1260        allocate(oldgrid(nsoilold+1))
1261        allocate(oldval(nsoilold+1))
1262        allocate(newval(nsoilmx))
1263        ! build vertical coordinate
1264        oldgrid(1)=0. ! ground
1265        oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
1266        do i=1,imold+1
1267         do j=1,jmold+1
1268           ! copy values
1269           oldval(1)=tsurfold(i,j)
1270!          oldval(1)=tsurfS(i,j)
1271           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1272          ! interpolate
1273          call interp_line(oldgrid,oldval,nsoilold+1,
1274     &                     mlayer,newval,nsoilmx)
1275         ! copy result in tsoilold
1276         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1277         enddo
1278        enddo
1279!       write(*,*)'tsoiloldnew(1,1,1):',tsoiloldnew(1,1,1)
1280        ! cleanup
1281        deallocate(oldgrid)
1282        deallocate(oldval)
1283        deallocate(newval)
1284       
1285       else
1286        tsoiloldnew(:,:,:)=tsoilold(:,:,:)
1287       endif ! of if (depthinterpol)
1288      endif ! of if (olddepthdef)
1289
1290      ! Do the horizontal interpolation
1291       call interp_horiz(tsoiloldnew,tsoilS,
1292     &                  imold,jmold,iim,jjm,nsoilmx,
1293     &                   rlonuold,rlatvold,rlonu,rlatv)
1294
1295      ! Reshape tsoilS to scalar grid as tsoil
1296       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS,tsoil)
1297
1298c-----------------------------------------------------------------------
1299c 6.3   Slab Ocean :
1300c-----------------------------------------------------------------------
1301      call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,nslay,
1302     &                   rlonuold,rlatvold,rlonu,rlatv)
1303      call gr_dyn_fi (nslay,iim+1,jjm+1,ngrid,tslabs,tslab)
1304
1305      call interp_horiz (rnatold,rnats,imold,jmold,iim,jjm,1,
1306     &                   rlonuold,rlatvold,rlonu,rlatv)
1307      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,rnats,rnat)
1308
1309      call interp_horiz (pctsrf_sicold,pctsrf_sics,imold,jmold,iim,
1310     &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
1311      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,pctsrf_sics,pctsrf_sic)
1312
1313      call interp_horiz (tsea_iceold,tsea_ices,imold,jmold,iim,jjm,1,
1314     &                   rlonuold,rlatvold,rlonu,rlatv)
1315      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsea_ices,tsea_ice)
1316
1317      call interp_horiz (sea_iceold,sea_ices,imold,jmold,iim,jjm,1,
1318     &                   rlonuold,rlatvold,rlonu,rlatv)
1319      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,sea_ices,sea_ice)
1320
1321c-----------------------------------------------------------------------
1322c 6.4 Variable 3d :
1323c-----------------------------------------------------------------------
1324     
1325c temperatures atmospheriques
1326      write (*,*) 'lect_start_archive: told ', told (1,jmold+1,1)  ! INFO
1327      call interp_vert
1328     &    (told,var,lmold,llm,apsold,bpsold,aps,bps,
1329     &     psold,(imold+1)*(jmold+1))
1330      write (*,*) 'lect_start_archive: var ', var (1,jmold+1,1)  ! INFO
1331      call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1332     &                   rlonuold,rlatvold,rlonu,rlatv)
1333      write (*,*) 'lect_start_archive: t ', t(1,jjp1,1)  ! INFO
1334
1335! Non-orographic GW
1336      call interp_vert
1337     &    (du_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
1338     &     psold,(imold+1)*(jmold+1))
1339      call interp_horiz(var,du_nonoro_gwdS,imold,jmold,iim,jjm,llm,
1340     &                   rlonuold,rlatvold,rlonu,rlatv)
1341      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,du_nonoro_gwdS,du_nonoro_gwd)
1342     
1343      call interp_vert
1344     &    (dv_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
1345     &     psold,(imold+1)*(jmold+1))
1346      call interp_horiz(var,dv_nonoro_gwdS,imold,jmold,iim,jjm,llm,
1347     &                   rlonuold,rlatvold,rlonu,rlatv)
1348      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,dv_nonoro_gwdS,dv_nonoro_gwd)
1349     
1350      call interp_vert
1351     &    (east_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
1352     &     psold,(imold+1)*(jmold+1))
1353      call interp_horiz(var,east_gwstressS,imold,jmold,iim,jjm,llm,
1354     &                   rlonuold,rlatvold,rlonu,rlatv)
1355      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,east_gwstressS,east_gwstress)
1356     
1357      call interp_vert
1358     &    (west_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
1359     &     psold,(imold+1)*(jmold+1))
1360      call interp_horiz(var,west_gwstressS,imold,jmold,iim,jjm,llm,
1361     &                   rlonuold,rlatvold,rlonu,rlatv)
1362      call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,west_gwstressS,west_gwstress)
1363
1364c q2 : pbl wind variance
1365      write (*,*) 'lect_start_archive: q2old ', q2old (1,2,1)  ! INFO
1366      call interp_vert (q2old,varp1,lmold+1,llm+1,
1367     &     apsold,bpsold,ap,bp,psold,(imold+1)*(jmold+1))
1368      write (*,*) 'lect_start_archive: varp1 ', varp1 (1,2,1)  ! INFO
1369      call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
1370     &                   rlonuold,rlatvold,rlonu,rlatv)
1371      write (*,*) 'lect_start_archive: q2s ', q2s (1,2,1)  ! INFO
1372      call gr_dyn_fi (llm+1,iim+1,jjm+1,ngrid,q2s,q2)
1373      write (*,*) 'lect_start_archive: q2 ', q2 (1,2)  ! INFO
1374c     write(47,*) 'q2',q2
1375
1376c calcul des champ de vent; passage en vent covariant
1377      write (*,*) 'lect_start_archive: uold ', uold (1,2,1)  ! INFO
1378      call interp_vert
1379     & (uold,var,lmold,llm,apsold,bpsold,aps,bps,
1380     &  psold,(imold+1)*(jmold+1))
1381      write (*,*) 'lect_start_archive: var ', var (1,2,1)  ! INFO
1382      call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1383     &                   rlonuold,rlatvold,rlonu,rlatv)
1384      write (*,*) 'lect_start_archive: us ', us (1,2,1)   ! INFO
1385
1386      call interp_vert
1387     & (vold,var,lmold,llm,
1388     &  apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1389      call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1390     &                   rlonuold,rlatvold,rlonu,rlatv)
1391      call scal_wind(us,vs,unat,vnat)
1392      write (*,*) 'lect_start_archive: unat ', unat (1,2,1)    ! INFO
1393      do l=1,llm
1394        do j = 1, jjp1
1395          do i=1,iip1
1396            ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1397c           ucov( i,j,l ) = 0
1398          end do
1399        end do
1400      end do
1401      write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1)  ! INFO
1402c     write(48,*) 'ucov',ucov
1403      do l=1,llm
1404        do j = 1, jjm
1405          do i=1,iim
1406            vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1407c           vcov( i,j,l ) = 0
1408          end do
1409          vcov( iip1,j,l ) = vcov( 1,j,l )
1410        end do
1411      end do
1412c     write(49,*) 'ucov',vcov
1413
1414      if (nqtot .gt. 0) then
1415c traceurs surface
1416      do iq = 1, nqtot
1417            call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
1418     &                  imold,jmold,iim,jjm,1,
1419     &                  rlonuold,rlatvold,rlonu,rlatv)
1420      enddo
1421
1422      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs,qsurf)
1423
1424c traceurs 3D
1425      do  iq = 1, nqtot
1426            call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1427     &        apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1428            call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1429     &                  rlonuold,rlatvold,rlonu,rlatv)
1430      enddo
1431cccccccccccccccccccccccccccccc     
1432c  make sure that sum of q = 1     
1433c dominent species is = 1 - sum(all other species)     
1434cccccccccccccccccccccccccccccc     
1435c      iqmax=1
1436c     
1437c      if (nqold.gt.10) then
1438c       do l=1,llm
1439c        do j=1,jjp1
1440c          do i=1,iip1
1441c           do iq=1,nqold
1442c            if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
1443c              iqmax=iq
1444c            endif
1445c           enddo
1446c           q(i,j,l,iqmax)=1.
1447c           qtot(i,j,l)=0
1448c           do iq=1,nqold
1449c            if (iq.ne.iqmax) then       
1450c              q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)       
1451c            endif
1452c           enddo !iq
1453c           do iq=1,nqold
1454c            qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
1455c            if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(*,*)' qtot(i,j,l)',
1456c     $    qtot(i,j,l)
1457c           enddo !iq
1458c          enddo !i   
1459c         enddo !j   
1460c       enddo !l 
1461c      endif
1462ccccccccccccccccccccccccccccccc
1463
1464c     Periodicite :
1465      do  iq = 1, nqtot
1466         do l=1, llm
1467            do j = 1, jjp1
1468               q(iip1,j,l,iq) = q(1,j,l,iq)
1469            end do
1470         end do
1471      enddo
1472     
1473!      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices,co2ice)
1474! no need to transfer "co2ice" any more; it is in qsurf(igcm_co2_ice)
1475
1476      endif !! if nqtot .ne. 0
1477
1478c-----------------------------------------------------------------------
1479c   Initialisation  h:  (passage de t -> h)
1480c-----------------------------------------------------------------------
1481
1482      DO l=1,llm
1483         DO j=1,jjp1
1484            DO i=1,iim
1485               h(i,j,l) = t(i,j,l)*((ps(i,j)/preff)**kappa)
1486            END DO
1487            h(iip1,j,l) =  h(1,j,l)
1488         END DO
1489      END DO
1490
1491
1492c***********************************************************************
1493c***********************************************************************
1494c     Fin subroutine lecture ini
1495c***********************************************************************
1496c***********************************************************************
1497
1498      deallocate(timelist)
1499      deallocate(rlonuold, rlatvold)
1500      deallocate(rlonvold, rlatuold)
1501      deallocate(apsold,bpsold)
1502      deallocate(uold)
1503      deallocate(vold)
1504      deallocate(told)
1505      deallocate(psold)
1506      deallocate(phisold)
1507      deallocate(qold)
1508      deallocate(co2iceold)
1509      deallocate(tsurfold)
1510      deallocate(emisold)
1511      deallocate(q2old)
1512      deallocate(tsoilold)
1513      deallocate(tsoiloldnew)
1514      deallocate(inertiedatold)
1515      deallocate(inertiedatoldnew)
1516      deallocate(surfithold)
1517      deallocate(mlayerold)
1518      deallocate(qsurfold)
1519      deallocate(var,varp1)
1520      deallocate(tslabold)
1521      deallocate(rnatold)
1522      deallocate(pctsrf_sicold)
1523      deallocate(tsea_iceold)
1524      deallocate(sea_iceold)
1525
1526      end
Note: See TracBrowser for help on using the repository browser.