source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/iniacademic.F90 @ 3110

Last change on this file since 3110 was 1529, checked in by Laurent Fairhead, 13 years ago

Necessary modifications to use the model in an aquaplanet or terraplanet configuration


Modifications nécessaires pour l'utilisation du moèle en configuration aquaplanète
ou terraplanète

Dans la dynamique:
Changement pour le paramètre iflag_phys:

auparavant : 0 pas de physique, 1 phyique, 2 rappel
ici on ajoute iflag_phys>=100 pour les aqua et terra

Du coup, dans leapfrog.F et gcm.F on appelle la physique si
iflag_phys=1.or.iflag_phys>=100
Dans iniacademic, on initialise si iflag_phys>=2 au lieu de =2
Dans gcm.F, on appelle en plus iniaqua (sous une clef CPP_EARTH)
Dans iniacademic, on met ps=108080 pour les aqua et terra pour répondre
à une specification CFMIP.

Dans la physique:
On ajoute phyaqua.F qui contient :
iniaqua : initialise les startphy.nc et limit.nc pour la phyique
zenang_an : calcule un ensolleillement moyen sur l'année en fonction de

la latiude (A. Campoy).

profil_sst : qui calcule différents profils latitudinaux de SSTs.
writelim : écrit le fichier limit.nc

Dans physiq.F
On ajoute la possibilité d'appeller un calcul d'ensoleillement moyen
sur l'année quand solarlong0=1000.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.8 KB
Line 
1!
2! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $
3!
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
9#ifdef CPP_IOIPSL
10  USE IOIPSL
11#else
12  ! if not using IOIPSL, we still need to use (a local version of) getin
13  USE ioipsl_getincom
14#endif
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  character(len=*),parameter :: modname="iniacademic"
75  character(len=80) :: abort_message
76
77  !-----------------------------------------------------------------------
78  ! 1. Initializations for Earth-like case
79  ! --------------------------------------
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     ttp=200.       ! Tropopause temperature (S.H. 200K)
139     CALL getin('ttp',ttp)
140     eps=0.         ! Deviation to N-S symmetry(~0-20K)
141     CALL getin('eps',eps)
142     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
143     CALL getin('delt_y',delt_y)
144     delt_z=10.     ! Vertical Gradient (S.H. 10K)
145     CALL getin('delt_z',delt_z)
146     ! Polar vortex
147     ok_pv=.false.
148     CALL getin('ok_pv',ok_pv)
149     phi_pv=-50.            ! Latitude of edge of vortex
150     CALL getin('phi_pv',phi_pv)
151     phi_pv=phi_pv*pi/180.
152     dphi_pv=5.             ! Width of the edge
153     CALL getin('dphi_pv',dphi_pv)
154     dphi_pv=dphi_pv*pi/180.
155     gam_pv=4.              ! -dT/dz vortex (in K/km)
156     CALL getin('gam_pv',gam_pv)
157
158     ! 2. Initialize fields towards which to relax
159     ! Friction
160     knewt_g=k_c_a
161     DO l=1,llm
162        zsig=presnivs(l)/preff
163        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
164        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
165     ENDDO
166     DO j=1,jjp1
167        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
168     ENDDO
169
170     ! Potential temperature
171     DO l=1,llm
172        zsig=presnivs(l)/preff
173        tetastrat=ttp*zsig**(-kappa)
174        tetapv=tetastrat
175        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
176           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
177        ENDIF
178        DO j=1,jjp1
179           ! Troposphere
180           ddsin=sin(rlatu(j))
181           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
182                -delt_z*(1.-ddsin*ddsin)*log(zsig)
183           if (planet_type=="giant") then
184             tetajl(j,l)=teta0+(delt_y*                   &
185                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
186                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
187                -delt_z*log(zsig)
188           endif
189           ! Profil stratospherique isotherme (+vortex)
190           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
191           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
192           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
193        ENDDO
194     ENDDO
195
196     !          CALL writefield('theta_eq',tetajl)
197
198     do l=1,llm
199        do j=1,jjp1
200           do i=1,iip1
201              ij=(j-1)*iip1+i
202              tetarappel(ij,l)=tetajl(j,l)
203           enddo
204        enddo
205     enddo
206
207     ! 3. Initialize fields (if necessary)
208     IF (.NOT. read_start) THEN
209        ! surface pressure
210        if (iflag_phys>2) then
211           ps(:)=preff
212        else
213           ps(:)=101080.
214        endif
215        ! ground geopotential
216        phis(:)=0.
217
218        CALL pression ( ip1jmp1, ap, bp, ps, p       )
219        if (disvert_type.eq.1) then
220          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
221        elseif (disvert_type.eq.2) then
222          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
223        else
224          write(abort_message,*) "Wrong value for disvert_type: ", &
225                              disvert_type
226          call abort_gcm(modname,abort_message,0)
227        endif
228        CALL massdair(p,masse)
229
230        ! bulk initialization of temperature
231        teta(:,:)=tetarappel(:,:)
232
233        ! geopotential
234        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
235
236        ! winds
237        if (ok_geost) then
238           call ugeostr(phi,ucov)
239        else
240           ucov(:,:)=0.
241        endif
242        vcov(:,:)=0.
243
244        ! bulk initialization of tracers
245        if (planet_type=="earth") then
246           ! Earth: first two tracers will be water
247           do i=1,nqtot
248              if (i == 1) q(:,:,i)=1.e-10
249              if (i == 2) q(:,:,i)=1.e-15
250              if (i.gt.2) q(:,:,i)=0.
251           enddo
252        else
253           q(:,:,:)=0
254        endif ! of if (planet_type=="earth")
255
256        ! add random perturbation to temperature
257        idum  = -1
258        zz = ran1(idum)
259        idum  = 0
260        do l=1,llm
261           do ij=iip2,ip1jm
262              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
263           enddo
264        enddo
265
266        ! maintain periodicity in longitude
267        do l=1,llm
268           do ij=1,ip1jmp1,iip1
269              teta(ij+iim,l)=teta(ij,l)
270           enddo
271        enddo
272
273     ENDIF ! of IF (.NOT. read_start)
274  endif academic_case
275
276END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.