source: TOOLS/MOZAIC/src/MOZAIC/tstmosaic.f90 @ 3363

Last change on this file since 3363 was 3326, checked in by omamce, 7 years ago

O.M. : Utility to generate interpolatio weights for OASIS-MCT
File size: 11.6 KB
Line 
1! -*- Mode: f90 -*-
2PROGRAM tstmosaic
3!!! -------------------------------------------------------------------------------
4!!! Tests des poids pour OASIS 2, Mozaic
5!!! Adapte a LMD5 et OPA
6!!! -------------------------------------------------------------------------------
7   !!
8   !! 1) Lit les fichiers de masques, grilles et surfaces des deux modeles
9   !! 2) Lit fichier de poids
10   !! 3) Lance une série de test
11   !!
12   !! --------------------------------------------------------------------------------
13   !! LMCE, CEA-Saclay
14   !! --------------------------------------------------------------------------------
15   !!
16   USE declare
17   USE modeles
18   USE mod_locio
19   USE fliocom
20   USE errioipsl
21   USE formula
22   USE mod_lbc
23   USE mod_inter
24   USE mod_lecdata
25   USE bords
26   USE mod_inipar
27   !!
28   IMPLICIT NONE
29   !!
30   INTEGER (kind=il), PARAMETER :: jpn = 16 ! Nombre de tests d'interpolation
31   INTEGER (kind=il) :: jpfirst = 8 ! Premier test avec bandes latitudes
32   !!
33   REAL (kind=rl)   , DIMENSION (:), ALLOCATABLE :: zoce, zo        ! Champ ocean
34   REAL (kind=rl)   , DIMENSION (:), ALLOCATABLE :: zatm, za        ! Champ atmosphere
35   REAL (kind=rl)   , DIMENSION (:), ALLOCATABLE :: zaoce           ! Fraction d'ocean dans la maille atmosphere
36   INTEGER (kind=il), DIMENSION (:), ALLOCATABLE :: iaoce
37   LOGICAL :: l_run = .FALSE.
38   LOGICAL :: l_land = .FALSE.
39   LOGICAL :: l_bord = .FALSE.
40   REAL (kind=rl) :: bila, bilo
41   REAL (kind=rd) :: zzd
42   CHARACTER (len=1) :: cw_wei = '2'
43   CHARACTER (LEN=3) :: lsens
44   CHARACTER (LEN=20) :: clarg
45   REAL (kind=rl) :: y_lat_min, y_lat_max
46   CHARACTER (LEN=1) :: c_lat_min, c_lat_max
47   INTEGER :: narg
48   !INTEGER, EXTERNAL :: iargc
49   !!
50   CHARACTER (len = 64) :: clout
51   INTEGER (kind = i_4) :: ncout
52   CHARACTER (len = 14) :: ctitle
53   !!
54   INTEGER (kind=il) :: ja, jo, jn, joi, joj, jai, jaj, ierr, jlat, jplat
55   REAL (kind=rl) :: dy, y_inf, y_sup
56   !! Initialisation
57   !CALL ipsldbg (.TRUE. )
58   ngrd = 11 ; nsrf = 12 ; nmsk = 13; nchk = 16 ; ndeb = 9 ; nbug1 = 10
59   nwei4o2a = 81 ; nwei8o2a = 82 ; nwei4a2o = 83 ; nwei8a2o = 84
60   !!
61   CALL inipar
62   CALL alloc_modeles
63   !!
64   ALLOCATE (zoce   (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'zoce', lreset=.TRUE., crout='tstmosaic')
65   ALLOCATE (zo     (jpon), STAT=ierr) ; CALL chk_allo (ierr, 'zo')
66   ALLOCATE (zatm   (jpan), STAT=ierr) ; CALL chk_allo (ierr, 'zatm')
67   ALLOCATE (za     (jpan), STAT=ierr) ; CALL chk_allo (ierr, 'za')
68   ALLOCATE (zaoce  (jpan), STAT=ierr) ; CALL chk_allo (ierr, 'zaoce')
69   ALLOCATE (iaoce  (jpan), STAT=ierr) ; CALL chk_allo (ierr, 'iaoce')
70   ALLOCATE (lacot  (jpan), STAT=ierr) ; CALL chk_allo (ierr, 'lacot')
71   ALLOCATE (laland (jpan), STAT=ierr) ; CALL chk_allo (ierr, 'laland')
72   !!
73   lwra2o = .FALSE. ; lwro2a = .TRUE.
74   !!
75   narg = iargc ()
76   DO ja = 1, narg
77      CALL getarg (ja, clarg)
78      SELECT CASE (TRIM(clarg))
79      CASE ('-r')
80         WRITE (UNIT=nout, FMT=*) 'Runoff case'
81         l_run = .TRUE.
82         cw_wei = '3'
83         jma2o = jma2or
84         lsens = 'a2o'
85      CASE ('-R')
86         l_run = .FALSE.
87      CASE ('-2') ; 
88         cw_wei = '2'
89      CASE ('-3') ; 
90         cw_wei = '3'
91         jma2o = jma2or
92      CASE ('-4') ;
93         cw_wei = '4'
94         jma2o = jma2or
95      CASE ('-5') ; 
96         cw_wei = '5'
97         jma2o = jma2oi
98      CASE ('-l') ;
99         l_land = .TRUE.
100      CASE ('-L') ;
101         l_land = .FALSE.
102      CASE ('-b')
103         l_bord = .TRUE.
104      CASE ('-B')
105         l_bord = .FALSE.
106      CASE ('-o')
107         lsens = 'o2a'
108      CASE ('-a')
109         lsens = 'a2o'
110      END SELECT
111   END DO
112   !!
113   !! Diag
114   !!
115   !! Fichiers de diagnostics
116   !!
117   OPEN (unit = nchk, file = 'tests.out', position = 'REWIND', action = 'write', status = 'replace' )
118   !!
119   !! Lecture des coordonnees
120   !!
121   CALL lecdata (lecweights = .TRUE., leco2amask = .TRUE., cw_num=cw_wei)
122   CALL bord_a
123   CALL bord_o
124   CALL bord_oa
125   !!
126   nva = COUNT (mask = (ko2a /= 0), DIM = 1)
127   nvo = COUNT (mask = (ka2o /= 0), DIM = 1)
128   WRITE (nout, *) 'MAX nva nvo ', MAXVAL(nva), MAXVAL (nvo)
129   WRITE (nout, *) 'MIN et MAX ko2a ', MINVAL (ko2a), MAXVAL (ko2a)
130   WRITE (nout, *) 'MIN et MAX ka2o ', MINVAL (ka2o), MAXVAL (ka2o)
131   !!
132   IF (MAXVAL (ko2a) .GT. jpon) WRITE (nout,*) 'Erreur :  MAXVAL (ko2a) .GT. jpon '
133   IF (MINVAL (ko2a) .LT. 0   ) WRITE (nout,*) 'Erreur :  MINVAL (ko2a) .LT. 0   '
134   IF (MAXVAL (ka2o) .GT. jpan) WRITE (nout,*) 'Erreur :  MAXVAL (ka2o) .GT. jpan '
135   IF (MINVAL (ka2o) .LT. 0   ) WRITE (nout,*) 'Erreur :  MINVAL (ka2o) .LT. 0   '
136   !!
137   !! Ecriture des champs dans un fichier NetCDF
138   clout = 'tstmosaic_' // TRIM (lsens)
139   CALL fliocrfd (TRIM(clout), (/ 'xa', 'ya', 'xo', 'yo'/), (/jpai, jpaj, jpoi, jpoj/), ncout, mode='REPLACE')
140   ! Quelques attributs globaux
141   CALL flioputa (ncout, '?', 'Conventions', 'COARDS')
142   !
143   CALL fliodefv (ncout, 'o2amask' , (/1_il, 2_il/), v_t=flio_r4)
144   CALL fliodefv (ncout, 'amask'   , (/1_il, 2_il/), v_t=flio_i1)
145   CALL fliodefv (ncout, 'tlmd_lon', (/1_il, 2_il/), v_t=flio_r4)
146   CALL fliodefv (ncout, 'tlmd_lat', (/1_il, 2_il/), v_t=flio_r4)
147   !
148   CALL fliodefv (ncout, 'omask'   , (/3_il, 4_il/), v_t=flio_i1)
149   CALL fliodefv (ncout, 'torc_lon', (/3_il, 4_il/), v_t=flio_r4)
150   CALL fliodefv (ncout, 'torc_lat', (/3_il, 4_il/), v_t=flio_r4)
151   !
152   CALL flioputv (ncout, 'omask'   , RESHAPE(iomskt ,(/jpoi, jpoj/)) )
153   CALL flioputv (ncout, 'torc_lon', RESHAPE(xolont ,(/jpoi, jpoj/)) )
154   CALL flioputv (ncout, 'torc_lat', RESHAPE(xolatt ,(/jpoi, jpoj/)) )
155   !
156   CALL flioputv (ncout, 'amask'   , RESHAPE(iamskt ,(/jpai, jpaj/)) )
157   CALL flioputv (ncout, 'o2amask' , RESHAPE(o2amask,(/jpai, jpaj/)) )
158   CALL flioputv (ncout, 'tlmd_lon', RESHAPE(xalont ,(/jpai, jpaj/)) )
159   CALL flioputv (ncout, 'tlmd_lat', RESHAPE(xalatt ,(/jpai, jpaj/)) )
160   !!
161   CALL inter_o2a (REAL (1-iomskt, kind = rl), zaoce) 
162   iaoce = 1
163   WHERE (zaoce .GT. 0) iaoce = 0
164   !!
165   !! Point cotier Atm
166   lacot = .FALSE.
167   WHERE (o2amask .GT. epsfrac .AND. o2amask .LT. (1.0_rl-epsfrac) ) lacot = .TRUE.
168   laland = .FALSE.
169   WHERE (o2amask .LT. (1.0_rl-epsfrac) ) laland = .TRUE.
170   !!
171   jplat = jpn-jpfirst
172   y_inf = MAX (-90.1_rl, MIN( yo_inf, ya_inf ) - 2.0_rl)
173   y_sup = MIN ( 90.1_rl, MAX (yo_sup, ya_sup ) + 2.0_rl)
174   dy = (y_sup-y_inf)/REAL(jplat+1,rl)
175   WRITE (nout,*) 'y_inf, y_suf, dy ; ', y_inf, y_sup, dy
176
177   !! Verifs
178   WRITE (UNIT=nout, FMT=*)   '    Nom               Int. Atm      Int. oce   Erreur'
179   ctitle = ' '
180   LoopCase: DO jn = 1, jpn
181      !WRITE (UNIT=nout, FMT=*) jn
182      IF ( lsens == 'o2a' ) THEN
183         SELECT CASE (jn)
184         CASE (1_il) ! Oceanic field = constant
185            zoce = 1.0_rl ; zatm = 0.0_rl
186            ctitle = 'Oce_constant'
187         CASE (2_il) ! Oceanic field = masked constant
188            zoce = REAL (1-iomskt, kind = rl)
189            ctitle = 'Oce_masked'
190         CASE (3_il) ! Oceanic field = latitude
191            zoce = ABS(xolatt) * REAL (1-iomskt, kind = rl)
192            ctitle = 'Oce_lat'
193         CASE (4_il) ! Oceanic field = longitude
194            zoce = 2.0_rl + COS (rad*xolont)**2 * REAL (1-iomskt, kind = rl)
195            ctitle = 'Oce_lon'
196         CASE (5_il) ! Oceanic field = COS(longitude) * latitude
197            zoce = 2.0_rl + ABS(xolatt) * COS (rad*xolont)**2 * REAL (1-iomskt, kind = rl)
198            ctitle = 'Oce_lonxlat'
199         CASE (6_il) ! Oce field = random
200            CALL RANDOM_NUMBER (zoce)
201            zoce = 1.0_rl + zoce * REAL (1-iomskt, kind = rl)
202            ctitle = 'Oce_random'
203         CASE (7_il) ! Oce field = 1 point
204            zoce = 0.0_rl 
205            zoce (jpon/2)  = 1.0_rl
206            ctitle = 'Oce_point'
207         CASE Default ! Oce field = bande de latitude
208            jlat  = jn-jpfirst
209            y_lat_min = REAL( jlat,rl)*dy + y_inf
210            y_lat_max = y_lat_min + dy
211            !WRITE (nout,*) jn, jlat, y_lat_min, y_lat_max
212            c_lat_min='S' ; c_lat_max='S'
213            IF ( y_lat_min >= 0.0_rl ) c_lat_min = 'N'
214            IF ( y_lat_max >= 0.0_rl ) c_lat_max = 'N'
215            zoce = 0.0_rl
216            WHERE ( xolatt .GE. y_lat_min .AND. xolatt .LE. y_lat_max ) zoce = 1.0_rl
217            zoce = zoce * REAL (1-iomskt, kind = rl)
218            WRITE (UNIT=ctitle, FMT='("Oce_", 1I2.2, 1A, "_", 1I2.2, 1A)') &
219               & INT(ABS(y_lat_min)), c_lat_min, INT(ABS(y_lat_max)), c_lat_max
220         END SELECT
221
222         CALL lbc (zoce, noperio, 'T', jpoi, jpoj)
223         zatm (:) = 0.0_rl
224         CALL inter_o2a (zoce, zatm)
225
226      END IF
227      IF ( lsens == 'a2o' ) THEN
228         SELECT CASE (jn)
229         CASE (1_il) ! Atm field = constant
230            zatm = 1.0_rl
231            ctitle = 'Atm_constant'
232         CASE (2_il) ! Atm field = masked constant
233            zatm = REAL (1-iaoce, kind = rl)
234            ctitle = 'Atm_masked'
235         CASE (3_il) ! Atm field = latitude
236            zatm = ABS(xalatt) * REAL (1-iaoce, kind = rl)
237            ctitle = 'Atm_lat'
238         CASE (4_il) ! Atm field = longitude * latitude
239            zatm= 2.0_rl + COS (rad*xalont)**2 * ABS(xalatt) * REAL (1-iaoce, kind = rl)
240            ctitle = 'Atm_lonxlat'
241         CASE (5_il) ! Atm field = longitude
242            zatm = 2.0_rl + COS (rad*xalont)**2 * REAL (1-iaoce, kind = rl)
243            ctitle = 'Atm_lon'
244         CASE (6_il) ! Atm field = random
245            CALL RANDOM_NUMBER (zatm)
246            zatm = 1.0_rl + zatm * REAL (1-iaoce, kind = rl)
247            ctitle = 'Atm_random'
248         CASE (7_il) ! Atm field = 1 point
249            zatm = 0.0_rl
250            zatm (1:jpan/2) = 1.0_rl
251            ctitle = 'Atm_north'
252         CASE Default ! Oce field = bande de latitude
253            jlat  = jn-jpfirst
254            y_lat_min = REAL( jlat,rl)*dy + y_inf
255            y_lat_max = y_lat_min + dy
256            c_lat_min='S' ; c_lat_max='S'
257            IF ( y_lat_min >= 0.0_rl ) c_lat_min = 'N'
258            IF ( y_lat_max >= 0.0_rl ) c_lat_max = 'N'
259            zatm = 0.0_rl
260            WHERE ( xalatt .GT. y_lat_min .AND. xalatt .LT. y_lat_max ) zatm = 1.0_rl
261            zatm = zatm * REAL (1-iaoce, kind = rl)
262            zoce = 0.0_rl
263            WRITE (UNIT=ctitle, FMT='("Atm_", 1I2.2, 1A, "_", 1I2.2, 1A)') &
264               & INT(ABS(y_lat_min)), c_lat_min, INT(ABS(y_lat_max)), c_lat_max
265         END SELECT
266
267         IF (l_land) WHERE (.NOT. laland) zatm = 0.0_rl
268         IF (l_bord) WHERE (.NOT. lacot ) zatm = 0.0_rl
269         CALL lbc (zatm, naperio, 'T', jpai, jpaj)
270         CALL inter_a2o (zatm, zoce)
271
272      END IF
273      !!
274      SELECT CASE (lsens) 
275      CASE ('a2o')
276         IF (l_run) THEN
277            bila = DOT_PRODUCT (zatm     , REAL((1-iamskp)*(1-iamskt),kind=rl))
278         ELSE
279            bila = DOT_PRODUCT (zatm*xasrft*o2amask, REAL((1-iamskp)*(1-iamskt),kind=rl))
280         END IF
281      CASE ('o2a')
282         bila = DOT_PRODUCT (zatm*xasrft*o2amask, REAL((1-iamskp)*(1-iamskt),kind=rl))
283      END SELECT
284      bilo = DOT_PRODUCT (zoce*xosrft, REAL((1-iomskp)*(1-iomskt),kind=rl))
285      !
286      IF ( bila+bilo .EQ. 0.0_rl ) THEN
287         WRITE (nout, FMT= '(1I3, 1A14, "  ", 2(1PE12.2, "  "), "  ", 1PE8.2)' ) &
288            & jn, TRIM (ctitle), bila, bilo
289      ELSE
290         WRITE (nout, FMT= '(1I3, 1A14, "  ", 2(1PE12.2, "  "), "  ", 1PE8.2)' ) &
291            & jn, TRIM (ctitle), bila, bilo, ABS((bila-bilo)/(bila+bilo))*2.0_rl
292      END IF
293      CALL fliodefv (ncout, TRIM (ctitle)//'_OCE', (/3, 4/), v_t=flio_r4)
294      !
295      CALL fliodefv (ncout, TRIM(ctitle)//'_ATM', (/1, 2/), v_t=flio_r4)
296      !
297      CALL flioputv (ncout, TRIM (ctitle)//'_OCE', RESHAPE (zoce, (/jpoi, jpoj/)) ) 
298      CALL flioputv (ncout, TRIM (ctitle)//'_ATM', RESHAPE (zatm, (/jpai, jpaj/)) ) 
299
300      !!
301   END DO LoopCase
302   !!
303   CALL flioclo (ncout) 
304   !!
305   STOP
306   !!
307END PROGRAM tstmosaic
308!!
Note: See TracBrowser for help on using the repository browser.