SUBROUTINE lect_start_archive(date,tsurf,tsoil,emis,q2, & t,ucov,vcov,ps,co2ice,h,phisold_newgrid, & q,qsurf,surfith,nid) c======================================================================= c c c Auteur: 05/1997 , 12/2003 : coord hybride FF c ------ c c c Objet: Lecture des variables d'un fichier "start_archive" c Plus besoin de régler ancienne valeurs grace c a l'allocation dynamique de memoire (Yann Wanherdrick) c c c c======================================================================= implicit none #include "dimensions.h" #include "dimphys.h" #include "surfdat.h" #include "comsoil.h" #include "dimradmars.h" #include "yomaer.h" #include "planete.h" #include "paramet.h" #include "comconst.h" #include "comvert.h" #include "comgeom2.h" #include "control.h" #include "logic.h" #include "description.h" #include "ener.h" #include "temps.h" #include "lmdstd.h" #include "netcdf.inc" !#include "tracer.h" #include"advtrac.h" c======================================================================= c Declarations c======================================================================= c Old variables dimensions (from file) c------------------------------------ INTEGER imold,jmold,lmold,nsoilold,nqold c et autres: c---------- INTEGER lnblnk EXTERNAL lnblnk c Variables pour les lectures des fichiers "ini" c-------------------------------------------------- ! INTEGER sizei, integer timelen,dimid ! INTEGER length ! parameter (length = 100) INTEGER tab0 INTEGER isoil,iq,iqmax CHARACTER*2 str2 ! REAL dimfirst(4) ! tableau contenant les 1ers elements des dimensions ! REAL dimlast(4) ! tableau contenant les derniers elements des dimensions ! REAL dimcycl(4) ! tableau contenant les periodes des dimensions ! CHARACTER*120 dimsource ! CHARACTER*16 dimname ! CHARACTER*80 dimtitle ! CHARACTER*40 dimunits ! INTEGER dimtype ! INTEGER dimord(4) ! tableau contenant l''ordre ! data dimord /1,2,3,4/ ! de sortie des dimensions ! INTEGER vardim(4) REAL date INTEGER memo ! character (len=50) :: tmpname c Variable histoire c------------------ REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants REAL h(iip1,jjp1,llm),ps(iip1,jjp1) REAL q(iip1,jjp1,llm,nqmx),qtot(iip1,jjp1,llm) c autre variables dynamique nouvelle grille c------------------------------------------ c!-*- ! integer klatdat,klongdat ! PARAMETER (klatdat=180,klongdat=360) c Physique sur grille scalaire c---------------------------- c variable physique c------------------ REAL tsurf(ngridmx) ! surface temperature REAL tsoil(ngridmx,nsoilmx) ! soil temperature REAL co2ice(ngridmx) ! CO2 ice layer REAL emis(ngridmx) REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx) c REAL phisfi(ngridmx) INTEGER i,j,l INTEGER nid,nvarid c REAL year_day,periheli,aphelie,peri_day c REAL obliquit,z0,emin_turb,lmixmin c REAL emissiv,emisice(2),albedice(2),tauvis c REAL iceradius(2) , dtemisice(2) ! EXTERNAL RAN1 ! REAL RAN1 ! EXTERNAL geopot,inigeom integer ierr ! integer ismin ! external ismin ! CHARACTER*80 datapath integer, dimension(4) :: start,count c Variable nouvelle grille naturelle au point scalaire c------------------------------------------------------ real us(iip1,jjp1,llm),vs(iip1,jjp1,llm) REAL phisold_newgrid(iip1,jjp1) REAL t(iip1,jjp1,llm) real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx) real inertiedatS(iip1,jjp1,nsoilmx) real co2iceS(iip1,jjp1),emisS(iip1,jjp1) REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqmx) real ptotal, co2icetotal c Var intermediaires : vent naturel, mais pas coord scalaire c----------------------------------------------------------- real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm) c Variable de l'ancienne grille c--------------------------------------------------------- real, dimension(:), allocatable :: timelist real, dimension(:), allocatable :: rlonuold, rlatvold real, dimension(:), allocatable :: rlonvold, rlatuold real, dimension(:), allocatable :: apsold,bpsold real, dimension(:), allocatable :: mlayerold real, dimension(:,:,:), allocatable :: uold,vold,told,q2old real, dimension(:,:,:), allocatable :: tsoilold,qsurfold real, dimension(:,:,:),allocatable :: tsoiloldnew ! tsoiloldnew: old soil values, but along new subterranean grid real, dimension(:,:,:), allocatable :: inertiedatold ! inertiedatoldnew: old inertia values, but along new subterranean grid real, dimension(:,:,:), allocatable :: inertiedatoldnew real, dimension(:,:), allocatable :: psold,phisold real, dimension(:,:), allocatable :: co2iceold,tsurfold real, dimension(:,:), allocatable :: emisold real, dimension(:,:,:,:), allocatable :: qold real tab_cntrl(100) real ptotalold, co2icetotalold logical :: olddepthdef=.false. ! flag ! olddepthdef=.true. if soil depths are in 'old' (unspecified) format logical :: depthinterpol=.false. ! flag ! depthinterpol=.true. if interpolation will be requiered logical :: therminertia_3D=.true. ! flag ! therminertia_3D=.true. if thermal inertia is 3D and read from datafile c Variable intermediaires iutilise pour l'extrapolation verticale c---------------------------------------------------------------- real, dimension(:,:,:), allocatable :: var,varp1 real, dimension(:), allocatable :: oldgrid, oldval real, dimension(:), allocatable :: newval ! real, dimension(:), allocatable :: oldmlayer real surfith(iip1,jjp1) ! surface thermal inertia ! real surfithfi(ngridmx) ! surface thermal inertia at old horizontal grid resolution real, dimension(:,:), allocatable :: surfithold ! flag which identifies if archive file is using old tracer names (qsurf01,...) logical :: oldtracernames=.false. integer :: counter character(len=30) :: txt ! to store some text real :: tmpval ! to store a temporary variable/value c======================================================================= ! 0. Preliminary stuff ! check if tracers follow old naming convention (q01, q02, q03, ...) counter=0 do iq=1,nqmx txt= " " write(txt,'(a1,i2.2)')'q',iq ierr=NF_INQ_VARID(nid,txt,nvarid) if (ierr.ne.NF_NOERR) then ! did not find old tracer name exit ! might as well stop here else ! found old tracer name counter=counter+1 endif enddo if (counter.eq.nqmx) then write(*,*) "lect_start_archive: tracers seem to follow old ", & "naming convention (q01, q02,...)" oldtracernames=.true. endif !----------------------------------------------------------------------- ! 1. Read data dimensions (i.e. size and length) !----------------------------------------------------------------------- ! 1.2 Read the various dimension lengths of data in file ierr= NF_INQ_DIMID(nid,"Time",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"temps",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,timelen) ierr= NF_INQ_DIMID(nid,"latitude",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"rlatu",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,jmold) jmold=jmold-1 ierr= NF_INQ_DIMID(nid,"longitude",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"rlonv",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,imold) imold=imold-1 ierr= NF_INQ_DIMID(nid,"altitude",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"sig_s",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,lmold) nqold=0 do write(str2,'(i2.2)') nqold+1 ierr= NF_INQ_VARID(nid,'q'//str2,dimid) ! write(*,*) 'q'//str2 if (ierr.eq.NF_NOERR) then nqold=nqold+1 else exit endif enddo ! 1.2.1 find out the # of subsurface_layers nsoilold=0 !dummy initialisation ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid) if (ierr.eq.NF_NOERR) then ierr=NF_INQ_DIMLEN(nid,dimid,nsoilold) if (ierr.ne.NF_NOERR) then write(*,*)'lec_start_archive: ', & 'Failed reading subsurface_layers length' endif else write(*,*)"lec_start_archive: did not find subsurface_layers" endif if (nsoilold.eq.0) then ! 'old' archive format; ! must use Tg//str2 fields to compute nsoilold write(*,*)"lec_start_archive: building nsoilold from Tg fields" do write(str2,'(i2.2)') nsoilold+1 ierr=NF_INQ_VARID(nid,'Tg'//str2,dimid) if (ierr.eq.NF_NOERR) then nsoilold=nsoilold+1 else exit endif enddo endif if (nsoilold.ne.nsoilmx) then ! interpolation will be required depthinterpol=.true. endif ! 1.3 Report dimensions write(*,*) "Start_archive dimensions:" write(*,*) "longitude: ",imold write(*,*) "latitude: ",jmold write(*,*) "altitude: ",lmold write(*,*) "tracers: ",nqold write(*,*) "subsurface_layers: ",nsoilold if (depthinterpol) then write(*,*) " => Warning, nsoilmx= ",nsoilmx write(*,*) ' which implies that you want subterranean interpola &tion.' write(*,*) ' Otherwise, set nsoilmx -in dimphys.h- to: ',nsoilold endif write(*,*) "time lenght: ",timelen write(*,*) !----------------------------------------------------------------------- ! 2. Allocate arrays to store datasets !----------------------------------------------------------------------- allocate(timelist(timelen)) allocate(rlonuold(imold+1), rlatvold(jmold)) allocate(rlonvold(imold+1), rlatuold(jmold+1)) allocate (apsold(lmold),bpsold(lmold)) allocate(uold(imold+1,jmold+1,lmold)) allocate(vold(imold+1,jmold+1,lmold)) allocate(told(imold+1,jmold+1,lmold)) allocate(psold(imold+1,jmold+1)) allocate(phisold(imold+1,jmold+1)) allocate(qold(imold+1,jmold+1,lmold,nqmx)) allocate(co2iceold(imold+1,jmold+1)) allocate(tsurfold(imold+1,jmold+1)) allocate(emisold(imold+1,jmold+1)) allocate(q2old(imold+1,jmold+1,lmold+1)) ! allocate(tsoilold(imold+1,jmold+1,nsoilmx)) allocate(tsoilold(imold+1,jmold+1,nsoilold)) allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx)) allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx)) ! surface thermal inertia at old horizontal grid resolution allocate(surfithold(imold+1,jmold+1)) allocate(mlayerold(nsoilold)) allocate(qsurfold(imold+1,jmold+1,nqmx)) allocate(var (imold+1,jmold+1,llm)) allocate(varp1 (imold+1,jmold+1,llm+1)) write(*,*) 'q2',ngridmx,nlayermx+1 write(*,*) 'q2S',iip1,jjp1,llm+1 write(*,*) 'q2old',imold+1,jmold+1,lmold+1 !----------------------------------------------------------------------- ! 3. Read time-independent data !----------------------------------------------------------------------- C----------------------------------------------------------------------- c 3.1. Lecture du tableau des parametres du run c (pour la lecture ulterieure de "ptotalold" et "co2icetotalold") c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "controle", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Lect_start_archive: is missing" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) #else ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort ENDIF c tab0 = 50 c----------------------------------------------------------------------- c 3.2 Lecture des longitudes et latitudes c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "rlonv", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: is missing" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort ENDIF c ierr = NF_INQ_VARID (nid, "rlatu", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: is missing" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort ENDIF c ierr = NF_INQ_VARID (nid, "rlonu", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: is missing" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort ENDIF c ierr = NF_INQ_VARID (nid, "rlatv", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: is missing" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort ENDIF c c----------------------------------------------------------------------- c 3.3. Lecture des niveaux verticaux c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "aps", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: is missing" apsold=0 PRINT*, " set to 0" ELSE #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, apsold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " ENDIF ENDIF c ierr = NF_INQ_VARID (nid, "bps", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: is missing" PRINT*, "It must be an old start_archive, lets look for sig_s" ierr = NF_INQ_VARID (nid, "sig_s", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Nothing to do..." CALL abort ENDIF ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort END IF c----------------------------------------------------------------------- c 3.4 Read Soil layers depths c----------------------------------------------------------------------- ierr=NF_INQ_VARID(nid,"soildepth",nvarid) if (ierr.ne.NF_NOERR) then write(*,*)'lect_start_archive: Could not find ' write(*,*)' => Assuming this is an archive in old format' olddepthdef=.true. depthinterpol=.true. ! this is how soil depth was defined in ye old days do isoil=1,nsoilold mlayerold(isoil)=sqrt(887.75/3.14)*((2.**(isoil-0.5))-1.) enddo else #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold) #else ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold) #endif if (ierr .NE. NF_NOERR) then PRINT*, "lect_start_archive: Failed reading " CALL abort endif endif !of if(ierr.ne.NF_NOERR) ! Read (or build) mlayer() if (depthinterpol) then ! Build (default) new soil depths (mlayer(:) is in comsoil.h), ! as in soil_settings.F write(*,*)' => Building default soil depths' do isoil=0,nsoilmx-1 mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) enddo write(*,*)' => mlayer: ',mlayer ! Also build (default) new soil interlayer depth layer(:) do isoil=1,nsoilmx layer(isoil)=sqrt(mlayer(0)*mlayer(1))* & ((mlayer(1)/mlayer(0))**(isoil-1)) enddo write(*,*)' => layer: ',layer else ! read mlayer() from file #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer) #else ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer) #endif if (ierr .NE. NF_NOERR) then PRINT*, "lect_start_archive: Failed reading " CALL abort endif ! Also build (default) soil interlayer depth layer(:) do isoil=1,nsoilmx layer(isoil)=sqrt(mlayer(0)*mlayer(1))* & ((mlayer(1)/mlayer(0))**(isoil-1)) enddo endif ! of if (depthinterpol) c----------------------------------------------------------------------- c 3.5 Read Soil thermal inertia c----------------------------------------------------------------------- ierr=NF_INQ_VARID(nid,"inertiedat",nvarid) if (ierr.ne.NF_NOERR) then write(*,*)'lect_start_archive: Could not find ' write(*,*)' => Assuming this is an archive in old format' therminertia_3D=.false. write(*,*)' => Thermal inertia will be read from reference file' volcapa=1.e6 write(*,*)' and soil volumetric heat capacity is set to ', & volcapa else #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold) #else ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold) #endif if (ierr .NE. NF_NOERR) then PRINT*, "lect_start_archive: Failed reading " CALL abort endif endif c----------------------------------------------------------------------- c 3.6 Lecture geopotentiel au sol c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "phisinit", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: is missing" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, phisold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort ENDIF C----------------------------------------------------------------------- c lecture de "ptotalold" et "co2icetotalold" c----------------------------------------------------------------------- ptotalold = tab_cntrl(tab0+49) co2icetotalold = tab_cntrl(tab0+50) c----------------------------------------------------------------------- c 4. Lecture du temps et choix c----------------------------------------------------------------------- c lecture du temps c ierr = NF_INQ_DIMID (nid, "Time", nvarid) IF (ierr .NE. NF_NOERR) THEN ierr = NF_INQ_DIMID (nid, "temps", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: