Changeset 1437 for LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/iniacademic.F
- Timestamp:
- Sep 30, 2010, 10:29:10 AM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/iniacademic.F
r1403 r1437 9 9 USE infotrac, ONLY : nqtot 10 10 USE control_mod 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 … … 84 97 if (planet_type=="earth") then 85 98 c 99 ! initialize planet radius, rotation rate,... 100 call conf_planete 101 86 102 time_0=0. 87 day_ref= 0103 day_ref=1 88 104 annee_ref=0 89 105 90 106 im = iim 91 107 jm = jjm 92 day_ini = 0 93 omeg = 4.*asin(1.)/86400. 94 rad = 6371229. 95 g = 9.8 96 daysec = 86400. 108 day_ini = 1 97 109 dtvr = daysec/REAL(day_step) 98 110 zdtvr=dtvr 99 kappa = 0.2857143100 cpp = 1004.70885101 preff = 101325.102 pa = 50000.103 111 etot0 = 0. 104 112 ptot0 = 0. … … 129 137 if (iflag_phys.eq.2) then 130 138 ! 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 139 140 ! 1. local parameters 141 ! by convention, winter is in the southern hemisphere 142 ! Geostrophic wind or no wind? 143 ok_geost=.TRUE. 144 CALL getin('ok_geost',ok_geost) 145 ! Constants for Newtonian relaxation and friction 146 k_f=1. !friction 147 CALL getin('k_j',k_f) 148 k_f=1./(daysec*k_f) 149 k_c_s=4. !cooling surface 150 CALL getin('k_c_s',k_c_s) 151 k_c_s=1./(daysec*k_c_s) 152 k_c_a=40. !cooling free atm 153 CALL getin('k_c_a',k_c_a) 154 k_c_a=1./(daysec*k_c_a) 155 ! Constants for Teta equilibrium profile 156 teta0=315. ! mean Teta (S.H. 315K) 157 CALL getin('teta0',teta0) 158 ttp=200. ! Tropopause temperature (S.H. 200K) 159 CALL getin('ttp',ttp) 160 eps=0. ! Deviation to N-S symmetry(~0-20K) 161 CALL getin('eps',eps) 162 delt_y=60. ! Merid Temp. Gradient (S.H. 60K) 163 CALL getin('delt_y',delt_y) 164 delt_z=10. ! Vertical Gradient (S.H. 10K) 165 CALL getin('delt_z',delt_z) 166 ! Polar vortex 167 ok_pv=.false. 168 CALL getin('ok_pv',ok_pv) 169 phi_pv=-50. ! Latitude of edge of vortex 170 CALL getin('phi_pv',phi_pv) 171 phi_pv=phi_pv*pi/180. 172 dphi_pv=5. ! Width of the edge 173 CALL getin('dphi_pv',dphi_pv) 174 dphi_pv=dphi_pv*pi/180. 175 gam_pv=4. ! -dT/dz vortex (in K/km) 176 CALL getin('gam_pv',gam_pv) 177 178 ! 2. Initialize fields towards which to relax 179 ! Friction 180 knewt_g=k_c_a 141 181 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) 182 zsig=presnivs(l)/preff 183 knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 184 kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 185 ENDDO 186 DO j=1,jjp1 187 clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 188 ENDDO 189 190 ! Potential temperature 191 DO l=1,llm 192 zsig=presnivs(l)/preff 193 tetastrat=ttp*zsig**(-kappa) 194 tetapv=tetastrat 195 IF ((ok_pv).AND.(zsig.LT.0.1)) THEN 196 tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) 197 ENDIF 198 DO j=1,jjp1 199 ! Troposphere 200 ddsin=sin(rlatu(j)) 201 tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin 202 & -delt_z*(1.-ddsin*ddsin)*log(zsig) 203 ! Profil stratospherique isotherme (+vortex) 204 w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. 205 tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv 206 tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 207 ENDDO 156 208 ENDDO ! of DO l=1,llm 209 210 ! CALL writefield('theta_eq',tetajl) 157 211 158 212 do l=1,llm … … 165 219 enddo 166 220 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)) 190 enddo 191 enddo 192 193 do l=1,llm 194 do ij=1,ip1jmp1,iip1 195 teta(ij+iim,l)=teta(ij,l) 196 enddo 197 enddo 198 199 221 ! 3. Initialize fields (if necessary) 222 IF (.NOT. read_start) THEN 223 ! surface pressure 224 ps(:)=preff 225 ! ground geopotential 226 phis(:)=0. 227 228 CALL pression ( ip1jmp1, ap, bp, ps, p ) 229 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 230 CALL massdair(p,masse) 231 232 ! bulk initialization of temperature 233 teta(:,:)=tetarappel(:,:) 234 235 ! geopotential 236 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 237 238 ! winds 239 if (ok_geost) then 240 call ugeostr(phi,ucov) 241 else 242 ucov(:,:)=0. 243 endif 244 vcov(:,:)=0. 245 246 ! bulk initialization of tracers 247 do i=1,nqtot 248 if (i.eq.1) q(:,:,i)=1.e-10 249 if (i.eq.2) q(:,:,i)=1.e-15 250 if (i.gt.2) q(:,:,i)=0. 251 enddo 252 253 ! add random perturbation to temperature 254 idum = -1 255 zz = ran1(idum) 256 idum = 0 257 do l=1,llm 258 do ij=iip2,ip1jm 259 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) 260 enddo 261 enddo 262 263 do l=1,llm 264 do ij=1,ip1jmp1,iip1 265 teta(ij+iim,l)=teta(ij,l) 266 enddo 267 enddo 200 268 201 269 c PRINT *,' Appel test_period avec tetarappel ' … … 204 272 c CALL test_period ( ucov,vcov,teta,q,p,phis ) 205 273 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. 274 ! initialize a traceur on one column 275 ! j=jjp1*3/4 276 ! i=iip1/2 277 ! ij=(j-1)*iip1+i 278 ! q(ij,:,3)=1. 279 280 ENDIF ! of IF (.NOT. read_start) 211 281 endif ! of if (iflag_phys.eq.2) 212 282
Note: See TracChangeset
for help on using the changeset viewer.