Changeset 66 for trunk/libf/dyn3dpar/iniacademic.F90
- Timestamp:
- Feb 16, 2011, 4:57:45 PM (14 years ago)
- File:
-
- 1 moved
-
trunk/libf/dyn3dpar/iniacademic.F90 (moved) (moved from trunk/libf/dyn3dpar/iniacademic.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/libf/dyn3dpar/iniacademic.F90
r64 r66 1 1 ! 2 ! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z emillour$2 ! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $ 3 3 ! 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 4 SUBROUTINE 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 11 9 #ifdef CPP_IOIPSL 12 10 USE IOIPSL 13 11 #else 14 ! if not using IOIPSL, we still need to use (a local version of) getin15 12 ! if not using IOIPSL, we still need to use (a local version of) getin 13 USE ioipsl_getincom 16 14 #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. 122 233 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 268 END SUBROUTINE iniacademic
Note: See TracChangeset
for help on using the changeset viewer.
