source: trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90

Last change on this file was 82, checked in by smasson, 10 years ago

good viscosity and diffusivity in trop075
File size: 10.6 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   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
10   !! $Id: ldfdyn_c2d.h90 3294 2012-01-28 16:44:18Z rblod $
11   !! Software governed by the CeCILL licence (NEMOGCM/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      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
35      !
36      INTEGER  ::   ji, jj
37      REAL(wp) ::   za00, zd_max, zetmax, zeumax, zefmax, zevmax
38      !!----------------------------------------------------------------------
39
40      IF(lwp) WRITE(numout,*)
41      IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient'
42      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
43
44      ! harmonic operator (ahm1, ahm2) : ( T- and F- points) (used for laplacian operators
45      ! ===============================                       whatever its orientation is)
46      IF( ln_dynldf_lap ) THEN
47         ! define ahm1 and ahm2 at the right grid point position
48         ! (USER: modify ahm1 and ahm2 following your desiderata)
49
50         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
51         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
52
53         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
54         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
55
56         za00 = ahm0 / zd_max
57         DO jj = 1, jpj
58            DO ji = 1, jpi
59               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
60               zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
61               ahm1(ji,jj) = za00 * zetmax
62               ahm2(ji,jj) = za00 * zefmax
63            END DO
64         END DO
65
66         IF( ln_dynldf_iso ) THEN
67            IF(lwp) WRITE(numout,*) '              Caution, as implemented now, the isopycnal part of momentum'
68            IF(lwp) WRITE(numout,*) '                 mixing use aht0 as eddy viscosity coefficient. Thus, it is'
69            IF(lwp) WRITE(numout,*) '                 uniform and you must be sure that your ahm is greater than'
70            IF(lwp) WRITE(numout,*) '                 aht0 everywhere in the model domain.'
71         ENDIF
72
73         ! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
74         ! ==============================================
75         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print, ld_cte = .true.  )
76         IF( cp_cfg == "orca" .AND.   jp_cfg == 1                    )   CALL ldf_dyn_c2d_orca( ld_print, ld_cte = .false. )
77         IF( cp_cfg == "trop" .AND.   jp_cfg == 075                  )   CALL ldf_dyn_c2d_orca( ld_print, ld_cte = .true.  )
78
79         ! Control print
80         IF( lwp .AND. ld_print ) THEN
81            WRITE(numout,*)
82            WRITE(numout,*) 'inildf: 2D ahm1 array'
83            CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
84            WRITE(numout,*)
85            WRITE(numout,*) 'inildf: 2D ahm2 array'
86            CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
87         ENDIF
88      ENDIF
89
90      ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator
91      ! =================================                      whatever its orientation is)
92      IF( ln_dynldf_bilap ) THEN
93         ! (USER: modify ahm3 and ahm4 following your desiderata)
94         ! Here: ahm is proportional to the cube of the maximum of the gridspacing
95         !       in the to horizontal direction
96
97         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
98         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
99
100         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
101         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
102
103         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
104         DO jj = 1, jpj
105            DO ji = 1, jpi
106               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
107               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
108               ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax
109               ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax
110            END DO
111         END DO
112
113         ! Control print
114         IF( lwp .AND. ld_print ) THEN
115            WRITE(numout,*)
116            WRITE(numout,*) 'inildf: ahm3 array'
117            CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
118            WRITE(numout,*)
119            WRITE(numout,*) 'inildf: ahm4 array'
120            CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
121         ENDIF
122      ENDIF
123      !
124   END SUBROUTINE ldf_dyn_c2d
125
126
127   SUBROUTINE ldf_dyn_c2d_orca( ld_print, ld_cte )
128      !!----------------------------------------------------------------------
129      !!                 ***  ROUTINE ldf_dyn_c2d  ***
130      !!
131      !!                   **** W A R N I N G ****
132      !!
133      !!                ORCA R2 and R4 configurations
134      !!                 
135      !!                   **** W A R N I N G ****
136      !!                 
137      !! ** Purpose :   initializations of the lateral viscosity for orca R2
138      !!
139      !! ** Method  :   blah blah blah...
140      !!
141      !!----------------------------------------------------------------------
142      USE ldftra_oce, ONLY:   aht0
143      !
144      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
145      LOGICAL, INTENT (in) ::   ld_cte     !
146      !
147      INTEGER  ::   ji, jj   ! dummy loop indices
148      INTEGER  ::   inum            ! local integers
149      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk,zam20s
150      REAL(wp), POINTER, DIMENSION(:,:)  :: zcof
151      !!----------------------------------------------------------------------
152      !                               
153      CALL wrk_alloc( jpi, jpj, zcof )
154      !
155      IF(lwp) WRITE(numout,*)
156      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
157      IF(lwp) WRITE(numout,*) '~~~~~~  --'
158      IF(lwp) WRITE(numout,*) '        orca ocean configuration'
159
160#if defined key_antarctic
161#     include "ldfdyn_antarctic.h90"
162#elif defined key_arctic
163#     include "ldfdyn_arctic.h90"
164#else
165      ! Read 2d integer array to specify western boundary increase in the
166      ! ===================== equatorial strip (20N-20S) defined at t-points
167     
168      CALL iom_open ( 'ahmcoef.nc', inum ) 
169      CALL iom_get  ( inum, jpdom_data, 'ahmcoef', zcof )
170      CALL iom_close( inum )
171
172      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
173      ! =================
174      ! define ahm1 and ahm2 at the right grid point position
175      ! (USER: modify ahm1 and ahm2 following your desiderata)
176     
177     
178      ! Decrease ahm to zahmeq m2/s in the tropics
179      ! (from 90 to 20 degre: ahm = constant
180      ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
181      ! from  2.5 to  0 degre: ahm = constant
182      ! symmetric in the south hemisphere)
183
184      zahmeq = aht0
185      IF (ld_cte) THEN   ;   zam20s = ahm0
186      ELSE               ;   zam20s = ahm0 * COS( rad * 20. )
187      END IF
188   
189      DO jj = 1, jpj
190         DO ji = 1, jpi
191            IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
192               IF (ld_cte) ahm2(ji,jj) = ahm0
193            ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
194               ahm2(ji,jj) = zahmeq
195            ELSE
196               ahm2(ji,jj) = zahmeq + (zam20s-zahmeq)/2.   &
197                  * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
198            ENDIF
199            IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
200               IF (ld_cte) ahm1(ji,jj) = ahm0
201            ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
202               ahm1(ji,jj) = zahmeq
203            ELSE
204               ahm1(ji,jj) = zahmeq + (zam20s-zahmeq)/2.   &
205                  * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
206            ENDIF
207         END DO
208      END DO
209
210      ! increase along western boundaries of equatorial strip
211      ! t-point
212      DO jj = 1, jpjm1
213         DO ji = 1, jpim1
214            zcoft = zcof(ji,jj) / 100.
215            ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
216         END DO
217      END DO
218      ! f-point
219      zcof(:,:) = zcof(:,:) * tmask(:,:,1)
220      DO jj = 1, jpjm1
221         DO ji = 1, jpim1   ! NO vector opt.
222            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
223            IF( zmsk == 0. ) THEN
224               zcoff = 1.
225            ELSE
226               zcoff = ( zcof(ji,jj+1) + zcof(ji+1,jj+1) + zcof(ji,jj) + zcof(ji,jj+1) )   &
227                     / (zmsk * 100.)
228            ENDIF
229            ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
230         END DO
231      END DO
232#endif
233     
234      ! Lateral boundary conditions on ( ahm1, ahm2 )
235      !                                ==============
236      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
237      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
238
239      ! Control print
240      IF( lwp .AND. ld_print ) THEN
241         WRITE(numout,*)
242         WRITE(numout,*) 'inildf: 2D ahm1 array'
243         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
244         WRITE(numout,*)
245         WRITE(numout,*) 'inildf: 2D ahm2 array'
246         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
247      ENDIF
248      !
249      CALL wrk_dealloc( jpi, jpj, zcof )
250      !
251   END SUBROUTINE ldf_dyn_c2d_orca
252
Note: See TracBrowser for help on using the repository browser.