source: trunk/procs/npc.pro @ 12

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

Initial import from ~/POST_IT/

File size: 3.8 KB
Line 
1;  Name: npc
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
50   FOR i = 0, jpi-1 DO BEGIN
51      FOR j = 0, jpj-1 DO BEGIN
52
53        ;  Indices des points T dans l ocean
54         i_ocean = where(tmask(i,j,*) EQ 1)
55
56         IF (i_ocean[0] NE -1) THEN BEGIN         ; on n'entre que si il y a des points ocean
57            ;
58            ; density profil
59                       s_z(*)= s3d(i,j,*)
60                       ;s_z(*) = rho(i,j,*)
61            ; 
62            ; 1. Static instability pointer
63            ; -----------------------------
64            ;
65              ds =(shift(s_z,-1)-s_z)(i_ocean(0:n_elements(i_ocean)-2))
66              ind_c = where(ds LT 0.)
67            ;
68            ; 2. Vertical mixing for each instable portion of the density profil
69            ;
70              IF ( ind_c(0) NE -1 ) THEN BEGIN
71                  ncompt=ncompt+1
72            ;      print, 'static instability at i,j=', i,j
73            ;
74            ; -->> the density profil is statically instable :
75            ; ikbot: last ocean level (just above the bottom)
76                  ikbot = n_elements(i_ocean)-1
77            ; vertical iteration
78                  jiter = 0
79                  WHILE ( (ind_c(0) NE -1) AND (jiter LT jpk-1) ) DO BEGIN             
80                    jiter = jiter+1                                                     
81                  ; ikup : the first static instability from the sea surface           
82                    ikup = ind_c(0)                                                                                                             
83                  ; the density profil is instable below ikup                           
84                  ; ikdown : bottom of the instable portion of the density profil       
85                  ; search of ikdown and vertical mixing from ikup to ikdown           
86                    ze3tot= e3t(ikup)                                                   
87                    zraua = s_z(ikup)                                                   
88                    jkdown = ikup+1
89;                                                                                       
90                    WHILE (jkdown LE ikbot AND zraua GT s_z(jkdown) ) DO BEGIN     
91                      ze3dwn =  e3t(jkdown)
92                      ze3tot =  ze3tot+ze3dwn
93                      zraua = ( zraua*(ze3tot-ze3dwn) + s_z(jkdown)*ze3dwn )/ze3tot
94                      jkdown = jkdown + 1
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.