New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
make_closea_masks.py in NEMO/branches/2019/ENHANCE-02_ISF_domcfg – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/make_closea_masks.py @ 11574

Last change on this file since 11574 was 9161, checked in by davestorkey, 6 years ago

Reformulation of closea module.
See ticket #2000
https://forge.ipsl.jussieu.fr/nemo/wiki/2017WP/ROBUST-14_Dave_Storkey-Closed_Seas_rewrite

  • Property svn:executable set to *
File size: 18.3 KB
Line 
1#!/usr/local/sci/bin/python2.7
2
3'''
4Routine to create closea mask fields based on old NEMO closea index definitions.
5Details of the grid and the bathymetry are read in from the domain_cfg.nc file and
6the closea_mask* fields are appended to the same domain_cfg.nc file.
7
8To use this routine:
9
10  1. Provide domain_cfg.nc file for your configuration.
11
12  2. Define closed seas for your configuration in Section 2
13     using indices in the old NEMO style. (Read the comments on
14     indexing in Section 2!). Examples are given for eORCA025
15     (UK version) for the three different options:
16        - just defining closed seas (and distribute fluxes over global ocean)
17        - defining closed seas with a RNF mapping for the American Great Lakes to the St Laurence Seaway
18        - defining closed seas with an EMPMR mapping for the American Great Lakes to the St Laurence Seaway
19
20  3. Choose whether to mask the closea_mask* fields. Not required
21     but makes the fields easier to check.
22
23  4. Module can be run in python or from linux command line if you
24     change the top line to point to your python installation. If
25     using from command line, type "make_closea_masks.py --help"
26     for usage.
27
28@author: Dave Storkey
29@date: Dec 2017
30'''
31import netCDF4 as nc
32import numpy as np
33import numpy.ma as ma
34
35def make_closea_masks(config=None,domcfg_file=None,mask=None):
36
37#=========================
38# 1. Read in domcfg file
39#=========================
40
41    if config is None:
42        raise Exception('configuration must be specified')
43
44    if domcfg_file is None:
45        raise Exception('domain_cfg file must be specified')
46
47    if mask is None:
48        mask=False
49
50    domcfg = nc.Dataset(domcfg_file,'r+')
51    lon = domcfg.variables['nav_lon'][:]
52    lat = domcfg.variables['nav_lat'][:]
53    top_level = domcfg.variables['top_level'][0][:]
54
55    nx = top_level.shape[1]
56    ny = top_level.shape[0]
57
58    # Generate 2D "i" and "j" fields for use in "where" statements.
59    # These are the Fortran indices, counting from 1, so we have to
60    # add 1 to np.arange because python counts from 0.
61
62    ones_2d = np.ones((ny,nx))
63    ii1d = np.arange(nx)+1
64    jj1d = np.arange(ny)+1
65    ii2d = ii1d * ones_2d
66    jj2d = np.transpose(jj1d*np.transpose(ones_2d)) 
67 
68#=====================================
69# 2. Closea definitions (old style)
70#=====================================
71
72    # NB. The model i and j indices defined here are Fortran indices,
73    #     ie. counting from 1 as in the NEMO code. Also the indices
74    #     of the arrays (ncsi1 etc) count from 1 in order to match
75    #     the Fortran code.
76    #     This means that you can cut and paste the definitions from
77    #     the NEMO code and change round brackets to square brackets.
78    #     But BEWARE: Fortran array(a:b) == Python array[a:b+1] !!!
79    #
80
81    # If use_runoff_box = True then specify runoff area as all sea points within
82    # a rectangular area. If use_runoff_box = False then specify a list of points
83    # as in the old NEMO code. Default to false.
84    use_runoff_box = False
85
86    #================================================================
87    if config == 'ORCA2':
88
89        num_closea = 4
90        max_runoff_points = 4
91
92        ncsnr = np.zeros(num_closea+1,dtype=np.int)                     ; ncstt = np.zeros(num_closea+1,dtype=np.int)
93        ncsi1 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj1 = np.zeros(num_closea+1,dtype=np.int)
94        ncsi2 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj2 = np.zeros(num_closea+1,dtype=np.int)
95        ncsir = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int) ; ncsjr = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int)
96
97        # Caspian Sea (spread over globe)
98        ncsnr[1]   =   1  ;  ncstt[1]   =   0   
99        ncsi1[1]   =  11  ;  ncsj1[1]   = 103
100        ncsi2[1]   =  17  ;  ncsj2[1]   = 112
101
102        # Great Lakes - North America - put at St Laurent mouth
103        ncsnr[2]   =   1  ;  ncstt[2]   =   2 
104        ncsi1[2]   =  97  ;  ncsj1[2]   = 107
105        ncsi2[2]   = 103  ;  ncsj2[2]   = 111
106        ncsir[2,1] = 110  ;  ncsjr[2,1] = 111           
107
108        # Black Sea (crossed by the cyclic boundary condition)
109        # put in Med Sea (north of Aegean Sea)
110        ncsnr[3:5] =   4  ;  ncstt[3:5] =   2           
111        ncsir[3:5,1] = 171;  ncsjr[3:5,1] = 106     
112        ncsir[3:5,2] = 170;  ncsjr[3:5,2] = 106 
113        ncsir[3:5,3] = 171;  ncsjr[3:5,3] = 105 
114        ncsir[3:5,4] = 170;  ncsjr[3:5,4] = 105 
115        # west part of the Black Sea     
116        ncsi1[3]   = 174  ;  ncsj1[3]   = 107     
117        ncsi2[3]   = 181  ;  ncsj2[3]   = 112     
118        # east part of the Black Sea     
119        ncsi1[4]   =   2  ;  ncsj1[4]   = 107     
120        ncsi2[4]   =   6  ;  ncsj2[4]   = 112     
121
122    #================================================================
123    elif config == 'eORCA1':
124
125        num_closea = 1
126        max_runoff_points = 1
127
128        ncsnr = np.zeros(num_closea+1,dtype=np.int)                     ; ncstt = np.zeros(num_closea+1,dtype=np.int)
129        ncsi1 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj1 = np.zeros(num_closea+1,dtype=np.int)
130        ncsi2 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj2 = np.zeros(num_closea+1,dtype=np.int)
131        ncsir = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int) ; ncsjr = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int)
132
133        # Caspian Sea  (spread over the globe)
134        ncsnr[1]   = 1    ; ncstt[1]   = 0           
135        ncsi1[1]   = 332  ; ncsj1[1]   = 243
136        ncsi2[1]   = 344  ; ncsj2[1]   = 275
137
138    #================================================================
139    elif config == 'eORCA025_UK':
140
141        num_closea = 10
142        max_runoff_points = 1
143
144        ncsnr = np.zeros(num_closea+1,dtype=np.int)                     ; ncstt = np.zeros(num_closea+1,dtype=np.int)
145        ncsi1 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj1 = np.zeros(num_closea+1,dtype=np.int)
146        ncsi2 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj2 = np.zeros(num_closea+1,dtype=np.int)
147        ncsir = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int) ; ncsjr = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int)
148
149        # Caspian Sea
150        ncsnr[1]   = 1    ; ncstt[1]   = 0     
151        ncsi1[1]   = 1330 ; ncsj1[1]   = 831
152        ncsi2[1]   = 1375 ; ncsj2[1]   = 981
153
154        # Aral Sea
155        ncsnr[2]   = 1    ; ncstt[2]   = 0     
156        ncsi1[2]   = 1376 ; ncsj1[2]   = 900
157        ncsi2[2]   = 1400 ; ncsj2[2]   = 981
158
159        # Azov Sea
160        ncsnr[3]   = 1    ; ncstt[3]   = 0     
161        ncsi1[3]   = 1284 ; ncsj1[3]   = 908
162        ncsi2[3]   = 1304 ; ncsj2[3]   = 933
163
164        # Lake Superior
165        ncsnr[4]   = 1    ; ncstt[4]   = 0     
166        ncsi1[4]   = 781  ; ncsj1[4]   = 905 
167        ncsi2[4]   = 815  ; ncsj2[4]   = 926 
168
169        # Lake Michigan
170        ncsnr[5]   = 1    ; ncstt[5]   = 0     
171        ncsi1[5]   = 795  ; ncsj1[5]   = 871             
172        ncsi2[5]   = 813  ; ncsj2[5]   = 905 
173
174        # Lake Huron part 1
175        ncsnr[6]   = 1    ; ncstt[6]   = 0     
176        ncsi1[6]   = 814  ; ncsj1[6]   = 882             
177        ncsi2[6]   = 825  ; ncsj2[6]   = 905 
178
179        # Lake Huron part 2
180        ncsnr[7]   = 1    ; ncstt[7]   = 0     
181        ncsi1[7]   = 826  ; ncsj1[7]   = 889             
182        ncsi2[7]   = 833  ; ncsj2[7]   = 905 
183
184        # Lake Erie
185        ncsnr[8]   = 1    ; ncstt[8]   = 0     
186        ncsi1[8]   = 816  ; ncsj1[8]   = 871             
187        ncsi2[8]   = 837  ; ncsj2[8]   = 881 
188
189        # Lake Ontario
190        ncsnr[9]   = 1    ; ncstt[9]   = 0     
191        ncsi1[9]   = 831  ; ncsj1[9]   = 882             
192        ncsi2[9]   = 847  ; ncsj2[9]   = 889 
193
194        # Lake Victoria
195        ncsnr[10]   = 1    ; ncstt[10]   = 0   
196        ncsi1[10]   = 1274 ; ncsj1[10]   = 672 
197        ncsi2[10]   = 1289 ; ncsj2[10]   = 687 
198
199    #================================================================
200    elif config == 'eORCA025_UK_rnf':
201
202        num_closea = 10
203        max_runoff_points = 1
204        use_runoff_box = True
205
206        ncsnr = np.zeros(num_closea+1,dtype=np.int)                     ; ncstt = np.zeros(num_closea+1,dtype=np.int)
207        ncsi1 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj1 = np.zeros(num_closea+1,dtype=np.int)
208        ncsi2 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj2 = np.zeros(num_closea+1,dtype=np.int)
209        ncsir1 = np.zeros(num_closea+1,dtype=np.int)                    ; ncsjr1 = np.zeros(num_closea+1,dtype=np.int)
210        ncsir2 = np.zeros(num_closea+1,dtype=np.int)                    ; ncsjr2 = np.zeros(num_closea+1,dtype=np.int)
211        ncsir = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int) ; ncsjr = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int)
212
213        # Caspian Sea
214        ncsnr[1]   = 1    ; ncstt[1]   = 0     
215        ncsi1[1]   = 1330 ; ncsj1[1]   = 831
216        ncsi2[1]   = 1375 ; ncsj2[1]   = 981
217
218        # Aral Sea
219        ncsnr[2]   = 1    ; ncstt[2]   = 0     
220        ncsi1[2]   = 1376 ; ncsj1[2]   = 900
221        ncsi2[2]   = 1400 ; ncsj2[2]   = 981
222
223        # Azov Sea
224        ncsnr[3]   = 1    ; ncstt[3]   = 0     
225        ncsi1[3]   = 1284 ; ncsj1[3]   = 908
226        ncsi2[3]   = 1304 ; ncsj2[3]   = 933
227
228        # Lake Superior
229        ncsnr[4]   = 1    ; ncstt[4]   = 1     
230        ncsi1[4]   = 781  ; ncsj1[4]   = 905 
231        ncsi2[4]   = 815  ; ncsj2[4]   = 926 
232        # runff points the St Laurence Seaway for all Great Lakes
233        ncsir1[4:10]   = 873 ; ncsjr1[4:10]   = 909 
234        ncsir2[4:10]   = 884 ; ncsjr2[4:10]   = 920 
235
236        # Lake Michigan
237        ncsnr[5]   = 1    ; ncstt[5]   = 1     
238        ncsi1[5]   = 795  ; ncsj1[5]   = 871             
239        ncsi2[5]   = 813  ; ncsj2[5]   = 905 
240
241        # Lake Huron part 1
242        ncsnr[6]   = 1    ; ncstt[6]   = 1     
243        ncsi1[6]   = 814  ; ncsj1[6]   = 882             
244        ncsi2[6]   = 825  ; ncsj2[6]   = 905 
245
246        # Lake Huron part 2
247        ncsnr[7]   = 1    ; ncstt[7]   = 1     
248        ncsi1[7]   = 826  ; ncsj1[7]   = 889             
249        ncsi2[7]   = 833  ; ncsj2[7]   = 905 
250
251        # Lake Erie
252        ncsnr[8]   = 1    ; ncstt[8]   = 1     
253        ncsi1[8]   = 816  ; ncsj1[8]   = 871             
254        ncsi2[8]   = 837  ; ncsj2[8]   = 881 
255
256        # Lake Ontario
257        ncsnr[9]   = 1    ; ncstt[9]   = 1     
258        ncsi1[9]   = 831  ; ncsj1[9]   = 882             
259        ncsi2[9]   = 847  ; ncsj2[9]   = 889 
260
261        # Lake Victoria
262        ncsnr[10]   = 1    ; ncstt[10]   = 0   
263        ncsi1[10]   = 1274 ; ncsj1[10]   = 672 
264        ncsi2[10]   = 1289 ; ncsj2[10]   = 687 
265
266    #================================================================
267    elif config == 'eORCA025_UK_empmr':
268
269        num_closea = 10
270        max_runoff_points = 1
271        use_runoff_box = True
272
273        ncsnr = np.zeros(num_closea+1,dtype=np.int)                     ; ncstt = np.zeros(num_closea+1,dtype=np.int)
274        ncsi1 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj1 = np.zeros(num_closea+1,dtype=np.int)
275        ncsi2 = np.zeros(num_closea+1,dtype=np.int)                     ; ncsj2 = np.zeros(num_closea+1,dtype=np.int)
276        ncsir1 = np.zeros(num_closea+1,dtype=np.int)                    ; ncsjr1 = np.zeros(num_closea+1,dtype=np.int)
277        ncsir2 = np.zeros(num_closea+1,dtype=np.int)                    ; ncsjr2 = np.zeros(num_closea+1,dtype=np.int)
278        ncsir = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int) ; ncsjr = np.zeros((num_closea+1,max_runoff_points+1),dtype=np.int)
279
280        # Caspian Sea
281        ncsnr[1]   = 1    ; ncstt[1]   = 0     
282        ncsi1[1]   = 1330 ; ncsj1[1]   = 831
283        ncsi2[1]   = 1375 ; ncsj2[1]   = 981
284
285        # Aral Sea
286        ncsnr[2]   = 1    ; ncstt[2]   = 0     
287        ncsi1[2]   = 1376 ; ncsj1[2]   = 900
288        ncsi2[2]   = 1400 ; ncsj2[2]   = 981
289
290        # Azov Sea
291        ncsnr[3]   = 1    ; ncstt[3]   = 0     
292        ncsi1[3]   = 1284 ; ncsj1[3]   = 908
293        ncsi2[3]   = 1304 ; ncsj2[3]   = 933
294
295        # Lake Superior
296        ncsnr[4]   = 1    ; ncstt[4]   = 2     
297        ncsi1[4]   = 781  ; ncsj1[4]   = 905
298        ncsi2[4]   = 815  ; ncsj2[4]   = 926 
299        # runff points the St Laurence Seaway for all Great Lakes
300        ncsir1[4:10]   = 873 ; ncsjr1[4:10]   = 909 
301        ncsir2[4:10]   = 884 ; ncsjr2[4:10]   = 920 
302
303        # Lake Michigan
304        ncsnr[5]   = 1    ; ncstt[5]   = 2     
305        ncsi1[5]   = 795  ; ncsj1[5]   = 871             
306        ncsi2[5]   = 813  ; ncsj2[5]   = 905 
307
308        # Lake Huron part 1
309        ncsnr[6]   = 1    ; ncstt[6]   = 2     
310        ncsi1[6]   = 814  ; ncsj1[6]   = 882             
311        ncsi2[6]   = 825  ; ncsj2[6]   = 905 
312
313        # Lake Huron part 2
314        ncsnr[7]   = 1    ; ncstt[7]   = 2     
315        ncsi1[7]   = 826  ; ncsj1[7]   = 889             
316        ncsi2[7]   = 833  ; ncsj2[7]   = 905 
317
318        # Lake Erie
319        ncsnr[8]   = 1    ; ncstt[8]   = 2     
320        ncsi1[8]   = 816  ; ncsj1[8]   = 871             
321        ncsi2[8]   = 837  ; ncsj2[8]   = 881 
322
323        # Lake Ontario
324        ncsnr[9]   = 1    ; ncstt[9]   = 2     
325        ncsi1[9]   = 831  ; ncsj1[9]   = 882             
326        ncsi2[9]   = 847  ; ncsj2[9]   = 889 
327
328        # Lake Victoria
329        ncsnr[10]   = 1    ; ncstt[10]   = 0   
330        ncsi1[10]   = 1274 ; ncsj1[10]   = 672 
331        ncsi2[10]   = 1289 ; ncsj2[10]   = 687 
332
333#=====================================
334# 3. Generate mask fields
335#=====================================
336
337    rnf_count = 0
338    empmr_count = 0
339
340    closea_mask = ma.zeros(top_level.shape,dtype=np.int)
341    temp_mask_rnf = ma.zeros(top_level.shape,dtype=np.int)
342    temp_mask_empmr = ma.zeros(top_level.shape,dtype=np.int)
343    closea_mask_rnf = ma.zeros(top_level.shape,dtype=np.int)
344    closea_mask_empmr = ma.zeros(top_level.shape,dtype=np.int)
345
346    for ics in range(num_closea):
347        closea_mask = ma.where( ( ii2d[:] >= ncsi1[ics+1] ) & ( ii2d[:] <= ncsi2[ics+1] ) &
348                                ( jj2d[:] >= ncsj1[ics+1] ) & ( jj2d[:] <= ncsj2[ics+1] ) &
349                                ( top_level == 1 ), ics+1, closea_mask)
350        if ncstt[ics+1] == 1:
351            rnf_count = rnf_count + 1
352            temp_mask_rnf[:] = 0
353            if use_runoff_box:
354                temp_mask_rnf = ma.where( ( ii2d[:] >= ncsir1[ics+1] ) & ( ii2d[:] <= ncsir2[ics+1] ) &
355                                       ( jj2d[:] >= ncsjr1[ics+1] ) & ( jj2d[:] <= ncsjr2[ics+1] ) &
356                                       ( top_level == 1 ), rnf_count, 0)
357            else:
358                for ir in range(ncsnr[ics+1]):
359                    temp_mask_rnf[ncsjr[ics+1],ncsjr[ics+1]] = rnf_count
360 
361            temp_mask_rnf = ma.where( closea_mask_rnf > 0, ma.minimum(temp_mask_rnf,closea_mask_rnf), temp_mask_rnf)
362            min_rnf = ma.amin(temp_mask_rnf[ma.where(temp_mask_rnf > 0)])
363            max_rnf = ma.amax(temp_mask_rnf[ma.where(temp_mask_rnf > 0)])
364            if min_rnf != max_rnf:
365                print 'min_rnf, max_rnf : ',min_rnf,max_rnf
366                raise Exception('Partially overlapping target rnf areas for two closed seas.')
367            else:
368                # source area:
369                closea_mask_rnf[ma.where(closea_mask==ics+1)] = min_rnf
370                # target area:
371                closea_mask_rnf[ma.where(temp_mask_rnf>0)] = min_rnf
372                # reset rnf_count:
373                rnf_count = min_rnf
374                   
375        if ncstt[ics+1] == 2:
376            empmr_count = empmr_count + 1
377            temp_mask_empmr[:] = 0
378            if use_runoff_box:
379                temp_mask_empmr = ma.where( ( ii2d[:] >= ncsir1[ics+1] ) & ( ii2d[:] <= ncsir2[ics+1] ) &
380                                          ( jj2d[:] >= ncsjr1[ics+1] ) & ( jj2d[:] <= ncsjr2[ics+1] ) &
381                                          ( top_level == 1 ), empmr_count, 0)
382            else:
383                for ir in range(ncsnr[ics+1]):
384                    temp_mask_empmr[ncsjr[ics+1],ncsjr[ics+1]] = empmr_count
385
386            temp_mask_empmr = ma.where( closea_mask_empmr > 0, ma.minimum(temp_mask_empmr,closea_mask_empmr), temp_mask_empmr)
387            min_empmr = ma.amin(temp_mask_empmr[ma.where(temp_mask_empmr > 0)])
388            max_empmr = ma.amax(temp_mask_empmr[ma.where(temp_mask_empmr > 0)])
389            if min_empmr != max_empmr:
390                raise Exception('Partially overlapping target empmr areas for two closed seas.')
391            else:
392                # source area:
393                closea_mask_empmr[ma.where(closea_mask==ics+1)] = min_empmr
394                # target area:
395                closea_mask_empmr[ma.where(temp_mask_empmr>0)] = min_empmr
396                # reset empmr_count:
397                empmr_count = min_empmr
398                   
399    if mask:
400        # apply land-sea mask if required
401        closea_mask.mask = np.where(top_level==0,True,False)
402        closea_mask_rnf.mask = np.where(top_level==0,True,False)
403        closea_mask_empmr.mask = np.where(top_level==0,True,False)
404
405#=====================================
406# 4. Append masks to domain_cfg file.
407#=====================================
408
409    domcfg.createVariable('closea_mask',datatype='i',dimensions=('y','x'),fill_value=closea_mask.fill_value,chunksizes=(1000,1000))
410    domcfg.variables['closea_mask'][:]=closea_mask
411    if rnf_count > 0:
412        domcfg.createVariable('closea_mask_rnf',datatype='i',dimensions=('y','x'),fill_value=closea_mask_rnf.fill_value,chunksizes=(1000,1000))
413        domcfg.variables['closea_mask_rnf'][:]=closea_mask_rnf
414    if empmr_count > 0:
415        domcfg.createVariable('closea_mask_empmr',datatype='i',dimensions=('y','x'),fill_value=closea_mask_empmr.fill_value,chunksizes=(1000,1000))
416        domcfg.variables['closea_mask_empmr'][:]=closea_mask_empmr
417
418    domcfg.close()
419
420
421if __name__=="__main__":
422    import argparse
423    parser = argparse.ArgumentParser()
424    parser.add_argument("-c", "--config", action="store",dest="config",default=None,
425                    help="configuration: eORCA1, eORCA025_UK")
426    parser.add_argument("-d", "--domcfg", action="store",dest="domcfg_file",default=None,
427                    help="domcfg file (input)")
428    parser.add_argument("-m", "--mask", action="store_true",dest="mask",default=False,
429                    help="mask output file based on top_level in domcfg file")
430
431    args = parser.parse_args()
432
433    make_closea_masks(config=args.config,domcfg_file=args.domcfg_file,mask=args.mask)
Note: See TracBrowser for help on using the repository browser.