source: LMDZ6/trunk/libf/phylmd/dyn1d/mod_1D_amma_read.f90 @ 5309

Last change on this file since 5309 was 5302, checked in by abarral, 8 days ago

Turn compar1d.h date_cas.h into module
Move fcg_racmo.h to obsolete

File size: 14.5 KB
Line 
1MODULE mod_1D_amma_read
2  USE netcdf, ONLY: nf90_get_var, nf90_open, nf90_noerr, nf90_open, nf90_nowrite, &
3          nf90_inq_dimid, nf90_inquire_dimension, nf90_strerror, nf90_inq_varid
4!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5!Declarations specifiques au cas AMMA
6        character*80 :: fich_amma
7! Option du cas AMMA ou on impose la discretisation verticale (Ap,Bp)
8        integer nlev_amma, nt_amma
9
10
11        integer year_ini_amma, day_ini_amma, mth_ini_amma
12        real heure_ini_amma
13        real day_ju_ini_amma   ! Julian day of amma first day
14        parameter (year_ini_amma=2006)
15        parameter (mth_ini_amma=7)
16        parameter (day_ini_amma=10)  ! 10 = 10Juil2006
17        parameter (heure_ini_amma=0.) !0h en secondes
18        real dt_amma
19        parameter (dt_amma=1800.)
20
21!profils initiaux:
22        real, allocatable::  plev_amma(:)
23
24        real, allocatable::  z_amma(:)
25        real, allocatable::  th_amma(:),q_amma(:)
26        real, allocatable::  u_amma(:)
27        real, allocatable::  v_amma(:)
28
29        real, allocatable::  th_ammai(:),q_ammai(:)
30        real, allocatable::  u_ammai(:)
31        real, allocatable::  v_ammai(:)
32        real, allocatable::  vitw_ammai(:)
33        real, allocatable::  ht_ammai(:)
34        real, allocatable::  hq_ammai(:)
35        real, allocatable::  vt_ammai(:)
36        real, allocatable::  vq_ammai(:)
37
38!forcings
39        real, allocatable::  ht_amma(:,:)
40        real, allocatable::  hq_amma(:,:)
41        real, allocatable::  vitw_amma(:,:)
42        real, allocatable::  lat_amma(:),sens_amma(:)
43
44!champs interpoles
45        real, allocatable::  vitw_profamma(:)
46        real, allocatable::  ht_profamma(:)
47        real, allocatable::  hq_profamma(:)
48        real lat_profamma,sens_profamma
49        real, allocatable::  vt_profamma(:)
50        real, allocatable::  vq_profamma(:)
51        real, allocatable::  th_profamma(:)
52        real, allocatable::  q_profamma(:)
53        real, allocatable::  u_profamma(:)
54        real, allocatable::  v_profamma(:)
55
56
57CONTAINS
58
59SUBROUTINE read_1D_cases
60      implicit none
61
62      INTEGER nid,rid,ierr
63
64      fich_amma='amma.nc'
65      print*,'fich_amma ',fich_amma
66      ierr = nf90_open(fich_amma,nf90_nowrite,nid)
67      print*,'fich_amma,nf90_nowrite,nid ',fich_amma,nf90_nowrite,nid
68      if (ierr.NE.nf90_noerr) then
69         write(*,*) 'ERROR: GROS Pb opening forcings nc file '
70         write(*,*) nf90_strerror(ierr)
71         stop ""
72      endif
73!.......................................................................
74      ierr=nf90_inq_dimid(nid,'lev',rid)
75      IF (ierr.NE.nf90_noerr) THEN
76         print*, 'Oh probleme lecture dimension zz'
77      ENDIF
78      ierr=nf90_inquire_dimension(nid,rid,len=nlev_amma)
79      print*,'OK nid,rid,nlev_amma',nid,rid,nlev_amma
80!.......................................................................
81      ierr=nf90_inq_dimid(nid,'time',rid)
82      print*,'nid,rid',nid,rid
83      nt_amma=0
84      IF (ierr.NE.nf90_noerr) THEN
85        stop 'probleme lecture dimension sens'
86      ENDIF
87      ierr=nf90_inquire_dimension(nid,rid,len=nt_amma)
88      print*,'nid,rid,nlev_amma',nid,rid,nt_amma
89
90!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
91!profils initiaux:
92        allocate(plev_amma(nlev_amma))
93       
94        allocate(z_amma(nlev_amma))
95        allocate(th_amma(nlev_amma),q_amma(nlev_amma))
96        allocate(u_amma(nlev_amma))
97        allocate(v_amma(nlev_amma))
98
99!forcings
100        allocate(ht_amma(nlev_amma,nt_amma))
101        allocate(hq_amma(nlev_amma,nt_amma))
102        allocate(vitw_amma(nlev_amma,nt_amma))
103        allocate(lat_amma(nt_amma),sens_amma(nt_amma))
104
105!profils initiaux:
106        allocate(th_ammai(nlev_amma),q_ammai(nlev_amma))
107        allocate(u_ammai(nlev_amma))
108        allocate(v_ammai(nlev_amma))
109        allocate(vitw_ammai(nlev_amma) )
110        allocate(ht_ammai(nlev_amma))
111        allocate(hq_ammai(nlev_amma))
112        allocate(vt_ammai(nlev_amma))
113        allocate(vq_ammai(nlev_amma))
114
115!champs interpoles
116        allocate(vitw_profamma(nlev_amma))
117        allocate(ht_profamma(nlev_amma))
118        allocate(hq_profamma(nlev_amma))
119        allocate(vt_profamma(nlev_amma))
120        allocate(vq_profamma(nlev_amma))
121        allocate(th_profamma(nlev_amma))
122        allocate(q_profamma(nlev_amma))
123        allocate(u_profamma(nlev_amma))
124        allocate(v_profamma(nlev_amma))
125
126        print*,'Allocations OK'
127        call read_amma(nid,nlev_amma,nt_amma                                  &
128     &     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma         &
129     &     ,ht_amma,hq_amma,sens_amma,lat_amma)
130
131END SUBROUTINE read_1D_cases
132
133
134
135!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
136SUBROUTINE deallocate_1D_cases
137!profils initiaux:
138        deallocate(plev_amma)
139       
140        deallocate(z_amma)
141        deallocate(th_amma,q_amma)
142        deallocate(u_amma)
143        deallocate(v_amma)
144
145        deallocate(th_ammai,q_ammai)
146        deallocate(u_ammai)
147        deallocate(v_ammai)
148        deallocate(vitw_ammai )
149        deallocate(ht_ammai)
150        deallocate(hq_ammai)
151        deallocate(vt_ammai)
152        deallocate(vq_ammai)
153       
154!forcings
155        deallocate(ht_amma)
156        deallocate(hq_amma)
157        deallocate(vitw_amma)
158        deallocate(lat_amma,sens_amma)
159
160!champs interpoles
161        deallocate(vitw_profamma)
162        deallocate(ht_profamma)
163        deallocate(hq_profamma)
164        deallocate(vt_profamma)
165        deallocate(vq_profamma)
166        deallocate(th_profamma)
167        deallocate(q_profamma)
168        deallocate(u_profamma)
169        deallocate(v_profamma)
170END SUBROUTINE deallocate_1D_cases
171
172
173END MODULE mod_1D_amma_read
174!=====================================================================
175      subroutine read_amma(nid,nlevel,ntime                          &
176     &     ,zz,pp,temp,qv,u,v,dw                   &
177     &     ,dt,dq,sens,flat)
178
179!program reading forcings of the AMMA case study
180      USE netcdf, ONLY: nf90_get_var, nf90_open, nf90_noerr, nf90_open, nf90_nowrite, &
181          nf90_inq_dimid, nf90_inquire_dimension, nf90_strerror, nf90_inq_varid
182      implicit none
183
184      integer ntime,nlevel
185
186      real zz(nlevel)
187      real temp(nlevel),pp(nlevel)
188      real qv(nlevel),u(nlevel)
189      real v(nlevel)
190      real dw(nlevel,ntime)
191      real dt(nlevel,ntime)
192      real dq(nlevel,ntime)
193      real flat(ntime),sens(ntime)
194
195
196      integer nid, ierr,rid
197      integer nbvar3d
198      parameter(nbvar3d=30)
199      integer var3didin(nbvar3d)
200
201       ierr=nf90_inq_varid(nid,"zz",var3didin(1))
202         if(ierr/=nf90_noerr) then
203           write(*,*) nf90_strerror(ierr)
204           stop 'lev'
205         endif
206
207
208      ierr=nf90_inq_varid(nid,"temp",var3didin(2))
209         if(ierr/=nf90_noerr) then
210           write(*,*) nf90_strerror(ierr)
211           stop 'temp'
212         endif
213
214      ierr=nf90_inq_varid(nid,"qv",var3didin(3))
215         if(ierr/=nf90_noerr) then
216           write(*,*) nf90_strerror(ierr)
217           stop 'qv'
218         endif
219
220      ierr=nf90_inq_varid(nid,"u",var3didin(4))
221         if(ierr/=nf90_noerr) then
222           write(*,*) nf90_strerror(ierr)
223           stop 'u'
224         endif
225
226      ierr=nf90_inq_varid(nid,"v",var3didin(5))
227         if(ierr/=nf90_noerr) then
228           write(*,*) nf90_strerror(ierr)
229           stop 'v'
230         endif
231
232      ierr=nf90_inq_varid(nid,"dw",var3didin(6))
233         if(ierr/=nf90_noerr) then
234           write(*,*) nf90_strerror(ierr)
235           stop 'dw'
236         endif
237
238      ierr=nf90_inq_varid(nid,"dt",var3didin(7))
239         if(ierr/=nf90_noerr) then
240           write(*,*) nf90_strerror(ierr)
241           stop 'dt'
242         endif
243
244      ierr=nf90_inq_varid(nid,"dq",var3didin(8))
245         if(ierr/=nf90_noerr) then
246           write(*,*) nf90_strerror(ierr)
247           stop 'dq'
248         endif
249     
250      ierr=nf90_inq_varid(nid,"sens",var3didin(9))
251         if(ierr/=nf90_noerr) then
252           write(*,*) nf90_strerror(ierr)
253           stop 'sens'
254         endif
255
256      ierr=nf90_inq_varid(nid,"flat",var3didin(10))
257         if(ierr/=nf90_noerr) then
258           write(*,*) nf90_strerror(ierr)
259           stop 'flat'
260         endif
261
262      ierr=nf90_inq_varid(nid,"pp",var3didin(11))
263         if(ierr/=nf90_noerr) then
264           write(*,*) nf90_strerror(ierr)
265      endif
266
267!dimensions lecture
268!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
269 
270         ierr = nf90_get_var(nid, var3didin(1), zz)
271         if(ierr/=nf90_noerr) then
272            write(*,*) nf90_strerror(ierr)
273            stop "getvarup"
274         endif
275!          write(*,*)'lecture z ok',zz
276
277         ierr = nf90_get_var(nid, var3didin(2), temp)
278         if(ierr/=nf90_noerr) then
279            write(*,*) nf90_strerror(ierr)
280            stop "getvarup"
281         endif
282!          write(*,*)'lecture th ok',temp
283
284         ierr = nf90_get_var(nid, var3didin(3), qv)
285         if(ierr/=nf90_noerr) then
286            write(*,*) nf90_strerror(ierr)
287            stop "getvarup"
288         endif
289!          write(*,*)'lecture qv ok',qv
290 
291         ierr = nf90_get_var(nid, var3didin(4), u)
292         if(ierr/=nf90_noerr) then
293            write(*,*) nf90_strerror(ierr)
294            stop "getvarup"
295         endif
296!          write(*,*)'lecture u ok',u
297
298         ierr = nf90_get_var(nid, var3didin(5), v)
299         if(ierr/=nf90_noerr) then
300            write(*,*) nf90_strerror(ierr)
301            stop "getvarup"
302         endif
303!          write(*,*)'lecture v ok',v
304
305         ierr = nf90_get_var(nid, var3didin(6), dw)
306         if(ierr/=nf90_noerr) then
307            write(*,*) nf90_strerror(ierr)
308            stop "getvarup"
309         endif
310!          write(*,*)'lecture w ok',dw
311
312         ierr = nf90_get_var(nid, var3didin(7), dt)
313         if(ierr/=nf90_noerr) then
314            write(*,*) nf90_strerror(ierr)
315            stop "getvarup"
316         endif
317!          write(*,*)'lecture dt ok',dt
318
319         ierr = nf90_get_var(nid, var3didin(8), dq)
320         if(ierr/=nf90_noerr) then
321            write(*,*) nf90_strerror(ierr)
322            stop "getvarup"
323         endif
324!          write(*,*)'lecture dq ok',dq
325
326         ierr = nf90_get_var(nid, var3didin(9), sens)
327         if(ierr/=nf90_noerr) then
328            write(*,*) nf90_strerror(ierr)
329            stop "getvarup"
330         endif
331!          write(*,*)'lecture sens ok',sens
332
333         ierr = nf90_get_var(nid, var3didin(10), flat)
334         if(ierr/=nf90_noerr) then
335            write(*,*) nf90_strerror(ierr)
336            stop "getvarup"
337         endif
338!          write(*,*)'lecture flat ok',flat
339
340         ierr = nf90_get_var(nid, var3didin(11), pp)
341         if(ierr/=nf90_noerr) then
342            write(*,*) nf90_strerror(ierr)
343            stop "getvarup"
344         endif
345!          write(*,*)'lecture pp ok',pp
346
347         return
348         end subroutine read_amma
349!======================================================================
350        SUBROUTINE interp_amma_time(day,day1,annee_ref                     &
351     &         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma       &
352     &         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma               &
353     &         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
354          USE compar1d_mod_h
355        implicit none
356
357!---------------------------------------------------------------------------------------
358! Time interpolation of a 2D field to the timestep corresponding to day
359!
360! day: current julian day (e.g. 717538.2)
361! day1: first day of the simulation
362! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
363! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
364!---------------------------------------------------------------------------------------
365
366! inputs:
367        integer annee_ref
368        integer nt_amma,nlev_amma
369        integer year_ini_amma
370        real day, day1,day_ini_amma,dt_amma
371        real vitw_amma(nlev_amma,nt_amma)
372        real ht_amma(nlev_amma,nt_amma)
373        real hq_amma(nlev_amma,nt_amma)
374        real lat_amma(nt_amma)
375        real sens_amma(nt_amma)
376! outputs:
377        real vitw_prof(nlev_amma)
378        real ht_prof(nlev_amma)
379        real hq_prof(nlev_amma)
380        real lat_prof,sens_prof
381! local:
382        integer it_amma1, it_amma2,k
383        real timeit,time_amma1,time_amma2,frac
384
385
386        if (forcing_type.eq.6) then
387! Check that initial day of the simulation consistent with AMMA case:
388       if (annee_ref.ne.2006) then
389        print*,'Pour AMMA, annee_ref doit etre 2006'
390        print*,'Changer annee_ref dans run.def'
391        stop
392       endif
393       if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
394        print*,'AMMA a d�but� le 10 juillet 2006',day1,day_ini_amma
395        print*,'Changer dayref dans run.def'
396        stop
397       endif
398       if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
399        print*,'AMMA a fini le 11 juillet'
400        print*,'Changer dayref ou nday dans run.def'
401        stop
402       endif
403       endif
404
405! Determine timestep relative to the 1st day of AMMA:
406!       timeit=(day-day1)*86400.
407!       if (annee_ref.eq.1992) then
408!        timeit=(day-day_ini_toga)*86400.
409!       else
410!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
411!       endif
412      timeit=(day-day_ini_amma)*86400
413
414! Determine the closest observation times:
415!       it_amma1=INT(timeit/dt_amma)+1
416!       it_amma2=it_amma1 + 1
417!       time_amma1=(it_amma1-1)*dt_amma
418!       time_amma2=(it_amma2-1)*dt_amma
419
420       it_amma1=INT(timeit/dt_amma)+1
421       IF (it_amma1 .EQ. nt_amma) THEN
422       it_amma2=it_amma1
423       ELSE
424       it_amma2=it_amma1 + 1
425       ENDIF
426       time_amma1=(it_amma1-1)*dt_amma
427       time_amma2=(it_amma2-1)*dt_amma
428
429       if (it_amma1 .gt. nt_amma) then
430        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
431     &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
432        stop
433       endif
434
435! time interpolation:
436       IF (it_amma1 .EQ. it_amma2) THEN
437          frac=0.
438       ELSE
439          frac=(time_amma2-timeit)/(time_amma2-time_amma1)
440          frac=max(frac,0.0)
441       ENDIF
442
443       lat_prof = lat_amma(it_amma2)                                       &
444     &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))
445       sens_prof = sens_amma(it_amma2)                                     &
446     &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
447
448       do k=1,nlev_amma
449        vitw_prof(k) = vitw_amma(k,it_amma2)                               &
450     &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
451        ht_prof(k) = ht_amma(k,it_amma2)                                   &
452     &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
453        hq_prof(k) = hq_amma(k,it_amma2)                                   &
454     &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
455        enddo
456
457        return
458        END
459
Note: See TracBrowser for help on using the repository browser.