| Line |
Branch |
Exec |
Source |
| 1 |
|
|
! |
| 2 |
|
|
! $Id: iniacademic.F90 3976 2021-08-06 14:30:26Z fhourdin $ |
| 3 |
|
|
! |
| 4 |
|
✗ |
SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) |
| 5 |
|
|
|
| 6 |
|
|
USE filtreg_mod, ONLY: inifilr |
| 7 |
|
|
USE infotrac |
| 8 |
|
|
USE control_mod, ONLY: day_step,planet_type |
| 9 |
|
|
USE IOIPSL, ONLY: getin |
| 10 |
|
|
USE Write_Field |
| 11 |
|
|
use exner_hyb_m, only: exner_hyb |
| 12 |
|
|
use exner_milieu_m, only: exner_milieu |
| 13 |
|
|
USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm |
| 14 |
|
|
USE logic_mod, ONLY: iflag_phys, read_start |
| 15 |
|
|
USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner |
| 16 |
|
|
USE temps_mod, ONLY: annee_ref, day_ini, day_ref |
| 17 |
|
|
USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 |
| 18 |
|
|
|
| 19 |
|
|
! Author: Frederic Hourdin original: 15/01/93 |
| 20 |
|
|
! The forcing defined here is from Held and Suarez, 1994, Bulletin |
| 21 |
|
|
! of the American Meteorological Society, 75, 1825. |
| 22 |
|
|
|
| 23 |
|
|
IMPLICIT NONE |
| 24 |
|
|
|
| 25 |
|
|
! Declararations: |
| 26 |
|
|
! --------------- |
| 27 |
|
|
|
| 28 |
|
|
include "dimensions.h" |
| 29 |
|
|
include "paramet.h" |
| 30 |
|
|
include "comgeom.h" |
| 31 |
|
|
include "academic.h" |
| 32 |
|
|
include "iniprint.h" |
| 33 |
|
|
|
| 34 |
|
|
! Arguments: |
| 35 |
|
|
! ---------- |
| 36 |
|
|
|
| 37 |
|
|
REAL,INTENT(OUT) :: time_0 |
| 38 |
|
|
|
| 39 |
|
|
! fields |
| 40 |
|
|
REAL,INTENT(OUT) :: vcov(ip1jm,llm) ! meridional covariant wind |
| 41 |
|
|
REAL,INTENT(OUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind |
| 42 |
|
|
REAL,INTENT(OUT) :: teta(ip1jmp1,llm) ! potential temperature (K) |
| 43 |
|
|
REAL,INTENT(OUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers (.../kg_of_air) |
| 44 |
|
|
REAL,INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa) |
| 45 |
|
|
REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass in grid cell (kg) |
| 46 |
|
|
REAL,INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential |
| 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 phi(ip1jmp1,llm) ! geopotentiel |
| 55 |
|
|
REAL ddsin,zsig,tetapv,w_pv ! variables auxiliaires |
| 56 |
|
|
real tetastrat ! potential temperature in the stratosphere, in K |
| 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,tetanoise ! Constantes pour polar vortex |
| 65 |
|
|
|
| 66 |
|
|
real zz,ran1 |
| 67 |
|
|
integer idum |
| 68 |
|
|
|
| 69 |
|
|
REAL zdtvr |
| 70 |
|
|
|
| 71 |
|
|
character(len=*),parameter :: modname="iniacademic" |
| 72 |
|
|
character(len=80) :: abort_message |
| 73 |
|
|
|
| 74 |
|
|
|
| 75 |
|
|
! Sanity check: verify that options selected by user are not incompatible |
| 76 |
|
✗ |
if ((iflag_phys==1).and. .not. read_start) then |
| 77 |
|
✗ |
write(lunout,*) trim(modname)," error: if read_start is set to ", & |
| 78 |
|
✗ |
" false then iflag_phys should not be 1" |
| 79 |
|
✗ |
write(lunout,*) "You most likely want an aquaplanet initialisation", & |
| 80 |
|
✗ |
" (iflag_phys >= 100)" |
| 81 |
|
✗ |
call abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.false.",1) |
| 82 |
|
|
endif |
| 83 |
|
|
|
| 84 |
|
|
!----------------------------------------------------------------------- |
| 85 |
|
|
! 1. Initializations for Earth-like case |
| 86 |
|
|
! -------------------------------------- |
| 87 |
|
|
! |
| 88 |
|
|
! initialize planet radius, rotation rate,... |
| 89 |
|
✗ |
call conf_planete |
| 90 |
|
|
|
| 91 |
|
✗ |
time_0=0. |
| 92 |
|
✗ |
day_ref=1 |
| 93 |
|
✗ |
annee_ref=0 |
| 94 |
|
|
|
| 95 |
|
✗ |
im = iim |
| 96 |
|
✗ |
jm = jjm |
| 97 |
|
✗ |
day_ini = 1 |
| 98 |
|
✗ |
dtvr = daysec/REAL(day_step) |
| 99 |
|
|
zdtvr=dtvr |
| 100 |
|
✗ |
etot0 = 0. |
| 101 |
|
✗ |
ptot0 = 0. |
| 102 |
|
✗ |
ztot0 = 0. |
| 103 |
|
✗ |
stot0 = 0. |
| 104 |
|
✗ |
ang0 = 0. |
| 105 |
|
|
|
| 106 |
|
|
if (llm == 1) then |
| 107 |
|
|
! specific initializations for the shallow water case |
| 108 |
|
|
kappa=1 |
| 109 |
|
|
endif |
| 110 |
|
|
|
| 111 |
|
✗ |
CALL iniconst |
| 112 |
|
✗ |
CALL inigeom |
| 113 |
|
✗ |
CALL inifilr |
| 114 |
|
|
|
| 115 |
|
|
! Initialize pressure and mass field if read_start=.false. |
| 116 |
|
✗ |
IF (.NOT. read_start) THEN |
| 117 |
|
|
! surface pressure |
| 118 |
|
✗ |
if (iflag_phys>2) then |
| 119 |
|
|
! specific value for CMIP5 aqua/terra planets |
| 120 |
|
|
! "Specify the initial dry mass to be equivalent to |
| 121 |
|
|
! a global mean surface pressure (101325 minus 245) Pa." |
| 122 |
|
✗ |
ps(:)=101080. |
| 123 |
|
|
else |
| 124 |
|
|
! use reference surface pressure |
| 125 |
|
✗ |
ps(:)=preff |
| 126 |
|
|
endif |
| 127 |
|
|
! ground geopotential |
| 128 |
|
✗ |
phis(:)=0. |
| 129 |
|
✗ |
CALL pression ( ip1jmp1, ap, bp, ps, p ) |
| 130 |
|
✗ |
if (pressure_exner) then |
| 131 |
|
✗ |
CALL exner_hyb( ip1jmp1, ps, p, pks, pk) |
| 132 |
|
|
else |
| 133 |
|
✗ |
call exner_milieu(ip1jmp1,ps,p,pks,pk) |
| 134 |
|
|
endif |
| 135 |
|
✗ |
CALL massdair(p,masse) |
| 136 |
|
|
ENDIF |
| 137 |
|
|
|
| 138 |
|
|
if (llm == 1) then |
| 139 |
|
|
! initialize fields for the shallow water case, if required |
| 140 |
|
|
if (.not.read_start) then |
| 141 |
|
|
phis(:)=0. |
| 142 |
|
|
q(:,:,:)=0 |
| 143 |
|
|
CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) |
| 144 |
|
|
endif |
| 145 |
|
|
endif |
| 146 |
|
|
|
| 147 |
|
✗ |
academic_case: if (iflag_phys >= 2) then |
| 148 |
|
|
! initializations |
| 149 |
|
|
|
| 150 |
|
|
! 1. local parameters |
| 151 |
|
|
! by convention, winter is in the southern hemisphere |
| 152 |
|
|
! Geostrophic wind or no wind? |
| 153 |
|
✗ |
ok_geost=.TRUE. |
| 154 |
|
✗ |
CALL getin('ok_geost',ok_geost) |
| 155 |
|
|
! Constants for Newtonian relaxation and friction |
| 156 |
|
✗ |
k_f=1. !friction |
| 157 |
|
✗ |
CALL getin('k_j',k_f) |
| 158 |
|
✗ |
k_f=1./(daysec*k_f) |
| 159 |
|
✗ |
k_c_s=4. !cooling surface |
| 160 |
|
✗ |
CALL getin('k_c_s',k_c_s) |
| 161 |
|
✗ |
k_c_s=1./(daysec*k_c_s) |
| 162 |
|
✗ |
k_c_a=40. !cooling free atm |
| 163 |
|
✗ |
CALL getin('k_c_a',k_c_a) |
| 164 |
|
✗ |
k_c_a=1./(daysec*k_c_a) |
| 165 |
|
|
! Constants for Teta equilibrium profile |
| 166 |
|
✗ |
teta0=315. ! mean Teta (S.H. 315K) |
| 167 |
|
✗ |
CALL getin('teta0',teta0) |
| 168 |
|
✗ |
ttp=200. ! Tropopause temperature (S.H. 200K) |
| 169 |
|
✗ |
CALL getin('ttp',ttp) |
| 170 |
|
✗ |
eps=0. ! Deviation to N-S symmetry(~0-20K) |
| 171 |
|
✗ |
CALL getin('eps',eps) |
| 172 |
|
✗ |
delt_y=60. ! Merid Temp. Gradient (S.H. 60K) |
| 173 |
|
✗ |
CALL getin('delt_y',delt_y) |
| 174 |
|
✗ |
delt_z=10. ! Vertical Gradient (S.H. 10K) |
| 175 |
|
✗ |
CALL getin('delt_z',delt_z) |
| 176 |
|
|
! Polar vortex |
| 177 |
|
✗ |
ok_pv=.false. |
| 178 |
|
✗ |
CALL getin('ok_pv',ok_pv) |
| 179 |
|
✗ |
phi_pv=-50. ! Latitude of edge of vortex |
| 180 |
|
✗ |
CALL getin('phi_pv',phi_pv) |
| 181 |
|
✗ |
phi_pv=phi_pv*pi/180. |
| 182 |
|
✗ |
dphi_pv=5. ! Width of the edge |
| 183 |
|
✗ |
CALL getin('dphi_pv',dphi_pv) |
| 184 |
|
✗ |
dphi_pv=dphi_pv*pi/180. |
| 185 |
|
✗ |
gam_pv=4. ! -dT/dz vortex (in K/km) |
| 186 |
|
✗ |
CALL getin('gam_pv',gam_pv) |
| 187 |
|
✗ |
tetanoise=0.005 |
| 188 |
|
✗ |
CALL getin('tetanoise',tetanoise) |
| 189 |
|
|
|
| 190 |
|
|
|
| 191 |
|
|
! 2. Initialize fields towards which to relax |
| 192 |
|
|
! Friction |
| 193 |
|
✗ |
knewt_g=k_c_a |
| 194 |
|
✗ |
DO l=1,llm |
| 195 |
|
✗ |
zsig=presnivs(l)/preff |
| 196 |
|
✗ |
knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) |
| 197 |
|
✗ |
kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) |
| 198 |
|
|
ENDDO |
| 199 |
|
✗ |
DO j=1,jjp1 |
| 200 |
|
✗ |
clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 |
| 201 |
|
|
ENDDO |
| 202 |
|
|
|
| 203 |
|
|
! Potential temperature |
| 204 |
|
✗ |
DO l=1,llm |
| 205 |
|
✗ |
zsig=presnivs(l)/preff |
| 206 |
|
✗ |
tetastrat=ttp*zsig**(-kappa) |
| 207 |
|
|
tetapv=tetastrat |
| 208 |
|
✗ |
IF ((ok_pv).AND.(zsig.LT.0.1)) THEN |
| 209 |
|
✗ |
tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) |
| 210 |
|
|
ENDIF |
| 211 |
|
✗ |
DO j=1,jjp1 |
| 212 |
|
|
! Troposphere |
| 213 |
|
✗ |
ddsin=sin(rlatu(j)) |
| 214 |
|
|
tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin & |
| 215 |
|
✗ |
-delt_z*(1.-ddsin*ddsin)*log(zsig) |
| 216 |
|
✗ |
if (planet_type=="giant") then |
| 217 |
|
|
tetajl(j,l)=teta0+(delt_y* & |
| 218 |
|
|
((sin(rlatu(j)*3.14159*eps+0.0001))**2) & |
| 219 |
|
|
/ ((rlatu(j)*3.14159*eps+0.0001)**2)) & |
| 220 |
|
✗ |
-delt_z*log(zsig) |
| 221 |
|
|
endif |
| 222 |
|
|
! Profil stratospherique isotherme (+vortex) |
| 223 |
|
✗ |
w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. |
| 224 |
|
✗ |
tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv |
| 225 |
|
✗ |
tetajl(j,l)=MAX(tetajl(j,l),tetastrat) |
| 226 |
|
|
ENDDO |
| 227 |
|
|
ENDDO |
| 228 |
|
|
|
| 229 |
|
|
! CALL writefield('theta_eq',tetajl) |
| 230 |
|
|
|
| 231 |
|
✗ |
do l=1,llm |
| 232 |
|
✗ |
do j=1,jjp1 |
| 233 |
|
✗ |
do i=1,iip1 |
| 234 |
|
✗ |
ij=(j-1)*iip1+i |
| 235 |
|
✗ |
tetarappel(ij,l)=tetajl(j,l) |
| 236 |
|
|
enddo |
| 237 |
|
|
enddo |
| 238 |
|
|
enddo |
| 239 |
|
|
|
| 240 |
|
|
! 3. Initialize fields (if necessary) |
| 241 |
|
✗ |
IF (.NOT. read_start) THEN |
| 242 |
|
|
! bulk initialization of temperature |
| 243 |
|
|
|
| 244 |
|
✗ |
IF (iflag_phys>10000) THEN |
| 245 |
|
|
! Particular case to impose a constant temperature T0=0.01*iflag_physx |
| 246 |
|
✗ |
teta(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp) |
| 247 |
|
|
ELSE |
| 248 |
|
✗ |
teta(:,:)=tetarappel(:,:) |
| 249 |
|
|
ENDIF |
| 250 |
|
|
|
| 251 |
|
|
! geopotential |
| 252 |
|
✗ |
CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) |
| 253 |
|
|
|
| 254 |
|
✗ |
DO l=1,llm |
| 255 |
|
✗ |
print*,'presnivs,play,l',presnivs(l),(pk(1,l)/cpp)**(1./kappa)*preff |
| 256 |
|
|
!pks(ij) = (cpp/preff) * ps(ij) |
| 257 |
|
|
!pk(ij,1) = .5*pks(ij) |
| 258 |
|
|
! pk = cpp * (p/preff)^kappa |
| 259 |
|
|
ENDDO |
| 260 |
|
|
|
| 261 |
|
|
! winds |
| 262 |
|
✗ |
if (ok_geost) then |
| 263 |
|
✗ |
call ugeostr(phi,ucov) |
| 264 |
|
|
else |
| 265 |
|
✗ |
ucov(:,:)=0. |
| 266 |
|
|
endif |
| 267 |
|
✗ |
vcov(:,:)=0. |
| 268 |
|
|
|
| 269 |
|
|
! bulk initialization of tracers |
| 270 |
|
✗ |
if (planet_type=="earth") then |
| 271 |
|
|
! Earth: first two tracers will be water |
| 272 |
|
✗ |
do i=1,nqtot |
| 273 |
|
✗ |
if (i == 1) q(:,:,i)=1.e-10 |
| 274 |
|
✗ |
if (i == 2) q(:,:,i)=1.e-15 |
| 275 |
|
✗ |
if (i.gt.2) q(:,:,i)=0. |
| 276 |
|
|
|
| 277 |
|
|
! CRisi: init des isotopes |
| 278 |
|
|
! distill de Rayleigh très simplifiée |
| 279 |
|
✗ |
if (ok_isotopes) then |
| 280 |
|
✗ |
if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then |
| 281 |
|
|
q(:,:,i)=q(:,:,iqpere(i)) & |
| 282 |
|
|
& *tnat(iso_num(i)) & |
| 283 |
|
|
& *(q(:,:,iqpere(i))/30.e-3) & |
| 284 |
|
✗ |
& **(alpha_ideal(iso_num(i))-1) |
| 285 |
|
|
endif |
| 286 |
|
✗ |
if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then |
| 287 |
|
✗ |
q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i))) |
| 288 |
|
|
endif |
| 289 |
|
|
endif !if (ok_isotopes) then |
| 290 |
|
|
|
| 291 |
|
|
enddo |
| 292 |
|
|
else |
| 293 |
|
✗ |
q(:,:,:)=0 |
| 294 |
|
|
endif ! of if (planet_type=="earth") |
| 295 |
|
|
|
| 296 |
|
✗ |
if (ok_iso_verif) then |
| 297 |
|
✗ |
call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc') |
| 298 |
|
|
endif !if (ok_iso_verif) then |
| 299 |
|
|
|
| 300 |
|
|
! add random perturbation to temperature |
| 301 |
|
✗ |
idum = -1 |
| 302 |
|
✗ |
zz = ran1(idum) |
| 303 |
|
✗ |
idum = 0 |
| 304 |
|
✗ |
do l=1,llm |
| 305 |
|
✗ |
do ij=iip2,ip1jm |
| 306 |
|
✗ |
teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum)) |
| 307 |
|
|
enddo |
| 308 |
|
|
enddo |
| 309 |
|
|
|
| 310 |
|
|
! maintain periodicity in longitude |
| 311 |
|
✗ |
do l=1,llm |
| 312 |
|
✗ |
do ij=1,ip1jmp1,iip1 |
| 313 |
|
✗ |
teta(ij+iim,l)=teta(ij,l) |
| 314 |
|
|
enddo |
| 315 |
|
|
enddo |
| 316 |
|
|
|
| 317 |
|
|
ENDIF ! of IF (.NOT. read_start) |
| 318 |
|
|
endif academic_case |
| 319 |
|
|
|
| 320 |
|
✗ |
END SUBROUTINE iniacademic |
| 321 |
|
|
|