source: CPL/oasis3/trunk/src/mod/oasis3/src/filling.f @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago

Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).

File size: 17.8 KB
Line 
1      SUBROUTINE filling (pfild, pworka, pworkb, pworkc,
2     $                    plon, plat, klon, klat, kmask, kmesh,
3     $                    kunit, cdfic, cdmet)
4C****
5C               *****************************
6C               * OASIS ROUTINE  -  LEVEL 3 *
7C               * -------------     ------- *
8C               *****************************
9C
10C**** *filling* - Fill up (and smooth or not) regional field with climatology 
11C
12C
13C     Purpose:
14C     -------
15C     Performs smooth or abrupt blending of a regional data set with a 
16C     global one for a Sea Surface Temperature or Sea Ice Extent field.
17C     The frequency of the global set can be interannual monthly, 
18C     climatological monthly or yearly. Sea-ice transitions are treated
19C     specifically: the ice is supposed to form or disappear on the first 
20C     day of the month.
21C
22C     N.B : The SST climatology must be in celsius. If not
23C           one has to change the value of local variable zmerg
24C
25C**   Interface:
26C     ---------
27C       *CALL*  *filling(pfild, pworka, pworkb, pworkc, klon, klat, 
28C                        kmask, kmesh, kunit, cdfic, cdmet)*
29C
30C     Input:
31C     -----
32C                pfild  : field to be completed (real 2D) 
33C                plon   : grid longitude array (real 2D)
34C                plat   : grid latitude array (real 2D)
35C                klon   : number of longitude (integer)
36C                klat   : number of latitude (integer)
37C                kmask  : mask array (integer 2D)
38C                kmesh  : overlapped neighbors array (integer 2D)
39C                kunit  : climatological data filename logical unit (integer)
40C                cdfic  : climatological data filename  (character)
41C                cdmet  : filling method (character)
42C
43C     Output:
44C     ------
45C                pfild  : field completed (real 2D) 
46C
47C     Workspace:
48C     ---------
49C     These are work arrays passed as arguments:
50C     pworka, pworkb, pworkc
51C
52C     Externals:
53C     ---------
54C     None
55C
56C     Reference:
57C     ---------
58C     See OASIS manual (1995)
59C
60C     History:
61C     -------
62C       Version   Programmer     Date      Description
63C       -------   ----------     ----      ----------- 
64C       1.0       L. Terray      94/01/01  created
65C       1.1       L. Terray      94/10/01  modified: lecture of climagrd
66C       2.0beta   L. Terray      95/12/10  modified: new structure
67C       2.0       L. Terray      96/02/01  modified: change on rewind()
68C       2.1       L. Terray      96/09/09  modified: flux correction term
69C       2.3       S. Valcke      99/04/30  added: printing levels
70C
71C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72C
73C* ----------------Include files and USE of modules---------------------------
74C
75      USE mod_kinds_oasis
76      USE mod_parameter
77      USE mod_analysis
78      USE mod_coast
79      USE mod_calendar
80      USE mod_unit
81      USE mod_smooth
82      USE mod_printing
83C
84C* ---------------------------- Argument declarations----------------------
85C
86      REAL (kind=ip_realwp_p) pfild(klon, klat), plon(klon, klat), 
87     $    plat(klon, klat)
88      REAL (kind=ip_realwp_p) pworka(klon, klat),  pworkb(klon, klat) 
89      REAL (kind=ip_realwp_p) pworkc(klon, klat)
90      INTEGER (kind=ip_intwp_p) kmask(klon,klat), kmesh(klon, klat)
91      CHARACTER*8 cdmet, cdfic
92C
93C* ---------------------------- Local declarations ----------------------
94C
95      CHARACTER*3 clsmooth, clfield
96      CHARACTER*2 clclim
97C
98C* ---------------------------- Poema verses ----------------------------
99C
100C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
101C
102C*    1. Initialization
103C        --------------
104C
105      IF (nlogprt .GE. 2) THEN
106          WRITE (UNIT = nulou,FMT = *) ' '
107          WRITE (UNIT = nulou,FMT = *) ' '
108          WRITE (UNIT = nulou,FMT = *) 
109     $    '           ROUTINE filling  -  Level 3'
110          WRITE (UNIT = nulou,FMT = *) 
111     $    '           ***************     *******'
112          WRITE (UNIT = nulou,FMT = *) ' '
113          WRITE (UNIT = nulou,FMT = *) 
114     $     ' Filling and smoothing of field with climatology'
115          WRITE (UNIT = nulou,FMT = *) ' '
116          WRITE (UNIT = nulou,FMT = *) ' '
117      ENDIF
118C
119C* initialize error flag for error routine
120C
121      iflag = 0
122C
123C* Get smoothing and time info
124C
125      clsmooth = cdmet(1:3)
126      clfield = cdmet(4:6)
127      clclim = cdmet(7:8)
128C
129C* Get congelation temperature for SST and SIE treatment.
130C
131      zmerg = -1.93
132C
133C
134C*    2. Getting climatological or interannual field
135C        -------------------------------------------
136C
137C* Monthly climatology
138C
139      IF (nlogprt .GE. 2) THEN
140          CALL prcout
141     $    ('Complementary data file is file =', cdfic, 1)
142      ENDIF
143      IF (clclim .EQ. 'SE' .OR. clclim .EQ. 'MO') THEN
144          IF (nlogprt .GE. 2) THEN
145              WRITE (UNIT = nulou,FMT = *) ' '
146          ENDIF
147C
148C* Monthly climatology
149C
150          IF (clclim .EQ. 'SE') THEN
151              IF (nlogprt .GE. 2) THEN
152                  CALL prtout
153     $            ('Monthly records skipped in climatology file
154     $            nsrec =', nsrec, 1)
155              ENDIF
156C
157C* Read the right records
158C
159              REWIND kunit 
160              IF (nsrec .eq. 1) THEN
161                  READ (kunit) pworka
162                  READ (kunit) pworkb
163                ELSEIF (nsrec .gt. 1 .and. nsrec .lt. 12) THEN
164                  DO 210 ji = 1, nsrec - 1
165                    READ (kunit) pworka
166 210              CONTINUE
167                  READ (kunit) pworka
168                  READ (kunit) pworkb
169                ELSEIF (nsrec .eq. 12 .or. nsrec .eq. 0) THEN
170                  DO 220 ji = 1, 11
171                    READ (kunit) pworka
172 220              CONTINUE
173                  READ (kunit) pworka
174                  REWIND kunit 
175                  READ (kunit) pworkb
176                ELSE
177                  WRITE (UNIT = nulou,FMT = *) '        ***WARNING***'
178                  WRITE (UNIT = nulou,FMT = *) 
179     $                ' ===>>> : wrong value for nsrec '
180                  WRITE (UNIT = nulou,FMT = *) 
181     $                ' ======                   ===== '
182                  WRITE (UNIT = nulou,FMT = *) ' '
183                  WRITE (UNIT = nulou,FMT = *) 
184     $                ' Variable nsrec = ',nsrec
185                  WRITE (UNIT = nulou,FMT = *) 
186     $                ' We STOP !!! Check value for variable nsrec '
187                  CALL HALTE ('STOP in filling')
188              ENDIF
189C
190C* Monthly interannual data
191C
192            ELSE IF(clclim .EQ. 'MO') THEN
193                IF (nlogprt .GE. 2) THEN
194                    CALL prtout
195     $                  ('Monthly records skipped in interannual file
196     $                  nmrec =', nmrec, 1)
197                ENDIF
198C
199C* Read the right records
200C
201              REWIND kunit
202              IF (nmrec .eq. 0) THEN
203                  READ (kunit) pworka
204                  READ (kunit) pworkb
205                ELSEIF (nmrec .ge. 1) THEN
206                  DO 230 ji = 1, nmrec
207                    READ (kunit) pworka
208 230              CONTINUE
209                  READ (kunit) pworka
210                  READ (kunit) pworkb
211              ENDIF
212          ENDIF 
213C
214C* Get all coefficients for time linear interpolation
215C
216          zpoid1 = FLOAT(njtwo - njnow) / FLOAT(njtwo - njone)
217          zpoid2 = 1.0 - zpoid1
218          zpoib1 = FLOAT(ndtwo - njnow) / FLOAT(ndtwo - ndone)
219          zpoib2 = 1.0 - zpoib1
220C
221C* Get climatological field for current day
222C
223C* Specific treatment for SST and SIE due to freezing temperature
224C  In both cases, one reads the temperature file and looks at 
225C  sea-ice transitions. One makes the assumption that the ice is
226C  formed or disappears on the first of the month.
227C
228          DO 240 jj = 1, klat
229            DO 250 ji = 1, klon
230C
231C* General case
232C
233              pworkc(ji,jj) = zpoid1 * pworka(ji,jj) +
234     $                        zpoid2 * pworkb(ji,jj)
235C
236C* If field is SST
237C
238              IF (clfield .EQ. 'SST') THEN
239C
240C* Sea Ice to Sea
241C
242                  IF (pworka(ji,jj) .LE. zmerg .AND. 
243     $                pworkb(ji,jj) .GT. zmerg) THEN
244                      IF (njnow .GT. 15) THEN
245                          CONTINUE
246                      ELSE
247                          pworkc(ji,jj) = zmerg * zpoib1 +
248     $                                    pworkb(ji,jj) * zpoib2
249                      ENDIF
250                  ENDIF
251C
252C* Sea to Sea Ice
253C
254                  IF (pworka(ji,jj) .GT. zmerg .AND. 
255     $                pworkb(ji,jj) .LE. zmerg) THEN
256                      IF (njnow .GT. 15) THEN
257                          pworkc(ji,jj) = pworka(ji,jj) * zpoib1 +
258     $                                    zmerg * zpoib2
259                      ENDIF
260                  ENDIF
261C
262C* Sea to Sea
263C
264                  IF (pworka(ji,jj) .GT. zmerg .AND. 
265     $                pworkb(ji,jj) .GT. zmerg) THEN
266                      pworkc(ji,jj) = zpoid1 * pworka(ji,jj) +
267     $                                zpoid2 * pworkb(ji,jj)
268                  ENDIF
269              ENDIF
270C
271C* If field is SIE
272C 
273              IF (clfield .EQ. 'SIE') THEN             
274C
275C* Sea Ice to Sea Ice
276C
277                  IF (pworka(ji,jj) .LE. zmerg .AND. 
278     $                pworkb(ji,jj) .LE. zmerg) THEN
279                      pworkc(ji,jj) = 1.0
280                  ENDIF
281C
282C* Sea Ice to Sea
283C
284                  IF (pworka(ji,jj) .LE. zmerg .AND. 
285     $                pworkb(ji,jj) .GT. zmerg) THEN
286                      IF (njnow .gt. 15) THEN
287                          pworkc(ji,jj) = 1.0
288                        ELSE
289                          pworkc(ji,jj) = 0.0
290                      ENDIF
291                  ENDIF
292C
293C* Sea to Sea Ice
294C
295                  IF (pworka(ji,jj) .GT. zmerg .AND. 
296     $                pworkb(ji,jj) .LE. zmerg) THEN
297                      IF (njnow .gt. 15) THEN
298                          pworkc(ji,jj) = 0.0
299                      ELSEIF (njnow .eq. 1) THEN
300                          pworkc(ji,jj) = 2.0
301                      ELSE
302                          pworkc(ji,jj) = 1.0
303                      ENDIF
304                  ENDIF
305C
306C* Sea to Sea
307C
308                  IF (pworka(ji,jj) .GT. zmerg .AND. 
309     $                pworkb(ji,jj) .GT. zmerg) THEN
310                      pworkc(ji,jj) = 0.0
311                  ENDIF
312              ENDIF
313 250        CONTINUE
314 240      CONTINUE
315C
316C* Annual climatology
317C
318      ELSE IF (clclim .EQ. 'AN') THEN
319          READ (kunit) pworkc
320      ELSE
321          CALL prcout ('This type of climatology cannot be used  -
322     $                 clim = ', clclim, 1)
323      ENDIF
324C
325C* Blending and smoothing
326C
327C* Storing of climatology field interpolated in intermediate array pworka
328C
329      DO 260 jj = 1, klat
330        DO 270 ji = 1, klon
331          pworka(ji,jj) = pworkc(ji,jj)
332 270    CONTINUE
333 260  CONTINUE
334C* SST case
335      IF (clfield .EQ. 'SST') THEN
336C
337C* Put interpolated values in intermediate array if:
338C     - the point is an atmosphere sea point
339C     - there is at least one underlying ocean sea point
340C
341          DO 280 jj = 1, klat
342            DO 290 ji = 1, klon
343              IF (kmesh(ji,jj) .NE. 0 .AND. kmask(ji,jj) .EQ. 0) THEN
344                  pworka(ji,jj) = pfild(ji,jj)
345              ENDIF
346 290        CONTINUE
347 280      CONTINUE
348C
349C* If needed, do coast mismatch correction
350C  --> Special treatment for coastal points where there is a
351C  --> a mismatch between ocean and atmosphere land sea masks
352C  --> The array npcoast is calculated in routine coasts.f
353C  --> icoor and jcoor are externally defined functions
354C
355          IF(nfcoast .EQ. 1) THEN
356C* If first iteration, initialize coast mismatch data
357              IF (lcoast) THEN
358                  itmsq = 1
359                  CALL coasts (plon, plat, kmask, itmsq,
360     $                         klon, klat, kmesh)
361C* Switch off flag
362                  lcoast = .FALSE.
363              ENDIF
364              IF (nlogprt .GE. 2) THEN
365                  WRITE (UNIT = nulan,FMT = *) ' '
366                  WRITE (UNIT = nulan,FMT = *) 
367     $            '               * Model run date * '
368                  WRITE (UNIT = nulan,FMT = *) 
369     $            '               * -------------- * '
370                  WRITE (UNIT = nulan,FMT = *) ' '
371                  WRITE (UNIT = nulan,FMT = *) 
372     $            '        Year  --->>> ',nanow
373                  WRITE (UNIT = nulan,FMT = *) 
374     $            '       Month  --->>> ',nmnow
375                  WRITE (UNIT = nulan,FMT = *) 
376     $            '         Day  --->>> ',njnow
377                  WRITE (UNIT = nulan,FMT = *) ' '
378                  WRITE (UNIT = nulan,FMT = *) 
379     $            '      Correct for coast mismatch '
380                  WRITE (UNIT = nulan,FMT = *) 
381     $            '      ************************** '
382                  WRITE (UNIT = nulan,FMT = *) ' '
383              ENDIF
384              DO 291 jc = 1, ncoast
385                iic = icoor(npcoast(jc,1), klon)
386                ijc = jcoor(npcoast(jc,1), klon)
387                ztmptot = 0.
388                DO 292 jn = 3, 2 + npcoast(jc,2)
389                  iiic = icoor(npcoast(jc,jn), klon)
390                  ijjc = jcoor(npcoast(jc,jn), klon)
391                  ztmptot = ztmptot + pfild(iiic,ijjc)
392 292            CONTINUE
393C
394C* Print old and new values
395C
396                IF (nlogprt .GE. 2) THEN
397                    WRITE (UNIT = nulan,FMT = *) 
398     $              ' Resetting SST at point (i,j)= ',iic,ijc
399                    WRITE (UNIT = nulan,FMT = *) 
400     $              ' Old SST value : SST = ', pfild(iic,ijc)
401                    WRITE (UNIT = nulan,FMT = *) 
402     $              ' Clim SST value : SST = ', pworka(iic,ijc)
403                ENDIF
404C
405C* Get new value
406C
407                pfild(iic,ijc) = ztmptot / float(npcoast(jc,2))
408                IF (nlogprt .GE. 2) THEN
409                    WRITE (UNIT = nulan,FMT = *) 
410     $              ' New SST value : SST = ', pfild(iic,ijc)
411                ENDIF
412C
413C* Copy new value into intermediate array
414C
415                pworka(iic,ijc) = pfild(iic,ijc)
416 291          CONTINUE 
417          ENDIF             
418C
419C* SIE case : we assign the new field to climatology i.e we do nothing
420C
421      ELSE IF (clfield .EQ. 'SIE') THEN
422          CONTINUE
423C
424C* General case: put interpolated value if there is an
425C                underlying grid point (whether land or sea)
426C
427      ELSE 
428          DO 293 jj = 1, klat
429            DO 294 ji = 1, klon
430              IF (kmesh(ji,jj) .NE.  0) THEN
431                  pworka(ji,jj) = pfild(ji,jj)
432              ENDIF
433 294        CONTINUE
434 293      CONTINUE
435      ENDIF
436C
437C
438C*    3. Do the smoothing on  borders (only for SST)
439C        ----------------------------
440C
441      IF (clfield .EQ. 'SST' .AND. clsmooth .EQ. 'SMO') THEN
442C
443C* First South and North borders
444C
445          DO 310 jk = nsltb, nslte
446            DO 320 ji = 1, klon
447              IF (kmesh(ji,jk) .NE. 0 .AND. kmask(ji,jk) .EQ. 0) THEN
448                  pworka(ji,jk) = qalfa * (jk - nsltb) * pfild(ji,jk)
449     $                + (1.-qalfa * (jk - nsltb)) * pworkc(ji,jk)
450              ENDIF
451 320        CONTINUE
452 310      CONTINUE
453          DO 330 jk = nnltb, nnlte, -1
454            DO 340 ji = 1, klon
455              IF (kmesh(ji,jk) .NE. 0 .AND. kmask(ji,jk) .EQ. 0) THEN
456                  pworka(ji,jk) = qalfa * (nnltb - jk) * pfild(ji,jk)
457     $                +  (1.-qalfa * (nnltb - jk)) * pworkc(ji,jk)
458              ENDIF
459 340        CONTINUE
460 330      CONTINUE
461C
462C* Western border
463C
464          DO 350 jj = nslte + 1, nnlte - 1
465            DO 360 ji = 1, nwlgmx
466              IF (kmesh(ji,jj) .EQ. 0 .AND. 
467     $            kmesh(ji + 1,jj) .NE. 0) THEN
468                  DO 370 jk = 1, nliss
469                    pworka(ji + jk - 1,jj) = qbeta*(jk - 1) *
470     $                                       pfild(ji + jk - 1,jj) +
471     $                                       pworkc(ji + jk - 1,jj) *
472     $                                       (1.-qbeta*(jk - 1))
473 370              CONTINUE
474              ENDIF
475 360        CONTINUE
476 350      CONTINUE
477C
478C* Eastern border
479C
480          DO 380 jj = nslte + 1, nnlte - 1
481            DO 390 ji = klon, nelgmx, -1
482              IF (kmesh(ji,jj) .EQ. 0 .AND. 
483     $            kmesh(ji - 1,jj) .NE. 0) THEN
484                  DO 391 jk = 1, nliss
485                    pworka(ji - jk + 1,jj) = qbeta*(jk - 1) *
486     $                                       pfild(ji - jk + 1,jj) +
487     $                                       pworkc(ji - jk + 1,jj) *
488     $                                       (1.-qbeta*(jk - 1))
489 391              CONTINUE
490              ENDIF
491 390        CONTINUE
492 380      CONTINUE
493C
494C* Calculate and write flux correction term cfldcor on unit nlucor
495C
496          CALL prcout('Calculate flux correction term', cfldcor, 1)
497          CALL szero(pworkb,klon*klat)
498          DO 392 jj = 1, klat
499            DO 393 ji = 1, klon
500C
501C* The correction term is estimated only for atmosphere sea points with 
502C  underlying ocean sea points
503C
504              IF (kmesh(ji,jj) .NE. 0 .AND. kmask(ji,jj) .EQ. 0) THEN
505                  pworkb(ji,jj) = pfild(ji,jj) - pworka(ji,jj)
506              ENDIF
507 393        CONTINUE
508 392      CONTINUE
509C
510C* Rewind unit nlucor and write
511C
512          IF (nlogprt .GE. 2) THEN
513              CALL prtout('Write it on logical unit ', nlucor, 2)
514          ENDIF
515          REWIND nlucor
516          CALL locwrite(cfldcor, pworkb, klon*klat, nlucor, iflag)
517          IF (iflag .NE. 0) THEN
518              CALL prcout
519     $            ('WARNING: problem in writing field',
520     $            cfldcor, 1)
521              CALL prtout
522     $            ('Error in writing logical unit', nlucor, 2)
523              CALL HALTE('STOP in filling')
524          ENDIF
525      ENDIF
526C
527C* Assign intermediate array pworka values to final array pfild
528C
529      DO 394 jj = 1, klat
530        DO 395 ji = 1, klon
531          pfild(ji,jj) = pworka(ji,jj)
532 395    CONTINUE
533 394  CONTINUE
534C
535C
536C*    4. End of routine
537C        --------------
538C
539      IF (nlogprt .GE. 2) THEN
540          WRITE (UNIT = nulou,FMT = *) ' '
541          WRITE (UNIT = nulou,FMT = *) 
542     $    '          --------- End of routine filling ---------'
543          CALL FLUSH (nulou)
544      ENDIF
545      RETURN
546      END
Note: See TracBrowser for help on using the repository browser.