- Timestamp:
- Nov 21, 2019, 4:43:45 PM (4 years ago)
- Location:
- LMDZ6/branches/Ocean_skin
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/dyn1d/1D_read_forc_cases.h (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Ocean_skin
-
LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1D_read_forc_cases.h
r2920 r3605 11 11 nq2=0 12 12 13 if (forcing_les .or. forcing_radconv &14 & .or. forcing_GCSSold .or. forcing_fire) then13 print*,'FORCING ,forcing_SCM',forcing_SCM 14 if (forcing_SCM) then 15 15 16 if (forcing_fire) then 17 !---------------------------------------------------------------------- 18 !read fire forcings from fire.nc 19 !---------------------------------------------------------------------- 20 fich_fire='fire.nc' 21 call read_fire(fich_fire,nlev_fire,nt_fire & 22 & ,height,tttprof,qtprof,uprof,vprof,e12prof & 23 & ,ugprof,vgprof,wfls,dqtdxls & 24 & ,dqtdyls,dqtdtls,thlpcar) 25 write(*,*) 'Forcing FIRE lu' 26 kmax=120 ! nombre de niveaux dans les profils et forcages 27 else 28 !---------------------------------------------------------------------- 29 ! Read profiles from files: prof.inp.001 and lscale.inp.001 30 ! (repris de readlesfiles) 31 !---------------------------------------------------------------------- 32 33 call readprofiles(nlev_max,kmax,nqtot,height, & 34 & tttprof,qtprof,uprof,vprof, & 35 & e12prof,ugprof,vgprof, & 36 & wfls,dqtdxls,dqtdyls,dqtdtls, & 37 & thlpcar,qprof,nq1,nq2) 38 endif 39 40 ! compute altitudes of play levels. 41 zlay(1) =zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) 42 do l = 2,llm 43 zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf) 44 enddo 45 46 !---------------------------------------------------------------------- 47 ! Interpolation of the profiles given on the input file to 48 ! model levels 49 !---------------------------------------------------------------------- 50 zlay(1) = zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) 51 do l=1,llm 52 ! Above the max altutide of the input file 53 54 if (zlay(l)<height(kmax)) mxcalc=l 55 56 frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1)) 57 ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1)) 58 if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile 59 temp(l) = ttt*(play(l)/pzero)**rkappa 60 teta(l) = ttt 61 else 62 temp(l) = ttt 63 teta(l) = ttt*(pzero/play(l))**rkappa 64 endif 65 print *,' temp,teta ',l,temp(l),teta(l) 66 q(l,1) = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1)) 67 u(l) = uprof(kmax)-frac*( uprof(kmax)- uprof(kmax-1)) 68 v(l) = vprof(kmax)-frac*( vprof(kmax)- vprof(kmax-1)) 69 ug(l) = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1)) 70 vg(l) = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1)) 71 IF (nq2>0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2) & 72 & -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2)) 73 omega(l)= wfls(kmax)-frac*( wfls(kmax)- wfls(kmax-1)) 74 75 dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1)) 76 dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1)) 77 do k=2,kmax 78 print *,'k l height(k) height(k-1) zlay(l) frac=',k,l,height(k),height(k-1),zlay(l),frac 79 frac = (height(k)-zlay(l))/(height(k)-height(k-1)) 80 if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k) 81 if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then 82 ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1)) 83 if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile 84 temp(l) = ttt*(play(l)/pzero)**rkappa 85 teta(l) = ttt 86 else 87 temp(l) = ttt 88 teta(l) = ttt*(pzero/play(l))**rkappa 89 endif 90 print *,' temp,teta ',l,temp(l),teta(l) 91 q(l,1) = qtprof(k)-frac*( qtprof(k)- qtprof(k-1)) 92 u(l) = uprof(k)-frac*( uprof(k)- uprof(k-1)) 93 v(l) = vprof(k)-frac*( vprof(k)- vprof(k-1)) 94 ug(l) = ugprof(k)-frac*( ugprof(k)- ugprof(k-1)) 95 vg(l) = vgprof(k)-frac*( vgprof(k)- vgprof(k-1)) 96 IF (nq2>0) q(l,nq1:nq2)=qprof(k,nq1:nq2) & 97 & -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2)) 98 omega(l)= wfls(k)-frac*( wfls(k)- wfls(k-1)) 99 dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1)) 100 dt_cooling(l)=thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1)) 101 elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1) 102 ttt =tttprof(1) 103 if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile 104 temp(l) = ttt*(play(l)/pzero)**rkappa 105 teta(l) = ttt 106 else 107 temp(l) = ttt 108 teta(l) = ttt*(pzero/play(l))**rkappa 109 endif 110 q(l,1) = qtprof(1) 111 u(l) = uprof(1) 112 v(l) = vprof(1) 113 ug(l) = ugprof(1) 114 vg(l) = vgprof(1) 115 omega(l)= wfls(1) 116 IF (nq2>0) q(l,nq1:nq2)=qprof(1,nq1:nq2) 117 dq_dyn(l,1) =dqtdtls(1) 118 dt_cooling(l)=thlpcar(1) 119 endif 120 enddo 121 122 temp(l)=max(min(temp(l),350.),150.) 123 rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 124 if (l .lt. llm) then 125 zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l)) 126 endif 127 omega2(l)=-rho(l)*omega(l) 128 omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s 129 if (l>1) then 130 if(zlay(l-1)>height(kmax)) then 131 omega(l)=0.0 132 omega2(l)=0.0 133 endif 134 endif 135 if(q(l,1)<0.) q(l,1)=0.0 136 q(l,2) = 0.0 137 enddo 138 139 endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire 140 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 141 !--------------------------------------------------------------------- 142 ! Forcing for GCSSold: 143 !--------------------------------------------------------------------- 144 if (forcing_GCSSold) then 145 fich_gcssold_ctl = './forcing.ctl' 146 fich_gcssold_dat = './forcing8.dat' 147 call copie(llm,play,psurf,fich_gcssold_ctl) 148 call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat, & 149 & ht_gcssold,hq_gcssold,hw_gcssold, & 150 & hu_gcssold,hv_gcssold, & 151 & hthturb_gcssold,hqturb_gcssold,Ts_gcssold, & 152 & imp_fcg_gcssold,ts_fcg_gcssold, & 153 & Tp_fcg_gcssold,Turb_fcg_gcssold) 154 print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold 155 endif ! forcing_GCSSold 156 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 157 !--------------------------------------------------------------------- 158 ! Forcing for RICO: 159 !--------------------------------------------------------------------- 160 if (forcing_rico) then 161 162 ! call writefield_phy('omega', omega,llm+1) 163 fich_rico = 'rico.txt' 164 call read_rico(fich_rico,nlev_rico,ps_rico,play & 165 & ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico & 166 & ,dth_rico,dqh_rico) 167 print*, ' on a lu et prepare RICO' 168 169 mxcalc=llm 170 print *, airefi, ' airefi ' 171 do l = 1, llm 172 rho(l) = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l))) 173 temp(l) = t_rico(l) 174 q(l,1) = q_rico(l) 175 q(l,2) = 0.0 176 u(l) = u_rico(l) 177 v(l) = v_rico(l) 178 ug(l)=u_rico(l) 179 vg(l)=v_rico(l) 180 omega(l) = -w_rico(l)*rg 181 omega2(l) = omega(l)/rg*airefi 182 enddo 183 endif 184 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 185 !--------------------------------------------------------------------- 186 ! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) : 187 !--------------------------------------------------------------------- 188 189 if (forcing_toga) then 190 191 ! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps): 192 fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt' 193 CALL read_togacoare(fich_toga,nlev_toga,nt_toga & 194 & ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga & 195 & ,ht_toga,vt_toga,hq_toga,vq_toga) 196 197 write(*,*) 'Forcing TOGA lu' 198 199 ! time interpolation for initial conditions: 200 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 201 CALL interp_toga_time(daytime,day1,annee_ref & 202 & ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga & 203 & ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga & 204 & ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga & 205 & ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof & 206 & ,ht_prof,vt_prof,hq_prof,vq_prof) 207 208 ! vertical interpolation: 209 CALL interp_toga_vertical(play,nlev_toga,plev_prof & 210 & ,t_prof,q_prof,u_prof,v_prof,w_prof & 211 & ,ht_prof,vt_prof,hq_prof,vq_prof & 212 & ,t_mod,q_mod,u_mod,v_mod,w_mod & 213 & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) 214 write(*,*) 'Profil initial forcing TOGA interpole' 215 216 ! initial and boundary conditions : 217 tsurf = ts_prof 218 write(*,*) 'SST initiale: ',tsurf 219 do l = 1, llm 220 temp(l) = t_mod(l) 221 q(l,1) = q_mod(l) 222 q(l,2) = 0.0 223 u(l) = u_mod(l) 224 v(l) = v_mod(l) 225 omega(l) = w_mod(l) 226 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 227 !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 228 !? omega2(l)=-rho(l)*omega(l) 229 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 230 d_t_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l)) 231 d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l)) 232 d_q_adv(l,2) = 0.0 233 enddo 234 235 endif ! forcing_toga 236 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 237 !--------------------------------------------------------------------- 238 ! Forcing from TWPICE experiment (Shaocheng et al. 2010) : 239 !--------------------------------------------------------------------- 240 241 if (forcing_twpice) then 242 !read TWP-ICE forcings 243 fich_twpice='d_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf' 244 call read_twpice(fich_twpice,nlev_twpi,nt_twpi & 245 & ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi & 246 & ,ht_twpi,vt_twpi,hq_twpi,vq_twpi) 247 248 write(*,*) 'Forcing TWP-ICE lu' 249 !Time interpolation for initial conditions using TOGA interpolation routine 250 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 251 CALL interp_toga_time(daytime,day1,annee_ref & 252 & ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi & 253 & ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi & 254 & ,ht_twpi,vt_twpi,hq_twpi,vq_twpi & 255 & ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp & 256 & ,u_proftwp,v_proftwp,w_proftwp & 257 & ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp) 258 259 ! vertical interpolation using TOGA interpolation routine: 260 ! write(*,*)'avant interp vert', t_proftwp 261 CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp & 262 & ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp & 263 & ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp & 264 & ,t_mod,q_mod,u_mod,v_mod,w_mod & 265 & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) 266 ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod 267 268 ! initial and boundary conditions : 269 ! tsurf = ts_proftwp 270 write(*,*) 'SST initiale: ',tsurf 271 do l = 1, llm 272 temp(l) = t_mod(l) 273 q(l,1) = q_mod(l) 274 q(l,2) = 0.0 275 u(l) = u_mod(l) 276 v(l) = v_mod(l) 277 omega(l) = w_mod(l) 278 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 279 280 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 281 !on applique le forcage total au premier pas de temps 282 !attention: signe different de toga 283 d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l)) 284 d_q_adv(l,1) = (hq_mod(l)+vq_mod(l)) 285 d_q_adv(l,2) = 0.0 286 enddo 287 288 endif !forcing_twpice 289 290 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 291 !--------------------------------------------------------------------- 292 ! Forcing from AMMA experiment (Couvreux et al. 2010) : 293 !--------------------------------------------------------------------- 294 295 if (forcing_amma) then 296 297 call read_1D_cases 298 299 write(*,*) 'Forcing AMMA lu' 300 301 !champs initiaux: 302 do k=1,nlev_amma 303 th_ammai(k)=th_amma(k) 304 q_ammai(k)=q_amma(k) 305 u_ammai(k)=u_amma(k) 306 v_ammai(k)=v_amma(k) 307 vitw_ammai(k)=vitw_amma(k,12) 308 ht_ammai(k)=ht_amma(k,12) 309 hq_ammai(k)=hq_amma(k,12) 310 vt_ammai(k)=0. 311 vq_ammai(k)=0. 312 enddo 313 omega(:)=0. 314 omega2(:)=0. 315 rho(:)=0. 316 ! vertical interpolation using TOGA interpolation routine: 317 ! write(*,*)'avant interp vert', t_proftwp 318 CALL interp_toga_vertical(play,nlev_amma,plev_amma & 319 & ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai & 320 & ,ht_ammai,vt_ammai,hq_ammai,vq_ammai & 321 & ,t_mod,q_mod,u_mod,v_mod,w_mod & 322 & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) 323 ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod 324 325 ! initial and boundary conditions : 326 ! tsurf = ts_proftwp 327 write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc 328 do l = 1, llm 329 ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp 330 ! temp(l) = t_mod(l)*(play(l)/pzero)**rkappa 331 temp(l) = t_mod(l) 332 q(l,1) = q_mod(l) 333 q(l,2) = 0.0 334 ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) 335 u(l) = u_mod(l) 336 v(l) = v_mod(l) 337 rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 338 omega(l) = w_mod(l)*(-rg*rho(l)) 339 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 340 341 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 342 !on applique le forcage total au premier pas de temps 343 !attention: signe different de toga 344 d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 345 !forcage en th 346 ! d_t_adv(l) = ht_mod(l) 347 d_q_adv(l,1) = hq_mod(l) 348 d_q_adv(l,2) = 0.0 349 dt_cooling(l)=0. 350 enddo 351 write(*,*) 'Prof initeforcing AMMA interpole temp39',temp(39) 352 353 354 ! ok_flux_surf=.false. 355 fsens=-1.*sens_amma(12) 356 flat=-1.*lat_amma(12) 357 358 endif !forcing_amma 359 360 361 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 362 !--------------------------------------------------------------------- 363 ! Forcing from DICE experiment (see file DICE_protocol_vn2-3.pdf) 364 !--------------------------------------------------------------------- 365 366 if (forcing_dice) then 367 !read DICE forcings 368 fich_dice='dice_driver.nc' 369 call read_dice(fich_dice,nlev_dice,nt_dice & 370 & ,zz_dice,plev_dice,t_dice,qv_dice,u_dice,v_dice,o3_dice & 371 & ,shf_dice,lhf_dice,lwup_dice,swup_dice,tg_dice,ustar_dice& 372 & ,psurf_dice,ug_dice,vg_dice,ht_dice,hq_dice & 373 & ,hu_dice,hv_dice,w_dice,omega_dice) 374 375 write(*,*) 'Forcing DICE lu' 376 377 !champs initiaux: 378 do k=1,nlev_dice 379 t_dicei(k)=t_dice(k) 380 qv_dicei(k)=qv_dice(k) 381 u_dicei(k)=u_dice(k) 382 v_dicei(k)=v_dice(k) 383 o3_dicei(k)=o3_dice(k) 384 ht_dicei(k)=ht_dice(k,1) 385 hq_dicei(k)=hq_dice(k,1) 386 hu_dicei(k)=hu_dice(k,1) 387 hv_dicei(k)=hv_dice(k,1) 388 w_dicei(k)=w_dice(k,1) 389 omega_dicei(k)=omega_dice(k,1) 390 enddo 391 omega(:)=0. 392 omega2(:)=0. 393 rho(:)=0. 394 ! vertical interpolation using TOGA interpolation routine: 395 ! write(*,*)'avant interp vert', t_proftwp 396 ! 397 ! CALL interp_dice_time(daytime,day1,annee_ref 398 ! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice 399 ! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice 400 ! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice 401 ! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice 402 ! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof 403 ! o ,ustar_prof,psurf_prof,ug_profd,vg_profd 404 ! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd 405 ! o ,omega_profd) 406 407 CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice & 408 & ,t_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei & 409 & ,ht_dicei,hq_dicei,hu_dicei,hv_dicei,w_dicei,omega_dicei& 410 & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & 411 & ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc) 412 413 ! Pour tester les advections horizontales de T et Q, on met w_mod et omega_mod ?? zero (MPL 20131108) 414 ! w_mod(:,:)=0. 415 ! omega_mod(:,:)=0. 416 417 ! write(*,*) 'Profil initial forcing DICE interpole',t_mod 418 ! Les forcages DICE sont donnes /jour et non /seconde ! 419 ht_mod(:)=ht_mod(:)/86400. 420 hq_mod(:)=hq_mod(:)/86400. 421 hu_mod(:)=hu_mod(:)/86400. 422 hv_mod(:)=hv_mod(:)/86400. 423 424 ! initial and boundary conditions : 425 write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc 426 do l = 1, llm 427 ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp 428 ! temp(l) = th_mod(l)*(play(l)/pzero)**rkappa 429 temp(l) = t_mod(l) 430 q(l,1) = qv_mod(l) 431 q(l,2) = 0.0 432 ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) 433 u(l) = u_mod(l) 434 v(l) = v_mod(l) 435 ug(l)=ug_dice(1) 436 vg(l)=vg_dice(1) 437 rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 438 ! omega(l) = w_mod(l)*(-rg*rho(l)) 439 omega(l) = omega_mod(l) 440 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 441 442 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 443 !on applique le forcage total au premier pas de temps 444 !attention: signe different de toga 445 d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 446 !forcage en th 447 ! d_t_adv(l) = ht_mod(l) 448 d_q_adv(l,1) = hq_mod(l) 449 d_q_adv(l,2) = 0.0 450 dt_cooling(l)=0. 451 enddo 452 write(*,*) 'Profil initial forcing DICE interpole temp39',temp(39) 453 454 455 ! ok_flux_surf=.false. 456 fsens=-1.*shf_dice(1) 457 flat=-1.*lhf_dice(1) 458 ! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par 459 ! le coefficient de trainee en surface cd**2=ustar*vent(k=1) 460 ! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface 461 ! MPL 05082013 462 ust=ustar_dice(1) 463 tg=tg_dice(1) 464 print *,'ust= ',ust 465 IF (tsurf .LE. 0.) THEN 466 tsurf= tg_dice(1) 467 ENDIF 468 psurf= psurf_dice(1) 469 solsw_in = (1.-albedo)/albedo*swup_dice(1) 470 sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1) 471 PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in 472 endif !forcing_dice 473 474 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 475 !--------------------------------------------------------------------- 476 ! Forcing from GABLS4 experiment 477 !--------------------------------------------------------------------- 478 479 !!!! Si la temperature de surface n'est pas impos??e: 480 481 if (forcing_gabls4) then 482 !read GABLS4 forcings 483 484 fich_gabls4='gabls4_driver.nc' 485 486 487 call read_gabls4(fich_gabls4,nlev_gabls4,nt_gabls4,nsol_gabls4,zz_gabls4,depth_sn_gabls4,ug_gabls4,vg_gabls4 & 488 & ,plev_gabls4,th_gabls4,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,ht_gabls4,hq_gabls4,tg_gabls4,tsnow_gabls4,snow_dens_gabls4) 489 490 write(*,*) 'Forcing GABLS4 lu' 491 492 !champs initiaux: 493 do k=1,nlev_gabls4 494 t_gabi(k)=t_gabls4(k) 495 qv_gabi(k)=qv_gabls4(k) 496 u_gabi(k)=u_gabls4(k) 497 v_gabi(k)=v_gabls4(k) 498 poub(k)=0. 499 ht_gabi(k)=ht_gabls4(k,1) 500 hq_gabi(k)=hq_gabls4(k,1) 501 ug_gabi(k)=ug_gabls4(k,1) 502 vg_gabi(k)=vg_gabls4(k,1) 503 enddo 504 505 omega(:)=0. 506 omega2(:)=0. 507 rho(:)=0. 508 ! vertical interpolation using TOGA interpolation routine: 509 ! write(*,*)'avant interp vert', t_proftwp 510 ! 511 ! CALL interp_dice_time(daytime,day1,annee_ref 512 ! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice 513 ! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice 514 ! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice 515 ! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice 516 ! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof 517 ! o ,ustar_prof,psurf_prof,ug_profd,vg_profd 518 ! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd 519 ! o ,omega_profd) 520 521 CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4 & 522 & ,t_gabi,qv_gabi,u_gabi,v_gabi,poub & 523 & ,ht_gabi,hq_gabi,ug_gabi,vg_gabi,poub,poub & 524 & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & 525 & ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc) 526 527 ! Les forcages GABLS4 ont l air d etre en K/S quoiqu en dise le fichier gabls4_driver.nc !? MPL 20141024 528 ! ht_mod(:)=ht_mod(:)/86400. 529 ! hq_mod(:)=hq_mod(:)/86400. 530 531 ! initial and boundary conditions : 532 write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc 533 do l = 1, llm 534 ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp 535 ! temp(l) = th_mod(l)*(play(l)/pzero)**rkappa 536 temp(l) = t_mod(l) 537 q(l,1) = qv_mod(l) 538 q(l,2) = 0.0 539 ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) 540 u(l) = u_mod(l) 541 v(l) = v_mod(l) 542 ug(l)=ug_mod(l) 543 vg(l)=vg_mod(l) 544 545 ! 546 ! tg=tsurf 547 ! 548 549 print *,'***** tsurf=',tsurf 550 rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 551 ! omega(l) = w_mod(l)*(-rg*rho(l)) 552 omega(l) = omega_mod(l) 553 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 554 555 556 557 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 558 !on applique le forcage total au premier pas de temps 559 !attention: signe different de toga 560 ! d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 561 !forcage en th 562 d_t_adv(l) = ht_mod(l) 563 d_q_adv(l,1) = hq_mod(l) 564 d_q_adv(l,2) = 0.0 565 dt_cooling(l)=0. 566 enddo 567 568 !--------------- Residus forcages du cas Dice (a supprimer) MPL 20141024--------------- 569 ! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par 570 ! le coefficient de trainee en surface cd**2=ustar*vent(k=1) 571 ! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface 572 ! MPL 05082013 573 ! ust=ustar_dice(1) 574 ! tg=tg_dice(1) 575 ! print *,'ust= ',ust 576 ! IF (tsurf .LE. 0.) THEN 577 ! tsurf= tg_dice(1) 578 ! ENDIF 579 ! psurf= psurf_dice(1) 580 ! solsw_in = (1.-albedo)/albedo*swup_dice(1) 581 ! sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1) 582 ! PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in 583 !-------------------------------------------------------------------------------------- 584 endif !forcing_gabls4 585 586 587 588 ! Forcing from Arm_Cu case 589 ! For this case, ifa_armcu.txt contains sensible, latent heat fluxes 590 ! large scale advective forcing,radiative forcing 591 ! and advective tendency of theta and qt to be applied 592 !--------------------------------------------------------------------- 593 594 if (forcing_armcu) then 595 ! read armcu forcing : 596 write(*,*) 'Avant lecture Forcing Arm_Cu' 597 fich_armcu = './ifa_armcu.txt' 598 CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu, & 599 & sens_armcu,flat_armcu,adv_theta_armcu, & 600 & rad_theta_armcu,adv_qt_armcu) 601 write(*,*) 'Forcing Arm_Cu lu' 602 603 !---------------------------------------------------------------------- 604 ! Read profiles from file: prof.inp.19 or prof.inp.40 605 ! For this case, profiles are given for two vertical resolution 606 ! 19 or 40 levels 607 ! 608 ! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html 609 ! Note that the initial profiles contain no liquid water! 610 ! (so potential temperature can be interpreted as liquid water 611 ! potential temperature and water vapor as total water) 612 ! profiles are given at full levels 613 !---------------------------------------------------------------------- 614 615 call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod, & 616 & v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp) 617 618 ! time interpolation for initial conditions: 619 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 620 621 print *,'Avant interp_armcu_time' 622 print *,'daytime=',daytime 623 print *,'day1=',day1 624 print *,'annee_ref=',annee_ref 625 print *,'year_ini_armcu=',year_ini_armcu 626 print *,'day_ju_ini_armcu=',day_ju_ini_armcu 627 print *,'nt_armcu=',nt_armcu 628 print *,'dt_armcu=',dt_armcu 629 print *,'nlev_armcu=',nlev_armcu 630 CALL interp_armcu_time(daytime,day1,annee_ref & 631 & ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu & 632 & ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu & 633 & ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof & 634 & ,adv_theta_prof,rad_theta_prof,adv_qt_prof) 635 write(*,*) 'Forcages interpoles dans temps' 636 637 ! No vertical interpolation if nlev imposed to 19 or 40 638 ! The vertical grid stops at 4000m # 600hPa 639 mxcalc=llm 640 641 ! initial and boundary conditions : 642 ! tsurf = ts_prof 643 ! tsurf read in lmdz1d.def 644 write(*,*) 'Tsurf initiale: ',tsurf 645 do l = 1, llm 646 play(l)=play_mod(l)*100. 647 presnivs(l)=play(l) 648 zlay(l)=height(l) 649 temp(l) = t_mod(l) 650 teta(l)=theta_mod(l) 651 q(l,1) = qv_mod(l)/1000. 652 ! No liquid water in the initial profil 653 q(l,2) = 0. 654 u(l) = u_mod(l) 655 ug(l)= u_mod(l) 656 v(l) = v_mod(l) 657 vg(l)= v_mod(l) 658 ! Advective forcings are given in K or g/kg ... per HOUR 659 ! IF(height(l).LT.1000) THEN 660 ! d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600. 661 ! d_q_adv(l,1) = adv_qt_prof/1000./3600. 662 ! d_q_adv(l,2) = 0.0 663 ! ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN 664 ! d_t_adv(l) = (adv_theta_prof + rad_theta_prof)* 665 ! : (1-(height(l)-1000.)/2000.) 666 ! d_t_adv(l) = d_t_adv(l)/3600. 667 ! d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.) 668 ! d_q_adv(l,1) = d_q_adv(l,1)/1000./3600. 669 ! d_q_adv(l,2) = 0.0 670 ! ELSE 671 ! d_t_adv(l) = 0.0 672 ! d_q_adv(l,1) = 0.0 673 ! d_q_adv(l,2) = 0.0 674 ! ENDIF 675 enddo 676 ! plev at half levels is given in proh.inp.19 or proh.inp.40 files 677 plev(1)= ap(llm+1)+bp(llm+1)*psurf 678 do l = 1, llm 679 plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf 680 print *,'Read_forc: l height play plev zlay temp', & 681 & l,height(l),play(l),plev(l),zlay(l),temp(l) 682 enddo 683 ! For this case, fluxes are imposed 684 fsens=-1*sens_prof 685 flat=-1*flat_prof 686 687 endif ! forcing_armcu 688 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 689 !--------------------------------------------------------------------- 690 ! Forcing from transition case of Irina Sandu 691 !--------------------------------------------------------------------- 692 693 if (forcing_sandu) then 694 write(*,*) 'Avant lecture Forcing SANDU' 695 696 ! read sanduref forcing : 697 fich_sandu = './ifa_sanduref.txt' 698 CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu) 699 700 write(*,*) 'Forcing SANDU lu' 701 702 !---------------------------------------------------------------------- 703 ! Read profiles from file: prof.inp.001 704 !---------------------------------------------------------------------- 705 706 call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs, & 707 & thl_profs,q_profs,u_profs,v_profs, & 708 & w_profs,omega_profs,o3mmr_profs) 709 710 ! time interpolation for initial conditions: 711 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 712 ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! 713 ! revoir 1DUTILS.h et les arguments 714 715 print *,'Avant interp_sandu_time' 716 print *,'daytime=',daytime 717 print *,'day1=',day1 718 print *,'annee_ref=',annee_ref 719 print *,'year_ini_sandu=',year_ini_sandu 720 print *,'day_ju_ini_sandu=',day_ju_ini_sandu 721 print *,'nt_sandu=',nt_sandu 722 print *,'dt_sandu=',dt_sandu 723 print *,'nlev_sandu=',nlev_sandu 724 CALL interp_sandu_time(daytime,day1,annee_ref & 725 & ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu & 726 & ,nlev_sandu & 727 & ,ts_sandu,ts_prof) 728 729 ! vertical interpolation: 730 print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu 731 CALL interp_sandu_vertical(play,nlev_sandu,plev_profs & 732 & ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs & 733 & ,omega_profs,o3mmr_profs & 734 & ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod & 735 & ,omega_mod,o3mmr_mod,mxcalc) 736 write(*,*) 'Profil initial forcing SANDU interpole' 737 738 ! initial and boundary conditions : 739 tsurf = ts_prof 740 write(*,*) 'SST initiale: ',tsurf 741 do l = 1, llm 742 temp(l) = t_mod(l) 743 tetal(l)=thl_mod(l) 744 q(l,1) = q_mod(l) 745 q(l,2) = 0.0 746 u(l) = u_mod(l) 747 v(l) = v_mod(l) 748 w(l) = w_mod(l) 749 omega(l) = omega_mod(l) 750 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 751 !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 752 !? omega2(l)=-rho(l)*omega(l) 753 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 754 ! d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) 755 ! d_q_adv(l,1) = vq_mod(l) 756 d_t_adv(l) = alpha*omega(l)/rcpd 757 d_q_adv(l,1) = 0.0 758 d_q_adv(l,2) = 0.0 759 enddo 760 761 endif ! forcing_sandu 762 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 763 !--------------------------------------------------------------------- 764 ! Forcing from Astex case 765 !--------------------------------------------------------------------- 766 767 if (forcing_astex) then 768 write(*,*) 'Avant lecture Forcing Astex' 769 770 ! read astex forcing : 771 fich_astex = './ifa_astex.txt' 772 CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex, & 773 & ug_astex,vg_astex,ufa_astex,vfa_astex) 774 775 write(*,*) 'Forcing Astex lu' 776 777 !---------------------------------------------------------------------- 778 ! Read profiles from file: prof.inp.001 779 !---------------------------------------------------------------------- 780 781 call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa, & 782 & thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa, & 783 & w_profa,tke_profa,o3mmr_profa) 784 785 ! time interpolation for initial conditions: 786 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 787 ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! 788 ! revoir 1DUTILS.h et les arguments 789 790 print *,'Avant interp_astex_time' 791 print *,'daytime=',daytime 792 print *,'day1=',day1 793 print *,'annee_ref=',annee_ref 794 print *,'year_ini_astex=',year_ini_astex 795 print *,'day_ju_ini_astex=',day_ju_ini_astex 796 print *,'nt_astex=',nt_astex 797 print *,'dt_astex=',dt_astex 798 print *,'nlev_astex=',nlev_astex 799 CALL interp_astex_time(daytime,day1,annee_ref & 800 & ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex & 801 & ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex & 802 & ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof & 803 & ,ufa_prof,vfa_prof) 804 805 ! vertical interpolation: 806 print *,'Avant interp_vertical: nlev_astex=',nlev_astex 807 CALL interp_astex_vertical(play,nlev_astex,plev_profa & 808 & ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa & 809 & ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa & 810 & ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod & 811 & ,tke_mod,o3mmr_mod,mxcalc) 812 write(*,*) 'Profil initial forcing Astex interpole' 813 814 ! initial and boundary conditions : 815 tsurf = ts_prof 816 write(*,*) 'SST initiale: ',tsurf 817 do l = 1, llm 818 temp(l) = t_mod(l) 819 tetal(l)=thl_mod(l) 820 q(l,1) = qv_mod(l) 821 q(l,2) = ql_mod(l) 822 u(l) = u_mod(l) 823 v(l) = v_mod(l) 824 w(l) = w_mod(l) 825 omega(l) = w_mod(l) 826 ! omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 827 ! rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 828 ! omega2(l)=-rho(l)*omega(l) 829 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 830 ! d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) 831 ! d_q_adv(l,1) = vq_mod(l) 832 d_t_adv(l) = alpha*omega(l)/rcpd 833 d_q_adv(l,1) = 0.0 834 d_q_adv(l,2) = 0.0 835 enddo 836 837 endif ! forcing_astex 838 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 839 !--------------------------------------------------------------------- 840 ! Forcing from standard case : 841 !--------------------------------------------------------------------- 842 843 if (forcing_case) then 844 845 write(*,*),'avant call read_1D_cas' 846 call read_1D_cas 847 write(*,*) 'Forcing read' 848 849 !Time interpolation for initial conditions using TOGA interpolation routine 850 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 851 CALL interp_case_time(day,day1,annee_ref & 852 ! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas & 853 & ,nt_cas,nlev_cas & 854 & ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas & 855 & ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas & 856 & ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas & 857 & ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas & 858 & ,uw_cas,vw_cas,q1_cas,q2_cas & 859 & ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas & 860 & ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas & 861 & ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas & 862 & ,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas,ustar_prof_cas & 863 & ,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas) 864 865 ! vertical interpolation using TOGA interpolation routine: 866 ! write(*,*)'avant interp vert', t_prof 867 CALL interp_case_vertical(play,nlev_cas,plev_prof_cas & 868 & ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas & 869 & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas & 870 & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & 871 & ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas & 872 & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas & 873 & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc) 874 ! write(*,*) 'Profil initial forcing case interpole',t_mod 875 876 ! initial and boundary conditions : 877 ! tsurf = ts_prof_cas 878 ts_cur = ts_prof_cas 879 psurf=plev_prof_cas(1) 880 write(*,*) 'SST initiale: ',tsurf 881 do l = 1, llm 882 temp(l) = t_mod_cas(l) 883 q(l,1) = q_mod_cas(l) 884 q(l,2) = 0.0 885 u(l) = u_mod_cas(l) 886 v(l) = v_mod_cas(l) 887 omega(l) = w_mod_cas(l) 888 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 889 890 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 891 !on applique le forcage total au premier pas de temps 892 !attention: signe different de toga 893 d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) 894 d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l)) 895 d_q_adv(l,2) = 0.0 896 d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l)) 897 ! correction bug d_u -> d_v (MM+MPL 20170310) 898 d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l)) 899 enddo 900 901 ! In case fluxes are imposed 902 IF (ok_flux_surf) THEN 903 fsens=sens_prof_cas 904 flat=lat_prof_cas 905 ENDIF 906 IF (ok_prescr_ust) THEN 907 ust=ustar_prof_cas 908 print *,'ust=',ust 909 ENDIF 910 911 endif !forcing_case 912 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 913 !--------------------------------------------------------------------- 914 ! Forcing from standard case : 915 !--------------------------------------------------------------------- 916 917 if (forcing_case2) then 918 919 write(*,*),'avant call read2_1D_cas' 920 call read2_1D_cas 921 write(*,*) 'Forcing read' 16 write(*,*),'avant call read_SCM' 17 call read_SCM_cas 18 write(*,*) 'Forcing read' 19 print*,'PS ps_cas',ps_cas 922 20 923 21 !Time interpolation for initial conditions using interpolation routine 924 22 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 925 CALL interp 2_case_time(daytime,day1,annee_ref &23 CALL interp_case_time_std(daytime,day1,annee_ref & 926 24 ! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas & 927 25 & ,nt_cas,nlev_cas & 928 26 & ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas & 929 & ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas & 27 & ,u_cas,v_cas,ug_cas,vg_cas & 28 & ,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas & 29 & ,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas & 930 30 & ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas & 931 31 & ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas & 932 32 & ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas & 933 33 ! 934 & ,ts_prof_cas,p lev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas &34 & ,ts_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas & 935 35 & ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas & 936 & ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & 36 & ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas & 37 & ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas & 38 & ,vitw_prof_cas,omega_prof_cas & 937 39 & ,du_prof_cas,hu_prof_cas,vu_prof_cas & 938 40 & ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas & … … 947 49 ! vertical interpolation using interpolation routine: 948 50 ! write(*,*)'avant interp vert', t_prof 949 CALL interp2_case_vertical (play,nlev_cas,plev_prof_cas &51 CALL interp2_case_vertical_std(play,nlev_cas,plev_prof_cas & 950 52 & ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas & 951 53 & ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas & 952 & ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & 54 & ,ug_prof_cas,vg_prof_cas & 55 & ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas & 56 57 & ,vitw_prof_cas,omega_prof_cas & 953 58 & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas & 954 59 & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & … … 956 61 ! 957 62 & ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas & 958 & ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas & 63 & ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas & 64 & ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas & 65 & ,w_mod_cas,omega_mod_cas & 959 66 & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas & 960 67 & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas & 961 68 & ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc) 962 69 963 ! write(*,*) 'Profil initial forcing case interpole',t_mod964 70 965 71 ! initial and boundary conditions : 966 72 ! tsurf = ts_prof_cas 73 psurf = ps_prof_cas 967 74 ts_cur = ts_prof_cas 968 psurf=plev_prof_cas(1)969 write(*,*) 'SST initiale: ',tsurf970 75 do l = 1, llm 971 76 temp(l) = t_mod_cas(l) … … 980 85 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 981 86 982 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 983 ! on applique le forcage total au premier pas de temps984 ! attention: signe different de toga985 d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) 986 d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) 987 ! d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))988 d_q_adv(l,1) = dq_mod_cas(l) 87 88 ! On effectue la somme du forcage total et de la decomposition 89 ! horizontal/vertical en supposant que soit l'un soit l'autre 90 ! sont remplis mais jamais les deux 91 92 d_t_adv(l) = dt_mod_cas(l)+ht_mod_cas(l)+vt_mod_cas(l) 93 d_q_adv(l,1) = dq_mod_cas(l)+hq_mod_cas(l)+vq_mod_cas(l) 989 94 d_q_adv(l,2) = 0.0 990 ! d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))991 d_ u_adv(l) = du_mod_cas(l)992 ! d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l)) 993 ! correction bug d_u -> d_v (MM+MPL 20170310)994 d_v_adv(l) = dv_mod_cas(l) 95 d_u_adv(l) = du_mod_cas(l)+hu_mod_cas(l)+vu_mod_cas(l) 96 d_v_adv(l) = dv_mod_cas(l)+hv_mod_cas(l)+vv_mod_cas(l) 97 98 !print*,'d_t_adv ',d_t_adv(1:20)*86400 99 995 100 enddo 996 101 … … 1006 111 ENDIF 1007 112 1008 endif !forcing_case2 1009 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1010 113 endif !forcing_SCM
Note: See TracChangeset
for help on using the changeset viewer.
