source: trunk/tools/density_binning/binning_neutral_and_co/npc.pro @ 2

Last change on this file since 2 was 2, checked in by post_it, 17 years ago

Initial import from ~/POST_IT/
File size: 3.8 KB
Line 
1;  Name: pnc
2;  -----
3;                     
4;  Purpose :
5;  --------
6;     Non penetrative convective adjustment scheme. solve the static
7;     instability of the water column.
8;       
9;  Method :
10;  -------
11;       The algorithm used converges in a maximium of jpk iterations.
12;      instabilities are treated when the vertical density gradient
13;      is less than 1.e-5.
14;
15;   Input : potential density (what ever the reference is) at t-point
16;   ------
17;
18;   Output :  adjusted potential density field
19;   -------
20;             
21;   References :
22;   -----------
23;       Madec et al., 1991, JPO, 21, 9, 1349-1371.
24;
25;   Modifications :
26;   --------------
27;      original : 90-09 (G. Madec)
28;      Additions : 01-06 (G. Madec) Idl version
29;----------------------------------------------------------------------
30FUNCTION npc, s3d
31@common
32;;
33;; Definition des variables
34;;
35;   Profils selon les niveaux du modele (suffixe _z)
36;
37   s_z = fltarr(jpk)     ; profil de la densite
38;
39; Tableau de sortie
40;
41   rhos = fltarr(jpi, jpj, jpk)
42   rhos = s3d
43;
44;  ====================================
45;  Loop over the horizontal domain (2D)
46;  ====================================
47;
48   ncompt = 0
49   FOR i = 0, jpi-1 DO BEGIN
50      FOR j = 0, jpj-1 DO BEGIN
51      ;
52      ;  Indices des points T dans l ocean
53         i_ocean = where(tmask(i,j,*) EQ 1)
54      ;
55         IF (i_ocean[0] NE -1) THEN BEGIN         ; on n'entre que si il y a des points ocean
56            ;
57            ; density profil
58                       s_z(*)= s3d(i,j,*)
59                       ;s_z(*) = rho(i,j,*)
60            ; 
61            ; 1. Static instability pointer
62            ; -----------------------------
63            ;
64              ds =(shift(s_z,-1)-s_z)(i_ocean(0:n_elements(i_ocean)-2))
65              ind_c = where(ds LT 0.)
66            ;
67            ; 2. Vertical mixing for each instable portion of the density profil
68            ;
69              IF ( ind_c(0) NE -1 ) THEN BEGIN
70                  ncompt=ncompt+1
71            ;      print, 'static instability at i,j=', i,j
72            ;
73            ; -->> the density profil is statically instable :
74            ; ikbot: last ocean level (just above the bottom)
75                  ikbot = n_elements(i_ocean)-1
76            ; vertical iteration
77                  jiter = 0
78                  WHILE ( (ind_c(0) NE -1) AND (jiter LT jpk-1) ) DO BEGIN              &$
79                    jiter = jiter+1                                                     &$
80                  ; ikup : the first static instability from the sea surface
81                    ikup = ind_c(0)                                                     &$                                                     
82                  ; the density profil is instable below ikup
83                  ; ikdown : bottom of the instable portion of the density profil
84                  ; search of ikdown and vertical mixing from ikup to ikdown
85                    ze3tot= e3t(ikup)                                                   &$
86                    zraua = s_z(ikup)                                                   &$
87                    jkdown = ikup+1                                                     &$
88;
89                    WHILE (jkdown LE ikbot AND zraua GT s_z(jkdown) ) DO BEGIN  &$
90                      ze3dwn =  e3t(jkdown)                                             &$
91                      ze3tot =  ze3tot+ze3dwn                                           &$
92                      zraua = ( zraua*(ze3tot-ze3dwn) + s_z(jkdown)*ze3dwn )/ze3tot     &$
93                      jkdown=jkdown+1                                                   &$
94;                      print, jkdown, zraua                                             &$
95                    ENDWHILE                                                            &$
96;
97                    FOR jkp = ikup,jkdown-1 DO BEGIN                                    &$
98                      s_z(jkp) = zraua                                                  &$
99                    ENDFOR                                                              &$
100                    ds =(shift(s_z, -1)-s_z)(i_ocean(0:n_elements(i_ocean)-2))          &$
101                    ind_c = where(ds LT 0.)                                             &$
102                  ENDWHILE
103              ENDIF
104            ; save the modifications
105              rhos(i,j,*) = s_z(*)
106;
107; <<-- no more static instability on slab jj
108;
109         ENDIF
110       ENDFOR
111     ENDFOR
112;
113     print, ' number of static instability treated : ', ncompt
114; sortie:
115   return,  rhos
116
117END
118
Note: See TracBrowser for help on using the repository browser.