source: LMDZ5/trunk/libf/phylmd/mod_1D_amma_read.F90 @ 2306

Last change on this file since 2306 was 2191, checked in by fhourdin, 10 years ago

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