!---------------------------------------------------------------------- ! forcing_les = .T. : Impose a constant cooling ! forcing_radconv = .T. : Pure radiative-convective equilibrium: !---------------------------------------------------------------------- nq1=0 nq2=0 if (forcing_les .or. forcing_radconv & & .or. forcing_GCSSold .or. forcing_fire) then if (forcing_fire) then !---------------------------------------------------------------------- !read fire forcings from fire.nc !---------------------------------------------------------------------- fich_fire='fire.nc' call read_fire(fich_fire,nlev_fire,nt_fire & & ,height,tttprof,qtprof,uprof,vprof,e12prof & & ,ugprof,vgprof,wfls,dqtdxls & & ,dqtdyls,dqtdtls,thlpcar) write(*,*) 'Forcing FIRE lu' kmax=120 ! nombre de niveaux dans les profils et forcages else !---------------------------------------------------------------------- ! Read profiles from files: prof.inp.001 and lscale.inp.001 ! (repris de readlesfiles) !---------------------------------------------------------------------- call readprofiles(nlev_max,kmax,nqtot,height, & & tttprof,qtprof,uprof,vprof, & & e12prof,ugprof,vgprof, & & wfls,dqtdxls,dqtdyls,dqtdtls, & & thlpcar,qprof,nq1,nq2) endif ! compute altitudes of play levels. zlay(1) =zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) do l = 2,llm zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf) enddo !---------------------------------------------------------------------- ! Interpolation of the profiles given on the input file to ! model levels !---------------------------------------------------------------------- zlay(1) = zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) do l=1,llm ! Above the max altutide of the input file if (zlay(l)0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2) & & -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2)) omega(l)= wfls(kmax)-frac*( wfls(kmax)- wfls(kmax-1)) dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1)) dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1)) do k=2,kmax frac = (height(k)-zlay(l))/(height(k)-height(k-1)) if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k) if(zlay(l)>height(k-1).and.zlay(l)0) q(l,nq1:nq2)=qprof(k,nq1:nq2) & & -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2)) omega(l)= wfls(k)-frac*( wfls(k)- wfls(k-1)) dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1)) dt_cooling(l)=thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1)) elseif(zlay(l)0) q(l,nq1:nq2)=qprof(1,nq1:nq2) dq_dyn(l,1) =dqtdtls(1) dt_cooling(l)=thlpcar(1) endif enddo temp(l)=max(min(temp(l),350.),150.) rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) if (l .lt. llm) then zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l)) endif omega2(l)=-rho(l)*omega(l) omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s if (l>1) then if(zlay(l-1)>height(kmax)) then omega(l)=0.0 omega2(l)=0.0 endif endif if(q(l,1)<0.) q(l,1)=0.0 q(l,2) = 0.0 enddo endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing for GCSSold: !--------------------------------------------------------------------- if (forcing_GCSSold) then fich_gcssold_ctl = './forcing.ctl' fich_gcssold_dat = './forcing8.dat' call copie(llm,play,psurf,fich_gcssold_ctl) call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat, & & ht_gcssold,hq_gcssold,hw_gcssold, & & hu_gcssold,hv_gcssold, & & hthturb_gcssold,hqturb_gcssold,Ts_gcssold, & & imp_fcg_gcssold,ts_fcg_gcssold, & & Tp_fcg_gcssold,Turb_fcg_gcssold) print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold endif ! forcing_GCSSold !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing for RICO: !--------------------------------------------------------------------- if (forcing_rico) then ! call writefield_phy('omega', omega,llm+1) fich_rico = 'rico.txt' call read_rico(fich_rico,nlev_rico,ps_rico,play & & ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico & & ,dth_rico,dqh_rico) print*, ' on a lu et prepare RICO' mxcalc=llm print *, airefi, ' airefi ' do l = 1, llm rho(l) = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l))) temp(l) = t_rico(l) q(l,1) = q_rico(l) q(l,2) = 0.0 u(l) = u_rico(l) v(l) = v_rico(l) ug(l)=u_rico(l) vg(l)=v_rico(l) omega(l) = -w_rico(l)*rg omega2(l) = omega(l)/rg*airefi enddo endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) : !--------------------------------------------------------------------- if (forcing_toga) then ! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps): fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt' CALL read_togacoare(fich_toga,nlev_toga,nt_toga & & ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga & & ,ht_toga,vt_toga,hq_toga,vq_toga) write(*,*) 'Forcing TOGA lu' ! time interpolation for initial conditions: write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 CALL interp_toga_time(daytime,day1,annee_ref & & ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga & & ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga & & ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga & & ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof & & ,ht_prof,vt_prof,hq_prof,vq_prof) ! vertical interpolation: CALL interp_toga_vertical(play,nlev_toga,plev_prof & & ,t_prof,q_prof,u_prof,v_prof,w_prof & & ,ht_prof,vt_prof,hq_prof,vq_prof & & ,t_mod,q_mod,u_mod,v_mod,w_mod & & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) write(*,*) 'Profil initial forcing TOGA interpole' ! initial and boundary conditions : tsurf = ts_prof write(*,*) 'SST initiale: ',tsurf do l = 1, llm temp(l) = t_mod(l) q(l,1) = q_mod(l) q(l,2) = 0.0 u(l) = u_mod(l) v(l) = v_mod(l) omega(l) = w_mod(l) omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) !? omega2(l)=-rho(l)*omega(l) alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l)) d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l)) d_q_adv(l,2) = 0.0 enddo endif ! forcing_toga !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing from TWPICE experiment (Shaocheng et al. 2010) : !--------------------------------------------------------------------- if (forcing_twpice) then !read TWP-ICE forcings fich_twpice='d_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf' call read_twpice(fich_twpice,nlev_twpi,nt_twpi & & ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi & & ,ht_twpi,vt_twpi,hq_twpi,vq_twpi) write(*,*) 'Forcing TWP-ICE lu' !Time interpolation for initial conditions using TOGA interpolation routine write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 CALL interp_toga_time(daytime,day1,annee_ref & & ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi & & ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi & & ,ht_twpi,vt_twpi,hq_twpi,vq_twpi & & ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp & & ,u_proftwp,v_proftwp,w_proftwp & & ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp) ! vertical interpolation using TOGA interpolation routine: ! write(*,*)'avant interp vert', t_proftwp CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp & & ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp & & ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp & & ,t_mod,q_mod,u_mod,v_mod,w_mod & & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod ! initial and boundary conditions : ! tsurf = ts_proftwp write(*,*) 'SST initiale: ',tsurf do l = 1, llm temp(l) = t_mod(l) q(l,1) = q_mod(l) q(l,2) = 0.0 u(l) = u_mod(l) v(l) = v_mod(l) omega(l) = w_mod(l) omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) !on applique le forcage total au premier pas de temps !attention: signe different de toga d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l)) d_q_adv(l,1) = (hq_mod(l)+vq_mod(l)) d_q_adv(l,2) = 0.0 enddo endif !forcing_twpice !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing from AMMA experiment (Couvreux et al. 2010) : !--------------------------------------------------------------------- if (forcing_amma) then call read_1D_cases write(*,*) 'Forcing AMMA lu' !champs initiaux: do k=1,nlev_amma th_ammai(k)=th_amma(k) q_ammai(k)=q_amma(k) u_ammai(k)=u_amma(k) v_ammai(k)=v_amma(k) vitw_ammai(k)=vitw_amma(k,12) ht_ammai(k)=ht_amma(k,12) hq_ammai(k)=hq_amma(k,12) vt_ammai(k)=0. vq_ammai(k)=0. enddo omega(:)=0. omega2(:)=0. rho(:)=0. ! vertical interpolation using TOGA interpolation routine: ! write(*,*)'avant interp vert', t_proftwp CALL interp_toga_vertical(play,nlev_amma,plev_amma & & ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai & & ,ht_ammai,vt_ammai,hq_ammai,vq_ammai & & ,t_mod,q_mod,u_mod,v_mod,w_mod & & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod ! initial and boundary conditions : ! tsurf = ts_proftwp write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc do l = 1, llm ! Ligne du dessous à decommenter si on lit theta au lieu de temp ! temp(l) = t_mod(l)*(play(l)/pzero)**rkappa temp(l) = t_mod(l) q(l,1) = q_mod(l) q(l,2) = 0.0 ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) u(l) = u_mod(l) v(l) = v_mod(l) rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) omega(l) = w_mod(l)*(-rg*rho(l)) omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) !on applique le forcage total au premier pas de temps !attention: signe different de toga d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) !forcage en th ! d_th_adv(l) = ht_mod(l) d_q_adv(l,1) = hq_mod(l) d_q_adv(l,2) = 0.0 dt_cooling(l)=0. enddo write(*,*) 'Prof initeforcing AMMA interpole temp39',temp(39) ! ok_flux_surf=.false. fsens=-1.*sens_amma(12) flat=-1.*lat_amma(12) endif !forcing_amma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing from DICE experiment (see file DICE_protocol_vn2-3.pdf) !--------------------------------------------------------------------- if (forcing_dice) then !read DICE forcings fich_dice='dice_driver.nc' call read_dice(fich_dice,nlev_dice,nt_dice & & ,zz_dice,plev_dice,th_dice,qv_dice,u_dice,v_dice,o3_dice & & ,shf_dice,lhf_dice,lwup_dice,swup_dice,tg_dice,ustar_dice& & ,psurf_dice,ug_dice,vg_dice,ht_dice,hq_dice & & ,hu_dice,hv_dice,w_dice,omega_dice) write(*,*) 'Forcing DICE lu' !champs initiaux: do k=1,nlev_dice th_dicei(k)=th_dice(k) qv_dicei(k)=qv_dice(k) u_dicei(k)=u_dice(k) v_dicei(k)=v_dice(k) o3_dicei(k)=o3_dice(k) ht_dicei(k)=ht_dice(k,1) hq_dicei(k)=hq_dice(k,1) hu_dicei(k)=hu_dice(k,1) hv_dicei(k)=hv_dice(k,1) w_dicei(k)=w_dice(k,1) omega_dicei(k)=omega_dice(k,1) enddo omega(:)=0. omega2(:)=0. rho(:)=0. ! vertical interpolation using TOGA interpolation routine: ! write(*,*)'avant interp vert', t_proftwp ! ! CALL interp_dice_time(daytime,day1,annee_ref ! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice ! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice ! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice ! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice ! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof ! o ,ustar_prof,psurf_prof,ug_profd,vg_profd ! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd ! o ,omega_profd) CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice & & ,th_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei & & ,ht_dicei,hq_dicei,hu_dicei,hv_dicei,w_dicei,omega_dicei& & ,th_mod,qv_mod,u_mod,v_mod,o3_mod & & ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc) ! Pour tester les advections horizontales de T et Q, on met w_mod et omega_mod à zero (MPL 20131108) ! w_mod(:,:)=0. ! omega_mod(:,:)=0. ! write(*,*) 'Profil initial forcing DICE interpole',t_mod ! Les forcages DICE sont donnes /jour et non /seconde ! ht_mod(:)=ht_mod(:)/86400. hq_mod(:)=hq_mod(:)/86400. hu_mod(:)=hu_mod(:)/86400. hv_mod(:)=hv_mod(:)/86400. ! initial and boundary conditions : write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc do l = 1, llm ! Ligne du dessous à decommenter si on lit theta au lieu de temp temp(l) = th_mod(l)*(play(l)/pzero)**rkappa ! temp(l) = t_mod(l) q(l,1) = qv_mod(l) q(l,2) = 0.0 ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) u(l) = u_mod(l) v(l) = v_mod(l) ug(l)=ug_dice(1) vg(l)=vg_dice(1) rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) ! omega(l) = w_mod(l)*(-rg*rho(l)) omega(l) = omega_mod(l) omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) !on applique le forcage total au premier pas de temps !attention: signe different de toga d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) !forcage en th ! d_th_adv(l) = ht_mod(l) d_q_adv(l,1) = hq_mod(l) d_q_adv(l,2) = 0.0 dt_cooling(l)=0. enddo write(*,*) 'Profil initial forcing DICE interpole temp39',temp(39) ! ok_flux_surf=.false. fsens=-1.*shf_dice(1) flat=-1.*lhf_dice(1) ! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par ! le coefficient de trainee en surface cd**2=ustar*vent(k=1) ! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface ! MPL 05082013 ust=ustar_dice(1) tg=tg_dice(1) print *,'ust= ',ust IF (tsurf .LE. 0.) THEN tsurf= tg_dice(1) ENDIF psurf= psurf_dice(1) solsw_in = (1.-albedo)/albedo*swup_dice(1) sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1) PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in endif !forcing_dice !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing from Arm_Cu case ! For this case, ifa_armcu.txt contains sensible, latent heat fluxes ! large scale advective forcing,radiative forcing ! and advective tendency of theta and qt to be applied !--------------------------------------------------------------------- if (forcing_armcu) then ! read armcu forcing : write(*,*) 'Avant lecture Forcing Arm_Cu' fich_armcu = './ifa_armcu.txt' CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu, & & sens_armcu,flat_armcu,adv_theta_armcu, & & rad_theta_armcu,adv_qt_armcu) write(*,*) 'Forcing Arm_Cu lu' !---------------------------------------------------------------------- ! Read profiles from file: prof.inp.19 or prof.inp.40 ! For this case, profiles are given for two vertical resolution ! 19 or 40 levels ! ! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html ! Note that the initial profiles contain no liquid water! ! (so potential temperature can be interpreted as liquid water ! potential temperature and water vapor as total water) ! profiles are given at full levels !---------------------------------------------------------------------- call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod, & & v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp) ! time interpolation for initial conditions: write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 print *,'Avant interp_armcu_time' print *,'daytime=',daytime print *,'day1=',day1 print *,'annee_ref=',annee_ref print *,'year_ini_armcu=',year_ini_armcu print *,'day_ju_ini_armcu=',day_ju_ini_armcu print *,'nt_armcu=',nt_armcu print *,'dt_armcu=',dt_armcu print *,'nlev_armcu=',nlev_armcu CALL interp_armcu_time(daytime,day1,annee_ref & & ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu & & ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu & & ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof & & ,adv_theta_prof,rad_theta_prof,adv_qt_prof) write(*,*) 'Forcages interpoles dans temps' ! No vertical interpolation if nlev imposed to 19 or 40 ! The vertical grid stops at 4000m # 600hPa mxcalc=llm ! initial and boundary conditions : ! tsurf = ts_prof ! tsurf read in lmdz1d.def write(*,*) 'Tsurf initiale: ',tsurf do l = 1, llm play(l)=play_mod(l)*100. presnivs(l)=play(l) zlay(l)=height(l) temp(l) = t_mod(l) teta(l)=theta_mod(l) q(l,1) = qv_mod(l)/1000. ! No liquid water in the initial profil q(l,2) = 0. u(l) = u_mod(l) ug(l)= u_mod(l) v(l) = v_mod(l) vg(l)= v_mod(l) ! Advective forcings are given in K or g/kg ... per HOUR ! IF(height(l).LT.1000) THEN ! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600. ! d_q_adv(l,1) = adv_qt_prof/1000./3600. ! d_q_adv(l,2) = 0.0 ! ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN ! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)* ! : (1-(height(l)-1000.)/2000.) ! d_th_adv(l) = d_th_adv(l)/3600. ! d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.) ! d_q_adv(l,1) = d_q_adv(l,1)/1000./3600. ! d_q_adv(l,2) = 0.0 ! ELSE ! d_th_adv(l) = 0.0 ! d_q_adv(l,1) = 0.0 ! d_q_adv(l,2) = 0.0 ! ENDIF enddo ! plev at half levels is given in proh.inp.19 or proh.inp.40 files plev(1)= ap(llm+1)+bp(llm+1)*psurf do l = 1, llm plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf print *,'Read_forc: l height play plev zlay temp', & & l,height(l),play(l),plev(l),zlay(l),temp(l) enddo ! For this case, fluxes are imposed fsens=-1*sens_prof flat=-1*flat_prof endif ! forcing_armcu !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing from transition case of Irina Sandu !--------------------------------------------------------------------- if (forcing_sandu) then write(*,*) 'Avant lecture Forcing SANDU' ! read sanduref forcing : fich_sandu = './ifa_sanduref.txt' CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu) write(*,*) 'Forcing SANDU lu' !---------------------------------------------------------------------- ! Read profiles from file: prof.inp.001 !---------------------------------------------------------------------- call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs, & & thl_profs,q_profs,u_profs,v_profs, & & w_profs,omega_profs,o3mmr_profs) ! time interpolation for initial conditions: write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! ! revoir 1DUTILS.h et les arguments print *,'Avant interp_sandu_time' print *,'daytime=',daytime print *,'day1=',day1 print *,'annee_ref=',annee_ref print *,'year_ini_sandu=',year_ini_sandu print *,'day_ju_ini_sandu=',day_ju_ini_sandu print *,'nt_sandu=',nt_sandu print *,'dt_sandu=',dt_sandu print *,'nlev_sandu=',nlev_sandu CALL interp_sandu_time(daytime,day1,annee_ref & & ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu & & ,nlev_sandu & & ,ts_sandu,ts_prof) ! vertical interpolation: print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu CALL interp_sandu_vertical(play,nlev_sandu,plev_profs & & ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs & & ,omega_profs,o3mmr_profs & & ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod & & ,omega_mod,o3mmr_mod,mxcalc) write(*,*) 'Profil initial forcing SANDU interpole' ! initial and boundary conditions : tsurf = ts_prof write(*,*) 'SST initiale: ',tsurf do l = 1, llm temp(l) = t_mod(l) tetal(l)=thl_mod(l) q(l,1) = q_mod(l) q(l,2) = 0.0 u(l) = u_mod(l) v(l) = v_mod(l) w(l) = w_mod(l) omega(l) = omega_mod(l) omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) !? omega2(l)=-rho(l)*omega(l) alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) ! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) ! d_q_adv(l,1) = vq_mod(l) d_th_adv(l) = alpha*omega(l)/rcpd d_q_adv(l,1) = 0.0 d_q_adv(l,2) = 0.0 enddo endif ! forcing_sandu !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- ! Forcing from Astex case !--------------------------------------------------------------------- if (forcing_astex) then write(*,*) 'Avant lecture Forcing Astex' ! read astex forcing : fich_astex = './ifa_astex.txt' CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex, & & ug_astex,vg_astex,ufa_astex,vfa_astex) write(*,*) 'Forcing Astex lu' !---------------------------------------------------------------------- ! Read profiles from file: prof.inp.001 !---------------------------------------------------------------------- call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa, & & thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa, & & w_profa,tke_profa,o3mmr_profa) ! time interpolation for initial conditions: write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! ! revoir 1DUTILS.h et les arguments print *,'Avant interp_astex_time' print *,'daytime=',daytime print *,'day1=',day1 print *,'annee_ref=',annee_ref print *,'year_ini_astex=',year_ini_astex print *,'day_ju_ini_astex=',day_ju_ini_astex print *,'nt_astex=',nt_astex print *,'dt_astex=',dt_astex print *,'nlev_astex=',nlev_astex CALL interp_astex_time(daytime,day1,annee_ref & & ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex & & ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex & & ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof & & ,ufa_prof,vfa_prof) ! vertical interpolation: print *,'Avant interp_vertical: nlev_astex=',nlev_astex CALL interp_astex_vertical(play,nlev_astex,plev_profa & & ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa & & ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa & & ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod & & ,tke_mod,o3mmr_mod,mxcalc) write(*,*) 'Profil initial forcing Astex interpole' ! initial and boundary conditions : tsurf = ts_prof write(*,*) 'SST initiale: ',tsurf do l = 1, llm temp(l) = t_mod(l) tetal(l)=thl_mod(l) q(l,1) = qv_mod(l) q(l,2) = ql_mod(l) u(l) = u_mod(l) v(l) = v_mod(l) w(l) = w_mod(l) omega(l) = w_mod(l) ! omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq ! rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) ! omega2(l)=-rho(l)*omega(l) alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) ! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) ! d_q_adv(l,1) = vq_mod(l) d_th_adv(l) = alpha*omega(l)/rcpd d_q_adv(l,1) = 0.0 d_q_adv(l,2) = 0.0 enddo endif ! forcing_astex