Index: LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90 (revision 3594) +++ LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_cases_read_std.F90 (revision 3595) @@ -393,5 +393,5 @@ ! Activation de quelques cles en fonction des variables disponibles !----------------------------------------------------------------------- -if ( 1 == 0 ) THEN +if ( 1 == 1 ) THEN if ( name_var(i) == 'temp_nudg' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp' if ( name_var(i) == 'qv_nudg' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv' Index: LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90 (revision 3594) +++ LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90 (revision 3595) @@ -59,4 +59,6 @@ ! DECLARATIONS !===================================================================== + +#undef OUTPUT_PHYS_SCM !--------------------------------------------------------------------- @@ -142,11 +144,4 @@ ! (cf read_tsurf1d.F) -real wwww -!vertical advection computation -! real d_t_z(llm), d_q_z(llm) -! real d_t_dyn_z(llm), dq_dyn_z(llm) -! real zz(llm) -! real zfact - !flag forcings logical :: nudge_wind=.true. @@ -191,5 +186,6 @@ real :: sfdt, cfdt real :: du_phys(llm),dv_phys(llm),dt_phys(llm) - real :: dt_dyn(llm) + real :: w_adv(llm),z_adv(llm) + real :: d_t_vert_adv(llm),d_u_vert_adv(llm),d_v_vert_adv(llm) real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm) real :: d_u_nudge(llm),d_v_nudge(llm) @@ -201,5 +197,5 @@ REAL, ALLOCATABLE, DIMENSION(:,:):: q REAL, ALLOCATABLE, DIMENSION(:,:):: dq - REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn + REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_vert_adv REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge @@ -269,5 +265,7 @@ dv_phys(:)=0. dt_phys(:)=0. - dt_dyn(:)=0. + d_t_vert_adv(:)=0. + d_u_vert_adv(:)=0. + d_v_vert_adv(:)=0. dt_cooling(:)=0. d_t_adv(:)=0. @@ -420,5 +418,5 @@ allocate(q(llm,nqtot)) ; q(:,:)=0. allocate(dq(llm,nqtot)) - allocate(dq_dyn(llm,nqtot)) + allocate(d_q_vert_adv(llm,nqtot)) allocate(d_q_adv(llm,nqtot)) allocate(d_q_nudge(llm,nqtot)) @@ -427,5 +425,5 @@ q(:,:) = 0. dq(:,:) = 0. - dq_dyn(:,:) = 0. + d_q_vert_adv(:,:) = 0. d_q_adv(:,:) = 0. d_q_nudge(:,:) = 0. @@ -808,5 +806,5 @@ ! raz for safety do l=1,llm - dq_dyn(l,1) = 0. + d_q_vert_adv(l,1) = 0. enddo endif @@ -819,5 +817,7 @@ ! !===================================================================== +#ifdef OUTPUT_PHYS_SCM CALL iophys_ini +#endif ! START OF THE TEMPORAL LOOP : !===================================================================== @@ -841,17 +841,9 @@ #include "1D_interp_cases.h" ! Vertical advection -! call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play) +! call lstendH(llm,nqtot,omega,d_t_vert_adv,d_q_vert_adv,q,temp,u,v,play) ! print*,'B d_t_adv ',d_t_adv(1:20)*86400 -! print*,'B dt_dyn ',dt_dyn(1:20)*86400 +! print*,'B d_t_vert_adv ',d_t_vert_adv(1:20)*86400 ! print*,'B dt omega ',omega - teta=temp*(pzero/play)**rkappa - do l=2,llm-1 - dt_dyn(l)=-(omega(l)*(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1)))/(pzero/play(l))**rkappa - dq_dyn(l,1)=-omega(l)*(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1)) - enddo - d_t_adv(:)=d_t_adv(:)+dt_dyn(:) - d_q_adv(:,1)=d_q_adv(:,1)+dq_dyn(:,1) - - print*,'OMEGA ',omega_mod_cas(10),omega(10) + !--------------------------------------------------------------------- ! Geopotential : @@ -863,4 +855,32 @@ & (play(l)-play(l+1))/(play(l)+play(l+1)) enddo + +!--------------------------------------------------------------------- +! Vertical advection +!--------------------------------------------------------------------- + + IF ( forc_w+forc_omega > 0 ) THEN + + IF ( forc_w == 1 ) THEN + w_adv=w_mod_cas + z_adv=phi/RG + ELSE + w_adv=omega + z_adv=play + ENDIF + + teta=temp*(pzero/play)**rkappa + do l=2,llm-1 + d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1)) + d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1)) + d_t_vert_adv(l)=-(w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)))/(pzero/play(l))**rkappa + d_q_vert_adv(l,1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1)) + enddo + d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:) + d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1) + + print*,'OMEGA ',w_adv(10),z_adv(10) + + ENDIF !--------------------------------------------------------------------- @@ -882,7 +902,4 @@ & presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm) endif - - CALL iophys_ecrit('dtadv',klev,'dtadv','K/day',86400*d_t_adv) - CALL iophys_ecrit('dtdyn',klev,'dtdyn','K/day',86400*dt_dyn) !--------------------------------------------------------------------- @@ -1021,4 +1038,10 @@ & d_q_nudge(l,1)=(qv_nudg_mod_cas(l)-q(l,1))/nudging_qv ENDDO + +#ifdef OUTPUT_PHYS_SCM + CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv) + CALL iophys_ecrit('z_adv',klev,'z_adv','K/day',z_adv) + CALL iophys_ecrit('dtadv',klev,'dtadv','K/day',86400*d_t_adv) + CALL iophys_ecrit('dtdyn',klev,'dtdyn','K/day',86400*d_t_vert_adv) CALL iophys_ecrit('qv',klev,'qv','g/kg',1000*q(:,1)) CALL iophys_ecrit('qvnud',klev,'qvnud','g/kg',1000*u_nudg_mod_cas) @@ -1027,6 +1050,9 @@ CALL iophys_ecrit('v',klev,'v','m/s',v) CALL iophys_ecrit('vnud',klev,'vnud','m/s',v_nudg_mod_cas) + CALL iophys_ecrit('temp',klev,'temp','K',temp) + CALL iophys_ecrit('tempnud',klev,'temp_nudg_mod_cas','K',temp_nudg_mod_cas) CALL iophys_ecrit('dtnud',klev,'dtnud','K/day',86400*d_t_nudge) CALL iophys_ecrit('dqnud',klev,'dqnud','K/day',1000*86400*d_q_nudge(:,1)) +#endif !