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

Last change on this file since 5270 was 5270, checked in by abarral, 2 weeks ago

Replace F77 netcdf library by F90 netcdf library

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        implicit none
355
356!---------------------------------------------------------------------------------------
357! Time interpolation of a 2D field to the timestep corresponding to day
358!
359! day: current julian day (e.g. 717538.2)
360! day1: first day of the simulation
361! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
362! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
363!---------------------------------------------------------------------------------------
364
365        INCLUDE "compar1d.h"
366
367! inputs:
368        integer annee_ref
369        integer nt_amma,nlev_amma
370        integer year_ini_amma
371        real day, day1,day_ini_amma,dt_amma
372        real vitw_amma(nlev_amma,nt_amma)
373        real ht_amma(nlev_amma,nt_amma)
374        real hq_amma(nlev_amma,nt_amma)
375        real lat_amma(nt_amma)
376        real sens_amma(nt_amma)
377! outputs:
378        real vitw_prof(nlev_amma)
379        real ht_prof(nlev_amma)
380        real hq_prof(nlev_amma)
381        real lat_prof,sens_prof
382! local:
383        integer it_amma1, it_amma2,k
384        real timeit,time_amma1,time_amma2,frac
385
386
387        if (forcing_type.eq.6) then
388! Check that initial day of the simulation consistent with AMMA case:
389       if (annee_ref.ne.2006) then
390        print*,'Pour AMMA, annee_ref doit etre 2006'
391        print*,'Changer annee_ref dans run.def'
392        stop
393       endif
394       if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
395        print*,'AMMA a d�but� le 10 juillet 2006',day1,day_ini_amma
396        print*,'Changer dayref dans run.def'
397        stop
398       endif
399       if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
400        print*,'AMMA a fini le 11 juillet'
401        print*,'Changer dayref ou nday dans run.def'
402        stop
403       endif
404       endif
405
406! Determine timestep relative to the 1st day of AMMA:
407!       timeit=(day-day1)*86400.
408!       if (annee_ref.eq.1992) then
409!        timeit=(day-day_ini_toga)*86400.
410!       else
411!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
412!       endif
413      timeit=(day-day_ini_amma)*86400
414
415! Determine the closest observation times:
416!       it_amma1=INT(timeit/dt_amma)+1
417!       it_amma2=it_amma1 + 1
418!       time_amma1=(it_amma1-1)*dt_amma
419!       time_amma2=(it_amma2-1)*dt_amma
420
421       it_amma1=INT(timeit/dt_amma)+1
422       IF (it_amma1 .EQ. nt_amma) THEN
423       it_amma2=it_amma1
424       ELSE
425       it_amma2=it_amma1 + 1
426       ENDIF
427       time_amma1=(it_amma1-1)*dt_amma
428       time_amma2=(it_amma2-1)*dt_amma
429
430       if (it_amma1 .gt. nt_amma) then
431        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
432     &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
433        stop
434       endif
435
436! time interpolation:
437       IF (it_amma1 .EQ. it_amma2) THEN
438          frac=0.
439       ELSE
440          frac=(time_amma2-timeit)/(time_amma2-time_amma1)
441          frac=max(frac,0.0)
442       ENDIF
443
444       lat_prof = lat_amma(it_amma2)                                       &
445     &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))
446       sens_prof = sens_amma(it_amma2)                                     &
447     &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
448
449       do k=1,nlev_amma
450        vitw_prof(k) = vitw_amma(k,it_amma2)                               &
451     &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
452        ht_prof(k) = ht_amma(k,it_amma2)                                   &
453     &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
454        hq_prof(k) = hq_amma(k,it_amma2)                                   &
455     &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
456        enddo
457
458        return
459        END
460
Note: See TracBrowser for help on using the repository browser.