Changeset 588


Ignore:
Timestamp:
Mar 19, 2012, 11:27:18 AM (13 years ago)
Author:
emillour
Message:

Generic GCM:
Some cleanup and bug fixing:

  • "cloudfrac" was not well written to restartfi (wrong size).
  • missing save attribute for "reffrad" in physiq.F90.
  • cleanup recomputation of surface pressure in newstart and change loop order in interp_horiz (which "fixes" an odd behaviour which fills some arrays with zeros, but only when using some versions of ifort!)

EM

Location:
trunk/LMDZ.GENERIC
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r543 r588  
    599599- Temperature grid for Planck calculations can now be refined through the parameter NTfac in radinc_h.
    600600  Default is NTfac = 1.0D-1, i.e. Delta T = 0.1 K
     601
     602== 16/03/2012 == JL
     603- Removed cpp3D and nonideal stuff.
     604
     605== 19/03/2012 == EM
     606Some cleanup and bug fixing:
     607- "cloudfrac" was not well written to restartfi (wrong size).
     608- missing save attribute for "reffrad" in physiq.F90.
     609- cleanup recomputation of surface pressure in newstart and change loop order
     610  in interp_horiz (which "fixes" an odd behaviour which fills some arrays with
     611  zeros, but only when using some versions of ifort!)

  • trunk/LMDZ.GENERIC/libf/dyn3d/interp_horiz.F

    r135 r588  
    1919c  """""""""
    2020       
    21        INTEGER imo, jmo ! dimensions ancienne grille (input)
    22        INTEGER imn,jmn  ! dimensions nouvelle grille (input)
     21       INTEGER,INTENT(IN) :: imo, jmo ! dimensions ancienne grille (input)
     22       INTEGER,INTENT(IN) :: imn,jmn  ! dimensions nouvelle grille (input)
    2323
    24        REAL rlonuo(imo+1)     !  Latitude et
    25        REAL rlatvo(jmo)       !  longitude des
    26        REAL rlonun(imn+1)     !  bord des
    27        REAL rlatvn(jmn)     !  cases "scalaires" (input)
     24       REAL,INTENT(IN) :: rlonuo(imo+1)     !  Latitude et
     25       REAL,INTENT(IN) :: rlatvo(jmo)       !  longitude des
     26       REAL,INTENT(IN) :: rlonun(imn+1)     !  bord des
     27       REAL,INTENT(IN) :: rlatvn(jmn)     !  cases "scalaires" (input)
    2828
    29        INTEGER lm ! dimension verticale (input)
    30        REAL varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
    31        REAL varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
     29       INTEGER,INTENT(IN) :: lm ! dimension verticale (input)
     30       REAL,INTENT(IN) :: varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
     31       REAL,INTENT(OUT) :: varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
    3232
    3333c Autres variables
    … …  
    3939       INTEGER ii,jj,l
    4040
    41        REAL airen ((imnmx2+1)*(jmnmx2+1)) ! aire dans la nouvelle grille
     41       REAL,SAVE :: airen ((imnmx2+1)*(jmnmx2+1)) ! aire dans la nouvelle grille
    4242       REAL airentotn   ! aire totale pole nord dans la nouvelle grille
    4343       REAL airentots   ! aire totale pole sud dans la nouvelle grille
    … …  
    4848c + des pouiemes (cas ou une maille est a cheval sur 2 ou 4 mailles)
    4949c  Il y a un test dans iniinterp_h pour s'assurer que ktotal < kmax
    50        INTEGER kmax, k, ktotal
     50       INTEGER kmax, k
     51       integer,save :: ktotal
    5152       parameter (kmax = 360*179 + 200000)
    5253c      parameter (kmax = 360*179 + 40000)
    5354
    54        INTEGER iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
    55        REAL intersec(kmax)
    56        REAL R
     55       INTEGER,SAVE :: iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
     56       REAL,SAVE :: intersec(kmax)
     57       REAL r
    5758       REAL totn, tots
    58        integer prev_sumdim
    59        save prev_sumdim
    60        data prev_sumdim /0/
    61        
     59       integer,save :: prev_sumdim=0
    6260
    63        logical firsttest, aire_ok
    64        save firsttest
    65        data firsttest /.true./
    66        data aire_ok /.true./
     61       logical,save :: firsttest=.true. , aire_ok=.true.
    6762
    68        integer imoS,jmoS,imnS,jmnS
    69        save imoS,jmoS,imnS,jmnS
    70        save ktotal,iik,jjk,jk,ik,intersec,airen
    71        REAL pi
     63       integer,save :: imoS,jmoS,imnS,jmnS
    7264
    7365c Test dimensions imnmx2 jmnmx2
    … …  
    115107      end if
    116108
    117       do l=1,lm
    118        do jj =1 , jmn+1
    119         do ii=1, imn+1
    120           varn(ii,jj,l) =0.
    121         end do
    122        end do
    123       end do
     109! initialize varn() to zero
     110      varn(1:imn+1,1:jmn+1,1:lm)=0.
    124111       
    125112c Interpolation horizontale
    … …  
    128115c de l'ancienne et la  nouvelle grille
    129116c
    130      
    131       do k=1,ktotal
    132         do l=1,lm
     117! Ehouarn 2012: for some strange reason, with ifort v12.x,
     118!               when the order of the loop below is changed
     119!               values of varn(:,:,l=2...) are then sometimes remain zero!   
     120      do l=1,lm
     121        do k=1,ktotal
    133122         varn(iik(k),jjk(k),l) = varn(iik(k),jjk(k),l)
    134123     &   + varo(ik(k), jk(k),l)*intersec(k)/airen(iik(k)
    … …  
    161150         end do
    162151
    163       ENDDO
     152      ENDDO ! of DO l=1, lm
    164153           
    165154
    … …  
    167156c  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST
    168157      if (firsttest) then
    169       pi=2.*asin(1.)
    170158      firsttest = .false.
    171159      write (*,*) 'INTERP. HORIZ. : TEST SUR LES AIRES:'
    … …  
    192180      end do
    193181      if (aire_ok) write(*,*) 'INTERP. HORIZ. : AIRES OK'
    194       endif
     182      endif ! of if (firsttest)
    195183
    196184c FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST
    … …  
    198186
    199187
    200         return
    201188        end

  • trunk/LMDZ.GENERIC/libf/dyn3d/lect_start_archive.F

    r253 r588  
    4646      INTEGER   imold,jmold,lmold,nsoilold,nqold
    4747
    48 c et autres:
    49 c----------
    50       INTEGER lnblnk
    51       EXTERNAL lnblnk
    5248
    5349c Variables pour les lectures des fichiers "ini"
    … …  
    109105c     REAL year_day,periheli,aphelie,peri_day
    110106c     REAL obliquit,z0,emin_turb,lmixmin
    111 c     REAL emissiv,emisice(2),albedice(2),tauvis
     107c     REAL emissiv,emisice(2),albedice(2)
    112108c     REAL iceradius(2) , dtemisice(2)
    113109
    … …  
    994990      co2icetotal = 0.
    995991      if (igcm_co2_ice.ne.0) then
    996       DO j=1,jjp1
     992        ! recast surface CO2 ice on new grid
     993        call interp_horiz(qsurfold(1,1,igcm_co2_ice),
     994     &                  qsurfs(1,1,igcm_co2_ice),
     995     &                  imold,jmold,iim,jjm,1,
     996     &                  rlonuold,rlatvold,rlonu,rlatv)
     997       DO j=1,jjp1
    997998         DO i=1,iim
    998999           !co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
    9991000           co2icetotal=co2icetotal+qsurfS(i,j,igcm_co2_ice)*aire(i,j)
    10001001         END DO
    1001       END DO
     1002       END DO
     1003      else
     1004        write(*,*) "Warning: No co2_ice surface tracer"
    10021005      endif
    10031006
    … …  
    10091012      write(*,*)'Ancienne grille: masse de la glace CO2:',co2icetotalold
    10101013      write(*,*)'Nouvelle grille: masse de la glace CO2:',co2icetotal
     1014      if (co2icetotalold.ne.0.) then
    10111015      write(*,*)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
     1016      endif
    10121017      write(*,*)
    10131018

  • trunk/LMDZ.GENERIC/libf/dyn3d/newstart.F

    r535 r588  
    307307      pa=tab_cntrl(17)  ! reference pressure at which coord is purely pressure
    308308      !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
    309       if (preff.eq.0) then
    310         preff=610
    311         pa=20
    312       endif
    313309      write(*,*) "Newstart: preff=",preff," pa=",pa
    314310      yes=' '
    … …  
    16711667            cloudfrac(ig,l)=0.5
    16721668         ENDDO
     1669! Ehouarn, march 2012: also add some initialisation for hice
     1670         hice(ig)=0.0
    16731671      ENDDO
    16741672     
    … …  
    16821680!      ENDDO
    16831681
    1684 
    16851682c=======================================================================
    16861683c   Correct pressure on the new grid (menu 0)
    16871684c=======================================================================
     1685
    16881686
    16891687      if ((choix_1.eq.0).and.(.not.flagps0)) then
    16901688        r = 1000.*8.31/mugaz
    16911689
    1692         phi0=0.0
    16931690        do j=1,jjp1
    16941691          do i=1,iip1
    1695              phi0 = phi0+phis(i,j)*aire(i,j)
    1696           end do
    1697         end do
    1698         phi0=phi0/airetot
    1699 
    1700         do j=1,jjp1
    1701           do i=1,iip1
    1702              ps(i,j) = ps(i,j) *
    1703 !     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
    1704      .            exp((phi0-phis(i,j)) /
     1692             ps(i,j) = ps(i,j) *
     1693     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
    17051694     .                                  (t(i,j,1) * r))
    17061695          end do
    17071696        end do
    1708  
    1709 !     we must renormalise pressures AGAIN, because exp(-phi) is nonlinear
    1710         ptot=0.0
     1697
     1698c periodicite de ps en longitude
    17111699        do j=1,jjp1
    1712            do i=1,iip1
    1713               ptot=ptot+ps(i,j)*aire(i,j)
    1714            enddo
    1715         enddo
    1716         do j=1,jjp1
    1717            do i=1,iip1
    1718               ps(i,j)=ps(i,j)*ptotn/ptot
    1719            enddo
    1720         enddo
    1721        
    1722 !     periodicity of surface ps in longitude
    1723         do j=1,jjp1
    1724            ps(1,j) = ps(iip1,j)
     1700          ps(1,j) = ps(iip1,j)
    17251701        end do
    1726        
    17271702      end if
    17281703
     1704         
    17291705c=======================================================================
    17301706c=======================================================================
    … …  
    17331709c    Initialisation de la physique / ecriture de newstartfi :
    17341710c=======================================================================
     1711
    17351712
    17361713      CALL inifilr
    … …  
    17641741      endif
    17651742
    1766 
    17671743C Calcul intermediaire
    17681744c
    … …  
    17981774      CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
    17991775     *                phi,w, pbaru,pbarv,day_ini+time )
    1800 c     CALL caldyn
    1801 c    $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
    1802 c    $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
    1803 
     1776
     1777         
    18041778      CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
    18051779      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)

  • trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90

    r538 r588  
    8080         R     = 8.314511E+0 *1000.E+0/mugaz
    8181         rcp   = R/cpp
    82       elseif((cpp.ne.cpp_c) .or. (mugaz.ne.mugaz_c))then
     82!      elseif((cpp.ne.cpp_c) .or. (mugaz.ne.mugaz_c))then
     83         ! Ehouarn: tolerate a small mismatch between computed/stored values
     84       elseif((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
     85              (abs(1.-mugaz/mugaz_c).gt.1.e-6)) then
    8386         if(check_cpp_match)then
    8487            print*,'Values do not match!'

  • trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

    r586 r588  
    522522               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
    523523
    524                tauaero(2*k+2,iaer)=max(temp*pweight,0.0)
    525                tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.0)
     524               tauaero(2*k+2,iaer)=max(temp*pweight,0.d0)
     525               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0)
    526526!
    527527            end do

  • trunk/LMDZ.GENERIC/libf/phystd/datareadnc.F

    r374 r588  
    4343
    4444      use datafile_mod, only: datadir
     45! to use  'getin'
     46      USE ioipsl_getincom
    4547      implicit none
    4648
    … …  
    7779      REAL        zthe(imdp1*jmdp1)
    7880
    79       INTEGER   lnblnk, ierr
    80       EXTERNAL    lnblnk
     81      INTEGER   ierr
    8182
    8283      INTEGER   unit,nvarid
    … …  
    110111c    Lecture NetCDF des donnees latitude et longitude
    111112c-----------------------------------------------------------------------
    112 
     113!    First get the corret value of dtadir, if not already done:
     114      ! default 'datadir' is set in "datadir_mod"
     115      call getin("datadir",datadir)
    113116      write(*,*) 'Choice of surface data is:'
    114117      filestring='ls '//trim(datadir)//'/*.nc'

  • trunk/LMDZ.GENERIC/libf/phystd/def_var.F90

    r135 r588  
    2727ierr=NF_REDEF(nid)
    2828
    29 print*,'in def_var.F90, dimids='
    30 print*,dimids
     29!print*,'in def_var.F90, dimids='
     30!print*,dimids
    3131
    3232! 2. Define the variable

  • trunk/LMDZ.GENERIC/libf/phystd/physdem1.F

    r253 r588  
    4848      REAL day_ini
    4949      INTEGER nsoil,nq
    50       integer ierr,idim1,idim2,idim3,idim4,idim5,nvarid
     50      integer ierr,idim1,idim2,idim3,idim4,idim5,idim6,nvarid
    5151
    5252      REAL phystep,time
    … …  
    8282      real pzsig(ngridmx),pzgam(ngridmx),pzthe(ngridmx)
    8383      integer ig
    84       INTEGER lnblnk
    85       EXTERNAL lnblnk
    8684
    8785! flag which identifies if we are using old tracer names (qsurf01,...)
    … …  
    144142      endif
    145143c
    146 !      ierr = NF_DEF_DIM (nid,"nlayer+1",llm+1,idim4)
    147144      ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
    148145      if (ierr.ne.NF_NOERR) then
    … …  
    157154        WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
    158155        write(6,*) NF_STRERROR(ierr)
     156      endif
     157
     158      ierr = NF_DEF_DIM (nid,"nlayer",llm,idim6)
     159      if (ierr.ne.NF_NOERR) then
     160        WRITE(6,*)'physdem1: Problem defining nlayer dimension'
     161        write(6,*) NF_STRERROR(ierr)
     162        call abort
    159163      endif
    160164
    … …  
    448452c Write the physical fields
    449453
    450 ! CO2 Ice Cover
    451 !
    452 !      ierr = NF_REDEF (nid)
    453 !#ifdef NC_DOUBLE
    454 !      ierr = NF_DEF_VAR (nid, "co2ice", NF_DOUBLE, 1, idim2,nvarid)
    455 !#else
    456 !      ierr = NF_DEF_VAR (nid, "co2ice", NF_FLOAT, 1, idim2,nvarid)
    457 !#endif
    458 !      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 13,
    459 !     .                        "CO2 ice cover")
    460 !      ierr = NF_ENDDEF(nid)
    461 !#ifdef NC_DOUBLE
    462 !      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,co2ice)
    463 !#else
    464 !      ierr = NF_PUT_VAR_REAL (nid,nvarid,co2ice)
    465 !#endif
    466 
    467 
    468 
    469 
    470454
    471455! Soil Thermal inertia
    … …  
    574558      ierr = NF_REDEF (nid)
    575559#ifdef NC_DOUBLE
    576       ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 1, idim2,nvarid)
    577 #else
    578       ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 1, idim2,nvarid)
     560      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 2,
     561     & (/idim2,idim6/),nvarid)
     562#else
     563      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 2,
     564     & (/idim2,idim6/),nvarid)
    579565#endif
    580566      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,

  • trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

    r586 r588  
    209209      real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
    210210
    211       real reffrad(ngridmx,nlayermx,naerkind)                ! aerosol effective radius (m)
     211      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
    212212
    213213!     Tendencies due to various processes:
    214214      real dqsurf(ngridmx,nqmx)
    215       real zdtlw(ngridmx,nlayermx)                            ! (K/s)
    216       real zdtsw(ngridmx,nlayermx)                            ! (K/s)
    217       save zdtlw, zdtsw
     215      real,save :: zdtlw(ngridmx,nlayermx)                    ! (K/s)
     216      real,save :: zdtsw(ngridmx,nlayermx)                    ! (K/s)
    218217
    219218      real cldtlw(ngridmx,nlayermx)                           ! (K/s) LW heating rate for clear areas
    … …  
    17241723            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
    17251724               ! to record global radiative balance
    1726                open(92,file="rad_bal.out",form='formatted',access='append')
     1725               open(92,file="rad_bal.out",form='formatted',position='append')
    17271726               write(92,*) zday,ISR/totarea,ASR/totarea,OLR/totarea
    17281727               close(92)
    1729                open(93,file="tem_bal.out",form='formatted',access='append')
     1728               open(93,file="tem_bal.out",form='formatted',position='append')
    17301729               write(93,*) zday,Ts1,Ts2,Ts3,TsS
    17311730               close(93)
    … …  
    17391738!     ------------------------------------------------------------------
    17401739      if(testradtimes)then
    1741          open(38,file="tau_phys.out",form='formatted',access='append')
     1740         open(38,file="tau_phys.out",form='formatted',position='append')
    17421741         ig=1
    17431742         do l=1,nlayer
    … …  
    18131812            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
    18141813               ! to record global water balance
    1815                open(98,file="h2o_bal.out",form='formatted',access='append')
     1814               open(98,file="h2o_bal.out",form='formatted',position='append')
    18161815               write(98,*) zday,icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
    18171816               close(98)

  • trunk/LMDZ.GENERIC/libf/phystd/tabfi.F

    r253 r588  
    176176      mugaz = tab_cntrl(tab0+8)
    177177      rcp = tab_cntrl(tab0+9)
     178      cpp=(8.314511/(mugaz/1000.0))/rcp
    178179      daysec = tab_cntrl(tab0+10)
    179180      dtphys = tab_cntrl(tab0+11)
    … …  
    502503          write(*,*)
    503504          write(*,*) ' mugaz (new value):',mugaz
    504           r=8.314/(mugaz/1000.0)
     505          r=8.314511/(mugaz/1000.0)
    505506          write(*,*) ' R (new value):',r
    506507
    … …  
    512513          write(*,*)
    513514          write(*,*) ' rcp (new value):',rcp
     515          r=8.314511/(mugaz/1000.0)
    514516          cpp=r/rcp
    515517          write(*,*) ' cpp (new value):',cpp
Note: See TracChangeset for help on using the changeset viewer.