source: LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_cases_read_std.f90 @ 5284

Last change on this file since 5284 was 5274, checked in by abarral, 7 days ago

Replace yomcst.h by existing module

File size: 59.3 KB
Line 
1!
2! $Id: mod_1D_cases_read.F90 2373 2015-10-13 17:28:01Z jyg $
3!
4MODULE mod_1D_cases_read_std
5  USE netcdf, ONLY: nf90_noerr, nf90_inq_varid, nf90_inq_dimid, nf90_inquire_dimension, nf90_open, nf90_nowrite, &
6          nf90_strerror, nf90_get_var
7
8!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9  !Declarations specifiques au cas standard
10  character*80 :: fich_cas
11  ! Discr?tisation
12  integer nlev_cas, nt_cas
13
14
15  !profils environnementaux
16  real, allocatable::  plev_cas(:,:),plevh_cas(:)
17  real, allocatable::  ap_cas(:),bp_cas(:)
18
19  real, allocatable::  z_cas(:,:),zh_cas(:)
20  real, allocatable::  t_cas(:,:),q_cas(:,:),qv_cas(:,:),ql_cas(:,:),qi_cas(:,:),rh_cas(:,:)
21  real, allocatable::  th_cas(:,:),thv_cas(:,:),thl_cas(:,:),rv_cas(:,:)
22  real, allocatable::  u_cas(:,:),v_cas(:,:),vitw_cas(:,:),omega_cas(:,:),tke_cas(:,:)
23
24  !forcing
25  real, allocatable::  ht_cas(:,:),vt_cas(:,:),dt_cas(:,:),dtrad_cas(:,:)
26  real, allocatable::  hth_cas(:,:),vth_cas(:,:),dth_cas(:,:)
27  real, allocatable::  hq_cas(:,:),vq_cas(:,:),dq_cas(:,:)
28  real, allocatable::  hr_cas(:,:),vr_cas(:,:),dr_cas(:,:)
29  real, allocatable::  hu_cas(:,:),vu_cas(:,:),du_cas(:,:)
30  real, allocatable::  hv_cas(:,:),vv_cas(:,:),dv_cas(:,:)
31  real, allocatable::  ug_cas(:,:),vg_cas(:,:)
32  real, allocatable::  temp_nudg_cas(:,:),qv_nudg_cas(:,:),u_nudg_cas(:,:),v_nudg_cas(:,:)
33  real, allocatable::  invtau_temp_nudg_cas(:,:),invtau_qv_nudg_cas(:,:),invtau_u_nudg_cas(:,:),invtau_v_nudg_cas(:,:)
34  real, allocatable::  lat_cas(:),sens_cas(:),tskin_cas(:),ts_cas(:),ps_cas(:),ustar_cas(:)
35  real, allocatable::  uw_cas(:,:),vw_cas(:,:),q1_cas(:,:),q2_cas(:,:),tkes_cas(:)
36
37  !champs interpoles
38  real, allocatable::  plev_prof_cas(:)
39  real, allocatable::  t_prof_cas(:)
40  real, allocatable::  theta_prof_cas(:)
41  real, allocatable::  thl_prof_cas(:)
42  real, allocatable::  thv_prof_cas(:)
43  real, allocatable::  q_prof_cas(:)
44  real, allocatable::  qv_prof_cas(:)
45  real, allocatable::  ql_prof_cas(:)
46  real, allocatable::  qi_prof_cas(:)
47  real, allocatable::  rh_prof_cas(:)
48  real, allocatable::  rv_prof_cas(:)
49  real, allocatable::  u_prof_cas(:)
50  real, allocatable::  v_prof_cas(:)       
51  real, allocatable::  vitw_prof_cas(:)
52  real, allocatable::  omega_prof_cas(:)
53  real, allocatable::  tke_prof_cas(:)
54  real, allocatable::  ug_prof_cas(:)
55  real, allocatable::  vg_prof_cas(:)
56  real, allocatable::  temp_nudg_prof_cas(:),qv_nudg_prof_cas(:),u_nudg_prof_cas(:),v_nudg_prof_cas(:)
57  real, allocatable::  invtau_temp_nudg_prof_cas(:),invtau_qv_nudg_prof_cas(:),invtau_u_nudg_prof_cas(:),invtau_v_nudg_prof_cas(:)
58
59  real, allocatable::  ht_prof_cas(:)
60  real, allocatable::  hth_prof_cas(:)
61  real, allocatable::  hq_prof_cas(:)
62  real, allocatable::  vt_prof_cas(:)
63  real, allocatable::  vth_prof_cas(:)
64  real, allocatable::  vq_prof_cas(:)
65  real, allocatable::  dt_prof_cas(:)
66  real, allocatable::  dth_prof_cas(:)
67  real, allocatable::  dtrad_prof_cas(:)
68  real, allocatable::  dq_prof_cas(:)
69  real, allocatable::  hu_prof_cas(:)
70  real, allocatable::  hv_prof_cas(:)
71  real, allocatable::  vu_prof_cas(:)
72  real, allocatable::  vv_prof_cas(:)
73  real, allocatable::  du_prof_cas(:)
74  real, allocatable::  dv_prof_cas(:)
75  real, allocatable::  uw_prof_cas(:)
76  real, allocatable::  vw_prof_cas(:)
77  real, allocatable::  q1_prof_cas(:)
78  real, allocatable::  q2_prof_cas(:)
79
80
81  real o3_cas,lat_prof_cas,sens_prof_cas,ts_prof_cas,tskin_prof_cas,ps_prof_cas,ustar_prof_cas,tkes_prof_cas
82  real orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough,rugos_cas,sand_cas,clay_cas
83
84
85
86CONTAINS
87
88
89  !**********************************************************************************************
90  SUBROUTINE read_SCM_cas
91    use netcdf, only: nf90_get_var
92    implicit none
93    INCLUDE "date_cas.h"
94
95    INTEGER nid,rid,ierr
96    INTEGER ii,jj,timeid
97    REAL, ALLOCATABLE :: time_val(:)
98
99    fich_cas='cas.nc'
100    print*,'fich_cas ',fich_cas
101    ierr = nf90_open(fich_cas,nf90_nowrite,nid)
102    print*,'fich_cas,nf90_nowrite,nid ',fich_cas,nf90_nowrite,nid
103    if (ierr.NE.nf90_noerr) then
104       write(*,*) 'ERROR: GROS Pb opening forcings nc file '
105       write(*,*) nf90_strerror(ierr)
106       stop ""
107    endif
108    !.......................................................................
109    ierr=nf90_inq_dimid(nid,'lat',rid)
110    IF (ierr.NE.nf90_noerr) THEN
111       print*, 'Oh probleme lecture dimension lat'
112    ENDIF
113    ierr=nf90_inquire_dimension(nid,rid,len=ii)
114    print*,'OK1 read_SCM_cas: nid,rid,lat',nid,rid,ii
115    !.......................................................................
116    ierr=nf90_inq_dimid(nid,'lon',rid)
117    IF (ierr.NE.nf90_noerr) THEN
118       print*, 'Oh probleme lecture dimension lon'
119    ENDIF
120    ierr=nf90_inquire_dimension(nid,rid,len=jj)
121    print*,'OK2 read_SCM_cas: nid,rid,lat',nid,rid,jj
122    !.......................................................................
123    ierr=nf90_inq_dimid(nid,'lev',rid)
124    IF (ierr.NE.nf90_noerr) THEN
125       print*, 'Oh probleme lecture dimension nlev'
126    ENDIF
127    ierr=nf90_inquire_dimension(nid,rid,len=nlev_cas)
128    print*,'OK3 read_SCM_cas: nid,rid,nlev_cas',nid,rid,nlev_cas
129    IF ( .NOT. ( nlev_cas > 10 .AND. nlev_cas < 200000 )) THEN
130       print*,'Valeur de nlev_cas peu probable'
131       STOP
132    ENDIF
133    !.......................................................................
134    ierr=nf90_inq_dimid(nid,'time',rid)
135    nt_cas=0
136    IF (ierr.NE.nf90_noerr) THEN
137       stop 'Oh probleme lecture dimension time'
138    ENDIF
139    ierr=nf90_inquire_dimension(nid,rid,len=nt_cas)
140    print*,'OK4 read_SCM_cas: nid,rid,nt_cas',nid,rid,nt_cas
141    ! Lecture de l'axe des temps
142    print*,'LECTURE DU TEMPS'
143    ierr=nf90_inq_varid(nid,'time',timeid)
144    if(ierr/=nf90_noerr) then
145       print *,'Variable time manquante dans cas.nc:'
146       ierr=nf90_noerr
147    else
148       allocate(time_val(nt_cas))
149       ierr = NF90_GET_VAR(nid,timeid,time_val)
150       if(ierr/=nf90_noerr) then
151          print *,'A Pb a la lecture de time cas.nc: '
152       endif
153    endif
154    IF (nt_cas>1) THEN
155       pdt_cas=time_val(2)-time_val(1)
156    ELSE
157       pdt_cas=0.
158    ENDIF
159
160
161!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
162    !profils moyens:
163    allocate(plev_cas(nlev_cas,nt_cas),plevh_cas(nlev_cas+1))       
164    allocate(z_cas(nlev_cas,nt_cas),zh_cas(nlev_cas+1))
165    allocate(ap_cas(nlev_cas+1),bp_cas(nlev_cas+1))
166    allocate(t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas),qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas), &
167         qi_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
168    allocate(th_cas(nlev_cas,nt_cas),thl_cas(nlev_cas,nt_cas),thv_cas(nlev_cas,nt_cas),rv_cas(nlev_cas,nt_cas))
169    allocate(u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas),vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas))
170    allocate(tke_cas(nlev_cas,nt_cas))
171    !forcing
172    allocate(ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas),dt_cas(nlev_cas,nt_cas),dtrad_cas(nlev_cas,nt_cas))
173    allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
174    allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
175    allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
176    allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
177    allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
178    allocate(ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas))
179    allocate(temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas))
180    allocate(u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas))
181    allocate(invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas))
182    allocate(invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_cas(nlev_cas,nt_cas))
183    allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas),tskin_cas(nt_cas),ps_cas(nt_cas),ustar_cas(nt_cas),tkes_cas(nt_cas))
184    allocate(uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas),q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas))
185
186
187
188    !champs interpoles
189    allocate(plev_prof_cas(nlev_cas))
190    allocate(t_prof_cas(nlev_cas))
191    allocate(theta_prof_cas(nlev_cas))
192    allocate(thl_prof_cas(nlev_cas))
193    allocate(thv_prof_cas(nlev_cas))
194    allocate(q_prof_cas(nlev_cas))
195    allocate(qv_prof_cas(nlev_cas))
196    allocate(ql_prof_cas(nlev_cas))
197    allocate(qi_prof_cas(nlev_cas))
198    allocate(rh_prof_cas(nlev_cas))
199    allocate(rv_prof_cas(nlev_cas))
200    allocate(u_prof_cas(nlev_cas))
201    allocate(v_prof_cas(nlev_cas))
202    allocate(vitw_prof_cas(nlev_cas))
203    allocate(omega_prof_cas(nlev_cas))
204    allocate(tke_prof_cas(nlev_cas))
205    allocate(ug_prof_cas(nlev_cas))
206    allocate(vg_prof_cas(nlev_cas))
207    allocate(temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas))
208    allocate(u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas))
209    allocate(invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas))
210    allocate(invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas))
211    allocate(ht_prof_cas(nlev_cas))
212    allocate(hth_prof_cas(nlev_cas))
213    allocate(hq_prof_cas(nlev_cas))
214    allocate(hu_prof_cas(nlev_cas))
215    allocate(hv_prof_cas(nlev_cas))
216    allocate(vt_prof_cas(nlev_cas))
217    allocate(vth_prof_cas(nlev_cas))
218    allocate(vq_prof_cas(nlev_cas))
219    allocate(vu_prof_cas(nlev_cas))
220    allocate(vv_prof_cas(nlev_cas))
221    allocate(dt_prof_cas(nlev_cas))
222    allocate(dth_prof_cas(nlev_cas))
223    allocate(dtrad_prof_cas(nlev_cas))
224    allocate(dq_prof_cas(nlev_cas))
225    allocate(du_prof_cas(nlev_cas))
226    allocate(dv_prof_cas(nlev_cas))
227    allocate(uw_prof_cas(nlev_cas))
228    allocate(vw_prof_cas(nlev_cas))
229    allocate(q1_prof_cas(nlev_cas))
230    allocate(q2_prof_cas(nlev_cas))
231
232    print*,'Allocations OK'
233    CALL read_SCM (nid,nlev_cas,nt_cas,                                                                     &
234         ap_cas,bp_cas,z_cas,plev_cas,zh_cas,plevh_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,                   &
235         ql_cas,qi_cas,rh_cas,rv_cas,u_cas,v_cas,vitw_cas,omega_cas,tke_cas,ug_cas,vg_cas,                            &
236         temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas,                                                     &
237         invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_cas,                         &
238         du_cas,hu_cas,vu_cas,                                                                                &
239         dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,               &
240         dr_cas,hr_cas,vr_cas,dtrad_cas,sens_cas,lat_cas,ts_cas,tskin_cas,ps_cas,ustar_cas,tkes_cas,            &
241         uw_cas,vw_cas,q1_cas,q2_cas,orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough, &
242         o3_cas,rugos_cas,clay_cas,sand_cas)
243    print*,'read_SCM cas OK'
244    do ii=1,nlev_cas
245       print*,'apres read_SCM_cas, plev_cas=',ii,plev_cas(ii,1)
246       !print*,'apres read_SCM, plev_cas=',ii,omega_cas(ii,nt_cas/2+1)
247    enddo
248
249
250  END SUBROUTINE read_SCM_cas
251
252
253!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
254  SUBROUTINE deallocate2_1D_cases
255    !profils environnementaux:
256    deallocate(plev_cas,plevh_cas)
257
258    deallocate(z_cas,zh_cas)
259    deallocate(ap_cas,bp_cas)
260    deallocate(t_cas,q_cas,qv_cas,ql_cas,qi_cas,rh_cas)
261    deallocate(th_cas,thl_cas,thv_cas,rv_cas)
262    deallocate(u_cas,v_cas,vitw_cas,omega_cas,tke_cas)
263
264    !forcing
265    deallocate(ht_cas,vt_cas,dt_cas,dtrad_cas)
266    deallocate(hq_cas,vq_cas,dq_cas)
267    deallocate(hth_cas,vth_cas,dth_cas)
268    deallocate(hr_cas,vr_cas,dr_cas)
269    deallocate(hu_cas,vu_cas,du_cas)
270    deallocate(hv_cas,vv_cas,dv_cas)
271    deallocate(ug_cas)
272    deallocate(vg_cas)
273    deallocate(lat_cas,sens_cas,tskin_cas,ts_cas,ps_cas,ustar_cas,tkes_cas,uw_cas,vw_cas,q1_cas,q2_cas)
274
275    !champs interpoles
276    deallocate(plev_prof_cas)
277    deallocate(t_prof_cas)
278    deallocate(theta_prof_cas)
279    deallocate(thl_prof_cas)
280    deallocate(thv_prof_cas)
281    deallocate(q_prof_cas)
282    deallocate(qv_prof_cas)
283    deallocate(ql_prof_cas)
284    deallocate(qi_prof_cas)
285    deallocate(rh_prof_cas)
286    deallocate(rv_prof_cas)
287    deallocate(u_prof_cas)
288    deallocate(v_prof_cas)
289    deallocate(vitw_prof_cas)
290    deallocate(omega_prof_cas)
291    deallocate(tke_prof_cas)
292    deallocate(ug_prof_cas)
293    deallocate(vg_prof_cas)
294    deallocate(temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas)
295    deallocate(invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas)
296    deallocate(ht_prof_cas)
297    deallocate(hq_prof_cas)
298    deallocate(hu_prof_cas)
299    deallocate(hv_prof_cas)
300    deallocate(vt_prof_cas)
301    deallocate(vq_prof_cas)
302    deallocate(vu_prof_cas)
303    deallocate(vv_prof_cas)
304    deallocate(dt_prof_cas)
305    deallocate(dtrad_prof_cas)
306    deallocate(dq_prof_cas)
307    deallocate(du_prof_cas)
308    deallocate(dv_prof_cas)
309    deallocate(t_prof_cas)
310    deallocate(u_prof_cas)
311    deallocate(v_prof_cas)
312    deallocate(uw_prof_cas)
313    deallocate(vw_prof_cas)
314    deallocate(q1_prof_cas)
315    deallocate(q2_prof_cas)
316
317  END SUBROUTINE deallocate2_1D_cases
318
319
320  !=====================================================================
321  SUBROUTINE read_SCM(nid,nlevel,ntime,                                       &
322       ap,bp,zz,pp,zzh,pph,temp,theta,thv,thl,qv,ql,qi,rh,rv,u,v,vitw,omega,tke,ug,vg,&
323       temp_nudg,qv_nudg,u_nudg,v_nudg,                                        &
324       invtau_temp_nudg,invtau_qv_nudg,invtau_u_nudg,invtau_v_nudg,             &
325       du,hu,vu,dv,hv,vv,dt,ht,vt,dq,hq,vq,                                    &
326       dth,hth,vth,dr,hr,vr,dtrad,sens,flat,ts,tskin,ps,ustar,tkes,uw,vw,q1,q2,       &
327       orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,          &
328       heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas)
329
330    !program reading forcing of the case study
331    use netcdf, only: nf90_get_var
332    implicit none
333    INCLUDE "compar1d.h"
334
335    integer ntime,nlevel,k,t
336
337    real ap(nlevel+1),bp(nlevel+1)
338    real zz(nlevel,ntime),zzh(nlevel+1)
339    real pp(nlevel,ntime),pph(nlevel+1)
340    !profils initiaux
341    real temp0(nlevel),qv0(nlevel),ql0(nlevel),qi0(nlevel),u0(nlevel),v0(nlevel),tke0(nlevel)
342    real pp0(nlevel)   
343    real temp(nlevel,ntime),qv(nlevel,ntime),ql(nlevel,ntime),qi(nlevel,ntime),rh(nlevel,ntime)
344    real theta(nlevel,ntime),thv(nlevel,ntime),thl(nlevel,ntime),rv(nlevel,ntime)
345    real u(nlevel,ntime),v(nlevel,ntime),tkes(ntime)
346    real temp_nudg(nlevel,ntime),qv_nudg(nlevel,ntime),u_nudg(nlevel,ntime),v_nudg(nlevel,ntime)
347    real invtau_temp_nudg(nlevel,ntime),invtau_qv_nudg(nlevel,ntime),invtau_u_nudg(nlevel,ntime),invtau_v_nudg(nlevel,ntime)
348    real ug(nlevel,ntime),vg(nlevel,ntime)
349    real vitw(nlevel,ntime),omega(nlevel,ntime),tke(nlevel,ntime)
350    real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
351    real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
352    real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
353    real dtrad(nlevel,ntime)
354    real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
355    real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime),hthl(nlevel,ntime)
356    real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
357    real flat(ntime),sens(ntime),ustar(ntime)
358    real uw(nlevel,ntime),vw(nlevel,ntime),q1(nlevel,ntime),q2(nlevel,ntime)
359    real ts(ntime),tskin(ntime),ps(ntime)
360    real orog_cas,albedo_cas,emiss_cas,q_skin_cas,mom_rough,heat_rough,o3_cas,rugos_cas,clay_cas,sand_cas
361    real apbp(nlevel+1),resul(nlevel,ntime),resul1(nlevel),resul2(ntime),resul3
362
363
364    integer nid, ierr,ierr1,ierr2,rid,i,int_test
365    integer nbvar3d
366    parameter(nbvar3d=78)
367    integer var3didin(nbvar3d),missing_var(nbvar3d)
368    character*13 name_var(1:nbvar3d)
369
370
371    !      data name_var/ &
372    !     ! coordonnees pression (n+1 niveaux) #4
373    !     & 'coor_par_a','coor_par_b','height_h','pressure_h',& ! #1-#4
374    !     ! coordonnees pression (n niveaux) #8
375    !     &'temp','qv','ql','qi','u','v','tke','pressure',& ! #5-#12
376    !     ! coordonnees pression + temps #42
377    !     &'w','omega','ug','vg','uadv','uadvh','uadvv','vadv','vadvh','vadvv','temp_adv','tadvh','tadvv',& !  #13 - #25
378    !     &'qv_adv','qvadvh','qvadvv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
379    !     & 'radv','radvh','radvv','radcool','q1','q2','ustress','vstress',                           & ! #33 - #40
380    !     & 'rh','temp_nudging','qv_nudging','u_nudging','v_nudging',                                       & ! #41-45
381    !     &'height_f','pressure_forc','tempt','theta','thv','thl','qvt','qlt','qit','rv','ut','vt',   & ! #46-58
382    !     ! coordonnees temps #12
383    !     &'tkes','sfc_sens_flx','sfc_lat_flx','ts','ps','ustar',&
384    !     &'orog','albedo','emiss','t_skin','q_skin','mom_rough','heat_rough',&
385    !     ! scalaires #4
386    !     &'o3','rugos','clay','sand'/
387
388
389
390    data name_var/ &
391                                ! coordonnees pression (n+1 niveaux) #4
392         'coor_par_a','coor_par_b','zf','pressure_h',& ! #1-#4
393                                ! coordonnees pression (n niveaux) #8
394         'ta','qv','ql','qi','ua','va','tke','pa',& ! #5-#12
395                                ! coordonnees pression + temps #46
396         'wa','wap','ug','vg','tnua_adv','tnua_advh','tnua_advv','tnva_adv','tnva_advh','tnva_advv','tnta_adv','tnta_advh',& !  #13 - #25
397         'tnta_advv','tnqv_adv','tnqv_advh','tnqv_advv','thadv','thadvh','thadvv','thladvh',                             & ! #26 - #32
398         'radv','radvh','radvv','tnta_rad','q1','q2','ustress','vstress',                           & ! #33 - #40
399         'rh','ta_nud','qv_nud','ua_nud','va_nud',                                       & ! #41-45
400         'zh_forc','pa_forc','tat','thetat','thetavt','thetalt','qvt','qlt','qit','rvt','uat','vat',   & ! #46-57
401         'nudging_constant_ta', 'nudging_constant_qv', 'nudging_constant_ua', 'nudging_constant_va',           & ! # 58-61
402                                ! coordonnees temps #12
403         'tkes','hfss','hfls','ts_forc','tskin','ps_forc','ustar', &                     ! 62-68
404                                  ! scalaires
405         'orog','albedo','emiss','q_skin','z0','z0h',       &                    ! 69-74
406         'O3','rugos','clay','sand'/                                                      ! 75-78
407
408
409    !-----------------------------------------------------------------------
410    ! First check that we are using a version > v2 of the 1D standard format
411    ! use the difference between 'temp' (old version) and 'ta' (new version)
412    !-----------------------------------------------------------------------
413
414
415    ierr=nf90_inq_varid(nid,'ta',int_test)
416    if(ierr/=nf90_noerr) then
417       print*, '++++++++++++++++++++++++++++++'
418       print*, 'variable ta missing in cas.nc '
419       print*, 'You are probably using an obsolete version of the 1D cases'
420       print*, 'please dowload the last version of the 1D archive from https://lmdz.lmd.jussieu.fr/pub/'
421       print*, '++++++++++++++++++++++++++++++'
422       CALL abort_gcm ('mod_1D_cases_read_std','bad version of 1D directory',0)
423    endif
424
425    !-----------------------------------------------------------------------
426    ! Checking availability of variable #i in the cas.nc file
427    !     missing_var=1 if the variable is missing
428    !-----------------------------------------------------------------------
429
430    do i=1,nbvar3d
431       missing_var(i)=0.
432       ierr=nf90_inq_varid(nid,name_var(i),var3didin(i))
433       print*, 'name_var(i)', name_var(i), var3didin(i)
434       if(ierr/=nf90_noerr) then
435          print *,'Variable manquante dans cas.nc:',i,name_var(i)
436          ierr=nf90_noerr
437          missing_var(i)=1
438       else
439
440          !-----------------------------------------------------------------------
441          ! Activating keys depending on the presence of specific variables in cas.nc
442          !-----------------------------------------------------------------------
443          if ( 1 == 1 ) THEN
444             ! A MODIFIER: il faudrait dire nudging_temp mais faut le declarer dans compar1d.h etc...       
445             !           if ( name_var(i) == 'temp_nudging' .and. nint(nudging_t)==0) stop 'Nudging inconsistency temp'
446             if ( name_var(i) == 'qv_nud' .and. nint(nudging_qv)==0) stop 'Nudging inconsistency qv'
447             if ( name_var(i) == 'ua_nud' .and. nint(nudging_u)==0) stop 'Nudging inconsistency u'
448             if ( name_var(i) == 'va_nud' .and. nint(nudging_v)==0) stop 'Nudging inconsistency v'
449          ELSE
450             print*,'GUIDAGE : CONSISTENCY CHECK DEACTIVATED FOR TESTS of SANDU/REF'
451          ENDIF
452
453          !-----------------------------------------------------------------------
454          ! Reading variables 1D (N+1) vertical variables (nlevelp1,lat,lon)
455          !-----------------------------------------------------------------------
456          if(i.LE.4) then
457             ierr = NF90_GET_VAR(nid,var3didin(i),apbp)
458             print *,'read_SCM(apbp), on a lu ',i,name_var(i)
459             if(ierr/=nf90_noerr) then
460                print *,'B Pb a la lecture de cas.nc: ',name_var(i)
461                stop "getvarup"
462             endif
463
464             !-----------------------------------------------------------------------
465             !  Reading 1D (N) vertical varialbes    (nlevel,lat,lon)   
466             !-----------------------------------------------------------------------
467          else if(i.gt.4.and.i.LE.12) then 
468             ierr = NF90_GET_VAR(nid,var3didin(i),resul1)
469             print *,'read_SCM(resul1), on a lu ',i,name_var(i)
470             if(ierr/=nf90_noerr) then
471                print *,'C Pb a la lecture de cas.nc: ',name_var(i)
472                stop "getvarup"
473             endif
474             print*,'Lecture de la variable #i ',i,name_var(i),minval(resul1),maxval(resul1)
475
476             !-----------------------------------------------------------------------
477             !  Reading 2D tim-vertical variables  (time,nlevel,lat,lon)
478             !  TBD : seems to be the same as above.
479             !-----------------------------------------------------------------------
480          else if(i.gt.12.and.i.LE.61) then
481             ierr = NF90_GET_VAR(nid,var3didin(i),resul)
482             print *,'read_SCM(resul), on a lu ',i,name_var(i)
483             if(ierr/=nf90_noerr) then
484                print *,'D Pb a la lecture de cas.nc: ',name_var(i)
485                stop "getvarup"
486             endif
487             print*,'Lecture de la variable #i ',i,name_var(i),minval(resul),maxval(resul)
488
489             !-----------------------------------------------------------------------
490             !  Reading 1D time variables (time,lat,lon)
491             !-----------------------------------------------------------------------
492          else if (i.gt.62.and.i.LE.75) then
493             ierr = NF90_GET_VAR(nid,var3didin(i),resul2)
494             print *,'read_SCM(resul2), on a lu ',i,name_var(i)
495             if(ierr/=nf90_noerr) then
496                print *,'E Pb a la lecture de cas.nc: ',name_var(i)
497                stop "getvarup"
498             endif
499             print*,'Lecture de la variable #i  ',i,name_var(i),minval(resul2),maxval(resul2)
500
501             !-----------------------------------------------------------------------
502             ! Reading scalar variables (lat,lon)
503             !-----------------------------------------------------------------------
504          else
505             ierr = NF90_GET_VAR(nid,var3didin(i),resul3)
506             print *,'read_SCM(resul3), on a lu ',i,name_var(i)
507             if(ierr/=nf90_noerr) then
508                print *,'F Pb a la lecture de cas.nc: ',name_var(i)
509                stop "getvarup"
510             endif
511             print*,'Lecture de la variable #i ',i,name_var(i),resul3
512          endif
513       endif
514
515       !-----------------------------------------------------------------------
516       ! Attributing variables
517       !-----------------------------------------------------------------------
518       select case(i)
519          !case(1) ; ap=apbp       ! donnees indexees en nlevel+1
520          ! case(2) ; bp=apbp
521       case(3) ; zzh=apbp
522       case(4) ; pph=apbp
523       case(5) ; temp0=resul1    ! donnees initiales
524       case(6) ; qv0=resul1
525       case(7) ; ql0=resul1
526       case(8) ; qi0=resul1
527       case(9) ; u0=resul1
528       case(10) ; v0=resul1
529       case(11) ; tke0=resul1
530       case(12) ; pp0=resul1
531       case(13) ; vitw=resul    ! donnees indexees en nlevel,time
532       case(14) ; omega=resul
533       case(15) ; ug=resul
534       case(16) ; vg=resul
535       case(17) ; du=resul
536       case(18) ; hu=resul
537       case(19) ; vu=resul
538       case(20) ; dv=resul
539       case(21) ; hv=resul
540       case(22) ; vv=resul
541       case(23) ; dt=resul
542       case(24) ; ht=resul
543       case(25) ; vt=resul
544       case(26) ; dq=resul
545       case(27) ; hq=resul
546       case(28) ; vq=resul
547       case(29) ; dth=resul
548       case(30) ; hth=resul
549       case(31) ; vth=resul
550       case(32) ; hthl=resul
551       case(33) ; dr=resul
552       case(34) ; hr=resul
553       case(35) ; vr=resul
554       case(36) ; dtrad=resul
555       case(37) ; q1=resul
556       case(38) ; q2=resul
557       case(39) ; uw=resul
558       case(40) ; vw=resul
559       case(41) ; rh=resul
560       case(42) ; temp_nudg=resul
561       case(43) ; qv_nudg=resul
562       case(44) ; u_nudg=resul
563       case(45) ; v_nudg=resul
564       case(46) ; zz=resul      ! donnees en time,nlevel pour profil initial
565       case(47) ; pp=resul
566       case(48) ; temp=resul
567       case(49) ; theta=resul
568       case(50) ; thv=resul
569       case(51) ; thl=resul
570       case(52) ; qv=resul
571       case(53) ; ql=resul
572       case(54) ; qi=resul
573       case(55) ; rv=resul
574       case(56) ; u=resul
575       case(57) ; v=resul
576       case(58) ; invtau_temp_nudg=resul
577       case(59) ; invtau_qv_nudg=resul
578       case(60) ; invtau_u_nudg=resul
579       case(61) ; invtau_v_nudg=resul
580       case(62) ; tkes=resul2   ! donnees indexees en time
581       case(63) ; sens=resul2
582       case(64) ; flat=resul2
583       case(65) ; ts=resul2
584       case(66) ; tskin=resul2       
585       case(67) ; ps=resul2
586       case(68) ; ustar=resul2
587       case(69) ; orog_cas=resul3      ! constantes
588       case(70) ; albedo_cas=resul3
589       case(71) ; emiss_cas=resul3
590       case(72) ; q_skin_cas=resul3
591       case(73) ; mom_rough=resul3
592       case(74) ; heat_rough=resul3
593       case(75) ; o3_cas=resul3       
594       case(76) ; rugos_cas=resul3
595       case(77) ; clay_cas=resul3
596       case(78) ; sand_cas=resul3
597       end select
598       resul=0.
599       resul1=0.
600       resul2=0.
601       resul3=0.
602    enddo
603    print*,'Lecture de la variable APRES ,sens ',minval(sens),maxval(sens)
604    print*,'Lecture de la variable APRES ,flat ',minval(flat),maxval(flat)
605
606    !CR:ATTENTION EN ATTENTE DE REGLER LA QUESTION DU PAS DE TEMPS INITIAL
607    do t=1,ntime
608       do k=1,nlevel
609          temp(k,t)=temp0(k)
610          qv(k,t)=qv0(k)
611          ql(k,t)=ql0(k)
612          qi(k,t)=qi0(k)
613          u(k,t)=u0(k)
614          v(k,t)=v0(k)
615          tke(k,t)=tke0(k)
616       enddo
617    enddo
618!!!! TRAVAIL : EN FONCTION DES DECISIONS SUR LA SPECIFICATION DE W
619!!!omega=-vitw*pres*rg/(rd*temp)
620    !-----------------------------------------------------------------------
621
622    return
623  END SUBROUTINE read_SCM
624  !======================================================================
625
626  !======================================================================
627
628  !**********************************************************************************************
629
630  !**********************************************************************************************
631  SUBROUTINE interp_case_time_std(day,day1,annee_ref                           &
632       !    &         ,year_cas,day_cas,nt_cas,pdt_forc,nlev_cas                         &
633       ,nt_cas,nlev_cas                                                   &
634       ,ts_cas,tskin_cas,ps_cas,plev_cas,t_cas,theta_cas,thv_cas,thl_cas            &
635       ,qv_cas,ql_cas,qi_cas,u_cas,v_cas                                  &
636       ,ug_cas,vg_cas,temp_nudg_cas,qv_nudg_cas,u_nudg_cas,v_nudg_cas     &
637       ,invtau_temp_nudg_cas,invtau_qv_nudg_cas,invtau_u_nudg_cas,invtau_v_nudg_cas     &
638       ,vitw_cas,omega_cas,tke_cas,du_cas,hu_cas,vu_cas             &
639       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas               &
640       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas                      &
641       ,lat_cas,sens_cas,ustar_cas                                        &
642       ,uw_cas,vw_cas,q1_cas,q2_cas,tkes_cas                               &
643       !
644       ,ts_prof_cas,tskin_prof_cas,ps_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas               &
645       ,thv_prof_cas,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas     &
646       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas                     &
647       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas     &
648       ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas     &     
649       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas  &
650       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas                   &
651       ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas                &
652       ,hq_prof_cas,vq_prof_cas,dth_prof_cas,hth_prof_cas,vth_prof_cas    &
653       ,lat_prof_cas,sens_prof_cas                                        &
654       ,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tkes_prof_cas)
655
656
657
658
659
660
661    implicit none
662
663    !---------------------------------------------------------------------------------------
664    ! Time interpolation of a 2D field to the timestep corresponding to day
665    !
666    ! day: current julian day (e.g. 717538.2)
667    ! day1: first day of the simulation
668    ! nt_cas: total nb of data in the forcing
669    ! pdt_cas: total time interval (in sec) between 2 forcing data
670    !---------------------------------------------------------------------------------------
671
672    INCLUDE "compar1d.h"
673    INCLUDE "date_cas.h"
674
675    ! inputs:
676    integer annee_ref
677    integer nt_cas,nlev_cas
678    real day, day1,day_cas
679    real ts_cas(nt_cas),tskin_cas(nt_cas),ps_cas(nt_cas)
680    real plev_cas(nlev_cas,nt_cas)
681    real t_cas(nlev_cas,nt_cas),theta_cas(nlev_cas,nt_cas)
682    real thv_cas(nlev_cas,nt_cas), thl_cas(nlev_cas,nt_cas)
683    real qv_cas(nlev_cas,nt_cas),ql_cas(nlev_cas,nt_cas),qi_cas(nlev_cas,nt_cas)
684    real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
685    real ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
686    real temp_nudg_cas(nlev_cas,nt_cas),qv_nudg_cas(nlev_cas,nt_cas)
687    real u_nudg_cas(nlev_cas,nt_cas),v_nudg_cas(nlev_cas,nt_cas)
688
689    real invtau_temp_nudg_cas(nlev_cas,nt_cas),invtau_qv_nudg_cas(nlev_cas,nt_cas)
690    real invtau_u_nudg_cas(nlev_cas,nt_cas),invtau_v_nudg_cas(nlev_cas,nt_cas)
691
692    real vitw_cas(nlev_cas,nt_cas),omega_cas(nlev_cas,nt_cas),tke_cas(nlev_cas,nt_cas)
693    real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
694    real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
695    real dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
696    real dth_cas(nlev_cas,nt_cas),hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas)
697    real dtrad_cas(nlev_cas,nt_cas)
698    real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
699    real lat_cas(nt_cas),sens_cas(nt_cas),tkes_cas(nt_cas)
700    real ustar_cas(nt_cas),uw_cas(nlev_cas,nt_cas),vw_cas(nlev_cas,nt_cas)
701    real q1_cas(nlev_cas,nt_cas),q2_cas(nlev_cas,nt_cas)
702
703    ! outputs:
704    real plev_prof_cas(nlev_cas)
705    real t_prof_cas(nlev_cas),theta_prof_cas(nlev_cas),thl_prof_cas(nlev_cas),thv_prof_cas(nlev_cas)
706    real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
707    real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
708    real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas)
709    real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
710    real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
711
712    real invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
713    real invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas)
714
715    real vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
716    real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
717    real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
718    real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
719    real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
720    real dtrad_prof_cas(nlev_cas)
721    real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
722    real lat_prof_cas,sens_prof_cas,tkes_prof_cas,ts_prof_cas,tskin_prof_cas,ps_prof_cas,ustar_prof_cas
723    real uw_prof_cas(nlev_cas),vw_prof_cas(nlev_cas),q1_prof_cas(nlev_cas),q2_prof_cas(nlev_cas)
724    ! local:
725    integer it_cas1, it_cas2,k
726    real timeit,time_cas1,time_cas2,frac
727
728
729    print*,'Check time',day1,day_ju_ini_cas,day_deb+1,pdt_cas
730    !       do k=1,nlev_cas
731    !       print*,'debut de interp2_case_time, plev_cas=',k,plev_cas(k,1)
732    !       enddo
733
734    ! On teste si la date du cas AMMA est correcte.
735    ! C est pour memoire car en fait les fichiers .def
736    ! sont censes etre corrects.
737    ! A supprimer a terme (MPL 20150623)
738    !     if ((forcing_type.eq.10).and.(1.eq.0)) then
739    ! Check that initial day of the simulation consistent with AMMA case:
740    !      if (annee_ref.ne.2006) then
741    !       print*,'Pour AMMA, annee_ref doit etre 2006'
742    !       print*,'Changer annee_ref dans run.def'
743    !       stop
744    !      endif
745    !      if (annee_ref.eq.2006 .and. day1.lt.day_cas) then
746    !       print*,'AMMA a debute le 10 juillet 2006',day1,day_cas
747    !       print*,'Changer dayref dans run.def'
748    !       stop
749    !      endif
750    !      if (annee_ref.eq.2006 .and. day1.gt.day_cas+1) then
751    !       print*,'AMMA a fini le 11 juillet'
752    !       print*,'Changer dayref ou nday dans run.def'
753    !       stop
754    !      endif
755    !      endif
756
757    ! Determine timestep relative to the 1st day:
758    !       timeit=(day-day1)*86400.
759    !       if (annee_ref.eq.1992) then
760    !        timeit=(day-day_cas)*86400.
761    !       else
762    !        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
763    !       endif
764    timeit=(day-day_ju_ini_cas)*86400
765    print *,'day=',day
766    print *,'day_ju_ini_cas=',day_ju_ini_cas
767    print *,'pdt_cas=',pdt_cas
768    print *,'timeit=',timeit
769    print *,'nt_cas=',nt_cas
770
771    ! Determine the closest observation times:
772    !       it_cas1=INT(timeit/pdt_cas)+1
773    !       it_cas2=it_cas1 + 1
774    !       time_cas1=(it_cas1-1)*pdt_cas
775    !       time_cas2=(it_cas2-1)*pdt_cas
776
777    it_cas1=INT(timeit/pdt_cas)+1
778    IF (it_cas1 .EQ. nt_cas) THEN
779       it_cas2=it_cas1
780    ELSE
781       it_cas2=it_cas1 + 1
782    ENDIF
783    time_cas1=(it_cas1-1)*pdt_cas
784    time_cas2=(it_cas2-1)*pdt_cas
785    !     print *,'timeit,pdt_cas,nt_cas=',timeit,pdt_cas,nt_cas
786    !     print *,'it_cas1,it_cas2,time_cas1,time_cas2=',it_cas1,it_cas2,time_cas1,time_cas2
787
788    if (it_cas1 .gt. nt_cas) then
789       write(*,*) 'PB-stop: day, day_ju_ini_cas,it_cas1, it_cas2, timeit: '            &
790            ,day,day_ju_ini_cas,it_cas1,it_cas2,timeit
791       stop
792    endif
793
794    ! time interpolation:
795    IF (it_cas1 .EQ. it_cas2) THEN
796       frac=0.
797    ELSE
798       frac=(time_cas2-timeit)/(time_cas2-time_cas1)
799       frac=max(frac,0.0)
800    ENDIF
801
802    lat_prof_cas = lat_cas(it_cas2)                                   &
803         -frac*(lat_cas(it_cas2)-lat_cas(it_cas1))
804    sens_prof_cas = sens_cas(it_cas2)                                 &
805         -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
806    tkes_prof_cas = tkes_cas(it_cas2)                                   &
807         -frac*(tkes_cas(it_cas2)-tkes_cas(it_cas1))
808    ts_prof_cas = ts_cas(it_cas2)                                     &
809         -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
810    tskin_prof_cas = tskin_cas(it_cas2)                                     &
811         -frac*(tskin_cas(it_cas2)-tskin_cas(it_cas1))
812    ps_prof_cas = ps_cas(it_cas2)                                     &
813         -frac*(ps_cas(it_cas2)-ps_cas(it_cas1))
814    ustar_prof_cas = ustar_cas(it_cas2)                               &
815         -frac*(ustar_cas(it_cas2)-ustar_cas(it_cas1))
816
817    do k=1,nlev_cas
818       plev_prof_cas(k) = plev_cas(k,it_cas2)                           &     
819            -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
820       t_prof_cas(k) = t_cas(k,it_cas2)                                 &       
821            -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
822       !print *,'k,frac,plev_cas1,plev_cas2=',k,frac,plev_cas(k,it_cas1),plev_cas(k,it_cas2)
823       theta_prof_cas(k) = theta_cas(k,it_cas2)                         &                     
824            -frac*(theta_cas(k,it_cas2)-theta_cas(k,it_cas1))
825       thv_prof_cas(k) = thv_cas(k,it_cas2)                             &         
826            -frac*(thv_cas(k,it_cas2)-thv_cas(k,it_cas1))
827       thl_prof_cas(k) = thl_cas(k,it_cas2)                             &             
828            -frac*(thl_cas(k,it_cas2)-thl_cas(k,it_cas1))
829       qv_prof_cas(k) = qv_cas(k,it_cas2)                               &
830            -frac*(qv_cas(k,it_cas2)-qv_cas(k,it_cas1))
831       ql_prof_cas(k) = ql_cas(k,it_cas2)                               &
832            -frac*(ql_cas(k,it_cas2)-ql_cas(k,it_cas1))
833       qi_prof_cas(k) = qi_cas(k,it_cas2)                               &
834            -frac*(qi_cas(k,it_cas2)-qi_cas(k,it_cas1))
835       u_prof_cas(k) = u_cas(k,it_cas2)                                 &
836            -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
837       v_prof_cas(k) = v_cas(k,it_cas2)                                 &
838            -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
839       ug_prof_cas(k) = ug_cas(k,it_cas2)                               &
840            -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
841       vg_prof_cas(k) = vg_cas(k,it_cas2)                               &
842            -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
843       temp_nudg_prof_cas(k) = temp_nudg_cas(k,it_cas2)                    &
844            -frac*(temp_nudg_cas(k,it_cas2)-temp_nudg_cas(k,it_cas1))
845       qv_nudg_prof_cas(k) = qv_nudg_cas(k,it_cas2)                        &
846            -frac*(qv_nudg_cas(k,it_cas2)-qv_nudg_cas(k,it_cas1))
847       u_nudg_prof_cas(k) = u_nudg_cas(k,it_cas2)                          &
848            -frac*(u_nudg_cas(k,it_cas2)-u_nudg_cas(k,it_cas1))
849       v_nudg_prof_cas(k) = v_nudg_cas(k,it_cas2)                          &
850            -frac*(v_nudg_cas(k,it_cas2)-v_nudg_cas(k,it_cas1))
851       invtau_temp_nudg_prof_cas(k) = invtau_temp_nudg_cas(k,it_cas2)                    &
852            -frac*(invtau_temp_nudg_cas(k,it_cas2)-invtau_temp_nudg_cas(k,it_cas1))
853       invtau_qv_nudg_prof_cas(k) = invtau_qv_nudg_cas(k,it_cas2)                        &
854            -frac*(invtau_qv_nudg_cas(k,it_cas2)-invtau_qv_nudg_cas(k,it_cas1))
855       invtau_u_nudg_prof_cas(k) = invtau_u_nudg_cas(k,it_cas2)                          &
856            -frac*(invtau_u_nudg_cas(k,it_cas2)-invtau_u_nudg_cas(k,it_cas1))
857       invtau_v_nudg_prof_cas(k) = invtau_v_nudg_cas(k,it_cas2)                          &
858            -frac*(invtau_v_nudg_cas(k,it_cas2)-invtau_v_nudg_cas(k,it_cas1))
859       vitw_prof_cas(k) = vitw_cas(k,it_cas2)                           &
860            -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
861       omega_prof_cas(k) = omega_cas(k,it_cas2)                         &
862            -frac*(omega_cas(k,it_cas2)-omega_cas(k,it_cas1))
863       tke_prof_cas(k) = tke_cas(k,it_cas2)                         &
864            -frac*(tke_cas(k,it_cas2)-tke_cas(k,it_cas1))
865       du_prof_cas(k) = du_cas(k,it_cas2)                               &
866            -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
867       hu_prof_cas(k) = hu_cas(k,it_cas2)                               &
868            -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
869       vu_prof_cas(k) = vu_cas(k,it_cas2)                               &
870            -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
871       dv_prof_cas(k) = dv_cas(k,it_cas2)                               &
872            -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
873       hv_prof_cas(k) = hv_cas(k,it_cas2)                               &
874            -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
875       vv_prof_cas(k) = vv_cas(k,it_cas2)                               &
876            -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
877       dt_prof_cas(k) = dt_cas(k,it_cas2)                               &
878            -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
879       ht_prof_cas(k) = ht_cas(k,it_cas2)                               &
880            -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
881       vt_prof_cas(k) = vt_cas(k,it_cas2)                               &
882            -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
883       dth_prof_cas(k) = dth_cas(k,it_cas2)                             &
884            -frac*(dth_cas(k,it_cas2)-dth_cas(k,it_cas1))
885       hth_prof_cas(k) = hth_cas(k,it_cas2)                             &
886            -frac*(hth_cas(k,it_cas2)-hth_cas(k,it_cas1))
887       vth_prof_cas(k) = vth_cas(k,it_cas2)                             &
888            -frac*(vth_cas(k,it_cas2)-vth_cas(k,it_cas1))
889       dtrad_prof_cas(k) = dtrad_cas(k,it_cas2)                         &
890            -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
891       dq_prof_cas(k) = dq_cas(k,it_cas2)                               &
892            -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
893       hq_prof_cas(k) = hq_cas(k,it_cas2)                               &
894            -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
895       vq_prof_cas(k) = vq_cas(k,it_cas2)                               &
896            -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
897       uw_prof_cas(k) = uw_cas(k,it_cas2)                                &
898            -frac*(uw_cas(k,it_cas2)-uw_cas(k,it_cas1))
899       vw_prof_cas(k) = vw_cas(k,it_cas2)                                &
900            -frac*(vw_cas(k,it_cas2)-vw_cas(k,it_cas1))
901       q1_prof_cas(k) = q1_cas(k,it_cas2)                                &
902            -frac*(q1_cas(k,it_cas2)-q1_cas(k,it_cas1))
903       q2_prof_cas(k) = q2_cas(k,it_cas2)                                &
904            -frac*(q2_cas(k,it_cas2)-q2_cas(k,it_cas1))
905    enddo
906
907    return
908  END SUBROUTINE interp_case_time_std
909
910  !**********************************************************************************************
911  !=====================================================================
912  SUBROUTINE interp2_case_vertical_std(play,plev,nlev_cas,plev_prof_cas                           &
913       ,t_prof_cas,th_prof_cas,thv_prof_cas,thl_prof_cas                                       &
914       ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                              &
915       ,ug_prof_cas,vg_prof_cas                                                                &
916       ,temp_nudg_prof_cas,qv_nudg_prof_cas,u_nudg_prof_cas,v_nudg_prof_cas                    &
917       ,invtau_temp_nudg_prof_cas,invtau_qv_nudg_prof_cas,invtau_u_nudg_prof_cas,invtau_v_nudg_prof_cas &     
918       ,vitw_prof_cas,omega_prof_cas,tke_prof_cas                                              &
919       ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                &
920       ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &
921       ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                 &
922       !
923       ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas                                        &
924       ,qv_mod_cas,ql_mod_cas,qi_mod_cas,u_mod_cas,v_mod_cas                                   &
925       ,ug_mod_cas,vg_mod_cas                                                                  &
926       ,temp_nudg_mod_cas,qv_nudg_mod_cas,u_nudg_mod_cas,v_nudg_mod_cas                        &
927       ,invtau_temp_nudg_mod_cas,invtau_qv_nudg_mod_cas,invtau_u_nudg_mod_cas,invtau_v_nudg_mod_cas                        &     
928       ,w_mod_cas,omega_mod_cas,tke_mod_cas                                                    &
929       ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                      &
930       ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
931       ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
932
933    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
934USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO                   &
935          , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA                  &
936          , R_ecc, R_peri, R_incl                                      &
937          , RA, RG, R1SA                                         &
938          , RSIGMA                                                     &
939          , R, RMD, RMV, RD, RV, RCPD                    &
940          , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12        &
941          , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w                    &
942          , RCW, RCS                                                 &
943          , RLVTT, RLSTT, RLMLT, RTT, RATM                           &
944          , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS            &
945          , RALPD, RBETD, RGAMD
946implicit none
947
948
949
950
951    !-------------------------------------------------------------------------
952    ! Vertical interpolation of generic case forcing data onto mod_casel levels
953    !-------------------------------------------------------------------------
954
955    integer nlevmax
956    parameter (nlevmax=41)
957    integer nlev_cas,mxcalc
958    !       real play(llm), plev_prof(nlevmax)
959    !       real t_prof(nlevmax),q_prof(nlevmax)
960    !       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
961    !       real ht_prof(nlevmax),vt_prof(nlevmax)
962    !       real hq_prof(nlevmax),vq_prof(nlevmax)
963
964    real play(llm), plev(llm+1), plev_prof_cas(nlev_cas)
965    real t_prof_cas(nlev_cas),th_prof_cas(nlev_cas),thv_prof_cas(nlev_cas),thl_prof_cas(nlev_cas)
966    real qv_prof_cas(nlev_cas),ql_prof_cas(nlev_cas),qi_prof_cas(nlev_cas)
967    real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
968    real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas),omega_prof_cas(nlev_cas),tke_prof_cas(nlev_cas)
969    real temp_nudg_prof_cas(nlev_cas),qv_nudg_prof_cas(nlev_cas)
970    real u_nudg_prof_cas(nlev_cas),v_nudg_prof_cas(nlev_cas)
971    real invtau_temp_nudg_prof_cas(nlev_cas),invtau_qv_nudg_prof_cas(nlev_cas)
972    real invtau_u_nudg_prof_cas(nlev_cas),invtau_v_nudg_prof_cas(nlev_cas)
973
974    real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
975    real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
976    real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
977    real dth_prof_cas(nlev_cas),hth_prof_cas(nlev_cas),vth_prof_cas(nlev_cas)
978    real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
979
980    real t_mod_cas(llm),theta_mod_cas(llm),thv_mod_cas(llm),thl_mod_cas(llm)
981    real qv_mod_cas(llm),ql_mod_cas(llm),qi_mod_cas(llm)
982    real u_mod_cas(llm),v_mod_cas(llm)
983    real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm),omega_mod_cas(llm),tke_mod_cas(llm+1)
984    real temp_nudg_mod_cas(llm),qv_nudg_mod_cas(llm)
985    real u_nudg_mod_cas(llm),v_nudg_mod_cas(llm)
986    real invtau_temp_nudg_mod_cas(llm),invtau_qv_nudg_mod_cas(llm)
987    real invtau_u_nudg_mod_cas(llm),invtau_v_nudg_mod_cas(llm)
988    real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
989    real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
990    real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
991    real dth_mod_cas(llm),hth_mod_cas(llm),vth_mod_cas(llm)
992    real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
993
994    integer l,k,k1,k2
995    real frac,frac1,frac2,fact
996
997
998
999    ! for variables defined at the middle of layers
1000
1001    do l = 1, llm
1002
1003       if (play(l).ge.plev_prof_cas(nlev_cas)) then
1004
1005          mxcalc=l
1006          !        print *,'debut interp2, mxcalc=',mxcalc
1007          k1=0
1008          k2=0
1009
1010          if (play(l).le.plev_prof_cas(1)) then
1011
1012             do k = 1, nlev_cas-1
1013                if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
1014                   k1=k
1015                   k2=k+1
1016                endif
1017             enddo
1018
1019             if (k1.eq.0 .or. k2.eq.0) then
1020                write(*,*) 'PB! k1, k2 = ',k1,k2
1021                write(*,*) 'l,play(l) = ',l,play(l)/100
1022                do k = 1, nlev_cas-1
1023                   write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
1024                enddo
1025             endif
1026
1027
1028
1029             frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
1030
1031             t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
1032             theta_mod_cas(l)= th_prof_cas(k2) - frac*(th_prof_cas(k2)-th_prof_cas(k1))
1033             if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
1034             thv_mod_cas(l)= thv_prof_cas(k2) - frac*(thv_prof_cas(k2)-thv_prof_cas(k1))
1035             thl_mod_cas(l)= thl_prof_cas(k2) - frac*(thl_prof_cas(k2)-thl_prof_cas(k1))
1036             qv_mod_cas(l)= qv_prof_cas(k2) - frac*(qv_prof_cas(k2)-qv_prof_cas(k1))
1037             ql_mod_cas(l)= ql_prof_cas(k2) - frac*(ql_prof_cas(k2)-ql_prof_cas(k1))
1038             qi_mod_cas(l)= qi_prof_cas(k2) - frac*(qi_prof_cas(k2)-qi_prof_cas(k1))
1039             u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
1040             v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
1041             ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
1042             vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
1043             temp_nudg_mod_cas(l)= temp_nudg_prof_cas(k2) - frac*(temp_nudg_prof_cas(k2)-temp_nudg_prof_cas(k1))
1044             qv_nudg_mod_cas(l)= qv_nudg_prof_cas(k2) - frac*(qv_nudg_prof_cas(k2)-qv_nudg_prof_cas(k1))
1045             u_nudg_mod_cas(l)= u_nudg_prof_cas(k2) - frac*(u_nudg_prof_cas(k2)-u_nudg_prof_cas(k1))
1046             v_nudg_mod_cas(l)= v_nudg_prof_cas(k2) - frac*(v_nudg_prof_cas(k2)-v_nudg_prof_cas(k1))
1047
1048             invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(k2) &
1049                  - frac*(invtau_temp_nudg_prof_cas(k2)-invtau_temp_nudg_prof_cas(k1))
1050             invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(k2) - frac*(invtau_qv_nudg_prof_cas(k2)-invtau_qv_nudg_prof_cas(k1))
1051             invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(k2) - frac*(invtau_u_nudg_prof_cas(k2)-invtau_u_nudg_prof_cas(k1))
1052             invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(k2) - frac*(invtau_v_nudg_prof_cas(k2)-invtau_v_nudg_prof_cas(k1))
1053
1054             w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
1055             omega_mod_cas(l)= omega_prof_cas(k2) - frac*(omega_prof_cas(k2)-omega_prof_cas(k1))
1056             du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
1057             hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
1058             vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
1059             dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
1060             hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
1061             vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
1062             dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
1063             ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
1064             vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
1065             dth_mod_cas(l)= dth_prof_cas(k2) - frac*(dth_prof_cas(k2)-dth_prof_cas(k1))
1066             hth_mod_cas(l)= hth_prof_cas(k2) - frac*(hth_prof_cas(k2)-hth_prof_cas(k1))
1067             vth_mod_cas(l)= vth_prof_cas(k2) - frac*(vth_prof_cas(k2)-vth_prof_cas(k1))
1068             dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
1069             hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
1070             vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
1071             dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
1072
1073          else !play>plev_prof_cas(1)
1074
1075             k1=1
1076             k2=2
1077             print *,'interp2_vert, k1,k2=',plev_prof_cas(k1),plev_prof_cas(k2)
1078             frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1079             frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1080             t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
1081             theta_mod_cas(l)= frac1*th_prof_cas(k1) - frac2*th_prof_cas(k2)
1082             if(theta_mod_cas(l).NE.0) t_mod_cas(l)= theta_mod_cas(l)*(play(l)/100000.)**(RD/RCPD)
1083             thv_mod_cas(l)= frac1*thv_prof_cas(k1) - frac2*thv_prof_cas(k2)
1084             thl_mod_cas(l)= frac1*thl_prof_cas(k1) - frac2*thl_prof_cas(k2)
1085             qv_mod_cas(l)= frac1*qv_prof_cas(k1) - frac2*qv_prof_cas(k2)
1086             ql_mod_cas(l)= frac1*ql_prof_cas(k1) - frac2*ql_prof_cas(k2)
1087             qi_mod_cas(l)= frac1*qi_prof_cas(k1) - frac2*qi_prof_cas(k2)
1088             u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
1089             v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
1090             ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
1091             vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
1092             temp_nudg_mod_cas(l)= frac1*temp_nudg_prof_cas(k1) - frac2*temp_nudg_prof_cas(k2)
1093             qv_nudg_mod_cas(l)= frac1*qv_nudg_prof_cas(k1) - frac2*qv_nudg_prof_cas(k2)
1094             u_nudg_mod_cas(l)= frac1*u_nudg_prof_cas(k1) - frac2*u_nudg_prof_cas(k2)
1095             v_nudg_mod_cas(l)= frac1*v_nudg_prof_cas(k1) - frac2*v_nudg_prof_cas(k2)
1096
1097             invtau_temp_nudg_mod_cas(l)= frac1*invtau_temp_nudg_prof_cas(k1) - frac2*invtau_temp_nudg_prof_cas(k2)
1098             invtau_qv_nudg_mod_cas(l)= frac1*invtau_qv_nudg_prof_cas(k1) - frac2*invtau_qv_nudg_prof_cas(k2)
1099             invtau_u_nudg_mod_cas(l)= frac1*invtau_u_nudg_prof_cas(k1) - frac2*invtau_u_nudg_prof_cas(k2)
1100             invtau_v_nudg_mod_cas(l)= frac1*invtau_v_nudg_prof_cas(k1) - frac2*invtau_v_nudg_prof_cas(k2)
1101
1102             w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
1103             omega_mod_cas(l)= frac1*omega_prof_cas(k1) - frac2*omega_prof_cas(k2)
1104             du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
1105             hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
1106             vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
1107             dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
1108             hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
1109             vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
1110             dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
1111             ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
1112             vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
1113             dth_mod_cas(l)= frac1*dth_prof_cas(k1) - frac2*dth_prof_cas(k2)
1114             hth_mod_cas(l)= frac1*hth_prof_cas(k1) - frac2*hth_prof_cas(k2)
1115             vth_mod_cas(l)= frac1*vth_prof_cas(k1) - frac2*vth_prof_cas(k2)
1116             dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
1117             hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
1118             vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
1119             dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
1120
1121          endif ! play.le.plev_prof_cas(1)
1122
1123       else ! above max altitude of forcing file
1124
1125          !jyg
1126          fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
1127          fact = max(fact,0.)                                           !jyg
1128          fact = exp(-fact)                                             !jyg
1129          t_mod_cas(l)= t_prof_cas(nlev_cas)                            !jyg
1130          theta_mod_cas(l)= th_prof_cas(nlev_cas)                       !jyg
1131          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
1132          thl_mod_cas(l)= thl_prof_cas(nlev_cas)                        !jyg
1133          qv_mod_cas(l)= qv_prof_cas(nlev_cas)*fact                     !jyg
1134          ql_mod_cas(l)= ql_prof_cas(nlev_cas)*fact                     !jyg
1135          qi_mod_cas(l)= qi_prof_cas(nlev_cas)*fact                     !jyg
1136          u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                       !jyg
1137          v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                       !jyg
1138          ug_mod_cas(l)= ug_prof_cas(nlev_cas)                          !jyg
1139          vg_mod_cas(l)= vg_prof_cas(nlev_cas)                          !jyg
1140          temp_nudg_mod_cas(l)= temp_nudg_prof_cas(nlev_cas)            !jyg
1141          qv_nudg_mod_cas(l)= qv_nudg_prof_cas(nlev_cas)                !jyg
1142          u_nudg_mod_cas(l)= u_nudg_prof_cas(nlev_cas)                  !jyg
1143          v_nudg_mod_cas(l)= v_nudg_prof_cas(nlev_cas)                  !jyg
1144
1145          invtau_temp_nudg_mod_cas(l)= invtau_temp_nudg_prof_cas(nlev_cas)            !jyg
1146          invtau_qv_nudg_mod_cas(l)= invtau_qv_nudg_prof_cas(nlev_cas)                !jyg
1147          invtau_u_nudg_mod_cas(l)= invtau_u_nudg_prof_cas(nlev_cas)                  !jyg
1148          invtau_v_nudg_mod_cas(l)= invtau_v_nudg_prof_cas(nlev_cas)                  !jyg
1149
1150          thv_mod_cas(l)= thv_prof_cas(nlev_cas)                        !jyg
1151          w_mod_cas(l)= 0.0                                             !jyg
1152          omega_mod_cas(l)= 0.0                                         !jyg
1153          du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
1154          hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                     !jyg
1155          vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                     !jyg
1156          dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
1157          hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                     !jyg
1158          vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                     !jyg
1159          dt_mod_cas(l)= dt_prof_cas(nlev_cas)
1160          ht_mod_cas(l)= ht_prof_cas(nlev_cas)                          !jyg
1161          vt_mod_cas(l)= vt_prof_cas(nlev_cas)                          !jyg
1162          dth_mod_cas(l)= dth_prof_cas(nlev_cas)
1163          hth_mod_cas(l)= hth_prof_cas(nlev_cas)                        !jyg
1164          vth_mod_cas(l)= vth_prof_cas(nlev_cas)                        !jyg
1165          dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
1166          hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
1167          vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
1168          dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
1169
1170       endif ! play
1171
1172    enddo ! l
1173
1174    ! for variables defined at layer interfaces (EV):
1175
1176
1177    do l = 1, llm+1
1178
1179       if (plev(l).ge.plev_prof_cas(nlev_cas)) then
1180
1181          mxcalc=l
1182          k1=0
1183          k2=0
1184
1185          if (plev(l).le.plev_prof_cas(1)) then
1186
1187             do k = 1, nlev_cas-1
1188                if (plev(l).le.plev_prof_cas(k).and. plev(l).gt.plev_prof_cas(k+1)) then
1189                   k1=k
1190                   k2=k+1
1191                endif
1192             enddo
1193
1194             if (k1.eq.0 .or. k2.eq.0) then
1195                write(*,*) 'PB! k1, k2 = ',k1,k2
1196                write(*,*) 'l,plev(l) = ',l,plev(l)/100
1197                do k = 1, nlev_cas-1
1198                   write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
1199                enddo
1200             endif
1201
1202             frac = (plev_prof_cas(k2)-plev(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
1203             tke_mod_cas(l)= tke_prof_cas(k2) - frac*(tke_prof_cas(k2)-tke_prof_cas(k1))
1204          else !play>plev_prof_cas(1)
1205             k1=1
1206             k2=2
1207             frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1208             frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
1209             tke_mod_cas(l)= frac1*tke_prof_cas(k1) - frac2*tke_prof_cas(k2)
1210
1211          endif ! plev.le.plev_prof_cas(1)
1212
1213       else ! above max altitude of forcing file
1214
1215          tke_mod_cas(l)=0.0
1216
1217       endif ! plev
1218
1219    enddo ! l
1220
1221
1222
1223    return
1224  end SUBROUTINE interp2_case_vertical_std
1225  !*****************************************************************************
1226
1227
1228
1229
1230
1231END MODULE mod_1D_cases_read_std
Note: See TracBrowser for help on using the repository browser.