Ignore:
Timestamp:
Dec 23, 2010, 5:38:42 PM (14 years ago)
Author:
lguez
Message:

Conversion to free source form for "disvert" and "iniacademic", no
other change. Bug fix in "fisrtilp": "fraca" could appear in an
expression while not defined.

Location:
LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d
Files:
2 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90

    r1463 r1472  
    1 !
    21! $Id$
    3 !
    4       SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    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, *)'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,*) 'ATTENTION discretisation z a ajuster'
     76           dsigmin=1.
     77        endif
     78        WRITE(LUNOUT,*) '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 = pi * (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                *(1. - tanh(2 * x / pi - 1.))**2 / 4.
     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, *)' BP '
     122  write(lunout, *) bp
     123  write(lunout, *)' 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, *) '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

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90

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