Changeset 1454 for LMDZ5/trunk/libf/dyn3d/iniacademic.F
- Timestamp:
- Nov 18, 2010, 1:01:24 PM (14 years ago)
- Location:
- LMDZ5/trunk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/dyn3d/iniacademic.F (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk
- Property svn:mergeinfo changed
/LMDZ5/branches/LMDZ5V1.0-dev (added) merged: 1436-1438,1441-1449,1452-1453
- Property svn:mergeinfo changed
-
LMDZ5/trunk/libf/dyn3d/iniacademic.F
r1403 r1454 8 8 USE filtreg_mod 9 9 USE infotrac, ONLY : nqtot 10 USE control_mod 10 USE control_mod, ONLY: day_step,planet_type 11 #ifdef CPP_IOIPSL 12 USE IOIPSL 13 #else 14 ! if not using IOIPSL, we still need to use (a local version of) getin 15 USE ioipsl_getincom 16 #endif 17 USE Write_Field 11 18 12 19 … … 70 77 REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches 71 78 REAL phi(ip1jmp1,llm) ! geopotentiel 72 REAL ddsin,teta rappelj,tetarappell,zsig79 REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires 73 80 real tetajl(jjp1,llm) 74 81 INTEGER i,j,l,lsup,ij 75 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 76 89 real zz,ran1 77 90 integer idum … … 82 95 ! 1. Initializations for Earth-like case 83 96 ! -------------------------------------- 84 if (planet_type=="earth") then 85 c 97 c 98 ! initialize planet radius, rotation rate,... 99 call conf_planete 100 86 101 time_0=0. 87 day_ref= 0102 day_ref=1 88 103 annee_ref=0 89 104 90 105 im = iim 91 106 jm = jjm 92 day_ini = 0 93 omeg = 4.*asin(1.)/86400. 94 rad = 6371229. 95 g = 9.8 96 daysec = 86400. 107 day_ini = 1 97 108 dtvr = daysec/REAL(day_step) 98 109 zdtvr=dtvr 99 kappa = 0.2857143100 cpp = 1004.70885101 preff = 101325.102 pa = 50000.103 110 etot0 = 0. 104 111 ptot0 = 0. … … 120 127 if (.not.read_start) then 121 128 phis(:)=0. 122 q(:,:,1)=1.e-10 123 q(:,:,2)=1.e-15 124 q(:,:,3:nqtot)=0. 129 q(:,:,:)=0 125 130 CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) 126 131 endif … … 129 134 if (iflag_phys.eq.2) then 130 135 ! initializations for the academic case 131 ps(:)=1.e5 132 phis(:)=0. 133 c--------------------------------------------------------------------- 134 135 taurappel=10.*daysec 136 137 c--------------------------------------------------------------------- 138 c Calcul de la temperature potentielle : 139 c -------------------------------------- 140 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 141 180 DO l=1,llm 142 zsig=ap(l)/preff+bp(l) 143 if (zsig.gt.0.3) then 144 lsup=l 145 tetarappell=1./8.*(-log(zsig)-.5) 146 DO j=1,jjp1 147 ddsin=sin(rlatu(j))-sin(pi/20.) 148 tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell) 149 ENDDO 150 else 151 c Choix isotherme au-dessus de 300 mbar 152 do j=1,jjp1 153 tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa 154 enddo 155 endif ! of if (zsig.gt.0.3) 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 156 207 ENDDO ! of DO l=1,llm 208 209 ! CALL writefield('theta_eq',tetajl) 157 210 158 211 do l=1,llm … … 165 218 enddo 166 219 167 c call dump2d(jjp1,llm,tetajl,'TEQ ') 168 169 CALL pression ( ip1jmp1, ap, bp, ps, p ) 170 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 171 CALL massdair(p,masse) 172 173 c intialisation du vent et de la temperature 174 teta(:,:)=tetarappel(:,:) 175 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 176 call ugeostr(phi,ucov) 177 vcov=0. 178 q(:,:,1 )=1.e-10 179 q(:,:,2 )=1.e-15 180 q(:,:,3:nqtot)=0. 181 182 183 c perturbation aleatoire sur la temperature 184 idum = -1 185 zz = ran1(idum) 186 idum = 0 187 do l=1,llm 188 do ij=iip2,ip1jm 189 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) 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 190 272 enddo 191 enddo 192 193 do l=1,llm 194 do ij=1,ip1jmp1,iip1 195 teta(ij+iim,l)=teta(ij,l) 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 196 279 enddo 197 enddo 198 199 200 201 c PRINT *,' Appel test_period avec tetarappel ' 202 c CALL test_period ( ucov,vcov,tetarappel,q,p,phis ) 203 c PRINT *,' Appel test_period avec teta ' 204 c CALL test_period ( ucov,vcov,teta,q,p,phis ) 205 206 c initialisation d'un traceur sur une colonne 207 j=jjp1*3/4 208 i=iip1/2 209 ij=(j-1)*iip1+i 210 q(ij,:,3)=1. 280 281 ENDIF ! of IF (.NOT. read_start) 211 282 endif ! of if (iflag_phys.eq.2) 212 283 213 else214 write(lunout,*)"iniacademic: planet types other than earth",215 & " not implemented (yet)."216 stop217 endif ! of if (planet_type=="earth")218 return219 284 END 220 285 c-----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
