Changeset 2720 for LMDZ5/branches/testing/libf/phylmd/dyn1d/1D_interp_cases.h

Timestamp:

Nov 30, 2016, 1:28:41 PM (8 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2664:2719 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/dyn1d/1D_interp_cases.h (modified) (3 diffs)

LMDZ5/branches/testing/libf/phylmd/dyn1d/1D_interp_cases.h

@@ -118 +118 @@
 ! vertical interpolation:
       CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice        &
-     &         ,th_dice,qv_dice,u_dice,v_dice,o3_dice                   &
+     &         ,t_dice,qv_dice,u_dice,v_dice,o3_dice                   &
      &         ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd,omega_profd &
-     &         ,th_mod,qv_mod,u_mod,v_mod,o3_mod                        &
+     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                        &
      &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
 !     do l = 1, llm
@@ -192 +192 @@
       endif ! forcing_dice
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Interpolation gabls4 forcing 
+!---------------------------------------------------------------------
+      if (forcing_gabls4 ) then
+
+       if (prt_level.ge.1) then
+        print*,'#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_gabls4=',&
+     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_gabls4
+       endif
+
+! time interpolation:
+      CALL interp_gabls4_time(daytime,day1,annee_ref                                     &
+     &             ,year_ini_gabls4,day_ju_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4  & 
+     &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                            &
+     &             ,ug_profg,vg_profg,ht_profg,hq_profg,tg_profg)
+
+        if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
+
+! vertical interpolation:
+! on re-utilise le programme interp_dice_vertical: les transformations sur 
+! plev_gabls4,th_gabls4,qv_gabls4,u_gabls4,v_gabls4 ne sont pas prises en compte.
+! seules celles sur ht_profg,hq_profg,ug_profg,vg_profg sont prises en compte.
+
+      CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4         &
+!    &         ,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,poub            &
+     &         ,poub,poub,poub,poub,poub                             &
+     &         ,ht_profg,hq_profg,ug_profg,vg_profg,poub,poub        &
+     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                      &
+     &         ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc)
+
+      do l = 1, llm
+       ug(l)= ug_mod(l)
+       vg(l)= vg_mod(l)
+       d_th_adv(l)=ht_mod(l)
+       d_q_adv(l,1)=hq_mod(l)
+      enddo
+
+      endif ! forcing_gabls4
+!---------------------------------------------------------------------
+
 !---------------------------------------------------------------------
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -766 +805 @@
       enddo
 
+! Faut-il multiplier par -1 ? (MPL 20160713)
+      IF(ok_flux_surf) THEN
+       fsens=sens_prof_cas
+       flat=lat_prof_cas
+      ENDIF
+!
+      IF (ok_prescr_ust) THEN
+       ust=ustar_prof_cas
+       print *,'ust=',ust
+      ENDIF
       endif ! forcing_case
 
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+!---------------------------------------------------------------------
+! Interpolation forcing standard case
+!---------------------------------------------------------------------
+      if (forcing_case2) then
+
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
+     &    daytime,day1,(daytime-day1)*86400.,                               &
+     &    (daytime-day1)*86400/pdt_cas
+
+! time interpolation:
+        CALL interp2_case_time(daytime,day1,annee_ref                                       &
+!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
+     &       ,nt_cas,nlev_cas                                                               &
+     &       ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas      &
+     &       ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas             &
+     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
+     &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
+     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
+!
+     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
+     &       ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas                              &
+     &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas    &
+     &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
+     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
+     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
+     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
+     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
+
+             ts_cur = ts_prof_cas 
+!            psurf=plev_prof_cas(1)
+             psurf=ps_prof_cas
+
+! vertical interpolation:
+      CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas                                              &
+     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                          &
+     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
+     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                      &
+     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
+     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
+     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                    &
+!
+     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas          &
+     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                         &
+     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
+     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
+     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
+
+
+      DO l=1,llm
+      teta(l)=temp(l)*(100000./play(l))**(rd/rcpd)
+      ENDDO
+!calcul de l'advection verticale a partir du omega
+!Calcul des gradients verticaux
+!initialisation
+      d_t_z(:)=0.
+      d_th_z(:)=0.
+      d_q_z(:)=0.
+      d_t_dyn_z(:)=0.
+      d_th_dyn_z(:)=0.
+      d_q_dyn_z(:)=0.
+      DO l=2,llm-1
+       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
+       d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
+       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
+      ENDDO
+      d_t_z(1)=d_t_z(2)
+      d_th_z(1)=d_th_z(2)
+      d_q_z(1)=d_q_z(2)
+      d_t_z(llm)=d_t_z(llm-1)
+      d_th_z(llm)=d_th_z(llm-1)
+      d_q_z(llm)=d_q_z(llm-1)
+
+!Calcul de l advection verticale
+      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
+      d_th_dyn_z(:)=w_mod_cas(:)*d_th_z(:)
+      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
+
+!wind nudging 
+      if (nudging_u.gt.0.) then
+        do l=1,llm
+           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
+        enddo
+      else
+        do l=1,llm
+        ug(l) = u_mod_cas(l) 
+        enddo
+      endif
+
+      if (nudging_v.gt.0.) then
+        do l=1,llm
+           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
+        enddo
+      else
+        do l=1,llm
+        vg(l) = v_mod_cas(l) 
+        enddo
+      endif
+
+      if (nudging_w.gt.0.) then
+        do l=1,llm
+           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
+        enddo
+      else
+        do l=1,llm
+        w(l) = w_mod_cas(l) 
+        enddo
+      endif
+
+!nudging of q and temp 
+      if (nudging_t.gt.0.) then
+        do l=1,llm
+           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
+        enddo
+      endif
+      if (nudging_q.gt.0.) then
+        do l=1,llm
+           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
+        enddo
+      endif
+
+      do l = 1, llm
+       omega(l) = w_mod_cas(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)
+
+!calcul advection
+!       if ((tend_u.eq.1).and.(tend_w.eq.0)) then
+!          d_u_adv(l)=du_mod_cas(l)
+!       else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
+!          d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
+!       endif
+!
+!       if ((tend_v.eq.1).and.(tend_w.eq.0)) then
+!          d_v_adv(l)=dv_mod_cas(l)
+!       else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
+!          d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
+!       endif
+!
+!-----------------------------------------------------
+        if (tadv.eq.1 .or. tadvh.eq.1) then
+           d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
+        else if (tadvv.eq.1) then
+! ATTENTION d_t_dyn_z pas calcule (voir twpice)
+           d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
+        endif
+        print *,'interp_case d_t_dyn_z=',d_t_dyn_z(l),d_q_dyn_z(l)
+
+! Verifier le signe !!
+        if (thadv.eq.1 .or. thadvh.eq.1) then
+           d_th_adv(l)=dth_mod_cas(l)
+           print *,'dthadv=',d_th_adv(l)*86400.
+        else if (thadvv.eq.1) then
+           d_th_adv(l)=hth_mod_cas(l)-d_th_dyn_z(l)
+        endif
+ 
+! Verifier le signe !!
+        if ((qadv.eq.1).and.(forc_w.eq.0)) then
+           d_q_adv(l,1)=dq_mod_cas(l)
+        else if ((qadvh.eq.1).and.(forc_w.eq.1)) then
+           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
+        endif
+         
+        if (trad.eq.1) then
+           tend_rayo=1
+           dt_cooling(l) = dtrad_mod_cas(l)
+!          print *,'dt_cooling=',dt_cooling(l)
+        else
+           dt_cooling(l) = 0.0
+        endif
+      enddo
+
+! Faut-il multiplier par -1 ? (MPL 20160713)
+      IF(ok_flux_surf) THEN
+       fsens=-1.*sens_prof_cas
+       flat=-1.*lat_prof_cas
+       print *,'1D_interp: sens,flat',fsens,flat
+      ENDIF
+!
+      IF (ok_prescr_ust) THEN
+       ust=ustar_prof_cas
+       print *,'ust=',ust
+      ENDIF
+      endif ! forcing_case2
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+