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.
ldfdyn_c2d.h90 in branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF – NEMO

source: branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 1839

Last change on this file since 1839 was 1839, checked in by edblockley, 14 years ago

1st commit for the mixed_ldfdyn branch
see ticket:659
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 11.1 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                      ***  ldfdyn_c2d.h90  ***
3   !!----------------------------------------------------------------------
4   !!   ldf_dyn_c2d  : set the lateral viscosity coefficients
5   !!   ldf_dyn_c2d_orca : specific case for orca r2 and r4
6   !!----------------------------------------------------------------------
7
8   !!----------------------------------------------------------------------
9   !!   OPA 9.0 , LOCEAN-IPSL (2005)
10   !! $Id$
11   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
12   !!----------------------------------------------------------------------
13
14   SUBROUTINE ldf_dyn_c2d( ld_print )
15      !!----------------------------------------------------------------------
16      !!                 ***  ROUTINE ldf_dyn_c2d  ***
17      !!                 
18      !! ** Purpose :   initializations of the horizontal ocean physics
19      !!
20      !! ** Method :
21      !!      2D eddy viscosity coefficients ( longitude, latitude )
22      !!
23      !!       harmonic operator   : ahm1 is defined at t-point
24      !!                             ahm2 is defined at f-point
25      !!           + isopycnal     : ahm3 is defined at u-point
26      !!           or geopotential   ahm4 is defined at v-point
27      !!           iso-model level : ahm3, ahm4 not used
28      !!
29      !!       biharmonic operator : ahm3 is defined at u-point
30      !!                             ahm4 is defined at v-point
31      !!                           : ahm1, ahm2 not used
32      !!
33      !!----------------------------------------------------------------------
34      !! * Arguments
35      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
36
37      !! * Local variables
38      INTEGER :: ji, jj
39      REAL(wp) ::   za00, zd_max, zetmax, zeumax, zefmax, zevmax
40      !!----------------------------------------------------------------------
41
42      IF(lwp) WRITE(numout,*)
43      IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient'
44      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
45      IF(lwp) WRITE(numout,*)
46
47      ! harmonic operator (ahm1, ahm2) : ( T- and F- points) (used for laplacian operators
48      ! ===============================                       whatever its orientation is)
49      IF( ln_dynldf_lap ) THEN
50         ! define ahm1 and ahm2 at the right grid point position
51         ! (USER: modify ahm1 and ahm2 following your desiderata)
52
53         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
54         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
55
56         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
57         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
58
59         za00 = ahm0 / zd_max
60         DO jj = 1, jpj
61            DO ji = 1, jpi
62               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
63               zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
64               ahm1(ji,jj) = za00 * zetmax
65               ahm2(ji,jj) = za00 * zefmax
66            END DO
67         END DO
68
69         IF( ln_dynldf_iso ) THEN
70            IF(lwp) WRITE(numout,*) '              Caution, as implemented now, the isopycnal part of momentum'
71            IF(lwp) WRITE(numout,*) '                 mixing use aht0 as eddy viscosity coefficient. Thus, it is'
72            IF(lwp) WRITE(numout,*) '                 uniform and you must be sure that your ahm is greater than'
73            IF(lwp) WRITE(numout,*) '                 aht0 everywhere in the model domain.'
74         ENDIF
75
76         ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahm1 ahm2)
77         ! ==============================================
78         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print )
79
80         ! Control print
81         IF( lwp .AND. ld_print ) THEN
82            WRITE(numout,*)
83            WRITE(numout,*) 'inildf: 2D ahm1 array'
84            CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
85            WRITE(numout,*)
86            WRITE(numout,*) 'inildf: 2D ahm2 array'
87            CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
88         ENDIF
89      ENDIF
90
91      ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator
92      ! =================================                      whatever its orientation is)
93      IF( ln_dynldf_bilap ) THEN
94         ! (USER: modify ahm3 and ahm4 following your desiderata)
95         ! Here: ahm is proportional to the cube of the maximum of the gridspacing
96         !       in the to horizontal direction
97
98         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
99         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
100
101         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
102         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
103
104         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
105         DO jj = 1, jpj
106            DO ji = 1, jpi
107               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
108               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
109               ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax
110               ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax
111            END DO
112         END DO
113
114         ! Control print
115         IF( lwp .AND. ld_print ) THEN
116            WRITE(numout,*)
117            WRITE(numout,*) 'inildf: ahm3 array'
118            CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
119            WRITE(numout,*)
120            WRITE(numout,*) 'inildf: ahm4 array'
121            CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
122         ENDIF
123      ENDIF
124
125
126   END SUBROUTINE ldf_dyn_c2d
127
128
129   SUBROUTINE ldf_dyn_c2d_orca( ld_print )
130      !!----------------------------------------------------------------------
131      !!                 ***  ROUTINE ldf_dyn_c2d  ***
132      !!
133      !!                   **** W A R N I N G ****
134      !!
135      !!                ORCA R2 and R4 configurations
136      !!                 
137      !!                   **** W A R N I N G ****
138      !!                 
139      !! ** Purpose :   initializations of the lateral viscosity for orca R2
140      !!
141      !! ** Method  :   blah blah blah...
142      !!
143      !!----------------------------------------------------------------------
144      !! * Modules used
145      USE ldftra_oce, ONLY : aht0
146
147      !! * Arguments
148      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
149
150      !! * Local variables
151      INTEGER ::   ji, jj, jn      ! dummy loop indices
152      INTEGER ::   inum            ! temporary logical unit
153      INTEGER ::   iim, ijm
154      INTEGER ::   ifreq, il1, il2, ij, ii
155      INTEGER, DIMENSION(jpidta,jpidta) ::   idata
156      INTEGER, DIMENSION(jpi   ,jpj   ) ::   icof
157
158      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
159
160      CHARACTER (len=15) ::   clexp
161      !!----------------------------------------------------------------------
162
163      IF(lwp) WRITE(numout,*)
164      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
165      IF(lwp) WRITE(numout,*) '~~~~~~  --'
166      IF(lwp) WRITE(numout,*)
167      IF(lwp) WRITE(numout,*) '        orca ocean model'
168      IF(lwp) WRITE(numout,*)
169
170#if defined key_antarctic
171#     include "ldfdyn_antarctic.h90"
172#elif defined key_arctic
173#     include "ldfdyn_arctic.h90"
174#else
175      ! Read 2d integer array to specify western boundary increase in the
176      ! ===================== equatorial strip (20N-20S) defined at t-points
177
178      CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
179      READ(inum,9101) clexp, iim, ijm
180      READ(inum,'(/)')
181      ifreq = 40
182      il1 = 1
183      DO jn = 1, jpidta/ifreq+1
184         READ(inum,'(/)')
185         il2 = MIN( jpidta, il1+ifreq-1 )
186         READ(inum,9201) ( ii, ji = il1, il2, 5 )
187         READ(inum,'(/)')
188         DO jj = jpjdta, 1, -1
189            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
190         END DO
191         il1 = il1 + ifreq
192      END DO
193     
194      DO jj = 1, nlcj
195         DO ji = 1, nlci
196            icof(ji,jj) = idata( mig(ji), mjg(jj) )
197         END DO
198      END DO
199      DO jj = nlcj+1, jpj
200         DO ji = 1, nlci
201            icof(ji,jj) = icof(ji,nlcj)
202         END DO
203      END DO
204      DO jj = 1, jpj
205         DO ji = nlci+1, jpi
206            icof(ji,jj) = icof(nlci,jj)
207         END DO
208      END DO
209
210 9101 FORMAT(1x,a15,2i8)
211 9201 FORMAT(3x,13(i3,12x))
212 9202 FORMAT(i3,41i3)
213
214
215      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
216      ! =================
217      ! define ahm1 and ahm2 at the right grid point position
218      ! (USER: modify ahm1 and ahm2 following your desiderata)
219     
220     
221      ! Decrease ahm to zahmeq m2/s in the tropics
222      ! (from 90 to 20 degre: ahm = constant
223      ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
224      ! from  2.5 to  0 degre: ahm = constant
225      ! symmetric in the south hemisphere)
226
227      zahmeq = aht0
228     
229      DO jj = 1, jpj
230         DO ji = 1, jpi
231            IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
232               ahm2(ji,jj) =  ahm0
233            ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
234               ahm2(ji,jj) =  zahmeq
235            ELSE
236               ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
237                  * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
238            ENDIF
239            IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
240               ahm1(ji,jj) =  ahm0
241            ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
242               ahm1(ji,jj) =  zahmeq
243            ELSE
244               ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
245                  * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
246            ENDIF
247         END DO
248      END DO
249
250      ! increase along western boundaries of equatorial strip
251      ! t-point
252      DO jj = 1, jpjm1
253         DO ji = 1, jpim1
254            zcoft = FLOAT( icof(ji,jj) ) / 100.
255            ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
256         END DO
257      END DO
258      ! f-point
259      icof(:,:) = icof(:,:) * tmask(:,:,1)
260      DO jj = 1, jpjm1
261         DO ji = 1, jpim1   ! NO vector opt.
262            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
263            IF( zmsk == 0. ) THEN
264               zcoff = 1.
265            ELSE
266               zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
267                     / (zmsk * 100.)
268            ENDIF
269            ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
270         END DO
271      END DO
272#endif
273     
274      ! Lateral boundary conditions on ( ahm1, ahm2 )
275      !                                ==============
276      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
277      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
278
279      ! Control print
280      IF( lwp .AND. ld_print ) THEN
281         WRITE(numout,*)
282         WRITE(numout,*) 'inildf: 2D ahm1 array'
283         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
284         WRITE(numout,*)
285         WRITE(numout,*) 'inildf: 2D ahm2 array'
286         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
287      ENDIF
288
289   END SUBROUTINE ldf_dyn_c2d_orca
Note: See TracBrowser for help on using the repository browser.