Changeset 2720 for LMDZ5/branches/testing/libf/phylmd/dyn1d/1D_read_forc_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_read_forc_cases.h (modified) (6 diffs)

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

@@ -367 +367 @@
       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 &
+     &     ,zz_dice,plev_dice,t_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              &
@@ -376 +376 @@
 !champs initiaux:
       do k=1,nlev_dice
-         th_dicei(k)=th_dice(k)
+         t_dicei(k)=t_dice(k)
          qv_dicei(k)=qv_dice(k)
          u_dicei(k)=u_dice(k)
@@ -405 +405 @@
 
       CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice       &
-     &         ,th_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei             &
+     &         ,t_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                       &
+     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                       &
      &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
 
@@ -425 +425 @@
       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) 
+!      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
@@ -473 +473 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !---------------------------------------------------------------------
+! Forcing from GABLS4 experiment
+!---------------------------------------------------------------------
+
+!!!! Si la temperature de surface n'est pas impos??e:
+  
+      if (forcing_gabls4) then
+!read GABLS4 forcings
+      
+      fich_gabls4='gabls4_driver.nc'
+      
+        
+      call read_gabls4(fich_gabls4,nlev_gabls4,nt_gabls4,nsol_gabls4,zz_gabls4,depth_sn_gabls4,ug_gabls4,vg_gabls4 &
+     & ,plev_gabls4,th_gabls4,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,ht_gabls4,hq_gabls4,tg_gabls4,tsnow_gabls4,snow_dens_gabls4)
+
+      write(*,*) 'Forcing GABLS4 lu'
+
+!champs initiaux:
+      do k=1,nlev_gabls4
+         t_gabi(k)=t_gabls4(k)
+         qv_gabi(k)=qv_gabls4(k)
+         u_gabi(k)=u_gabls4(k)
+         v_gabi(k)=v_gabls4(k)
+         poub(k)=0.
+         ht_gabi(k)=ht_gabls4(k,1)
+         hq_gabi(k)=hq_gabls4(k,1)
+         ug_gabi(k)=ug_gabls4(k,1)
+         vg_gabi(k)=vg_gabls4(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_gabls4,nt_gabls4,plev_gabls4       &
+     &         ,t_gabi,qv_gabi,u_gabi,v_gabi,poub                  &
+     &         ,ht_gabi,hq_gabi,ug_gabi,vg_gabi,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)
+
+! Les forcages GABLS4 ont l air d etre en K/S quoiqu en dise le fichier gabls4_driver.nc !? MPL 20141024
+!     ht_mod(:)=ht_mod(:)/86400.
+!     hq_mod(:)=hq_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_mod(l)
+       vg(l)=vg_mod(l)
+       
+!
+!       tg=tsurf
+!       
+
+       print *,'***** tsurf=',tsurf
+       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     
+
+!--------------- Residus forcages du cas Dice (a supprimer) MPL 20141024---------------
+! 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_gabls4
+
+
+
 ! Forcing from Arm_Cu case                   
 ! For this case, ifa_armcu.txt contains sensible, latent heat fluxes
@@ -797 +909 @@
       endif !forcing_case
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Forcing from standard case :
+!---------------------------------------------------------------------
+
+      if (forcing_case2) then
+
+         write(*,*),'avant call read2_1D_cas'
+         call read2_1D_cas
+         write(*,*) 'Forcing read'
+
+!Time interpolation for initial conditions using interpolation routine
+         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
+        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)
+
+      do l = 1, nlev_cas
+      print *,'apres 1ere interp: plev_cas, plev_prof_cas=',l,plev_cas(l,1),plev_prof_cas(l)
+      enddo
+
+! vertical interpolation using interpolation routine:
+!      write(*,*)'avant interp vert', t_prof
+      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)
+
+!       write(*,*) 'Profil initial forcing case interpole',t_mod
+
+! initial and boundary conditions :
+!      tsurf = ts_prof_cas
+      ts_cur = ts_prof_cas 
+      psurf=plev_prof_cas(1)
+      write(*,*) 'SST initiale: ',tsurf
+      do l = 1, llm
+       temp(l) = t_mod_cas(l)
+       q(l,1) = qv_mod_cas(l)
+       q(l,2) = ql_mod_cas(l)
+       u(l) = u_mod_cas(l)
+       ug(l)= u_mod_cas(l)
+       v(l) = v_mod_cas(l)
+       vg(l)= v_mod_cas(l)
+       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)
+!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_cas(l)+vt_mod_cas(l))
+       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
+!      d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
+       d_q_adv(l,1) = dq_mod_cas(l)
+       d_q_adv(l,2) = 0.0
+!      d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
+       d_u_adv(l) = du_mod_cas(l)
+!      d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
+       d_u_adv(l) = dv_mod_cas(l)
+      enddo     
+
+! Faut-il multiplier par -1 ? (MPL 20160713)
+       IF (ok_flux_surf) THEN
+       fsens=-1.*sens_prof_cas
+       flat=-1.*lat_prof_cas
+       ENDIF
+!
+       IF (ok_prescr_ust) THEN
+       ust=ustar_prof_cas
+       print *,'ust=',ust
+       ENDIF
+
+      endif !forcing_case2
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+