Ignore:
Timestamp:
Aug 2, 2011, 11:13:07 AM (13 years ago)
Author:
emillour
Message:

Generic GCM

  • Massive update to version 0.7

EM+RW

Location:
trunk/LMDZ.GENERIC/libf/dyn3d
Files:
2 deleted
11 edited

Legend:

Unmodified
Added
Removed

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

    r135 r253  
    135135c     ENDDO
    136136
     137! increment tracers
    137138      DO iq = 1, nq
    138139         DO k = 1,llm
    139140            DO j = 1,ip1jmp1
     141!               pq(j,k,iq)=  pdqfi(j,k,iq) * pdt
    140142               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
    141                pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt )
     143               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt ) ! forbid negative tracer values
    142144            ENDDO
    143145         ENDDO
    144146      ENDDO
     147!      print*,'MAJOR TRACER MOD in addfi for test'
    145148                                       
    146149

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

    r135 r253  
    11      SUBROUTINE disvert
     2
     3! to use  'getin'
     4      USE ioipsl_getincom
    25
    36c    Auteur :  F. Forget Y. Wanherdrick, P. Levan
    … …  
    1114#include "comconst.h"
    1215#include "logic.h"
     16
     17#include "callkeys.h"
     18
     19
    1320c
    1421c=======================================================================
    … …  
    2633      REAL zsig(llm)
    2734      INTEGER np,ierr
    28       integer :: ierr1,ierr2,ierr3,ierr4
     35      integer :: ierr1,ierr2,ierr3,ierr4 
    2936      REAL x
    3037
    … …  
    3946      real z, ps,p
    4047
     48
     49      real psurf, hmax ! for autozlevs
     50
    4151c
    4252c-----------------------------------------------------------------------
    … …  
    4454      pi=2.*ASIN(1.)
    4555
    46 ! Ouverture possible de fichiers typiquement E.T.
    47 
    48          open(99,file="esasig.def",status='old',form='formatted',
    49      s   iostat=ierr2)
    50          if(ierr2.ne.0) then
    51               close(99)
    52               open(99,file="z2sig.def",status='old',form='formatted',
    53      s        iostat=ierr4)
    54          endif
    55 
    56 c-----------------------------------------------------------------------
    57 c   cas 1 on lit les options dans esasig.def:
    58 c   ----------------------------------------
    59 
    60       IF(ierr2.eq.0) then
    61 
    62 c        Lecture de esasig.def :
    63 c        Systeme peu souple, mais qui respecte en theorie
    64 c        La conservation de l'energie (conversion Energie potentielle
    65 c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
    66 
    67          PRINT*,'*****************************'
    68          PRINT*,'WARNING Lecture de esasig.def'
    69          PRINT*,'*****************************'
    70          READ(99,*) h
    71          READ(99,*) dz0
    72          READ(99,*) dz1
    73          READ(99,*) nhaut
    74          CLOSE(99)
    75 
    76          dz0=dz0/h
    77          dz1=dz1/h
    78 
    79          sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
    80 
    81          esig=1.
    82 
    83          do l=1,20
    84             esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
    85          enddo
    86          csig=(1./sig1-1.)/(exp(esig)-1.)
    87 
    88          DO L = 2, llm
    89             zz=csig*(exp(esig*(l-1.))-1.)
    90             sig(l) =1./(1.+zz)
    91      &      * tanh(.5*(llm+1-l)/nhaut)
    92          ENDDO
    93          sig(1)=1.
    94          sig(llm+1)=0.
    95          quoi      = 1. + 2.* kappa
    96          s( llm )  = 1.
    97          s(llm-1) = quoi
    98          IF( llm.gt.2 )  THEN
    99             DO  ll = 2, llm-1
    100                l         = llm+1 - ll
    101                quand     = sig(l+1)/ sig(l)
    102                s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
    103             ENDDO
    104          END IF
    105 c
    106          snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
    107          DO l = 1, llm
    108             s(l)    = s(l)/ snorm
    109          ENDDO
    110 
    111 c-----------------------------------------------------------------------
    112 c   cas 2 on lit les options dans z2sig.def:
    113 c   ----------------------------------------
    114 
    115       ELSE IF(ierr4.eq.0) then
     56      open(99,file="z2sig.def",status='old',form='formatted',
     57     s     iostat=ierr4)
     58
     59
     60      autozlevs=.false.
     61      PRINT *,'Auto-discretise vertical levels ?'
     62      call getin("autozlevs",autozlevs)
     63      write(*,*) " autozlevs = ", autozlevs
     64
     65      write(*,*)"Operate in kastprof mode?"
     66      kastprof=.false.
     67      call getin("kastprof",kastprof)
     68      write(*,*)" kastprof = ",kastprof
     69
     70      print*,'kast=',kastprof
     71
     72      pceil=100.0 ! Pascals
     73      PRINT *,'Ceiling pressure (Pa) ?'
     74      call getin("pceil",pceil)
     75      write(*,*) " pceil = ", pceil
     76
     77      if(autozlevs.and.iim.gt.1)then
     78         print*,'autozlevs no good in 3D...'
     79         call abort
     80      endif
     81
     82      if(kastprof.and.iim.gt.1)then
     83         print*,'kastprof no good in 3D...'
     84         call abort
     85      endif
     86
     87      psurf=610. ! default value for psurf
     88      PRINT *,'Surface pressure (Pa) ?'
     89      call getin("psurf",psurf)
     90      write(*,*) " psurf = ",psurf
     91
     92      if(kastprof)then
     93
     94        sig(1)=1
     95        do l=2,llm
     96                                !sig(l)=1. - real(l-1)/real(llm) ! uses linear sigma spacing
     97                                !sig(l)=exp(-real(l-1)*h/real(llm)) ! uses log sigma spacing
     98                                !sig(l)=exp(-real(l-1)*Hmax/real(llm)) ! uses log sigma spacing
     99           sig(l)=(pceil/psurf)**(real(l-1)/real(llm)) ! uses log sigma spacing
     100           
     101        end do
     102        sig(llm+1)=0
     103
     104      elseIF(ierr4.eq.0)then
    116105         PRINT*,'****************************'
    117106         PRINT*,'Lecture de z2sig.def'
    … …  
    124113         CLOSE(99)
    125114
    126          sig(1) =1
    127          do l=2,llm
     115
     116         if(autozlevs)then
     117            open(91,file="z2sig.def",form='formatted')
     118            read(91,*) h
     119            DO l=1,llm-2
     120               read(91,*) Hmax
     121            enddo
     122            close(91)
     123
     124            print*,'Hmax = ',Hmax,' km'
     125            print*,'Auto-shifting h in disvert.F to:'
     126!            h = Hmax / log(psurf/100.0)
     127            h = Hmax / log(psurf/pceil)
     128            print*,'h = ',h,' km'
     129        endif
     130       
     131        sig(1)=1
     132        do l=2,llm
    128133           sig(l) = 0.5 * ( exp(-zsig(l)/h) + exp(-zsig(l-1)/h) )
    129          end do
    130          sig(llm+1) =0
     134        end do
     135        sig(llm+1)=0
    131136
    132137c-----------------------------------------------------------------------
    133138      ELSE
    134          write(*,*) 'didn t you forget something ??? '
    135          write(*,*) 'We need the file  z2sig.def ! (OR esasig.def) '
     139         write(*,*) 'didn_t you forget something?'
     140         write(*,*) 'We need the file  z2sig.def!'! (OR esasig.def) '
    136141         stop
    137142      ENDIF
    … …  
    206211      else
    207212         bps(llm) = bps(llm-1)**2 / bps(llm-2)
    208          aps(llm) = 0.
     213         aps(llm) = 0. ! what the hell is this???
    209214      end if
    210215
    … …  
    213218      PRINT *,' APs'
    214219      PRINT *,  aps
     220
    215221
    216222      DO l = 1, llm

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

    r135 r253  
    5151      INTEGER ierr, nid, nvarid, nqold
    5252      CHARACTER  str3*3,yes*1
     53
     54
     55!     added by RW for test
     56      real pmean,airetot
     57      integer ij
    5358
    5459c-----------------------------------------------------------------------
    … …  
    421426     *rrage est differente de la valeur  dtinteg =',i4//)
    422427
     428
     429
     430
     431
     432
    423433      RETURN
    424434      END

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

    r135 r253  
    132132! flag to set/remove calls to groupeun
    133133      logical callgroupeun
    134       parameter (callgroupeun = .true.)
     134      parameter (callgroupeun = .false.)
    135135
    136136!     added by RDW for tests without dynamics
    … …  
    138138      parameter (calldyn = .true.)
    139139
     140!     added by RW for test
     141!      real pmean,airetot
     142
     143!     added by FF for dissipation / energy conservation tests
     144      logical dissip_conservative
     145      parameter (dissip_conservative = .false.)
     146      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
     147!     d (theta)/ d t (pot. temperature) due to the creation
     148!     of thermal energy by dissipation
     149      REAL dtetaecdt(ip1jmp1,llm)
     150      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
     151      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
     152
    140153
    141154c-----------------------------------------------------------------------
    … …  
    151164      CALL iniadvtrac(nq,numvanle)
    152165
     166
    153167      CALL dynetat0("start.nc",nqmx,vcov,ucov,
    154168     .              teta,q,masse,ps,phis, time_0)
    155169
    156170      CALL defrun_new( 99, .TRUE. )
    157 
    158171
    159172c  on recalcule eventuellement le pas de temps
    … …  
    176189      CALL iniconst
    177190      CALL inigeom
    178 
    179191      CALL inifilr
    180192
    … …  
    192204      CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
    193205
    194 c
     206
    195207c  numero de stockage pour les fichiers de redemarrage:
    196208
    … …  
    303315      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    304316
    305 
    306317      if(calldyn)then
    307318         CALL caldyn
    … …  
    314325c   -------------------------------------------------------------
    315326
     327
     328
    316329      if (tracer) then
    317330       IF( forward. OR . leapf )  THEN
    … …  
    360373     $            dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,
    361374     $            finvmaold)
     375          else
     376             print*,'Currently no dynamics in this GCM...'
    362377          endif
    363378
    … …  
    407422c      ajout des tendances physiques:
    408423c      ------------------------------
     424!        if(1.eq.2)then
    409425          CALL addfi( nqmx, dtphys, leapf, forward   ,
    410426     $                  ucov, vcov, teta , q   ,ps , masse,
    411427     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
     428!       else
     429!          print*,'Currently no physics in this GCM...'
     430!       endif
    412431c
    413432       ENDIF
    … …  
    436455c        Dissipation horizontale
    437456c        ~~~~~~~~~~~~~~~~~~~~~~~
     457
     458
     459         if(dissip_conservative) then
     460!     calculate kinetic energy before dissipation
     461            call covcont(llm,ucov,vcov,ucont,vcont)
     462            call enercin(vcov,ucov,vcont,ucont,ecin0)
     463         end if
     464
    438465         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
    439466
    440467         CALL addit( ijp1llm,ucov ,dudis,ucov )
    441468         CALL addit( ijmllm ,vcov ,dvdis,vcov )
     469
     470
     471         if (dissip_conservative) then
     472C           On rajoute la tendance due a la transform. Ec -> E therm. cree
     473C           lors de la dissipation
     474            call covcont(llm,ucov,vcov,ucont,vcont)
     475            call enercin(vcov,ucov,vcont,ucont,ecin)
     476            dtetaecdt(:,:) = (ecin0(:,:)-ecin(:,:))/ pk(:,:)
     477            dhdis(:,:)     = dhdis(:,:) + dtetaecdt(:,:)
     478         endif
     479
    442480         CALL addit( ijp1llm,teta ,dhdis,teta )
    443481

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

    r135 r253  
    4848
    4949
    50 !     MODIFICATION TO TEST WITHOUT TRACER ADVECTION BY RDW
    51         !do iq=1,nq
    52         !  iadv(iq)=0
    53         !enddo
    54         !print*,'TRACERS ARE NO LONGER ADVECTED IN THE GCM!!!!!!'
    55 
     50!     MODIFICATION TO TEST OTHER SCHEMES BY RDW
     51!        do iq=1,nq
     52!           iadv(iq)=1
     53!        enddo
     54!        print*,'IADV SET TO 1 IN iniadvtrac!!!!'
    5655
    5756        do iq=1,nq

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

    r135 r253  
    5151c         FF 2004       
    5252         if (pseudoalt(llm).lt.160.) then
     53!     currently disabled for the universal model!!!
    5354             fac_mid=2  ! coeff  pour dissipation aux basses/moyennes altitudes
    5455             fac_up=10 ! coeff multiplicateur pour dissipation hautes altitudes
    5556             startalt=90. ! altitude en Km de la transition mid / up
    5657             delta=20.! Intervalle (km) pour le changement mid / up
     58
     59!             fac_mid=1  ! coeff  pour dissipation aux basses/moyennes altitudes
     60!             fac_up=1 ! coeff multiplicateur pour dissipation hautes altitudes
     61!             startalt=100. ! altitude en Km de la transition mid / up
     62!             delta=20.! Intervalle (km) pour le changement mid / up
     63
    5764         else ! thermosphere model
    5865             fac_mid=2 ! coeff pour dissipation aux basses/moyennes altitudes

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

    r135 r253  
    22     &     t,ucov,vcov,ps,h,phisold_newgrid,
    33     &     q,qsurf,surfith,nid)
     4
    45c=======================================================================
    56c
    … …  
    11021103c-----------------------------------------------------------------------
    11031104!      write(*,*) 'Soil'
     1105
     1106      !print*,'Problem in lect_start_archive interpolating'
     1107      !print*,'to new resolution!!'
     1108
    11041109      ! Recast temperatures along soil depth, if necessary
    11051110      if (olddepthdef) then
    … …  
    11091114        do i=1,imold+1
    11101115         do j=1,jmold+1
     1116
     1117            !if(i.gt.iip1 .or. j.gt.jjp1)then
     1118               !print*,'Problem in lect_start_archive interpolating'
     1119               !print*,'to new resolution!!'
     1120               !call abort
     1121            !endif
     1122
    11111123           ! copy values
    1112            oldval(1)=tsurfS(i,j)
     1124           oldval(1)=tsurfold(i,j)
     1125!          oldval(1)=tsurfS(i,j)
    11131126           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
    11141127           ! build vertical coordinate
    … …  
    11421155         do j=1,jmold+1
    11431156           ! copy values
    1144            oldval(1)=tsurfS(i,j)
     1157           oldval(1)=tsurfold(i,j)
     1158!          oldval(1)=tsurfS(i,j)
    11451159           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
    11461160          ! interpolate

  • trunk/LMDZ.GENERIC/libf/dyn3d/logic.h

    r135 r253  
    33
    44      COMMON/logic/ purmats,physic,forward,leapf,apphys,grireg,
    5      *  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,hybrid
     5     *  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,hybrid,autozlevs
    66
    77      LOGICAL purmats,physic,forward,leapf,apphys,grireg,statcl,conser,
    8      * apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,hybrid
     8     * apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,hybrid,autozlevs
    99
    1010c-----------------------------------------------------------------------

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

    r135 r253  
    154154      real :: MONS_coeffN ! coeff for northern hemisphere
    155155!      real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
     156
     157!     added by RW for test
     158      real pmean, phi0
     159
     160!     added by BC for equilibrium temperature startup
     161      real teque
     162
     163!     added by BC for cloud fraction setup
     164      REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
     165      REAL totalfrac(ngridmx)
     166
     167!     added by RW for nuketharsis
     168      real fact1
     169      real fact2
     170
    156171
    157172c sortie visu pour les champs dynamiques
    … …  
    425440        CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
    426441     .        day_ini,time,
    427      .        tsurf,tsoil,emis,q2,qsurf)
    428        
     442     .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
     443     .        cloudfrac,totalfrac,hice)
     444
    429445        ! copy albedo and soil thermal inertia
    430446        do i=1,ngridmx
    … …  
    457473      write(*,*)
    458474      write(*,*) 'flat : no topography ("aquaplanet")'
     475      write(*,*) 'nuketharsis : no Tharsis bulge'
    459476      write(*,*) 'bilball : uniform albedo and thermal inertia'
    460477      write(*,*) 'coldspole : cold subsurface and high albedo at S.pole'
    … …  
    471488      write(*,*) 'watercapn : H20 ice on permanent N polar cap '
    472489      write(*,*) 'watercaps : H20 ice on permanent S polar cap '
     490      write(*,*) 'noacglac  : H2O ice across Noachis Terra'
    473491      write(*,*) 'oborealis : H2O ice across Vastitas Borealis'
     492      write(*,*) 'iceball   : Thick ice layer all over surface'
    474493      write(*,*) 'wetstart  : start with a wet atmosphere'
    475       write(*,*) 'isotherm : Isothermal Temperatures, wind set to zero'
     494      write(*,*) 'isotherm  : Isothermal Temperatures, wind set to zero'
     495      write(*,*) 'radequi   : Earth-like radiative equilibrium temperature
     496     $ profile (lat-alt) and winds set to zero'
    476497      write(*,*) 'coldstart : Start X K above the CO2 frost point and
    477498     $set wind to zero (assumes 100% CO2)'
    … …  
    479500      write(*,*) 'ptot : change total pressure'
    480501      write(*,*) 'emis : change surface emissivity'
    481       write(*,*) 'therm_ini_s: Set soil thermal inertia to reference sur
     502      write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur
    482503     &face values'
    483       write(*,*) 'subsoilice_n: Put deep underground ice layer in northe
    484      &rn hemisphere'
    485       write(*,*) 'subsoilice_s: Put deep underground ice layer in southe
    486      &rn hemisphere'
    487       write(*,*) 'mons_ice: Put underground ice layer according to MONS-
    488      &derived data'
     504!      write(*,*) 'subsoilice_n : Put deep underground ice layer in northe
     505!     &rn hemisphere'
     506!      write(*,*) 'subsoilice_s : Put deep underground ice layer in southe
     507!     &rn hemisphere'
     508!      write(*,*) 'mons_ice : Put underground ice layer according to MONS-
     509!     &derived data'
    489510
    490511        write(*,*)
    … …  
    510531     &    ' were not set to zero ! => set calllott to F'                   
    511532
    512 c        Choice for surface pressure
     533c        Choice of surface pressure
    513534         yes=' '
    514535         do while ((yes.ne.'y').and.(yes.ne.'n'))
    … …  
    533554         end if
    534555
    535 c       bilball : albedo, inertie thermique uniforme
    536 c       --------------------------------------------
     556c       'nuketharsis : no tharsis bulge for Early Mars'
     557c       ---------------------------------------------
     558        else if (modif(1:len_trim(modif)) .eq. 'nuketharsis') then
     559
     560           DO j=1,jjp1       
     561              DO i=1,iim
     562                 ig=1+(j-2)*iim +i
     563                 if(j.eq.1) ig=1
     564                 if(j.eq.jjp1) ig=ngridmx
     565
     566                 fact1=(((rlonv(i)*180./pi)+100)**2 +
     567     &                (rlatu(j)*180./pi)**2)/65**2
     568                 fact2=exp( -fact1**2.5 )
     569
     570                 phis(i,j) = phis(i,j) - (phis(i,j)+4000.*g)*fact2
     571
     572!                 if(phis(i,j).gt.2500.*g)then
     573!                    if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
     574!                       phis(i,j)=2500.*g
     575!                    endif
     576!                 endif
     577
     578              ENDDO
     579           ENDDO
     580          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     581
     582
     583c       bilball : uniform albedo, thermal inertia
     584c       -----------------------------------------
    537585        else if (modif(1:len_trim(modif)) .eq. 'bilball') then
    538586          write(*,*) 'constante albedo and iner.therm:'
    … …  
    654702             enddo
    655703          enddo
     704
     705
    656706
    657707c        Correction pour la conservation des traceurs
    … …  
    908958c      --------------------------------------------------
    909959        else if (modif(1:len_trim(modif)) .eq. 'noglacier') then
     960           if (igcm_h2o_ice.eq.0) then
     961             write(*,*) "No water ice tracer! Can't use this option"
     962             stop
     963           endif
    910964           do ig=1,ngridmx
    911965             j=(ig-2)/iim +2
    … …  
    921975c      --------------------------------------------------
    922976        else if (modif(1:len_trim(modif)) .eq. 'watercapn') then
    923            do ig=1,ngridmx
    924              j=(ig-2)/iim +2
    925               if(ig.eq.1) j=1
    926               if (rlatu(j)*180./pi.gt.80.) then
    927                    qsurf(ig,igcm_h2o_ice)=1.e5
    928                    write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
    929      &              qsurf(ig,igcm_h2o_ice)
    930                    write(*,*)'     ==> Ice mesh South boundary (deg)= ',
    931      &              rlatv(j)*180./pi
    932                 end if
    933            enddo
     977           if (igcm_h2o_ice.eq.0) then
     978             write(*,*) "No water ice tracer! Can't use this option"
     979             stop
     980           endif
     981
     982           DO j=1,jjp1       
     983              DO i=1,iim
     984                 ig=1+(j-2)*iim +i
     985                 if(j.eq.1) ig=1
     986                 if(j.eq.jjp1) ig=ngridmx
     987
     988                 if (rlatu(j)*180./pi.gt.80.) then
     989                    do isoil=1,nsoilmx
     990                       ith(i,j,isoil) = 18000. ! thermal inertia
     991                    enddo
     992                   write(*,*)'     ==> Ice mesh North boundary (deg)= ',
     993     &                   rlatv(j-1)*180./pi
     994                 end if
     995              ENDDO
     996           ENDDO
     997           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     998
     999c$$$           do ig=1,ngridmx
     1000c$$$             j=(ig-2)/iim +2
     1001c$$$              if(ig.eq.1) j=1
     1002c$$$              if (rlatu(j)*180./pi.gt.80.) then
     1003c$$$
     1004c$$$                   qsurf(ig,igcm_h2o_ice)=1.e5
     1005c$$$                   qsurf(ig,igcm_h2o_vap)=0.0!1.e5
     1006c$$$
     1007c$$$                   write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
     1008c$$$     &              qsurf(ig,igcm_h2o_ice)
     1009c$$$
     1010c$$$                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
     1011c$$$     &              rlatv(j)*180./pi
     1012c$$$                end if
     1013c$$$           enddo
    9341014
    9351015c      watercaps : H20 ice on permanent southern cap
    9361016c      -------------------------------------------------
    9371017        else if (modif(1:len_trim(modif)) .eq. 'watercaps') then
    938            do ig=1,ngridmx
    939                j=(ig-2)/iim +2
    940                if(ig.eq.1) j=1
    941                if (rlatu(j)*180./pi.lt.-80.) then
    942                    qsurf(ig,igcm_h2o_ice)=1.e5
    943                    write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
    944      &              qsurf(ig,igcm_h2o_ice)
    945                    write(*,*)'     ==> Ice mesh North boundary (deg)= ',
    946      &              rlatv(j-1)*180./pi
    947                end if
    948            enddo
    949 
    950 c       oborealis : H2O ice across Vastitas Borealis
     1018           if (igcm_h2o_ice.eq.0) then
     1019              write(*,*) "No water ice tracer! Can't use this option"
     1020              stop
     1021           endif
     1022
     1023           DO j=1,jjp1       
     1024              DO i=1,iim
     1025                 ig=1+(j-2)*iim +i
     1026                 if(j.eq.1) ig=1
     1027                 if(j.eq.jjp1) ig=ngridmx
     1028
     1029                 if (rlatu(j)*180./pi.lt.-80.) then
     1030                    do isoil=1,nsoilmx
     1031                       ith(i,j,isoil) = 18000. ! thermal inertia
     1032                    enddo
     1033                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
     1034     &                   rlatv(j-1)*180./pi
     1035                 end if
     1036              ENDDO
     1037           ENDDO
     1038           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1039
     1040c$$$           do ig=1,ngridmx
     1041c$$$               j=(ig-2)/iim +2
     1042c$$$               if(ig.eq.1) j=1
     1043c$$$               if (rlatu(j)*180./pi.lt.-80.) then
     1044c$$$                  qsurf(ig,igcm_h2o_ice)=1.e5
     1045c$$$                  qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
     1046c$$$
     1047c$$$                  write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
     1048c$$$     &                 qsurf(ig,igcm_h2o_ice)
     1049c$$$                  write(*,*)'     ==> Ice mesh North boundary (deg)= ',
     1050c$$$     &                 rlatv(j-1)*180./pi
     1051c$$$               end if
     1052c$$$           enddo
     1053
     1054
     1055c       noacglac : H2O ice across highest terrain
    9511056c       --------------------------------------------
    952         else if (modif(1:len_trim(modif)) .eq. 'oborealis') then
     1057        else if (modif(1:len_trim(modif)) .eq. 'noacglac') then
     1058           if (igcm_h2o_ice.eq.0) then
     1059             write(*,*) "No water ice tracer! Can't use this option"
     1060             stop
     1061           endif
    9531062          DO j=1,jjp1       
    9541063             DO i=1,iim
    … …  
    9571066                if(j.eq.jjp1) ig=ngridmx
    9581067
    959                 if(phis(i,j).lt.-4000.)then
    960                     qsurf(ig,igcm_h2o_ice)=1.e5
    961                     phis(i,j)=-4000.0
     1068                if(phis(i,j).gt.1000.*g)then
     1069                    alb(i,j) = 0.5 ! snow value
     1070                    do isoil=1,nsoilmx
     1071                       ith(i,j,isoil) = 18000. ! thermal inertia
     1072                       ! this leads to rnat set to 'ocean' in physiq.F90
     1073                       ! actually no, because it is soil not surface
     1074                    enddo
    9621075                endif
    9631076             ENDDO
    9641077          ENDDO
     1078          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1079          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
     1080          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     1081
     1082
     1083
     1084c       oborealis : H2O oceans across Vastitas Borealis
     1085c       -----------------------------------------------
     1086        else if (modif(1:len_trim(modif)) .eq. 'oborealis') then
     1087           if (igcm_h2o_ice.eq.0) then
     1088             write(*,*) "No water ice tracer! Can't use this option"
     1089             stop
     1090           endif
     1091          DO j=1,jjp1       
     1092             DO i=1,iim
     1093                ig=1+(j-2)*iim +i
     1094                if(j.eq.1) ig=1
     1095                if(j.eq.jjp1) ig=ngridmx
     1096
     1097                if(phis(i,j).lt.-4000.*g)then
     1098!                if( (phis(i,j).lt.-4000.*g)
     1099!     &               .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
     1100!                    phis(i,j)=-8000.0*g ! proper ocean
     1101                    phis(i,j)=-4000.0*g ! scanty ocean
     1102
     1103                    alb(i,j) = 0.07 ! oceanic value
     1104                    do isoil=1,nsoilmx
     1105                       ith(i,j,isoil) = 18000. ! thermal inertia
     1106                       ! this leads to rnat set to 'ocean' in physiq.F90
     1107                       ! actually no, because it is soil not surface
     1108                    enddo
     1109                endif
     1110             ENDDO
     1111          ENDDO
     1112          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1113          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
     1114          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     1115
     1116c       iborealis : H2O ice in Northern plains
     1117c       --------------------------------------
     1118        else if (modif(1:len_trim(modif)) .eq. 'iborealis') then
     1119           if (igcm_h2o_ice.eq.0) then
     1120             write(*,*) "No water ice tracer! Can't use this option"
     1121             stop
     1122           endif
     1123          DO j=1,jjp1       
     1124             DO i=1,iim
     1125                ig=1+(j-2)*iim +i
     1126                if(j.eq.1) ig=1
     1127                if(j.eq.jjp1) ig=ngridmx
     1128
     1129                if(phis(i,j).lt.-4000.*g)then
     1130                   qsurf(ig,igcm_h2o_ice)=1.e3
     1131                endif
     1132             ENDDO
     1133          ENDDO
     1134          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1135          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
     1136          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     1137
     1138
     1139c       oceanball : H2O liquid everywhere
     1140c       ----------------------------
     1141        else if (modif(1:len_trim(modif)) .eq. 'oceanball') then
     1142           if (igcm_h2o_ice.eq.0) then
     1143             write(*,*) "No water ice tracer! Can't use this option"
     1144             stop
     1145           endif
     1146          DO j=1,jjp1       
     1147             DO i=1,iim
     1148                ig=1+(j-2)*iim +i
     1149                if(j.eq.1) ig=1
     1150                if(j.eq.jjp1) ig=ngridmx
     1151
     1152                qsurf(ig,igcm_h2o_vap)=0.0    ! for ocean, this is infinite source
     1153                qsurf(ig,igcm_h2o_ice)=0.0
     1154                alb(i,j) = 0.07 ! ocean value
     1155
     1156                do isoil=1,nsoilmx
     1157                   ith(i,j,isoil) = 18000. ! thermal inertia
     1158                   !ith(i,j,isoil) = 50. ! extremely low for test
     1159                   ! this leads to rnat set to 'ocean' in physiq.F90
     1160                enddo
     1161
     1162             ENDDO
     1163          ENDDO
     1164          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1165          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
     1166          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
     1167
     1168c       iceball : H2O ice everywhere
     1169c       ----------------------------
     1170        else if (modif(1:len_trim(modif)) .eq. 'iceball') then
     1171           if (igcm_h2o_ice.eq.0) then
     1172             write(*,*) "No water ice tracer! Can't use this option"
     1173             stop
     1174           endif
     1175          DO j=1,jjp1       
     1176             DO i=1,iim
     1177                ig=1+(j-2)*iim +i
     1178                if(j.eq.1) ig=1
     1179                if(j.eq.jjp1) ig=ngridmx
     1180
     1181                qsurf(ig,igcm_h2o_vap)=-50.    ! for ocean, this is infinite source
     1182                qsurf(ig,igcm_h2o_ice)=50.     ! == to 0.05 m of oceanic ice
     1183                alb(i,j) = 0.6 ! ice albedo value
     1184
     1185                do isoil=1,nsoilmx
     1186                   !ith(i,j,isoil) = 18000. ! thermal inertia
     1187                   ! this leads to rnat set to 'ocean' in physiq.F90
     1188                enddo
     1189
     1190             ENDDO
     1191          ENDDO
     1192          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1193          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
    9651194
    9661195c       isotherm : Isothermal temperatures and no winds
    967 c       ------------------------------------------------
     1196c       -----------------------------------------------
    9681197        else if (modif(1:len_trim(modif)) .eq. 'isotherm') then
    9691198
    … …  
    9941223          call initial0(ngridmx*(llm+1),q2)
    9951224
     1225c       radequi : Radiative equilibrium profile of temperatures and no winds
     1226c       --------------------------------------------------------------------
     1227        else if (modif(1:len_trim(modif)) .eq. 'radequi') then
     1228
     1229          write(*,*)'radiative equilibrium temperature profile'       
     1230
     1231          do ig=1, ngridmx
     1232             teque= 335.0-60.0*sin(latfi(ig))*sin(latfi(ig))-
     1233     &            10.0*cos(latfi(ig))*cos(latfi(ig))
     1234             tsurf(ig) = MAX(220.0,teque)
     1235          end do
     1236          do l=2,nsoilmx
     1237             do ig=1, ngridmx
     1238                tsoil(ig,l) = tsurf(ig)
     1239             end do
     1240          end do
     1241          DO j=1,jjp1
     1242             DO i=1,iim
     1243                Do l=1,llm
     1244                   teque=335.-60.0*sin(rlatu(j))*sin(rlatu(j))-
     1245     &                  10.0*cos(rlatu(j))*cos(rlatu(j))
     1246                   Tset(i,j,l)=MAX(220.0,teque)
     1247                end do
     1248             end do
     1249          end do
     1250          flagtset=.true.
     1251          call initial0(llm*ip1jmp1,ucov)
     1252          call initial0(llm*ip1jm,vcov)
     1253          call initial0(ngridmx*(llm+1),q2)
    9961254
    9971255c       coldstart : T set 1K above CO2 frost point and no winds
    … …  
    10691327!        enddo
    10701328
    1071 !       subsoilice_n: Put deep ice layer in northern hemisphere soil
    1072 !       ------------------------------------------------------------
    1073 
    1074         else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
    1075 
    1076          write(*,*)'From which latitude (in deg.), up to the north pole,
    1077      &should we put subterranean ice?'
    1078          ierr=1
    1079          do while (ierr.ne.0)
    1080           read(*,*,iostat=ierr) val
    1081           if (ierr.eq.0) then ! got a value
    1082             ! do a sanity check
    1083             if((val.lt.0.).or.(val.gt.90)) then
    1084               write(*,*)'Latitude should be between 0 and 90 deg. !!!'
    1085               ierr=1
    1086             else ! find corresponding jref (nearest latitude)
    1087               ! note: rlatu(:) contains decreasing values of latitude
    1088               !       starting from PI/2 to -PI/2
    1089               do j=1,jjp1
    1090                 if ((rlatu(j)*180./pi.ge.val).and.
    1091      &              (rlatu(j+1)*180./pi.le.val)) then
    1092                   ! find which grid point is nearest to val:
    1093                   if (abs(rlatu(j)*180./pi-val).le.
    1094      &                abs((rlatu(j+1)*180./pi-val))) then
    1095                    jref=j
    1096                   else
    1097                    jref=j+1
    1098                   endif
    1099                  
    1100                  write(*,*)'Will use nearest grid latitude which is:',
    1101      &                     rlatu(jref)*180./pi
    1102                 endif
    1103               enddo ! of do j=1,jjp1
    1104             endif ! of if((val.lt.0.).or.(val.gt.90))
    1105           endif !of if (ierr.eq.0)
    1106          enddo ! of do while
    1107 
    1108          ! Build layers() (as in soil_settings.F)
    1109          val2=sqrt(mlayer(0)*mlayer(1))
    1110          val3=mlayer(1)/mlayer(0)
    1111          do isoil=1,nsoilmx
    1112            layer(isoil)=val2*(val3**(isoil-1))
    1113          enddo
    1114 
    1115          write(*,*)'At which depth (in m.) does the ice layer begin?'
    1116          write(*,*)'(currently, the deepest soil layer extends down to:'
    1117      &              ,layer(nsoilmx),')'
    1118          ierr=1
    1119          do while (ierr.ne.0)
    1120           read(*,*,iostat=ierr) val2
    1121 !         write(*,*)'val2:',val2,'ierr=',ierr
    1122           if (ierr.eq.0) then ! got a value, but do a sanity check
    1123             if(val2.gt.layer(nsoilmx)) then
    1124               write(*,*)'Depth should be less than ',layer(nsoilmx)
    1125               ierr=1
    1126             endif
    1127             if(val2.lt.layer(1)) then
    1128               write(*,*)'Depth should be more than ',layer(1)
    1129               ierr=1
    1130             endif
    1131           endif
    1132          enddo ! of do while
    1133          
    1134          ! find the reference index iref the depth corresponds to
    1135 !        if (val2.lt.layer(1)) then
    1136 !         iref=1
    1137 !        else
    1138           do isoil=1,nsoilmx-1
    1139            if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
    1140      &       then
    1141              iref=isoil
    1142              exit
    1143            endif
    1144           enddo
    1145 !        endif
    1146          
    1147 !        write(*,*)'iref:',iref,'  jref:',jref
    1148 !        write(*,*)'layer',layer
    1149 !        write(*,*)'mlayer',mlayer
    1150          
    1151          ! thermal inertia of the ice:
    1152          ierr=1
    1153          do while (ierr.ne.0)
    1154           write(*,*)'What is the value of subterranean ice thermal inert
    1155      &ia? (e.g.: 2000)'
    1156           read(*,*,iostat=ierr)iceith
    1157          enddo ! of do while
    1158          
    1159          ! recast ithfi() onto ith()
    1160          call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    1161          
    1162          do j=1,jref
    1163 !           write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
    1164             do i=1,iip1 ! loop on longitudes
    1165              ! Build "equivalent" thermal inertia for the mixed layer
    1166              ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
    1167      &                     (((val2-layer(iref))/(ith(i,j,iref)**2))+
    1168      &                      ((layer(iref+1)-val2)/(iceith)**2)))
    1169              ! Set thermal inertia of lower layers
    1170              do isoil=iref+2,nsoilmx
    1171               ith(i,j,isoil)=iceith ! ice
    1172              enddo
    1173             enddo ! of do i=1,iip1
    1174          enddo ! of do j=1,jjp1
    1175          
    1176 
    1177          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1178 
    1179 !         do i=1,nsoilmx
    1180 !         write(*,*)'i:',i,'ithfi(1,i):',ithfi(1,i)
    1181 !        enddo
     1329
     1330
     1331
     1332c$$$!       subsoilice_n: Put deep ice layer in northern hemisphere soil
     1333c$$$!       ------------------------------------------------------------
     1334c$$$
     1335c$$$    else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
     1336c$$$
     1337c$$$         write(*,*)'From which latitude (in deg.), up to the north pole,
     1338c$$$     &should we put subterranean ice?'
     1339c$$$     ierr=1
     1340c$$$     do while (ierr.ne.0)
     1341c$$$      read(*,*,iostat=ierr) val
     1342c$$$      if (ierr.eq.0) then ! got a value
     1343c$$$        ! do a sanity check
     1344c$$$        if((val.lt.0.).or.(val.gt.90)) then
     1345c$$$          write(*,*)'Latitude should be between 0 and 90 deg. !!!'
     1346c$$$          ierr=1
     1347c$$$        else ! find corresponding jref (nearest latitude)
     1348c$$$          ! note: rlatu(:) contains decreasing values of latitude
     1349c$$$          !       starting from PI/2 to -PI/2
     1350c$$$          do j=1,jjp1
     1351c$$$            if ((rlatu(j)*180./pi.ge.val).and.
     1352c$$$     &              (rlatu(j+1)*180./pi.le.val)) then
     1353c$$$              ! find which grid point is nearest to val:
     1354c$$$              if (abs(rlatu(j)*180./pi-val).le.
     1355c$$$     &                abs((rlatu(j+1)*180./pi-val))) then
     1356c$$$               jref=j
     1357c$$$              else
     1358c$$$               jref=j+1
     1359c$$$              endif
     1360c$$$             
     1361c$$$             write(*,*)'Will use nearest grid latitude which is:',
     1362c$$$     &                     rlatu(jref)*180./pi
     1363c$$$            endif
     1364c$$$          enddo ! of do j=1,jjp1
     1365c$$$        endif ! of if((val.lt.0.).or.(val.gt.90))
     1366c$$$      endif !of if (ierr.eq.0)
     1367c$$$     enddo ! of do while
     1368c$$$
     1369c$$$         ! Build layers() (as in soil_settings.F)
     1370c$$$     val2=sqrt(mlayer(0)*mlayer(1))
     1371c$$$     val3=mlayer(1)/mlayer(0)
     1372c$$$     do isoil=1,nsoilmx
     1373c$$$       layer(isoil)=val2*(val3**(isoil-1))
     1374c$$$     enddo
     1375c$$$
     1376c$$$         write(*,*)'At which depth (in m.) does the ice layer begin?'
     1377c$$$         write(*,*)'(currently, the deepest soil layer extends down to:'
     1378c$$$     &              ,layer(nsoilmx),')'
     1379c$$$     ierr=1
     1380c$$$     do while (ierr.ne.0)
     1381c$$$      read(*,*,iostat=ierr) val2
     1382c$$$!     write(*,*)'val2:',val2,'ierr=',ierr
     1383c$$$      if (ierr.eq.0) then ! got a value, but do a sanity check
     1384c$$$        if(val2.gt.layer(nsoilmx)) then
     1385c$$$          write(*,*)'Depth should be less than ',layer(nsoilmx)
     1386c$$$          ierr=1
     1387c$$$        endif
     1388c$$$        if(val2.lt.layer(1)) then
     1389c$$$          write(*,*)'Depth should be more than ',layer(1)
     1390c$$$          ierr=1
     1391c$$$        endif
     1392c$$$      endif
     1393c$$$         enddo ! of do while
     1394c$$$     
     1395c$$$     ! find the reference index iref the depth corresponds to
     1396c$$$!    if (val2.lt.layer(1)) then
     1397c$$$!     iref=1
     1398c$$$!    else
     1399c$$$      do isoil=1,nsoilmx-1
     1400c$$$       if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
     1401c$$$     &       then
     1402c$$$         iref=isoil
     1403c$$$         exit
     1404c$$$       endif
     1405c$$$      enddo
     1406c$$$!    endif
     1407c$$$     
     1408c$$$!    write(*,*)'iref:',iref,'  jref:',jref
     1409c$$$!    write(*,*)'layer',layer
     1410c$$$!    write(*,*)'mlayer',mlayer
     1411c$$$         
     1412c$$$     ! thermal inertia of the ice:
     1413c$$$     ierr=1
     1414c$$$     do while (ierr.ne.0)
     1415c$$$          write(*,*)'What is the value of subterranean ice thermal inert
     1416c$$$     &ia? (e.g.: 2000)'
     1417c$$$          read(*,*,iostat=ierr)iceith
     1418c$$$     enddo ! of do while
     1419c$$$     
     1420c$$$     ! recast ithfi() onto ith()
     1421c$$$     call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     1422c$$$     
     1423c$$$     do j=1,jref
     1424c$$$!       write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
     1425c$$$        do i=1,iip1 ! loop on longitudes
     1426c$$$         ! Build "equivalent" thermal inertia for the mixed layer
     1427c$$$         ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
     1428c$$$     &                     (((val2-layer(iref))/(ith(i,j,iref)**2))+
     1429c$$$     &                      ((layer(iref+1)-val2)/(iceith)**2)))
     1430c$$$         ! Set thermal inertia of lower layers
     1431c$$$         do isoil=iref+2,nsoilmx
     1432c$$$          ith(i,j,isoil)=iceith ! ice
     1433c$$$         enddo
     1434c$$$        enddo ! of do i=1,iip1
     1435c$$$     enddo ! of do j=1,jjp1
     1436c$$$     
     1437c$$$
     1438c$$$     CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1439c$$$
     1440c$$$!         do i=1,nsoilmx
     1441c$$$!     write(*,*)'i:',i,'ithfi(1,i):',ithfi(1,i)
     1442c$$$!    enddo
    11821443
    11831444       
    1184 !       subsoilice_s: Put deep ice layer in southern hemisphere soil
    1185 !       ------------------------------------------------------------
    1186 
    1187         else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then
    1188 
    1189          write(*,*)'From which latitude (in deg.), down to the south pol
    1190      &e, should we put subterranean ice?'
    1191          ierr=1
    1192          do while (ierr.ne.0)
    1193           read(*,*,iostat=ierr) val
    1194           if (ierr.eq.0) then ! got a value
    1195             ! do a sanity check
    1196             if((val.gt.0.).or.(val.lt.-90)) then
    1197               write(*,*)'Latitude should be between 0 and -90 deg. !!!'
    1198               ierr=1
    1199             else ! find corresponding jref (nearest latitude)
    1200               ! note: rlatu(:) contains decreasing values of latitude
    1201               !       starting from PI/2 to -PI/2
    1202               do j=1,jjp1
    1203                 if ((rlatu(j)*180./pi.ge.val).and.
    1204      &              (rlatu(j+1)*180./pi.le.val)) then
    1205                   ! find which grid point is nearest to val:
    1206                   if (abs(rlatu(j)*180./pi-val).le.
    1207      &                abs((rlatu(j+1)*180./pi-val))) then
    1208                    jref=j
    1209                   else
    1210                    jref=j+1
    1211                   endif
    1212                  
    1213                  write(*,*)'Will use nearest grid latitude which is:',
    1214      &                     rlatu(jref)*180./pi
    1215                 endif
    1216               enddo ! of do j=1,jjp1
    1217             endif ! of if((val.lt.0.).or.(val.gt.90))
    1218           endif !of if (ierr.eq.0)
    1219          enddo ! of do while
    1220 
    1221          ! Build layers() (as in soil_settings.F)
    1222          val2=sqrt(mlayer(0)*mlayer(1))
    1223          val3=mlayer(1)/mlayer(0)
    1224          do isoil=1,nsoilmx
    1225            layer(isoil)=val2*(val3**(isoil-1))
    1226          enddo
    1227 
    1228          write(*,*)'At which depth (in m.) does the ice layer begin?'
    1229          write(*,*)'(currently, the deepest soil layer extends down to:'
    1230      &              ,layer(nsoilmx),')'
    1231          ierr=1
    1232          do while (ierr.ne.0)
    1233           read(*,*,iostat=ierr) val2
    1234 !         write(*,*)'val2:',val2,'ierr=',ierr
    1235           if (ierr.eq.0) then ! got a value, but do a sanity check
    1236             if(val2.gt.layer(nsoilmx)) then
    1237               write(*,*)'Depth should be less than ',layer(nsoilmx)
    1238               ierr=1
    1239             endif
    1240             if(val2.lt.layer(1)) then
    1241               write(*,*)'Depth should be more than ',layer(1)
    1242               ierr=1
    1243             endif
    1244           endif
    1245          enddo ! of do while
    1246          
    1247          ! find the reference index iref the depth corresponds to
    1248           do isoil=1,nsoilmx-1
    1249            if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
    1250      &       then
    1251              iref=isoil
    1252              exit
    1253            endif
    1254           enddo
    1255          
    1256 !        write(*,*)'iref:',iref,'  jref:',jref
    1257          
    1258          ! thermal inertia of the ice:
    1259          ierr=1
    1260          do while (ierr.ne.0)
    1261           write(*,*)'What is the value of subterranean ice thermal inert
    1262      &ia? (e.g.: 2000)'
    1263           read(*,*,iostat=ierr)iceith
    1264          enddo ! of do while
    1265          
    1266          ! recast ithfi() onto ith()
    1267          call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    1268          
    1269          do j=jref,jjp1
    1270 !           write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
    1271             do i=1,iip1 ! loop on longitudes
    1272              ! Build "equivalent" thermal inertia for the mixed layer
    1273              ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
    1274      &                     (((val2-layer(iref))/(ith(i,j,iref)**2))+
    1275      &                      ((layer(iref+1)-val2)/(iceith)**2)))
    1276              ! Set thermal inertia of lower layers
    1277              do isoil=iref+2,nsoilmx
    1278               ith(i,j,isoil)=iceith ! ice
    1279              enddo
    1280             enddo ! of do i=1,iip1
    1281          enddo ! of do j=jref,jjp1
    1282          
    1283 
    1284          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1285 
    1286 c       'mons_ice' : use MONS data to build subsurface ice table
    1287 c       --------------------------------------------------------
    1288         else if (modif(1:len_trim(modif)).eq.'mons_ice') then
    1289        
    1290        ! 1. Load MONS data
    1291         call load_MONS_data(MONS_Hdn,MONS_d21)
    1292        
    1293         ! 2. Get parameters from user
    1294         ierr=1
    1295         do while (ierr.ne.0)
    1296           write(*,*) "Coefficient to apply to MONS 'depth' in Northern",
    1297      &               " Hemisphere?"
    1298           write(*,*) " (should be somewhere between 3.2e-4 and 1.3e-3)"
    1299           read(*,*,iostat=ierr) MONS_coeffN
    1300         enddo
    1301         ierr=1
    1302         do while (ierr.ne.0)
    1303           write(*,*) "Coefficient to apply to MONS 'depth' in Southern",
    1304      &               " Hemisphere?"
    1305           write(*,*) " (should be somewhere between 3.2e-4 and 1.3e-3)"
    1306           read(*,*,iostat=ierr) MONS_coeffS
    1307         enddo
    1308         ierr=1
    1309         do while (ierr.ne.0)
    1310           write(*,*) "Value of subterranean ice thermal inertia?"
    1311           write(*,*) " (e.g.: 2000, or perhaps 2290)"
    1312           read(*,*,iostat=ierr) iceith
    1313         enddo
    1314        
    1315         ! 3. Build subterranean thermal inertia
    1316        
    1317         ! initialise subsurface inertia with reference surface values
    1318         do isoil=1,nsoilmx
    1319           ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx)
    1320         enddo
    1321         ! recast ithfi() onto ith()
    1322         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
    1323        
    1324         do i=1,iip1 ! loop on longitudes
    1325           do j=1,jjp1 ! loop on latitudes
    1326             ! set MONS_coeff
    1327             if (rlatu(j).ge.0) then ! northern hemisphere
    1328               ! N.B: rlatu(:) contains decreasing values of latitude
    1329               !       starting from PI/2 to -PI/2
    1330               MONS_coeff=MONS_coeffN
    1331             else ! southern hemisphere
    1332               MONS_coeff=MONS_coeffS
    1333             endif
    1334             ! check if we should put subterranean ice
    1335             if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14%
    1336               ! compute depth at which ice lies:
    1337               val=MONS_d21(i,j)*MONS_coeff
    1338               ! compute val2= the diurnal skin depth of surface inertia
    1339               ! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1
    1340               val2=ith(i,j,1)*1.e-6*sqrt(88775./3.14159)
    1341               if (val.lt.val2) then
    1342                 ! ice must be below the (surface inertia) diurnal skin depth
    1343                 val=val2
    1344               endif
    1345               if (val.lt.layer(nsoilmx)) then ! subterranean ice
    1346                 ! find the reference index iref that depth corresponds to
    1347                 iref=0
    1348                 do isoil=1,nsoilmx-1
    1349                  if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1)))
    1350      &             then
    1351                    iref=isoil
    1352                    exit
    1353                  endif
    1354                 enddo
    1355                 ! Build "equivalent" thermal inertia for the mixed layer
    1356                 ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
    1357      &                     (((val-layer(iref))/(ith(i,j,iref+1)**2))+
    1358      &                      ((layer(iref+1)-val)/(iceith)**2)))
    1359                 ! Set thermal inertia of lower layers
    1360                 do isoil=iref+2,nsoilmx
    1361                   ith(i,j,isoil)=iceith
    1362                 enddo
    1363               endif ! of if (val.lt.layer(nsoilmx))
    1364             endif ! of if (MONS_Hdn(i,j).lt.14.0)
    1365           enddo ! do j=1,jjp1
    1366         enddo ! do i=1,iip1
    1367        
    1368 ! Check:
    1369 !         do i=1,iip1
    1370 !           do j=1,jjp1
    1371 !             do isoil=1,nsoilmx
    1372 !               write(77,*) i,j,isoil,"  ",ith(i,j,isoil)
    1373 !             enddo
    1374 !           enddo
    1375 !        enddo
    1376 
    1377         ! recast ith() into ithfi()
    1378         CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
    1379        
    1380         else
    1381           write(*,*) '       Unknown (misspelled?) option!!!'
     1445c$$$!       subsoilice_s: Put deep ice layer in southern hemisphere soil
     1446c$$$!       ------------------------------------------------------------
     1447c$$$
     1448c$$$    else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then
     1449c$$$
     1450c$$$         write(*,*)'From which latitude (in deg.), down to the south pol
     1451c$$$     &e, should we put subterranean ice?'
     1452c$$$     ierr=1
     1453c$$$     do while (ierr.ne.0)
     1454c$$$      read(*,*,iostat=ierr) val
     1455c$$$      if (ierr.eq.0) then ! got a value
     1456c$$$        ! do a sanity check
     1457c$$$        if((val.gt.0.).or.(val.lt.-90)) then
     1458c$$$          write(*,*)'Latitude should be between 0 and -90 deg. !!!'
     1459c$$$          ierr=1
     1460c$$$        else ! find corresponding jref (nearest latitude)
     1461c$$$          ! note: rlatu(:) contains decreasing values of latitude
     1462c$$$          !       starting from PI/2 to -PI/2
     1463c$$$          do j=1,jjp1
     1464c$$$            if ((rlatu(j)*180./pi.ge.val).and.
     1465c$$$     &              (rlatu(j+1)*180./pi.le.val)) then
     1466c$$$              ! find which grid point is nearest to val:
     1467c$$$              if (abs(rlatu(j)*180./pi-val).le.
     1468c$$$     &                abs((rlatu(j+1)*180./pi-val))) then
     1469c$$$               jref=j
     1470c$$$              else
     1471c$$$               jref=j+1
     1472c$$$              endif
     1473c$$$             
     1474c$$$             write(*,*)'Will use nearest grid latitude which is:',
     1475c$$$     &                     rlatu(jref)*180./pi
     1476c$$$            endif
     1477c$$$          enddo ! of do j=1,jjp1
     1478c$$$        endif ! of if((val.lt.0.).or.(val.gt.90))
     1479c$$$      endif !of if (ierr.eq.0)
     1480c$$$     enddo ! of do while
     1481c$$$
     1482c$$$         ! Build layers() (as in soil_settings.F)
     1483c$$$     val2=sqrt(mlayer(0)*mlayer(1))
     1484c$$$     val3=mlayer(1)/mlayer(0)
     1485c$$$     do isoil=1,nsoilmx
     1486c$$$       layer(isoil)=val2*(val3**(isoil-1))
     1487c$$$     enddo
     1488c$$$
     1489c$$$         write(*,*)'At which depth (in m.) does the ice layer begin?'
     1490c$$$         write(*,*)'(currently, the deepest soil layer extends down to:'
     1491c$$$     &              ,layer(nsoilmx),')'
     1492c$$$     ierr=1
     1493c$$$     do while (ierr.ne.0)
     1494c$$$      read(*,*,iostat=ierr) val2
     1495c$$$!     write(*,*)'val2:',val2,'ierr=',ierr
     1496c$$$      if (ierr.eq.0) then ! got a value, but do a sanity check
     1497c$$$        if(val2.gt.layer(nsoilmx)) then
     1498c$$$          write(*,*)'Depth should be less than ',layer(nsoilmx)
     1499c$$$          ierr=1
     1500c$$$        endif
     1501c$$$        if(val2.lt.layer(1)) then
     1502c$$$          write(*,*)'Depth should be more than ',layer(1)
     1503c$$$          ierr=1
     1504c$$$        endif
     1505c$$$      endif
     1506c$$$         enddo ! of do while
     1507c$$$     
     1508c$$$     ! find the reference index iref the depth corresponds to
     1509c$$$      do isoil=1,nsoilmx-1
     1510c$$$       if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
     1511c$$$     &       then
     1512c$$$         iref=isoil
     1513c$$$         exit
     1514c$$$       endif
     1515c$$$      enddo
     1516c$$$     
     1517c$$$!    write(*,*)'iref:',iref,'  jref:',jref
     1518c$$$         
     1519c$$$     ! thermal inertia of the ice:
     1520c$$$     ierr=1
     1521c$$$     do while (ierr.ne.0)
     1522c$$$          write(*,*)'What is the value of subterranean ice thermal inert
     1523c$$$     &ia? (e.g.: 2000)'
     1524c$$$          read(*,*,iostat=ierr)iceith
     1525c$$$     enddo ! of do while
     1526c$$$     
     1527c$$$     ! recast ithfi() onto ith()
     1528c$$$     call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     1529c$$$     
     1530c$$$     do j=jref,jjp1
     1531c$$$!       write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi
     1532c$$$        do i=1,iip1 ! loop on longitudes
     1533c$$$         ! Build "equivalent" thermal inertia for the mixed layer
     1534c$$$         ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
     1535c$$$     &                     (((val2-layer(iref))/(ith(i,j,iref)**2))+
     1536c$$$     &                      ((layer(iref+1)-val2)/(iceith)**2)))
     1537c$$$         ! Set thermal inertia of lower layers
     1538c$$$         do isoil=iref+2,nsoilmx
     1539c$$$          ith(i,j,isoil)=iceith ! ice
     1540c$$$         enddo
     1541c$$$        enddo ! of do i=1,iip1
     1542c$$$     enddo ! of do j=jref,jjp1
     1543c$$$     
     1544c$$$
     1545c$$$     CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1546
     1547
     1548c$$$c       'mons_ice' : use MONS data to build subsurface ice table
     1549c$$$c       --------------------------------------------------------
     1550c$$$        else if (modif(1:len_trim(modif)).eq.'mons_ice') then
     1551c$$$       
     1552c$$$       ! 1. Load MONS data
     1553c$$$        call load_MONS_data(MONS_Hdn,MONS_d21)
     1554c$$$       
     1555c$$$        ! 2. Get parameters from user
     1556c$$$        ierr=1
     1557c$$$    do while (ierr.ne.0)
     1558c$$$          write(*,*) "Coefficient to apply to MONS 'depth' in Northern",
     1559c$$$     &               " Hemisphere?"
     1560c$$$          write(*,*) " (should be somewhere between 3.2e-4 and 1.3e-3)"
     1561c$$$          read(*,*,iostat=ierr) MONS_coeffN
     1562c$$$        enddo
     1563c$$$        ierr=1
     1564c$$$    do while (ierr.ne.0)
     1565c$$$          write(*,*) "Coefficient to apply to MONS 'depth' in Southern",
     1566c$$$     &               " Hemisphere?"
     1567c$$$          write(*,*) " (should be somewhere between 3.2e-4 and 1.3e-3)"
     1568c$$$          read(*,*,iostat=ierr) MONS_coeffS
     1569c$$$        enddo
     1570c$$$        ierr=1
     1571c$$$        do while (ierr.ne.0)
     1572c$$$          write(*,*) "Value of subterranean ice thermal inertia?"
     1573c$$$          write(*,*) " (e.g.: 2000, or perhaps 2290)"
     1574c$$$          read(*,*,iostat=ierr) iceith
     1575c$$$        enddo
     1576c$$$       
     1577c$$$        ! 3. Build subterranean thermal inertia
     1578c$$$       
     1579c$$$        ! initialise subsurface inertia with reference surface values
     1580c$$$        do isoil=1,nsoilmx
     1581c$$$          ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx)
     1582c$$$        enddo
     1583c$$$        ! recast ithfi() onto ith()
     1584c$$$    call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
     1585c$$$       
     1586c$$$        do i=1,iip1 ! loop on longitudes
     1587c$$$          do j=1,jjp1 ! loop on latitudes
     1588c$$$            ! set MONS_coeff
     1589c$$$            if (rlatu(j).ge.0) then ! northern hemisphere
     1590c$$$              ! N.B: rlatu(:) contains decreasing values of latitude
     1591c$$$          !       starting from PI/2 to -PI/2
     1592c$$$              MONS_coeff=MONS_coeffN
     1593c$$$            else ! southern hemisphere
     1594c$$$              MONS_coeff=MONS_coeffS
     1595c$$$            endif
     1596c$$$            ! check if we should put subterranean ice
     1597c$$$            if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14%
     1598c$$$              ! compute depth at which ice lies:
     1599c$$$              val=MONS_d21(i,j)*MONS_coeff
     1600c$$$              ! compute val2= the diurnal skin depth of surface inertia
     1601c$$$              ! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1
     1602c$$$              val2=ith(i,j,1)*1.e-6*sqrt(88775./3.14159)
     1603c$$$              if (val.lt.val2) then
     1604c$$$                ! ice must be below the (surface inertia) diurnal skin depth
     1605c$$$                val=val2
     1606c$$$              endif
     1607c$$$              if (val.lt.layer(nsoilmx)) then ! subterranean ice
     1608c$$$                ! find the reference index iref that depth corresponds to
     1609c$$$                iref=0
     1610c$$$                do isoil=1,nsoilmx-1
     1611c$$$                 if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1)))
     1612c$$$     &             then
     1613c$$$               iref=isoil
     1614c$$$               exit
     1615c$$$             endif
     1616c$$$                enddo
     1617c$$$                ! Build "equivalent" thermal inertia for the mixed layer
     1618c$$$                ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
     1619c$$$     &                     (((val-layer(iref))/(ith(i,j,iref+1)**2))+
     1620c$$$     &                      ((layer(iref+1)-val)/(iceith)**2)))
     1621c$$$            ! Set thermal inertia of lower layers
     1622c$$$                do isoil=iref+2,nsoilmx
     1623c$$$                  ith(i,j,isoil)=iceith
     1624c$$$                enddo
     1625c$$$              endif ! of if (val.lt.layer(nsoilmx))
     1626c$$$            endif ! of if (MONS_Hdn(i,j).lt.14.0)
     1627c$$$          enddo ! do j=1,jjp1
     1628c$$$        enddo ! do i=1,iip1
     1629c$$$       
     1630c$$$! Check:
     1631c$$$!         do i=1,iip1
     1632c$$$!           do j=1,jjp1
     1633c$$$!             do isoil=1,nsoilmx
     1634c$$$!               write(77,*) i,j,isoil,"  ",ith(i,j,isoil)
     1635c$$$!             enddo
     1636c$$$!           enddo
     1637c$$$!    enddo
     1638c$$$
     1639c$$$        ! recast ith() into ithfi()
     1640c$$$        CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
     1641c$$$       
     1642c$$$    else
     1643c$$$          write(*,*) '       Unknown (misspelled?) option!!!'
    13821644        end if ! of if (modif(1:len_trim(modif)) .eq. '...') elseif ...
    13831645       
     1646
     1647
     1648
     1649
     1650
     1651
     1652
     1653
     1654
     1655
    13841656       enddo ! of do ! infinite loop on liste of changes
    13851657
    … …  
    13871659
    13881660 
     1661c=======================================================================
     1662c   Initialisation for cloud fraction and oceanic ice (added by BC 2010)
     1663c=======================================================================
     1664      DO ig=1, ngridmx
     1665         totalfrac(ig)=0.5
     1666         DO l=1,llm
     1667            cloudfrac(ig,l)=0.5
     1668         ENDDO
     1669      ENDDO
     1670     
     1671c========================================================
     1672
     1673!      DO ig=1,ngridmx
     1674!         IF(tsurf(ig) .LT. 273.16-1.8) THEN
     1675!            hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
     1676!            hice(ig)=min(hice(ig),1.0)
     1677!         ENDIF
     1678!      ENDDO
     1679
     1680
    13891681c=======================================================================
    13901682c   Correct pressure on the new grid (menu 0)
    … …  
    13931685      if ((choix_1.eq.0).and.(.not.flagps0)) then
    13941686        r = 1000.*8.31/mugaz
     1687
     1688        phi0=0.0
     1689        do j=1,jjp1
     1690          do i=1,iip1
     1691             phi0 = phi0+phis(i,j)*aire(i,j)
     1692          end do
     1693        end do
     1694        phi0=phi0/airetot
    13951695
    13961696        do j=1,jjp1
    … …  
    13981698             ps(i,j) = ps(i,j) *
    13991699!     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
    1400      .            exp((-phis(i,j)) /
     1700     .            exp((phi0-phis(i,j)) /
    14011701     .                                  (t(i,j,1) * r))
    14021702          end do
    14031703        end do
    14041704 
    1405 c periodicity of surface ps in longitude
     1705!     we must renormalise pressures AGAIN, because exp(-phi) is nonlinear
     1706        ptot=0.0
    14061707        do j=1,jjp1
    1407           ps(1,j) = ps(iip1,j)
     1708           do i=1,iip1
     1709              ptot=ptot+ps(i,j)*aire(i,j)
     1710           enddo
     1711        enddo
     1712        do j=1,jjp1
     1713           do i=1,iip1
     1714              ps(i,j)=ps(i,j)*ptotn/ptot
     1715           enddo
     1716        enddo
     1717       
     1718!     periodicity of surface ps in longitude
     1719        do j=1,jjp1
     1720           ps(1,j) = ps(iip1,j)
    14081721        end do
     1722       
    14091723      end if
    14101724
     1725
    14111726c=======================================================================
    14121727c=======================================================================
    … …  
    14151730c    Initialisation de la physique / ecriture de newstartfi :
    14161731c=======================================================================
    1417 
    14181732
    14191733      CALL inifilr
    … …  
    14471761      endif
    14481762
     1763
    14491764C Calcul intermediaire
    14501765c
    … …  
    14901805C
    14911806
     1807!      do ig=1,ngridmx
     1808!         print*,'surface ice in newstart=',qsurf(ig,igcm_h2o_ice)
     1809!         print*,'surface liq in newstart=',qsurf(ig,igcm_h2o_vap)
     1810!      enddo
     1811!      stop
     1812
     1813
    14921814      call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
    14931815     .              dtphys,float(day_ini),
    14941816     .              time,tsurf,tsoil,emis,q2,qsurf,
    1495      .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
     1817     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
     1818     .              cloudfrac,totalfrac,hice)
    14961819
    14971820c=======================================================================

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

    r135 r253  
    6868      INTEGER*4 day_ini_fi
    6969
     70!     added by FF for cloud fraction setup
     71      REAL hice(ngridmx)
     72      REAL cloudfrac(ngridmx,nlayermx),totalcloudfrac(ngridmx)
     73
     74
    7075c Variable naturelle / grille scalaire
    7176c ------------------------------------
    … …  
    7883      REAL emisS(ip1jmp1)
    7984
     85!     added by FF for cloud fraction setup
     86      REAL hiceS(ip1jmp1)
     87      REAL cloudfracS(ip1jmp1,llm),totalcloudfracS(ip1jmp1)
     88
    8089c Variables intermediaires : vent naturel, mais pas coord scalaire
    8190c----------------------------------------------------------------
    … …  
    127136
    128137      CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,nqmx,day_ini_fi,timefi,
    129      .      tsurf,tsoil,emis,q2,qsurf)
     138     .      tsurf,tsoil,emis,q2,qsurf,
     139!       change FF 05/2011
     140     .       cloudfrac,totalcloudfrac,hice)
     141
     142
    130143
    131144       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
    … …  
    233246      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
    234247      call gr_fi_dyn(nqmx,ngridmx,iip1,jjp1,qsurf,qsurfS)
     248      call gr_fi_dyn(llm,ngridmx,iip1,jjp1,cloudfrac,cloudfracS)
     249      call gr_fi_dyn(1,ngridmx,iip1,jjp1,hice,hiceS)
     250      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totalcloudfrac,totalcloudfracS)
    235251
    236252c=======================================================================
    … …  
    376392! Note: no need to write volcapa, it is stored in "controle" table
    377393
     394c-----------------------------------------------------------------------
     395c Ecriture du champs  cloudfrac,hice,totalcloudfrac
     396c-----------------------------------------------------------------------
     397      call write_archive(nid,ntime,'hice',
     398     &         'Height of oceanic ice','m',2,hiceS)
     399      call write_archive(nid,ntime,'totalcloudfrac',
     400     &        'Total cloud Fraction','',2,totalcloudfracS)
     401      call write_archive(nid,ntime,'cloudfrac'
     402     &        ,'Cloud fraction','',3,cloudfracS)
     403c-----------------------------------------------------------------------
     404c Fin
     405c-----------------------------------------------------------------------
    378406      ierr=NF_CLOSE(nid)
    379 c-----------------------------------------------------------------------
    380 c Fin
    381 c-----------------------------------------------------------------------
    382407
    383408      end

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

    r135 r253  
    7272c   Test pour savoir si on advecte a ce pas de temps
    7373      IF ( iadvtr.EQ.iapp_tracvl ) THEN
     74!      print*,'WARNING: ALL TRACER ADVECTION HALTED IN TRACVL.F'
     75!         if(2.eq.1)then
     76
    7477
    7578cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
Note: See TracChangeset for help on using the changeset viewer.