Changeset 66 for trunk/libf/dyn3d


Ignore:
Timestamp:
Feb 16, 2011, 4:57:45 PM (14 years ago)
Author:
emillour
Message:

EM: Mise a niveau par rapport a la version terrestre (LMDZ5V2.0-dev, rev 1487)

  • Mise a jour des scripts (terrestres) 'makegcm' et 'create_make_gcm'
  • Ajout du script 'makelmdz' (version "amelioree", en Bash, de makegcm)
  • Mise a jour des routines dans phylmd (sauf regr_lat_time_climoz_m.F)
  • disvert (dans dyn3d et dyn3dpar): passage au Fortran 90
  • parallel.F90 (dyn3dpar): correction bug
  • etat0_netcdf.F90 (dyn3d et dyn3dpar) : mise a jour mineure
  • ce0l.F90 (dyn3dpar) : correction bug
  • abort_gcm.F (dyn3dpar) : correction bug
  • ugeostr.F90 (dyn3d et dyn3dpar) : passage au Fortran 90
  • fluxstokenc_p.F (dyn3dpar) : correction bug
  • iniacademic.F90 (dyn3d et dyn3dpar) : passage au Fortran 90
  • friction_p.F (dyn3dpar) : correction bug
  • infotrac.F90 (dyn3d et dyn3dpar) : correction bug mineur sur lecture traceurs
  • caladvtrac.F (dyn3d) : modifications cosmetiques
Location:
trunk/libf/dyn3d
Files:
3 edited
3 moved

Legend:

Unmodified
Added
Removed

  • trunk/libf/dyn3d/caladvtrac.F

    r7 r66  
    88     *                   flxw, pk)
    99c
    10       USE infotrac
     10      USE infotrac, ONLY : nqtot
    1111      USE control_mod, ONLY : iapp_tracvl,planet_type
    1212 

  • trunk/libf/dyn3d/disvert.F90

    r64 r66  
    1 !
    2 ! $Id: disvert.F 1403 2010-07-01 09:02:53Z fairhead $
    3 !
    4       SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
     1! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
    52
    6 c    Auteur :  P. Le Van .
    7 c
    8       IMPLICIT NONE
     3SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    94
    10 #include "dimensions.h"
    11 #include "paramet.h"
    12 #include "iniprint.h"
    13 #include "logic.h"
    14 c
    15 c=======================================================================
    16 c
    17 c
    18 c    s = sigma ** kappa   :  coordonnee  verticale
    19 c    dsig(l)            : epaisseur de la couche l ds la coord.  s
    20 c    sig(l)             : sigma a l'interface des couches l et l-1
    21 c    ds(l)              : distance entre les couches l et l-1 en coord.s
    22 c
    23 c=======================================================================
    24 c
    25       REAL pa,preff
    26       REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
    27       REAL presnivs(llm)
    28 c
    29 c   declarations:
    30 c   -------------
    31 c
    32       REAL sig(llm+1),dsig(llm)
    33        real zzz(1:llm+1)
    34        real dzz(1:llm)
    35       real zk,zkm1,dzk1,dzk2,k0,k1
    36 c
    37       INTEGER l
    38       REAL snorm,dsigmin
    39       REAL alpha,beta,gama,delta,deltaz,h
    40       INTEGER np,ierr
    41       REAL pi,x
     5  ! Auteur : P. Le Van
    426
    43       REAL SSUM
    44 c
    45 c-----------------------------------------------------------------------
    46 c
    47       pi=2.*ASIN(1.)
     7  IMPLICIT NONE
    488
    49       OPEN(99,file='sigma.def',status='old',form='formatted',
    50      s   iostat=ierr)
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "iniprint.h"
     12  include "logic.h"
    5113
    52 c-----------------------------------------------------------------------
    53 c   cas 1 on lit les options dans sigma.def:
    54 c   ----------------------------------------
     14  ! s = sigma ** kappa : coordonnee verticale
     15  ! dsig(l) : epaisseur de la couche l ds la coord. s
     16  ! sig(l) : sigma a l'interface des couches l et l-1
     17  ! ds(l) : distance entre les couches l et l-1 en coord.s
    5518
    56       IF (ierr.eq.0) THEN
     19  REAL pa, preff
     20  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
     21  REAL presnivs(llm)
    5722
    58       READ(99,*) h           ! hauteur d'echelle 8.
    59       READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
    60       READ(99,*) beta        ! facteur d'acroissement en haut 1.3
    61       READ(99,*) k0          ! nombre de couches dans la transition surf
    62       READ(99,*) k1          ! nombre de couches dans la transition haute
    63       CLOSE(99)
    64       alpha=deltaz/(llm*h)
    65       write(lunout,*)'h,alpha,k0,k1,beta'
     23  REAL sig(llm+1), dsig(llm)
     24  real zk, zkm1, dzk1, dzk2, k0, k1
    6625
    67 c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
     26  INTEGER l
     27  REAL dsigmin
     28  REAL alpha, beta, deltaz, h
     29  INTEGER iostat
     30  REAL pi, x
    6831
    69       alpha=deltaz/tanh(1./k0)*2.
    70       zkm1=0.
    71       sig(1)=1.
    72       do l=1,llm
    73         sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
    74      + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
     32  !-----------------------------------------------------------------------
     33
     34  pi = 2 * ASIN(1.)
     35
     36  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     37
     38  IF (iostat == 0) THEN
     39     ! cas 1 on lit les options dans sigma.def:
     40     READ(99, *) h ! hauteur d'echelle 8.
     41     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
     42     READ(99, *) beta ! facteur d'acroissement en haut 1.3
     43     READ(99, *) k0 ! nombre de couches dans la transition surf
     44     READ(99, *) k1 ! nombre de couches dans la transition haute
     45     CLOSE(99)
     46     alpha=deltaz/(llm*h)
     47     write(lunout, *)'disvert: h, alpha, k0, k1, beta'
     48
     49     alpha=deltaz/tanh(1./k0)*2.
     50     zkm1=0.
     51     sig(1)=1.
     52     do l=1, llm
     53        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
     54             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
    7555        zk=-h*log(sig(l+1))
    7656
    7757        dzk1=alpha*tanh(l/k0)
    7858        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
    79         write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
     59        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
    8060        zkm1=zk
    81       enddo
     61     enddo
    8262
    83       sig(llm+1)=0.
     63     sig(llm+1)=0.
    8464
    85 c
    86        DO 2  l = 1, llm
    87        dsig(l) = sig(l)-sig(l+1)
    88    2   CONTINUE
    89 c
     65     DO l = 1, llm
     66        dsig(l) = sig(l)-sig(l+1)
     67     end DO
     68  ELSE
     69     if (ok_strato) then
     70        if (llm==39) then
     71           dsigmin=0.3
     72        else if (llm==50) then
     73           dsigmin=1.
     74        else
     75           WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster'
     76           dsigmin=1.
     77        endif
     78        WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin
     79     endif
    9080
    91       ELSE
    92 c-----------------------------------------------------------------------
    93 c   cas 2 ancienne discretisation (LMD5...):
    94 c   ----------------------------------------
     81     DO l = 1, llm
     82        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
    9583
    96       WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
     84        IF (ok_strato) THEN
     85           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     86                *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     87        ELSE
     88           dsig(l) = 1.0 + 7.0 * SIN(x)**2
     89        ENDIF
     90     ENDDO
     91     dsig = dsig / sum(dsig)
     92     sig(llm+1) = 0.
     93     DO l = llm, 1, -1
     94        sig(l) = sig(l+1) + dsig(l)
     95     ENDDO
     96  ENDIF
    9797
    98       if (ok_strato) then
    99          if (llm==39) then
    100             dsigmin=0.3
    101          else if (llm==50) then
    102             dsigmin=1.
    103          else
    104             WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
    105             dsigmin=1.
    106          endif
    107          WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
    108       endif
     98  DO l=1, llm
     99     nivsigs(l) = REAL(l)
     100  ENDDO
    109101
    110       h=7.
    111       snorm  = 0.
    112       DO l = 1, llm
    113          x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
     102  DO l=1, llmp1
     103     nivsig(l)= REAL(l)
     104  ENDDO
    114105
    115          IF (ok_strato) THEN
    116            dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
    117      &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
    118          ELSE
    119            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    120          ENDIF
     106  ! .... Calculs de ap(l) et de bp(l) ....
     107  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    121108
    122          snorm = snorm + dsig(l)
    123       ENDDO
    124       snorm = 1./snorm
    125       DO l = 1, llm
    126          dsig(l) = dsig(l)*snorm
    127       ENDDO
    128       sig(llm+1) = 0.
    129       DO l = llm, 1, -1
    130          sig(l) = sig(l+1) + dsig(l)
    131       ENDDO
     109  bp(llmp1) = 0.
    132110
    133       ENDIF
     111  DO l = 1, llm
     112     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
     113     ap(l) = pa * ( sig(l) - bp(l) )
     114  ENDDO
    134115
     116  bp(1)=1.
     117  ap(1)=0.
    135118
    136       DO l=1,llm
    137         nivsigs(l) = REAL(l)
    138       ENDDO
     119  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    139120
    140       DO l=1,llmp1
    141         nivsig(l)= REAL(l)
    142       ENDDO
     121  write(lunout, *)'disvert: BP '
     122  write(lunout, *) bp
     123  write(lunout, *)'disvert: AP '
     124  write(lunout, *) ap
    143125
    144 c
    145 c    ....  Calculs  de ap(l) et de bp(l)  ....
    146 c    .........................................
    147 c
    148 c
    149 c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
    150 c
     126  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
     127  write(lunout, *)'couches calcules pour une pression de surface =', preff
     128  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
     129  write(lunout, *)'8km'
     130  DO l = 1, llm
     131     dpres(l) = bp(l) - bp(l+1)
     132     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
     133     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
     134          log(preff/presnivs(l))*8. &
     135          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
     136          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
     137  ENDDO
    151138
    152       bp(llmp1) =   0.
     139  write(lunout, *) 'disvert: PRESNIVS '
     140  write(lunout, *) presnivs
    153141
    154       DO l = 1, llm
    155 cc
    156 ccc    ap(l) = 0.
    157 ccc    bp(l) = sig(l)
    158 
    159       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    160       ap(l) = pa * ( sig(l) - bp(l) )
    161 c
    162       ENDDO
    163 
    164       bp(1)=1.
    165       ap(1)=0.
    166 
    167       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    168 
    169       write(lunout,*)' BP '
    170       write(lunout,*)  bp
    171       write(lunout,*)' AP '
    172       write(lunout,*)  ap
    173 
    174       write(lunout,*)
    175      .'Niveaux de pressions approximatifs aux centres des'
    176       write(lunout,*)'couches calcules pour une pression de surface =',
    177      .                 preff
    178       write(lunout,*)
    179      .     'et altitudes equivalentes pour une hauteur d echelle de'
    180       write(lunout,*)'8km'
    181       DO l = 1, llm
    182        dpres(l) = bp(l) - bp(l+1)
    183        presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
    184        write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
    185      .        log(preff/presnivs(l))*8.
    186      .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
    187      .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
    188       ENDDO
    189 
    190       write(lunout,*)' PRESNIVS '
    191       write(lunout,*)presnivs
    192 
    193       RETURN
    194       END
     142END SUBROUTINE disvert

  • trunk/libf/dyn3d/etat0_netcdf.F90

    r1 r66  
    11!
    2 ! $Id: etat0_netcdf.F90 1425 2010-09-02 13:45:23Z lguez $
     2! $Id: etat0_netcdf.F90 1486 2011-02-11 12:07:39Z fairhead $
    33!
    44!-------------------------------------------------------------------------------
    … …  
    9898  REAL    :: dummy
    9999  LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
    100   LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod
     100  LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats
    101101  INTEGER :: iflag_radia, flag_aerosol
    102102  REAL    :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
    … …  
    130130!--- CONSTRUCT A GRID
    131131  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
     132                   callstats,                                           &
    132133                   solarlong0,seuil_inversion,                          &
    133134                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &

  • trunk/libf/dyn3d/infotrac.F90

    r37 r66  
    177177          ! Continue to read tracer.def
    178178          DO iq=1,nqtrue
    179              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     179             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    180180          END DO
    181181          CLOSE(90) 
    … …  
    234234          ! Continue to read tracer.def
    235235          DO iq=1,nqtrue
    236              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     236             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    237237          END DO
    238238          CLOSE(90) 
    … …  
    316316       tname(new_iq)= tnom_0(iq)
    317317       IF (iadv(new_iq)==0) THEN
    318           ttext(new_iq)=str1(1:len_trim(str1))
     318          ttext(new_iq)=trim(str1)
    319319       ELSE
    320           ttext(new_iq)=str1(1:len_trim(str1))//descrq(iadv(new_iq))
     320          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
    321321       END IF
    322322
    … …  
    381381    DEALLOCATE(tracnam)
    382382
    383 999 FORMAT (i2,1x,i2,1x,a15)
    384 
    385383  END SUBROUTINE infotrac_init
    386384

  • trunk/libf/dyn3d/iniacademic.F90

    r64 r66  
    11!
    2 ! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z emillour $
     2! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $
    33!
    4 c
    5 c
    6       SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    7 
    8       USE filtreg_mod
    9       USE infotrac, ONLY : nqtot
    10       USE control_mod, ONLY: day_step,planet_type
     4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     5
     6  USE filtreg_mod
     7  USE infotrac, ONLY : nqtot
     8  USE control_mod, ONLY: day_step,planet_type
    119#ifdef CPP_IOIPSL
    12       USE IOIPSL
     10  USE IOIPSL
    1311#else
    14 ! if not using IOIPSL, we still need to use (a local version of) getin
    15       USE ioipsl_getincom
     12  ! if not using IOIPSL, we still need to use (a local version of) getin
     13  USE ioipsl_getincom
    1614#endif
    17       USE Write_Field
    18  
    19 
    20 c%W%    %G%
    21 c=======================================================================
    22 c
    23 c   Author:    Frederic Hourdin      original: 15/01/93
    24 c   -------
    25 c
    26 c   Subject:
    27 c   ------
    28 c
    29 c   Method:
    30 c   --------
    31 c
    32 c   Interface:
    33 c   ----------
    34 c
    35 c      Input:
    36 c      ------
    37 c
    38 c      Output:
    39 c      -------
    40 c
    41 c=======================================================================
    42       IMPLICIT NONE
    43 c-----------------------------------------------------------------------
    44 c   Declararations:
    45 c   ---------------
    46 
    47 #include "dimensions.h"
    48 #include "paramet.h"
    49 #include "comvert.h"
    50 #include "comconst.h"
    51 #include "comgeom.h"
    52 #include "academic.h"
    53 #include "ener.h"
    54 #include "temps.h"
    55 #include "iniprint.h"
    56 #include "logic.h"
    57 
    58 c   Arguments:
    59 c   ----------
    60 
    61       real time_0
    62 
    63 c   variables dynamiques
    64       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    65       REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    66       REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
    67       REAL ps(ip1jmp1)                       ! pression  au sol
    68       REAL masse(ip1jmp1,llm)                ! masse d'air
    69       REAL phis(ip1jmp1)                     ! geopotentiel au sol
    70 
    71 c   Local:
    72 c   ------
    73 
    74       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    75       REAL pks(ip1jmp1)                      ! exner au  sol
    76       REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    77       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    78       REAL phi(ip1jmp1,llm)                  ! geopotentiel
    79       REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
    80       real tetajl(jjp1,llm)
    81       INTEGER i,j,l,lsup,ij
    82 
    83       REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
    84       REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
    85       LOGICAL ok_geost             ! Initialisation vent geost. ou nul
    86       LOGICAL ok_pv                ! Polar Vortex
    87       REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
    88      
    89       real zz,ran1
    90       integer idum
    91 
    92       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
    93 
    94 c-----------------------------------------------------------------------
    95 ! 1. Initializations for Earth-like case
    96 ! --------------------------------------
    97 c
    98 
    99         print *, 'This is iniacademic'
    100 
    101         ! initialize planet radius, rotation rate,...
    102         call conf_planete
    103 
    104         time_0=0.
    105         day_ref=1
    106         annee_ref=0
    107 
    108         im         = iim
    109         jm         = jjm
    110         day_ini    = 1
    111         dtvr    = daysec/REAL(day_step)
    112         zdtvr=dtvr
    113         etot0      = 0.
    114         ptot0      = 0.
    115         ztot0      = 0.
    116         stot0      = 0.
    117         ang0       = 0.
    118 
    119         if (llm.eq.1) then
    120           ! specific initializations for the shallow water case
    121           kappa=1
     15  USE Write_Field
     16
     17  !   Author:    Frederic Hourdin      original: 15/01/93
     18  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
     19  ! of the American Meteorological Society, 75, 1825.
     20
     21  IMPLICIT NONE
     22
     23  !   Declararations:
     24  !   ---------------
     25
     26  include "dimensions.h"
     27  include "paramet.h"
     28  include "comvert.h"
     29  include "comconst.h"
     30  include "comgeom.h"
     31  include "academic.h"
     32  include "ener.h"
     33  include "temps.h"
     34  include "iniprint.h"
     35  include "logic.h"
     36
     37  !   Arguments:
     38  !   ----------
     39
     40  real time_0
     41
     42  !   variables dynamiques
     43  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
     44  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
     45  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
     46  REAL ps(ip1jmp1)                       ! pression  au sol
     47  REAL masse(ip1jmp1,llm)                ! masse d'air
     48  REAL phis(ip1jmp1)                     ! geopotentiel au sol
     49
     50  !   Local:
     51  !   ------
     52
     53  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     54  REAL pks(ip1jmp1)                      ! exner au  sol
     55  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
     56  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
     57  REAL phi(ip1jmp1,llm)                  ! geopotentiel
     58  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
     59  real tetastrat ! potential temperature in the stratosphere, in K
     60  real tetajl(jjp1,llm)
     61  INTEGER i,j,l,lsup,ij
     62
     63  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
     64  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     65  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
     66  LOGICAL ok_pv                ! Polar Vortex
     67  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
     68
     69  real zz,ran1
     70  integer idum
     71
     72  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     73
     74  !-----------------------------------------------------------------------
     75  ! 1. Initializations for Earth-like case
     76  ! --------------------------------------
     77  !
     78 
     79  write(lunout,*) 'Iniacademic'
     80 
     81  ! initialize planet radius, rotation rate,...
     82  call conf_planete
     83
     84  time_0=0.
     85  day_ref=1
     86  annee_ref=0
     87
     88  im         = iim
     89  jm         = jjm
     90  day_ini    = 1
     91  dtvr    = daysec/REAL(day_step)
     92  zdtvr=dtvr
     93  etot0      = 0.
     94  ptot0      = 0.
     95  ztot0      = 0.
     96  stot0      = 0.
     97  ang0       = 0.
     98
     99  if (llm == 1) then
     100     ! specific initializations for the shallow water case
     101     kappa=1
     102  endif
     103
     104  CALL iniconst
     105  CALL inigeom
     106  CALL inifilr
     107
     108  if (llm == 1) then
     109     ! initialize fields for the shallow water case, if required
     110     if (.not.read_start) then
     111        phis(:)=0.
     112        q(:,:,:)=0
     113        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
     114     endif
     115  endif
     116
     117  academic_case: if (iflag_phys == 2) then
     118     ! initializations
     119
     120     ! 1. local parameters
     121     ! by convention, winter is in the southern hemisphere
     122     ! Geostrophic wind or no wind?
     123     ok_geost=.TRUE.
     124     CALL getin('ok_geost',ok_geost)
     125     ! Constants for Newtonian relaxation and friction
     126     k_f=1.                !friction
     127     CALL getin('k_j',k_f)
     128     k_f=1./(daysec*k_f)
     129     k_c_s=4.  !cooling surface
     130     CALL getin('k_c_s',k_c_s)
     131     k_c_s=1./(daysec*k_c_s)
     132     k_c_a=40. !cooling free atm
     133     CALL getin('k_c_a',k_c_a)
     134     k_c_a=1./(daysec*k_c_a)
     135     ! Constants for Teta equilibrium profile
     136     teta0=315.     ! mean Teta (S.H. 315K)
     137     CALL getin('teta0',teta0)
     138     write(lunout,*) 'Iniacademic - teta0 ',teta0
     139     write(lunout,*) 'Iniacademic - rad ',rad
     140     ttp=200.       ! Tropopause temperature (S.H. 200K)
     141     CALL getin('ttp',ttp)
     142     eps=0.         ! Deviation to N-S symmetry(~0-20K)
     143     CALL getin('eps',eps)
     144     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
     145     CALL getin('delt_y',delt_y)
     146     delt_z=10.     ! Vertical Gradient (S.H. 10K)
     147     CALL getin('delt_z',delt_z)
     148     ! Polar vortex
     149     ok_pv=.false.
     150     CALL getin('ok_pv',ok_pv)
     151     phi_pv=-50.            ! Latitude of edge of vortex
     152     CALL getin('phi_pv',phi_pv)
     153     phi_pv=phi_pv*pi/180.
     154     dphi_pv=5.             ! Width of the edge
     155     CALL getin('dphi_pv',dphi_pv)
     156     dphi_pv=dphi_pv*pi/180.
     157     gam_pv=4.              ! -dT/dz vortex (in K/km)
     158     CALL getin('gam_pv',gam_pv)
     159
     160     ! 2. Initialize fields towards which to relax
     161     ! Friction
     162     knewt_g=k_c_a
     163     DO l=1,llm
     164        zsig=presnivs(l)/preff
     165        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
     166        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
     167     ENDDO
     168     DO j=1,jjp1
     169        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
     170     ENDDO
     171
     172     ! Potential temperature
     173     DO l=1,llm
     174        zsig=presnivs(l)/preff
     175        tetastrat=ttp*zsig**(-kappa)
     176        tetapv=tetastrat
     177        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
     178           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
     179        ENDIF
     180        DO j=1,jjp1
     181           ! Troposphere
     182           ddsin=sin(rlatu(j))
     183           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
     184                -delt_z*(1.-ddsin*ddsin)*log(zsig)
     185           !! Aymeric -- tests particuliers
     186           if (planet_type=="giant") then
     187             tetajl(j,l)=teta0+(delt_y*                   &
     188                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
     189                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
     190                -delt_z*log(zsig)
     191           endif
     192           ! Profil stratospherique isotherme (+vortex)
     193           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
     194           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
     195           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
     196        ENDDO
     197     ENDDO
     198
     199     !          CALL writefield('theta_eq',tetajl)
     200
     201     do l=1,llm
     202        do j=1,jjp1
     203           do i=1,iip1
     204              ij=(j-1)*iip1+i
     205              tetarappel(ij,l)=tetajl(j,l)
     206           enddo
     207        enddo
     208     enddo
     209     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
     210
     211     ! 3. Initialize fields (if necessary)
     212     IF (.NOT. read_start) THEN
     213        ! surface pressure
     214        ps(:)=preff
     215        ! ground geopotential
     216        phis(:)=0.
     217
     218        CALL pression ( ip1jmp1, ap, bp, ps, p       )
     219        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     220        CALL massdair(p,masse)
     221
     222        ! bulk initialization of temperature
     223        teta(:,:)=tetarappel(:,:)
     224
     225        ! geopotential
     226        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     227
     228        ! winds
     229        if (ok_geost) then
     230           call ugeostr(phi,ucov)
     231        else
     232           ucov(:,:)=0.
    122233        endif
    123        
    124         CALL iniconst
    125         CALL inigeom
    126         CALL inifilr
    127 
    128         if (llm.eq.1) then
    129           ! initialize fields for the shallow water case, if required
    130           if (.not.read_start) then
    131             phis(:)=0.
    132             q(:,:,:)=0
    133             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    134           endif
    135         endif
    136 
    137         if (iflag_phys.eq.2) then
    138           ! initializations for the academic case
    139          
    140 !         if (planet_type=="earth") then
    141 
    142           ! 1. local parameters
    143           ! by convention, winter is in the southern hemisphere
    144           ! Geostrophic wind or no wind?
    145           ok_geost=.TRUE.
    146           CALL getin('ok_geost',ok_geost)
    147           ! Constants for Newtonian relaxation and friction
    148           k_f=1.                !friction
    149           CALL getin('k_j',k_f)
    150           k_f=1./(daysec*k_f)
    151           k_c_s=4.  !cooling surface
    152           CALL getin('k_c_s',k_c_s)
    153           k_c_s=1./(daysec*k_c_s)
    154           k_c_a=40. !cooling free atm
    155           CALL getin('k_c_a',k_c_a)
    156           k_c_a=1./(daysec*k_c_a)
    157           ! Constants for Teta equilibrium profile
    158           teta0=315.     ! mean Teta (S.H. 315K)
    159           CALL getin('teta0',teta0)
    160           print *, 'iniacademic - teta0 ', teta0
    161           print *, 'iniacademic - rad ', rad
    162           ttp=200.       ! Tropopause temperature (S.H. 200K)
    163           CALL getin('ttp',ttp)
    164           eps=0.         ! Deviation to N-S symmetry(~0-20K)
    165           CALL getin('eps',eps)
    166           delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
    167           CALL getin('delt_y',delt_y)
    168           delt_z=10.     ! Vertical Gradient (S.H. 10K)
    169           CALL getin('delt_z',delt_z)
    170           ! Polar vortex
    171           ok_pv=.false.
    172           CALL getin('ok_pv',ok_pv)
    173           phi_pv=-50.            ! Latitude of edge of vortex
    174           CALL getin('phi_pv',phi_pv)
    175           phi_pv=phi_pv*pi/180.
    176           dphi_pv=5.             ! Width of the edge
    177           CALL getin('dphi_pv',dphi_pv)
    178           dphi_pv=dphi_pv*pi/180.
    179           gam_pv=4.              ! -dT/dz vortex (in K/km)
    180           CALL getin('gam_pv',gam_pv)
    181          
    182           ! 2. Initialize fields towards which to relax
    183           ! Friction
    184           knewt_g=k_c_a
    185           DO l=1,llm
    186            zsig=presnivs(l)/preff
    187            knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
    188            kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
    189           ENDDO
    190           DO j=1,jjp1
    191             clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
    192           ENDDO
    193          
    194           ! Potential temperature
    195           DO l=1,llm
    196            zsig=presnivs(l)/preff
    197            tetastrat=ttp*zsig**(-kappa)
    198            tetapv=tetastrat
    199            IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
    200                tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
    201            ENDIF
    202            DO j=1,jjp1
    203              ! Troposphere
    204              ddsin=sin(rlatu(j))
    205              tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
    206      &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    207              !! Aymeric -- tests particuliers
    208              if (planet_type=="giant") then
    209              tetajl(j,l)=teta0+(delt_y*
    210      &          ((sin(rlatu(j)*3.14159*eps+0.0001))**2)
    211      &          / ((rlatu(j)*3.14159*eps+0.0001)**2))
    212      &          -delt_z*log(zsig)
    213 !!!             ddsin=sin(2.5*3.14159*rlatu(j))
    214 !!!             tetajl(j,l)=teta0-delt_y*ddsin*ddsin
    215 !!!!     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    216              endif
    217              ! Profil stratospherique isotherme (+vortex)
    218              w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
    219              tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    220              tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
    221            ENDDO
    222           ENDDO ! of DO l=1,llm
    223  
    224 !          CALL writefield('theta_eq',tetajl)
    225 
    226           do l=1,llm
    227             do j=1,jjp1
    228               do i=1,iip1
    229                  ij=(j-1)*iip1+i
    230                  tetarappel(ij,l)=tetajl(j,l)
    231               enddo
    232             enddo
    233           enddo
    234           PRINT *, 'iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
    235 
    236 
    237 !         else
    238 !          write(lunout,*)"iniacademic: planet types other than earth",
    239 !     &                   " not implemented (yet)."
    240 !          stop
    241 !         endif ! of if (planet_type=="earth")
    242 
    243           ! 3. Initialize fields (if necessary)
    244           IF (.NOT. read_start) THEN
    245             ! surface pressure
    246             ps(:)=preff
    247             ! ground geopotential
    248             phis(:)=0.
    249            
    250             CALL pression ( ip1jmp1, ap, bp, ps, p       )
    251             CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    252             CALL massdair(p,masse)
    253 
    254             ! bulk initialization of temperature
    255             teta(:,:)=tetarappel(:,:)
    256            
    257             ! geopotential
    258             CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    259            
    260             ! winds
    261             if (ok_geost) then
    262               call ugeostr(phi,ucov)
    263             else
    264               ucov(:,:)=0.
    265             endif
    266             vcov(:,:)=0.
    267            
    268             ! bulk initialization of tracers
    269             if (planet_type=="earth") then
    270               ! Earth: first two tracers will be water
    271               do i=1,nqtot
    272                 if (i.eq.1) q(:,:,i)=1.e-10
    273                 if (i.eq.2) q(:,:,i)=1.e-15
    274                 if (i.gt.2) q(:,:,i)=0.
    275               enddo
    276             else
    277               q(:,:,:)=0
    278             endif ! of if (planet_type=="earth")
    279 
    280             ! add random perturbation to temperature
    281             idum  = -1
    282             zz = ran1(idum)
    283             idum  = 0
    284             do l=1,llm
    285               do ij=iip2,ip1jm
    286                 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
    287               enddo
    288             enddo
    289 
    290             ! maintain periodicity in longitude
    291             do l=1,llm
    292               do ij=1,ip1jmp1,iip1
    293                 teta(ij+iim,l)=teta(ij,l)
    294               enddo
    295             enddo
    296 
    297           ENDIF ! of IF (.NOT. read_start)
    298         endif ! of if (iflag_phys.eq.2)
    299        
    300       END
    301 c-----------------------------------------------------------------------
     234        vcov(:,:)=0.
     235
     236        ! bulk initialization of tracers
     237        if (planet_type=="earth") then
     238           ! Earth: first two tracers will be water
     239           do i=1,nqtot
     240              if (i == 1) q(:,:,i)=1.e-10
     241              if (i == 2) q(:,:,i)=1.e-15
     242              if (i.gt.2) q(:,:,i)=0.
     243           enddo
     244        else
     245           q(:,:,:)=0
     246        endif ! of if (planet_type=="earth")
     247
     248        ! add random perturbation to temperature
     249        idum  = -1
     250        zz = ran1(idum)
     251        idum  = 0
     252        do l=1,llm
     253           do ij=iip2,ip1jm
     254              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
     255           enddo
     256        enddo
     257
     258        ! maintain periodicity in longitude
     259        do l=1,llm
     260           do ij=1,ip1jmp1,iip1
     261              teta(ij+iim,l)=teta(ij,l)
     262           enddo
     263        enddo
     264
     265     ENDIF ! of IF (.NOT. read_start)
     266  endif academic_case
     267
     268END SUBROUTINE iniacademic

  • trunk/libf/dyn3d/ugeostr.F90

    r64 r66  
    11!
    2 ! $Id: ugeostr.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: ugeostr.F90 1474 2011-01-14 11:04:45Z lguez $
    33!
    4       subroutine ugeostr(phi,ucov)
     4subroutine ugeostr(phi,ucov)
    55
     6  ! Calcul du vent covariant geostrophique a partir du champ de
     7  ! geopotentiel.
     8  ! We actually compute: (1 - cos^8 \phi) u_g
     9  ! to have a wind going smoothly to 0 at the equator.
     10  ! We assume that the surface pressure is uniform so that model
     11  ! levels are pressure levels.
    612
    7 c  Calcul du vent covariant geostrophique a partir du champs de
    8 c  geopotentiel en supposant que le vent au sol est nul.
     13  implicit none
    914
    10       implicit none
     15  include "dimensions.h"
     16  include "paramet.h"
     17  include "comconst.h"
     18  include "comgeom2.h"
    1119
    12 #include "dimensions.h"
    13 #include "paramet.h"
    14 #include "comconst.h"
    15 #include "comgeom2.h"
     20  real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm)
     21  real um(jjm,llm),fact,u(iip1,jjm,llm)
     22  integer i,j,l
    1623
    17       real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm)
    18       real um(jjm,llm),fact,u(iip1,jjm,llm)
    19       integer i,j,l
     24  real zlat
    2025
    21       real zlat
     26  um(:,:)=0 ! initialize um()
    2227
    23       um(:,:)=0 ! initialize um()
     28  DO j=1,jjm
    2429
    25       DO j=1,jjm
     30     if (abs(sin(rlatv(j))).lt.1.e-4) then
     31        zlat=1.e-4
     32     else
     33        zlat=rlatv(j)
     34     endif
     35     fact=cos(zlat)
     36     fact=fact*fact
     37     fact=fact*fact
     38     fact=fact*fact
     39     fact=(1.-fact)/ &
     40          (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j)))
     41     fact=-fact/rad
     42     DO l=1,llm
     43        DO i=1,iim
     44           u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l))
     45           um(j,l)=um(j,l)+u(i,j,l)/REAL(iim)
     46        ENDDO
     47     ENDDO
     48  ENDDO
     49  call dump2d(jjm,llm,um,'Vent-u geostrophique')
    2650
    27          if (abs(sin(rlatv(j))).lt.1.e-4) then
    28              zlat=1.e-4
    29          else
    30              zlat=rlatv(j)
    31          endif
    32          fact=cos(zlat)
    33          fact=fact*fact
    34          fact=fact*fact
    35          fact=fact*fact
    36          fact=(1.-fact)/
    37      s    (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j)))
    38          fact=-fact/rad
    39          DO l=1,llm
    40             DO i=1,iim
    41                u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l))
    42                um(j,l)=um(j,l)+u(i,j,l)/REAL(iim)
    43             ENDDO
    44          ENDDO
    45       ENDDO
    46       call dump2d(jjm,llm,um,'Vent-u geostrophique')
     51  !   calcul des champ de vent:
    4752
    48 c
    49 c-----------------------------------------------------------------------
    50 c   calcul des champ de vent:
    51 c   -------------------------
     53  DO l=1,llm
     54     DO i=1,iip1
     55        ucov(i,1,l)=0.
     56        ucov(i,jjp1,l)=0.
     57     end DO
     58     DO  j=2,jjm
     59        DO  i=1,iim
     60           ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j)
     61        end DO
     62        ucov(iip1,j,l)=ucov(1,j,l)
     63     end DO
     64  end DO
    5265
    53       DO 301 l=1,llm
    54          DO 302 i=1,iip1
    55             ucov(i,1,l)=0.
    56             ucov(i,jjp1,l)=0.
    57 302      CONTINUE
    58          DO 304 j=2,jjm
    59             DO 305 i=1,iim
    60                ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j)
    61 305         CONTINUE
    62             ucov(iip1,j,l)=ucov(1,j,l)
    63 304      CONTINUE
    64 301   CONTINUE
     66  print *, 301
    6567
    66       print*,301
    67 
    68       return
    69       end
     68end subroutine ugeostr
Note: See TracChangeset for help on using the changeset viewer.