source: trunk/src/tropflux.pro

Last change on this file was 205, checked in by pinsard, 10 years ago

website not any more hosted by LOCEAN

  • Property svn:executable set to *
  • Property svn:keywords set to Id URL
File size: 21.2 KB
Line 
1;+
2;
3; ============
4; tropflux.pro
5; ============
6;
7; .. function:: tropflux(yyyymmddb,yyyymmdde)
8;
9; DESCRIPTION
10; ===========
11;
12; This program computes net heat flux components on the 1° oaflux grid.
13;
14; all input variables are corrected for mean bias and variability.
15;
16; gustiness correction is applied for wind speed based on cronin's climatological
17; gustiness values.
18;
19; :file:`${PROJECT_ID}/mask_oaflux_30N30S.nc`
20; containing
21; OAFLUX grid
22; has been produced by
23; :func:`oaflux_mask_30n30s`.
24;
25; :file:`${PROJECT_OD}/TropFlux_sst_{yyyymmdd}_{yyyymmdd}.nc`
26; containing
27; sst corrected on OAFLUX grid
28; has been produced by
29; :func:`sst_correction_ncdf`.
30;
31; :file:`${PROJECT_OD}/TropFlux_ws_{yyyymmdd}_{yyyymmdd}.nc`
32; containing
33; ws corrected on OAFLUX grid
34; has been produced by
35; :func:`ws_correction_ncdf`.
36;
37; :file:`${PROJECT_OD}/TropFlux_gustiness_{yyyymmdd}_{yyyymmdd}.nc`
38; containing
39; ++
40; has been produced by
41; :func:`cronin_gustiness_ncdf`.
42;
43; :file:`${PROJECT_OD}/TropFlux_swr_{yyyymmdd}_{yyyymmdd}_BLND.nc`
44; containing
45; ws corrected on OAFLUX grid
46; has been produced by
47; :func:`tropflux_swr_blnd`.
48;
49; :file:`${PROJECT_OD}/TropFlux_lwr_{yyyymmdd}_{yyyymmdd}.nc`
50; containing
51; lwr corrected on OAFLUX grid
52; has been produced by
53; :func:`lwr_correction_ncdf`.
54;
55; containing
56; t2m corrected on OAFLUX grid
57; has been produced by
58; :func:`t2m_correction_ncdf`.
59;
60; :file:`${PROJECT_OD}/TropFlux_q2m_{yyyymmdd}_{yyyymmdd}.nc`
61; containing
62; q2m corrected on OAFLUX grid
63; has been produced by
64; :func:`d2m_to_q2m_erai`.
65;
66; net heat flux components are written
67; in :file:`${PROJECT_OD}/TropFlux_{yyyymmdd}_{yyyymmdd}_coarev3.nc`
68; if this file not already exists.
69;
70; This output file
71; :file:`${PROJECT_OD}/TropFlux_{yyyymmdd}_{yyyymmdd}_coarev3.nc`
72; will be used by
73; :func:`tropflux_nrt_ncdf`.
74;
75;     .. graphviz::
76;
77;        digraph tropflux {
78;
79;           mask [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/mask_oaflux_30N30S.nc"];
80;           file_sst [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_sst_{yyyymmdd}_{yyyymmdd}.nc"];
81;           file_ws [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_ws_{yyyymmdd}_{yyyymmdd}.nc"];
82;           file_wg [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_gustiness_{yyyymmdd}_{yyyymmdd}.nc"];
83;           file_swr [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_swr_{yyyymmdd}_{yyyymmdd}_BLND.nc"];
84;           file_lwr [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_lwr_{yyyymmdd}_{yyyymmdd}.nc"];
85;           file_t2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_t2m_{yyyymmdd}_{yyyymmdd}.nc"];
86;           file_q2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_q2m_{yyyymmdd}_{yyyymmdd}.nc"];
87;
88;           file_out[shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_{yyyymmdd}_{yyyymmdd}_coarev3.nc"];
89;
90;           tropflux [shape=box,
91;           fontname=Courier,
92;           color=blue,
93;           URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/tropflux.pro",
94;           label="${PROJECT}/src/tropflux.pro"];
95;
96;           {mask file_sst file_ws file_wg file_swr file_lwr file_t2m file_q2m} -> {tropflux} -> {file_out}
97;
98;        }
99;
100; SEE ALSO
101; ========
102;
103; :ref:`project_profile.sh`
104;
105; Used by :ref:`tropflux.sh`
106;
107; :func:`report <saxo:report>`
108; :func:`isadirectory <saxo:isadirectory>`
109; :func:`isafile <saxo:isafile>`
110; :func:`initncdf <saxo:initncdf>`
111; :func:`ncdf_lec <saxo:ncdf_lec>`
112; :func:`read_ncdf <saxo:read_ncdf>`
113; :func:`caldat <saxo:caldat>`
114; :func:`julday <saxo:julday>`
115; :func:`ncdf_quickwrite <saxo:ncdf_quickwrite>`
116;
117; :func:`qsee`
118; :func:`lwnet_clark`
119; :func:`calc_claud_vlat`
120; :func:`cor30a`
121;
122; .. note::
123;
124;    :func:`lwdown_clark` is not used here but available.
125;
126; EXAMPLES
127; ========
128;
129; .. code-block:: idl
130;
131;    yyyymmddb = 20000101L
132;    yyyymmdde = 20001231L
133;    result = tropflux(yyyymmddb, yyyymmdde)
134;    print, result
135;
136; TODO
137; ====
138;
139; describe usage of tau
140;
141; compare time between all inputs after reading
142;
143; make it work :
144;
145; .. parsed-literal::
146;
147;     % Compiled module: EXTRAPOLATE.
148;     % Attempt to subscript LAND with OK is out of range.
149;     % Execution halted at: EXTRAPOLATE       176
150;       /usr/home/fplod/SAXO_DIR/SRC/Interpolation/extrapolate.pro
151;     %                      TROPFLUX_19890101_20091231  450
152;       /usr/home/fplod/incas/tropflux/tropflux_ws/src/TropFlux_19890101_20091231.p
153;       ro
154;     %                      $MAIN$
155;     % Program caused arithmetic error: Floating overflow
156;     % Program caused arithmetic error: Floating illegal operand
157;
158; the incriminated line is :
159;
160; .. code-block:: idl
161;
162;    tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y(*,1)) & tab(ocean)=x & m=fi    nite(tab) & tab=extrapolate(tab,m) & lwr(*,*,jt)=tab*msk+valmask*(1-msk)
163;
164; to make progress I DO NOT extrapolate anymore lwr, sen and lat
165; !!!!! must be reactivated
166;
167; when I DO NOT extrapolate lwr, sen and lat , program end with :
168;
169; .. parsed-literal::
170;
171;    Written to
172;    !/usr/work/incas/fplod/tropflux_d/TropFlux_19890101_20091231_coarev3.nc
173;    -------------------------
174;    % Program caused arithmetic error: Floating divide by 0
175;    % Program caused arithmetic error: Floating underflow
176;    % Program caused arithmetic error: Floating overflow
177;    % Program caused arithmetic error: Floating illegal operand
178;
179; explain why :func:`lwdown_clark` is not used
180;
181; check if K or °C in input
182;
183; get rid of:
184;
185; .. parsed-literal::
186;
187;    % date 1: 19880101 is not found in the time axis.
188;
189; why :
190;
191; .. code-block:: idl
192;
193;    da1=19880101
194;
195; 1st date 19890101 no ?
196;
197; why da1-1 (with da1=19880101) when reading gustiness file :
198;
199; .. code-block:: idl
200;
201;    wg=read_ncdf('wg',da1-1,da2,file=fullfilename_wg,/nostr)
202;
203; coding rules
204;
205; KNOWN ISSUES
206; ============
207;
208; test of existence of fullfilename_msk not very efficient because
209; MUST_EXIST keyword of :func:`isafile <saxo:isafile>` not yet implemented
210;
211; EVOLUTIONS
212; ==========
213;
214; $Id$
215;
216; $URL$
217;
218; - fplod 20120322
219;
220;   * pro -> func
221;   * rename with lower case TropFlux_19890101_20091231 become tropflux
222;   * taking project_overwrite into account
223;   * try to add compile_opt
224;   * hard coded da1 and da2 replaced by yyyymmddb and yyyymmdde parameters
225;   * get rid of timegen
226;
227; - fplod 20110830T135832Z cratos (Linux)
228;
229;   * replace tt by time
230;
231; - fplod 20110830T084129Z aedon.locean-ipsl.upmc.fr (Darwin)
232;
233;   * add tau in output file
234;
235; - fplod 20110822T090838Z aedon.locean-ipsl.upmc.fr (Darwin)
236;
237;    * split some multiple lines statement to investigate
238;      pb of extrapolation for lwr
239;
240; - fplod 20110819T144332Z aedon.locean-ipsl.upmc.fr (Darwin)
241;
242;   * add _coarev3 to filename output
243;   * check if filename output exists
244;
245; - fplod 20110809T110911Z aedon.locean-ipsl.upmc.fr (Darwin)
246;
247;   * complete description
248;   * remove v* from filenames (in and out)
249;   * usage of ${PROJECT_OD}
250;   * remove return statement
251;   * add test on IO files
252;
253; - fplod 20101217T140745Z aedon.locean-ipsl.upmc.fr (Darwin)
254;
255;   * remove hard coded directory for mask_oaflux_30N30S.nc
256;
257; - fplod 20101214T112131Z aedon.locean-ipsl.upmc.fr (Darwin)
258;
259;   * add graph
260;
261; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin)
262;
263;   * minimal header
264;
265; - pbk 2008
266;
267;   * creation
268;
269;-
270;
271function tropflux $
272, yyyymmddb $
273, yyyymmdde
274;
275compile_opt idl2, strictarrsubs, logical_predicate
276;
277@cm_4cal
278@cm_4data
279@cm_4mesh
280@cm_4data
281@cm_project
282;
283; Return to caller if errors
284ON_ERROR, 2
285;
286result = -1
287;
288usage = 'result = tropflux_swr_nrt(yyyymmddb, yyyymmdde)'
289nparam = N_PARAMS()
290IF (nparam NE 2) THEN BEGIN
291    ras = report(['Incorrect number of arguments.' $
292    + '!C' $
293    + 'Usage : ' + usage])
294    return, result
295ENDIF
296;
297; check for input directory
298;
299; test if ${PROJECT_ID} defined
300CASE project_id_env OF
301    '' :  BEGIN
302        msg = 'eee : ${PROJECT_ID} is not defined'
303        ras = report(msg)
304        return, result
305    END
306    ELSE : BEGIN
307        msg = 'iii : ${PROJECT_ID} is ' + project_id_env
308        ras = report(msg)
309    END
310ENDCASE
311;
312iodirin = isadirectory(project_id_env)
313;
314; existence and protection of ${PROJECT_ID}
315IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
316    msg = 'eee : the directory' + iodirin  + ' is not accessible.'
317    ras = report(msg)
318    return, result
319ENDIF
320;
321; build mask filename
322filename_msk='mask_oaflux_30N30S.nc'
323;
324; check if this file exists
325msg='iii : looking for ' + filename_msk
326ras = report(msg)
327fullfilename_msk = isafile(iodirin + filename_msk, NEW=0, /MUST_EXIST)
328IF fullfilename_msk[0] EQ '' THEN BEGIN
329    msg = 'eee : the file ' + fullfilename_msk + ' was not found.'
330    ras = report(msg)
331    return, result
332ENDIF
333;
334; test if ${PROJECT_OD} defined
335CASE project_od_env OF
336    '' : BEGIN
337        msg = 'eee : ${PROJECT_OD} is not defined'
338        ras = report(msg)
339        return, result
340    END
341    ELSE : BEGIN
342        msg = 'iii : ${PROJECT_OD} is ' + project_od_env
343        ras = report(msg)
344    END
345ENDCASE
346;
347; check if output data will be possible
348iodirout = isadirectory(project_od_env)
349;
350; existence and protection for reading
351IF (FILE_TEST(iodirout, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
352    msg = 'eee : the directory' + iodirout  + ' is not accessible.'
353    ras = report(msg)
354    return, result
355ENDIF
356;
357; existence and protection for writing
358IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN
359    msg = 'eee : the directory' + iodirout  + ' was not found.'
360    ras = report(msg)
361    return, result
362ENDIF
363;
364; build sst data filename
365filename_sst='TropFlux_sst_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '.nc'
366;
367; check if this file exists
368msg='iii : looking for ' + filename_sst
369ras = report(msg)
370fullfilename_sst = isafile(iodirout + filename_sst, NEW=0, /MUST_EXIST)
371IF fullfilename_sst[0] EQ '' THEN BEGIN
372    msg = 'eee : the file ' + fullfilename_sst + ' was not found.'
373    ras = report(msg)
374    return, result
375ENDIF
376;
377; build ws data filename
378filename_ws='TropFlux_ws_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '.nc'
379;
380; check if this file exists
381msg='iii : looking for ' + filename_ws
382ras = report(msg)
383fullfilename_ws = isafile(iodirout + filename_ws, NEW=0, /MUST_EXIST)
384IF fullfilename_ws[0] EQ '' THEN BEGIN
385    msg = 'eee : the file ' + fullfilename_ws + ' was not found.'
386    ras = report(msg)
387    return, result
388ENDIF
389;
390; build swr data filename
391filename_swr='TropFlux_swr_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '_BLND.nc'
392;
393; check if this file exists
394msg='iii : looking for ' + filename_swr
395ras = report(msg)
396fullfilename_swr = isafile(iodirout + filename_swr, NEW=0, /MUST_EXIST)
397IF fullfilename_swr[0] EQ '' THEN BEGIN
398    msg = 'eee : the file ' + fullfilename_swr + ' was not found.'
399    ras = report(msg)
400    return, result
401ENDIF
402;
403; build lwr data filename
404filename_lwr='TropFlux_lwr_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '.nc'
405;
406; check if this file exists
407msg='iii : looking for ' + filename_lwr
408ras = report(msg)
409fullfilename_lwr = isafile(iodirout + filename_lwr, NEW=0, /MUST_EXIST)
410IF fullfilename_lwr[0] EQ '' THEN BEGIN
411    msg = 'eee : the file ' + fullfilename_lwr + ' was not found.'
412    ras = report(msg)
413    return, result
414ENDIF
415;
416; build t2m data filename
417filename_t2m='TropFlux_t2m_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '.nc'
418;
419; check if this file exists
420msg='iii : looking for ' + filename_t2m
421ras = report(msg)
422fullfilename_t2m = isafile(iodirout + filename_t2m, NEW=0, /MUST_EXIST)
423IF fullfilename_t2m[0] EQ '' THEN BEGIN
424    msg = 'eee : the file ' + fullfilename_t2m + ' was not found.'
425    ras = report(msg)
426    return, result
427ENDIF
428;
429; build q2m data filename
430filename_q2m='TropFlux_q2m_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '.nc'
431;
432; check if this file exists
433msg='iii : looking for ' + filename_q2m
434ras = report(msg)
435fullfilename_q2m = isafile(iodirout + filename_q2m, NEW=0, /MUST_EXIST)
436IF fullfilename_q2m[0] EQ '' THEN BEGIN
437    msg = 'eee : the file ' + fullfilename_q2m + ' was not found.'
438    ras = report(msg)
439    return, result
440ENDIF
441;
442; build wg data filename
443filename_wg='TropFlux_gustiness_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '.nc'
444;
445; check if this file exists
446msg='iii : looking for ' + filename_wg
447ras = report(msg)
448fullfilename_wg = isafile(iodirout + filename_wg, NEW=0, /MUST_EXIST)
449IF fullfilename_wg[0] EQ '' THEN BEGIN
450    msg = 'eee : the file ' + fullfilename_wg + ' was not found.'
451    ras = report(msg)
452    return, result
453ENDIF
454;
455; build output filename
456filename_out = 'TropFlux_' +  string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '_coarev3.nc'
457fullfilename_out = iodirout + filename_out
458; in order to avoid unexpected overwritten
459IF (FILE_TEST(fullfilename_out) EQ 1) THEN BEGIN
460    msg = 'eee : the file ' + fullfilename_out  + ' already exists.'
461    ras = report(msg)
462    return, result
463ENDIF
464;
465initncdf, fullfilename_msk
466msk=ncdf_lec(fullfilename_msk,var='msk')
467;
468initncdf, fullfilename_sst
469;
470ws=read_ncdf('ws',yyyymmddb-.5d,yyyymmdde,file=fullfilename_ws,/nostr)
471timein=24.d*(time-julday(1,1,1957,0,0,0))
472jpt=n_elements(timein)
473da=jul2date(time[0])
474cda0=string(da,format='(i8.8)')
475da=jul2date(time[jpt-1])
476cda1=string(da,format='(i8.8)')
477print, 'ws in tropflux first date ', cda0
478print, 'ws in ws_correction_ncdf last date ' , cda1
479;
480wg=read_ncdf('wg',yyyymmddb-.5d-1,yyyymmdde,file=fullfilename_wg,/nostr)
481timein=24.d*(time-julday(1,1,1957,0,0,0))
482jpt=n_elements(timein)
483da=jul2date(time[0])
484cda0=string(da,format='(i8.8)')
485da=jul2date(time[jpt-1])
486cda1=string(da,format='(i8.8)')
487print, 'wg in tropflux first date ', cda0
488print, 'wg in ws_correction_ncdf last date ' , cda1
489tt=time
490jpt=n_elements(time)
491sst=read_ncdf('sst',yyyymmddb-.5d,yyyymmdde,file=fullfilename_sst,/nostr)
492timein=24.d*(time-julday(1,1,1957,0,0,0))
493jpt=n_elements(timein)
494da=jul2date(time[0])
495cda0=string(da,format='(i8.8)')
496da=jul2date(time[jpt-1])
497cda1=string(da,format='(i8.8)')
498print, 'sst in tropflux first date ', cda0
499print, 'sst in ws_correction_ncdf last date ' , cda1
500sst=reform(sst-273.15)
501swd=read_ncdf('swr',yyyymmddb-.5d,yyyymmdde,file=fullfilename_swr,/nostr)
502timein=24.d*(time-julday(1,1,1957,0,0,0))
503jpt=n_elements(timein)
504da=jul2date(time[0])
505cda0=string(da,format='(i8.8)')
506da=jul2date(time[jpt-1])
507cda1=string(da,format='(i8.8)')
508print, 'swr in tropflux first date ', cda0
509print, 'swr in ws_correction_ncdf last date ' , cda1
510lw=read_ncdf('lwr',yyyymmddb-.5d,yyyymmdde,file=fullfilename_lwr,/nostr)
511timein=24.d*(time-julday(1,1,1957,0,0,0))
512jpt=n_elements(timein)
513da=jul2date(time[0])
514cda0=string(da,format='(i8.8)')
515da=jul2date(time[jpt-1])
516cda1=string(da,format='(i8.8)')
517print, 'lwr in tropflux first date ', cda0
518print, 'lwr in ws_correction_ncdf last date ' , cda1
519swd=swd/0.94  ; converting from net swr to downward swr
520;
521t2m=read_ncdf('t2m',yyyymmddb-.5d,yyyymmdde,file=fullfilename_t2m,/nostr)-273.15   ; in C
522timein=24.d*(time-julday(1,1,1957,0,0,0))
523jpt=n_elements(timein)
524da=jul2date(time[0])
525cda0=string(da,format='(i8.8)')
526da=jul2date(time[jpt-1])
527cda1=string(da,format='(i8.8)')
528print, 't2m in tropflux first date ', cda0
529print, 't2m in ws_correction_ncdf last date ' , cda1
530q2m=read_ncdf('q2m',yyyymmddb-.5d,yyyymmdde,file=fullfilename_q2m,/nostr)    ; in g/kg
531timein=24.d*(time-julday(1,1,1957,0,0,0))
532jpt=n_elements(timein)
533da=jul2date(time[0])
534cda0=string(da,format='(i8.8)')
535da=jul2date(time[jpt-1])
536cda1=string(da,format='(i8.8)')
537print, 'q2m in tropflux first date ', cda0
538print, 'q2m in ws_correction_ncdf last date ' , cda1
539;
540w=sqrt(ws*ws+wg*wg)  ; wind corrected for gustiness
541w=ws
542tmask=msk
543help, ws,w,wg,u,sst,swd,t2m,q2m
544ocean=where(msk eq 1,compl=land)
545valmask=1.e20
546time=tt
547jpt=n_elements(time)
548;
549; Constants for flux computation
550;
551zu=10.                          ; height of wind speed measurement (m)
552us=0.                           ; surf current (m/s)
553zt=2.                           ; Height of air T measurement (m)
554zq=2.                           ; height of humidity measurement (m)
555P=1008.                         ; Pressure
556zi=600.                         ; Inversion height (m)
557jcool=0                         ; Compute cool-skin
558jwave=0                         ; No waves
559twave=5.
560hwave=1.
561;
562caldat, time,mon,day,yea
563swr=fltarr(jpi,jpj,jpt)+1.e20
564lwr=fltarr(jpi,jpj,jpt)+1.e20
565lat=fltarr(jpi,jpj,jpt)+1.e20
566sen=fltarr(jpi,jpj,jpt)+1.e20
567lwnet_clrk=fltarr(jpi,jpj,jpt)+1.e20
568;Ch=fltarr(jpi,jpj,jpt)+1.e20
569;Ce=fltarr(jpi,jpj,jpt)+1.e20
570junk=fltarr(jpi,jpj,jpt)+1.e20
571;
572for jt=0,jpt-1 do begin
573    jday=time[jt]-julday(1,1,yea[jt])
574    print, 'Computing Fluxes ',jt,' / ',jpt-1
575    ;
576    ;  P=msl(*,*,jt)
577    ;  P=P(ocean)
578    ;
579    ; wind speed (m/s)
580    wn=w[*,*,jt]
581    wn=wn[ocean]
582    ; Bulk sst (°C)
583    ts=sst[*,*,jt]
584    ts=ts[ocean]
585    ; 2m Air T (°C)
586    t=t2m[*,*,jt]
587    t=t[ocean]
588    ; Sea surface sat. spec. humidity (g/kg)
589    qs=qsee(ts,P)
590    ; 2m AIr specific humidity  (g/kg)
591    q=q2m[*,*,jt]
592    q=q[ocean]
593    ; Downward solar flux (W/m2)
594    Rs=swd[*,*,jt]
595    Rs=Rs[ocean]
596    ylat=gphit[ocean]
597    ;  cld=calc_cloud(jday,Rs,ylat)
598    cld=calc_cloud_vlat(jday,Rs,ylat)
599    ; Downward IR flux (W/m2)
600    ;  Rl=lwdown_clark(ts,q,cld,t,P)
601    Rl=lw[*,*,jt]
602    Rl=Rl[ocean]
603    rain=0.
604    lw_clrk=-lwnet_clark(ts,q,cld,t,P)
605    ;  junk[*,*,jt]=lw_clrk
606    ;
607;stop
608    y=cor30a(wn,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,ylat,jcool,jwave,twave,hwave)
609    ; A few punctual missing values (coare does not converge): filled by spatial extrapolation
610    tab=fltarr(jpi,jpj)+!values.f_nan
611    x=reform(y[*,0])
612    tab[ocean]=x
613    m=finite(tab)
614    help, tab, m, swr
615    tab=extrapolate(tab,m)
616    swr[*,*,jt]=tab*msk+valmask*(1-msk)
617    ;
618    tab=fltarr(jpi,jpj)+!values.f_nan
619    x=reform(y[*,1])
620    tab[ocean]=x
621    m=finite(tab)
622    help, tab, m, lwr
623    ;!!!!!!!!!!!!+++++++++++++++  tab=extrapolate(tab,m)
624    print, 'www : extrapolation for lwr is not done'
625    tab=tab
626    lwr[*,*,jt]=tab*msk+valmask*(1-msk)
627    ;
628    tab=fltarr(jpi,jpj)+!values.f_nan
629    x=reform(y[*,2])
630    tab[ocean]=x
631    m=finite(tab)
632    help, tab, m, lat
633    ;!!!!!!!!!!!!+++++++++++++++  tab=extrapolate(tab,m)
634    print, 'www : extrapolation for lat is not done'
635    tab=tab
636    lat[*,*,jt]=tab*msk+valmask*(1-msk)
637    tab=fltarr(jpi,jpj)+!values.f_nan
638    x=reform(y[*,3])
639    tab[ocean]=x
640    m=finite(tab)
641    help, tab, m, sen
642    ;!!!!!!!!!!!!+++++++++++++++  tab=extrapolate(tab,m)
643    print, 'www : extrapolation for sen is not done'
644    tab=tab
645    sen[*,*,jt]=tab*msk+valmask*(1-msk)
646    ;  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(lw_clrk) & tab[ocean]=x & m=finite(tab) & tab=extrapolate(tab,m) & lwnet_clrk[*,*,jt]=tab*msk+valmask*(1-msk)
647    tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y[*,5]) & tab[ocean]=x & m=finite(tab) & tab=extrapolate(tab,m) & tau[*,*,jt]=tab*msk+valmask*(1-msk)
648    ;  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y[*,6]) & tab[ocean]=x & m=finite(tab) & tab=extrapolate(tab,m) & Ch[*,*,jt]=tab*msk+valmask*(1-msk)
649    ;  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y[*,7]) & tab[ocean]=x & m=finite(tab) & tab=extrapolate(tab,m) & Ce[*,*,jt]=tab*msk+valmask*(1-msk)
650    tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y[*,8]) & tab[ocean]=x & m=finite(tab) & tab=extrapolate(tab,m) & wg[*,*,jt]=tab*msk+valmask*(1-msk)
651endfor
652jpt=n_elements(time)
653;
654cda0=string(jul2date(time[0]),format='(i8.8)')
655cda1=string(jul2date(time[jpt-1]),format='(i8.8)')
656time=time-julday(1,1,1950,00,00,00)
657xlon=reform(glamt[*,0])
658ylat=reform(gphit[0,*])
659;
660ncfile='!' + fullfilename_out
661lon_attr={units:'degrees_east',long_name:'Longitude'}
662lat_attr={units:'degrees_north',long_name:'Latitude'}
663time_attr={units:'days since 1950-01-01 00:00:00',long_name:'Time axis',time_origin:'1950-JAN-01 00:00:00'}
664globattr={source:'Fluxes for the Global Tropical Ocean - TropFlux',timerange:cda0+' - '+cda1}
665swr_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net shortwave radiation',short_name:'swr',axis:'TYX'}
666lwr_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net longwave radiation',short_name:'lwr',axis:'TYX'}
667lwr_clrk_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net longwave radiation from clark',short_name:'lwr',axis:'TYX'}
668lhf_attr={units:'W/m2',missing_value:valmask,long_name:'Surface latent flux',short_name:'lhf',axis:'TYX'}
669shf_attr={units:'W/m2',missing_value:valmask,long_name:'Surface sensible flux',short_name:'shf',axis:'TYX'}
670wg_attr={units:'m/s',missing_value:valmask,long_name:'COARE convective gustiness',short_name:'wg',axis:'TYX'}
671tau_attr={units:'N/m2',missing_value:valmask,long_name:'Wind stress magnitude',short_name:'tau',axis:'TYX'}
672;
673help, swr,lwr,lat,sen,time,xlon,ylat
674;
675ncfields = 'tau[longitude,latitude,time]=tau:tau_attr; ' $
676+ 'wg[longitude,latitude,time]=wg:wg_attr; ' $
677+ 'swr[longitude,latitude,time]=swr:swr_attr; ' $
678;         +'lwr_coare[longitude,latitude,time]=lwr:lwr_attr; ' $
679+'lwr[longitude,latitude,time]=lwnet_clrk:lwr_clrk_attr; ' $
680+'lhf[longitude,latitude,time]=lat:lhf_attr; ' $
681+'shf[longitude,latitude,time]=sen:shf_attr; ' $
682;         +'Ch[longitude,latitude,time]=Ch:Ch_attr; ' $
683;         +'Ce[longitude,latitude,time]=Ce:Ce_attr; ' $
684+ 'longitude[]=xlon:lon_attr; ' $
685+ 'latitude[]=ylat:lat_attr; ' $
686+ 'time[*time]=time:time_attr ' $
687+ ' @ globattr'
688;
689@ncdf_quickwrite
690;
691result = 0
692return, result
693;
694end
Note: See TracBrowser for help on using the repository browser.