source: LMDZ6/trunk/libf/phylmd/dyn1d/1D_interp_cases.h @ 5456

Last change on this file since 5456 was 4297, checked in by evignon, 2 years ago

technical corrections (+cleaning) to make standard BOMEX case work

  • Property svn:keywords set to Id
File size: 8.7 KB
RevLine 
[2017]1
[3780]2        print*,'FORCING CASE forcing_case2'
[3594]3!       print*,                                                             &
4!    & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
5!    &    daytime,day1,(daytime-day1)*86400.,                               &
6!    &    (daytime-day1)*86400/pdt_cas
[2191]7
8! time interpolation:
[3686]9        CALL interp_case_time_std(daytime,day1,annee_ref                                       &
10!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
11     &       ,nt_cas,nlev_cas                                                               &
[4297]12     &       ,ts_cas,tskin_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas      &
[3686]13     &       ,u_cas,v_cas,ug_cas,vg_cas                                                     &
14     &       ,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas                               &
[4104]15     &       ,invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_cas    &
[3781]16     &       ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas                                       &
[3686]17     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
18     &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
[3781]19     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                                           &
[2716]20!
[4297]21     &       ,ts_prof_cas,tskin_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
[3686]22     &       ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas                              &
23     &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                                 &
24     &       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas           &
[4104]25     &       ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas  &     
[3781]26     &       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                     &
[3686]27     &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
28     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
29     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
30     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
[3781]31     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
[3780]32! EV tg instead of ts_cur
[4297]33             tg = ts_prof_cas
34             if ((tg .eq. 0.) .and. (tskin_prof_cas .ne. 0.)) then
35                tg=tskin_prof_cas
36              endif
37
[2716]38!            psurf=plev_prof_cas(1)
39             psurf=ps_prof_cas
40
41! vertical interpolation:
[3781]42      CALL interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                                              &
[3686]43     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                       &
44     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
45     &         ,ug_prof_cas,vg_prof_cas                                                                   &
46     &         ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                       &
[4104]47     &         ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas &
48
[3781]49     &         ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                                 &
[3686]50     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
51     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
52     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                    &
[2716]53!
[3686]54     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas          &
55     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas                                                 &
56     &         ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                           &
[4104]57     &         ,invtau_temp_nudg_mod_cas,invtau_qv_nudg_mod_cas,invtau_u_nudg_mod_cas,invtau_v_nudg_mod_cas   &     
[3781]58     &         ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                       &
[3686]59     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
60     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
61     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
[2716]62
63
64      DO l=1,llm
65      teta(l)=temp(l)*(100000./play(l))**(rd/rcpd)
66      ENDDO
67!calcul de l'advection verticale a partir du omega
68!Calcul des gradients verticaux
69!initialisation
70      d_t_z(:)=0.
71      d_th_z(:)=0.
72      d_q_z(:)=0.
[2920]73      d_u_z(:)=0.
74      d_v_z(:)=0.
[2716]75      d_t_dyn_z(:)=0.
76      d_th_dyn_z(:)=0.
77      d_q_dyn_z(:)=0.
[2920]78      d_u_dyn_z(:)=0.
79      d_v_dyn_z(:)=0.
[3541]80      if (1==0) then
81         DO l=2,llm-1
82          d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
83          d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
84          d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
85          d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
86          d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
87         ENDDO
88      else
89         DO l=2,llm-1
[3686]90            IF (omega(l)>0.) THEN
[3541]91             d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
92             d_th_z(l)=(teta(l+1)-teta(l))/(play(l+1)-play(l))
93             d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
94             d_u_z(l)=(u(l+1)-u(l))/(play(l+1)-play(l))
95             d_v_z(l)=(v(l+1)-v(l))/(play(l+1)-play(l))
[3686]96            ELSE
[3541]97             d_t_z(l)=(temp(l-1)-temp(l))/(play(l-1)-play(l))
98             d_th_z(l)=(teta(l-1)-teta(l))/(play(l-1)-play(l))
99             d_q_z(l)=(q(l-1,1)-q(l,1))/(play(l-1)-play(l))
100             d_u_z(l)=(u(l-1)-u(l))/(play(l-1)-play(l))
101             d_v_z(l)=(v(l-1)-v(l))/(play(l-1)-play(l))
[3686]102            ENDIF
103         ENDDO
[3541]104      endif
[2716]105      d_t_z(1)=d_t_z(2)
[3541]106      d_t_z(1)=d_t_z(2)
[2716]107      d_th_z(1)=d_th_z(2)
108      d_q_z(1)=d_q_z(2)
[2920]109      d_u_z(1)=d_u_z(2)
110      d_v_z(1)=d_v_z(2)
[2716]111      d_t_z(llm)=d_t_z(llm-1)
112      d_th_z(llm)=d_th_z(llm-1)
113      d_q_z(llm)=d_q_z(llm-1)
[2920]114      d_u_z(llm)=d_u_z(llm-1)
115      d_v_z(llm)=d_v_z(llm-1)
[2716]116
[3541]117! TRAVAIL : PRENDRE DES NOTATIONS COHERENTES POUR W
118      do l = 1, llm
119! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
[3781]120        print*, l, llm
121        print*, play(l), temp(l)
[3541]122       omega(l) = -w_mod_cas(l)*play(l)*rg/(rd*temp(l))
123      enddo
124
[2716]125!Calcul de l advection verticale
[2920]126! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170310)
[3541]127      d_t_dyn_z(:)=omega(:)*d_t_z(:)
128      d_th_dyn_z(:)=omega(:)*d_th_z(:)
129      d_q_dyn_z(:)=omega(:)*d_q_z(:)
130      d_u_dyn_z(:)=omega(:)*d_u_z(:)
131      d_v_dyn_z(:)=omega(:)*d_v_z(:)
[2716]132
[2920]133!geostrophic wind
134      if (forc_geo.eq.1) then
135        do l=1,llm
136        ug(l) = ug_mod_cas(l) 
137        vg(l) = vg_mod_cas(l) 
138        enddo
139      endif
[2716]140
141      do l = 1, llm
[3592]142
143!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2920]144! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
[3592]145       !!! omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
146       omega(l) = omega_mod_cas(l)
147       omega2(l)= omega_mod_cas(l)/rg*airefi ! flxmass_w calcule comme ds physiq
148
[3593]149! On effectue la somme du forcage total et de la decomposition
150! horizontal/vertical en supposant que soit l'un soit l'autre
151! sont remplis mais jamais les deux
[2716]152
[3593]153       d_t_adv(l) = dt_mod_cas(l)+ht_mod_cas(l)+vt_mod_cas(l)
154       d_q_adv(l,1) = dq_mod_cas(l)+hq_mod_cas(l)+vq_mod_cas(l)
155       d_q_adv(l,2) = 0.0
156       d_u_adv(l) = du_mod_cas(l)+hu_mod_cas(l)+vu_mod_cas(l)
157       d_v_adv(l) = dv_mod_cas(l)+hv_mod_cas(l)+vv_mod_cas(l)
[2716]158
[3592]159!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
160!! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
161!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
162!if (forc_w==1) then
163!          d_q_adv(l,1)=d_q_adv(l,1)-d_q_dyn_z(l)
164!          d_t_adv(l)=d_t_adv(l)-d_t_dyn_z(l)
165!          d_v_adv(l)=d_v_adv(l)-d_v_dyn_z(l)
166!          d_u_adv(l)=d_u_adv(l)-d_u_dyn_z(l)
167!       endif
168!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2716]169         
[3686]170        if (trad.eq.1) then
[2716]171           tend_rayo=1
172           dt_cooling(l) = dtrad_mod_cas(l)
173!          print *,'dt_cooling=',dt_cooling(l)
174        else
175           dt_cooling(l) = 0.0
176        endif
177      enddo
178
179! Faut-il multiplier par -1 ? (MPL 20160713)
180      IF(ok_flux_surf) THEN
181       fsens=-1.*sens_prof_cas
182       flat=-1.*lat_prof_cas
183       print *,'1D_interp: sens,flat',fsens,flat
184      ENDIF
185!
186      IF (ok_prescr_ust) THEN
187       ust=ustar_prof_cas
188       print *,'ust=',ust
189      ENDIF
Note: See TracBrowser for help on using the repository browser.