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.
crsfld.F90 in branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90 @ 6101

Last change on this file since 6101 was 6101, checked in by cbricaud, 8 years ago

correction of bugs from last update and improvments for CRS
  • Property svn:keywords set to Id
File size: 24.4 KB
Line 
1MODULE crsfld
2   !!======================================================================
3   !!                     ***  MODULE  crsdfld  ***
4   !!  Ocean coarsening :  coarse ocean fields
5   !!=====================================================================
6   !!   2012-07  (J. Simeon, C. Calone, G. Madec, C. Ethe)
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   crs_fld       : create the standard output files for coarse grid and prep
11   !!                       other variables needed to be passed to TOP
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and tracers
14   USE dom_oce         ! ocean space and time domain
15   USE ldftra_oce      ! ocean active tracers: lateral physics
16   USE sbc_oce         ! Surface boundary condition: ocean fields
17   USE zdf_oce         ! vertical  physics: ocean fields
18   USE zdfddm          ! vertical  physics: double diffusion
19   USe zdfmxl
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE in_out_manager  ! I/O manager
22   USE timing          ! preformance summary
23   USE wrk_nemo        ! working array
24   USE crs
25   USE crsdom
26   USE crslbclnk
27   USE iom
28   USE zdfmxl_crs
29   USE eosbn2
30   USE zdfevd_crs
31   USE zdftke
32   USE zdftke_crs
33
34!   USE ieee_arithmetic
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   crs_fld                 ! routines called by step.F90
40
41
42   !! * Substitutions
43#  include "zdfddm_substitute.h90"
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id $
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE crs_fld( kt )
54      !!---------------------------------------------------------------------
55      !!                  ***  ROUTINE crs_fld  ***
56      !!                   
57      !! ** Purpose :   Basic output of coarsened dynamics and tracer fields
58      !!      NETCDF format is used by default
59      !!      1. Accumulate in time the dimensionally-weighted fields
60      !!      2. At time of output, rescale [1] by dimension and time
61      !!         to yield the spatial and temporal average.
62      !!  See. diawri_dimg.h90, sbcmod.F90
63      !!
64      !! ** Method  : 
65      !!----------------------------------------------------------------------
66      !!
67     
68      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
69      !!
70      INTEGER               ::   ji, jj, jk              ! dummy loop indices
71      !!
72      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3t, zfse3u, zfse3v, zfse3w ! 3D workspace for e3
73      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zs , ztmp
74      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d,z2d_crs
75      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_crs, zs_crs, zerr_crs,zmax_crs
76      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmp_crs
77      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: avte_crs
78      REAL(wp)       :: z2dcrsu, z2dcrsv
79      REAL(wp)       :: zmin,zmax,icnt1,icnt2
80      INTEGER :: i,j,ijis,ijie,ijjs,ijje
81      REAL(wp)       :: zw,zwp1,zum1,zu,zvm1,zv,zsnm,zsm,z
82      REAL(wp)       :: zerr, zerr0, zerr1, zmean
83      INTEGER,DIMENSION(4,3) :: ind
84      REAL(wp),DIMENSION(4) :: zwgt
85      INTEGER ::  iji,ijj
86      INTEGER :: jl,jm,jn
87      !!
88      !!----------------------------------------------------------------------
89
90      IF( nn_timing == 1 )   CALL timing_start('crs_fld')
91
92      !  Initialize arrays
93      CALL wrk_alloc( jpi, jpj, jpk, zfse3t, zfse3w )
94      CALL wrk_alloc( jpi, jpj, jpk, zfse3u, zfse3v )
95      CALL wrk_alloc( jpi, jpj, jpk, zt, zs , ztmp        )
96      CALL wrk_alloc( jpi, jpj,      z2d            )
97      !
98      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs, zerr_crs,zmax_crs )
99      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, ztmp_crs )
100      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, 4, avte_crs )
101      CALL wrk_alloc( jpi_crs, jpj_crs, z2d_crs     )
102
103      ! Depth work arrrays
104      zfse3t(:,:,:) = fse3t(:,:,:)
105      zfse3u(:,:,:) = fse3u(:,:,:)
106      zfse3v(:,:,:) = fse3v(:,:,:)
107      zfse3w(:,:,:) = fse3w(:,:,:)
108
109      IF( kt == nit000  ) THEN
110         tsn_crs  (:,:,:,:) = 0._wp    ! temp/sal  array, now
111         un_crs   (:,:,:  ) = 0._wp    ! u-velocity
112         vn_crs   (:,:,:  ) = 0._wp    ! v-velocity
113         wn_crs   (:,:,:  ) = 0._wp    ! w
114         avt_crs  (:,:,:  ) = 0._wp    ! avt
115         hdivb_crs(:,:,:  ) = 0._wp    ! hdiv
116         hdivn_crs(:,:,:  ) = 0._wp    ! hdiv
117         rke_crs  (:,:,:  ) = 0._wp    ! rke
118         sshn_crs (:,:    ) = 0._wp    ! ssh
119         utau_crs (:,:    ) = 0._wp    ! taux
120         vtau_crs (:,:    ) = 0._wp    ! tauy
121         wndm_crs (:,:    ) = 0._wp    ! wind speed
122         qsr_crs  (:,:    ) = 0._wp    ! qsr
123         emp_crs  (:,:    ) = 0._wp    ! emp
124         emp_b_crs(:,:    ) = 0._wp    ! emp
125         rnf_crs  (:,:    ) = 0._wp    ! runoff
126         fr_i_crs (:,:    ) = 0._wp    ! ice cover
127      ENDIF
128
129      CALL iom_swap( "nemo_crs" )    ! swap on the coarse grid
130
131      ! 2. Coarsen fields at each time step
132      ! --------------------------------------------------------
133
134      !  Temperature
135      zt(:,:,:) = tsb(:,:,:,jp_tem)  ;      zt_crs(:,:,:) = 0._wp
136      CALL crs_dom_ope( zt, 'VOL', 'T', tmask, zt_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
137      tsb_crs(:,:,:,jp_tem) = zt_crs(:,:,:)
138      zt(:,:,:) = tsn(:,:,:,jp_tem)  ;      zt_crs(:,:,:) = 0._wp
139      CALL crs_dom_ope( zt, 'VOL', 'T', tmask, zt_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
140      tsn_crs(:,:,:,jp_tem) = zt_crs(:,:,:)
141
142      CALL iom_put( "toce", tsn_crs(:,:,:,jp_tem) )    ! temp
143      CALL iom_put( "sst" , tsn_crs(:,:,1,jp_tem) )    ! sst
144
145      !n2 before
146      zt(:,:,:) = rn2b(:,:,:)  ;      zt_crs(:,:,:) = 0._wp
147      CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
148      rb2_crs(:,:,:) = zt_crs(:,:,:)
149      CALL iom_put("rb2_crs",rb2_crs)
150
151      !  Salinity
152      zs(:,:,:) = tsb(:,:,:,jp_sal)  ;      zs_crs(:,:,:) = 0._wp
153      CALL crs_dom_ope( zs, 'VOL', 'T', tmask, zs_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
154      tsb_crs(:,:,:,jp_sal) = zs_crs(:,:,:)
155      zs(:,:,:) = tsn(:,:,:,jp_sal)  ;      zs_crs(:,:,:) = 0._wp
156      CALL crs_dom_ope( zs, 'VOL', 'T', tmask, zs_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
157      tsn_crs(:,:,:,jp_sal) = zs_crs(:,:,:)
158
159      CALL iom_put( "soce" , tsn_crs(:,:,:,jp_sal) )    ! sal
160      CALL iom_put( "sss"  , tsn_crs(:,:,1,jp_sal) )    ! sss
161
162      !  U-velocity
163      CALL crs_dom_ope( ub, 'SUM', 'U', umask, ub_crs, p_e12=e2u, p_e3=zfse3u, p_surf_crs=e2e3u_msk, psgn=-1.0 )
164      CALL crs_dom_ope( un, 'SUM', 'U', umask, un_crs, p_e12=e2u, p_e3=zfse3u, p_surf_crs=e2e3u_msk, psgn=-1.0 )
165      !cbr
166      ub_crs(:,:,:) = ub_crs(:,:,:)*umask_crs(:,:,:)
167      un_crs(:,:,:) = un_crs(:,:,:)*umask_crs(:,:,:)
168      !
169      zt(:,:,:) = 0._wp     ;    zs(:,:,:) = 0._wp  ;   zt_crs(:,:,:) = 0._wp   ;    zs_crs(:,:,:) = 0._wp
170      DO jk = 1, jpkm1
171         DO jj = 2, jpjm1
172            DO ji = 2, jpim1   
173               zt(ji,jj,jk)  = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
174               zs(ji,jj,jk)  = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 
175            END DO
176         END DO
177      END DO
178      CALL crs_dom_ope( zt, 'SUM', 'U', umask, zt_crs, p_e12=e2u, p_e3=zfse3u, p_surf_crs=e2e3u_msk, psgn=-1.0 )
179      CALL crs_dom_ope( zs, 'SUM', 'U', umask, zs_crs, p_e12=e2u, p_e3=zfse3u, p_surf_crs=e2e3u_msk, psgn=-1.0 )
180
181      CALL iom_put( "uoce"  , un_crs )   ! i-current
182      CALL iom_put( "uocet" , zt_crs )   ! uT
183      CALL iom_put( "uoces" , zs_crs )   ! uS
184
185      !  V-velocity
186      CALL crs_dom_ope( vb, 'SUM', 'V', vmask, vb_crs, p_e12=e1v, p_e3=zfse3v, p_surf_crs=e1e3v_msk, psgn=-1.0 )
187      CALL crs_dom_ope( vn, 'SUM', 'V', vmask, vn_crs, p_e12=e1v, p_e3=zfse3v, p_surf_crs=e1e3v_msk, psgn=-1.0 )
188      vb_crs(:,:,:) = vb_crs(:,:,:)*vmask_crs(:,:,:)
189      vn_crs(:,:,:) = vn_crs(:,:,:)*vmask_crs(:,:,:)
190      !                                                                                 
191      zt(:,:,:) = 0._wp     ;    zs(:,:,:) = 0._wp  ;   zt_crs(:,:,:) = 0._wp   ;    zs_crs(:,:,:) = 0._wp
192      DO jk = 1, jpkm1
193         DO jj = 2, jpjm1
194            DO ji = 2, jpim1   
195               zt(ji,jj,jk)  = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
196               zs(ji,jj,jk)  = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) ) 
197            END DO
198         END DO
199      END DO
200      CALL crs_dom_ope( zt, 'SUM', 'V', vmask, zt_crs, p_e12=e1v, p_e3=zfse3v, p_surf_crs=e1e3v_msk, psgn=-1.0 )
201      CALL crs_dom_ope( zs, 'SUM', 'V', vmask, zs_crs, p_e12=e1v, p_e3=zfse3v, p_surf_crs=e1e3v_msk, psgn=-1.0 )
202 
203      CALL iom_put( "voce"  , vn_crs )   ! i-current
204      CALL iom_put( "vocet" , zt_crs )   ! vT
205      CALL iom_put( "voces" , zs_crs )   ! vS
206
207     
208      !  Kinetic energy
209      CALL crs_dom_ope( rke, 'VOL', 'T', tmask, rke_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
210      CALL iom_put( "eken", rke_crs )
211
212      !  Horizontal divergence ( following OPA_SRC/DYN/divcur.F90 )
213      DO jk = 1, jpkm1
214         DO ji = 2, jpi_crsm1
215            DO jj = 2, jpj_crsm1
216               IF( tmask_crs(ji,jj,jk ) > 0 ) THEN
217                  !z2dcrsu =  ( un_crs(ji  ,jj  ,jk) * crs_surfu_wgt(ji  ,jj  ,jk) ) &
218                  !   &     - ( un_crs(ji-1,jj  ,jk) * crs_surfu_wgt(ji-1,jj  ,jk) )
219                  !z2dcrsv =  ( vn_crs(ji  ,jj  ,jk) * crs_surfv_wgt(ji  ,jj  ,jk) ) &
220                  !   &     - ( vn_crs(ji  ,jj-1,jk) * crs_surfv_wgt(ji  ,jj-1,jk) )
221                  !
222                  !IF( crs_volt_wgt(ji,jj,jk) .NE. 0._wp ) hdivn_crs(ji,jj,jk) = ( z2dcrsu + z2dcrsv ) / crs_volt_wgt(ji,jj,jk)
223                  z2dcrsu =  ( un_crs(ji  ,jj  ,jk) * e2e3u_msk(ji  ,jj  ,jk) ) &
224                     &     - ( un_crs(ji-1,jj  ,jk) * e2e3u_msk(ji-1,jj  ,jk) )
225                  z2dcrsv =  ( vn_crs(ji  ,jj  ,jk) * e1e3v_msk(ji  ,jj  ,jk) ) &
226                     &     - ( vn_crs(ji  ,jj-1,jk) * e1e3v_msk(ji  ,jj-1,jk) )
227                  !
228                  !cbr
229                  !bug1: il manquait le facvol_t(ji,jj,jk) ds la division ; ca creait des grosses erreurs de Wcrs ( vu en recalculant la divergence 3D )
230                  !bug2: mm test que bug1: on n'obtient tjs pas zero
231                  !on a la div calculée via ocean_volume_crs_t puis w via  e3t_crs ; or ,e1t_crs(ji,jj)*e2t_crs(ji,jj)*e3t_crs(ji,jj,jk) NE ocean_volume_crs_t*crs_volt_wgt(ji,jj,jk)
232                  !exp (117,211,74) : e1*e2*e3=235206030060.005 / ocean_volume_crs_t * facvol = 235205585307.810
233                  !                   e1*e2*e3-cean_volume_crs_t * facvol/(cean_volume_crs_t * facvol) ~1.e-6) 
234                  IF( ocean_volume_crs_t(ji,jj,jk) .NE. 0._wp ) hdivn_crs(ji,jj,jk) = ( z2dcrsu + z2dcrsv )
235
236                  z2dcrsu =  ( ub_crs(ji  ,jj  ,jk) * e2e3u_msk(ji  ,jj  ,jk) ) &
237                     &     - ( ub_crs(ji-1,jj  ,jk) * e2e3u_msk(ji-1,jj  ,jk) )
238                  z2dcrsv =  ( vb_crs(ji  ,jj  ,jk) * e1e3v_msk(ji  ,jj  ,jk) ) &
239                     &     - ( vb_crs(ji  ,jj-1,jk) * e1e3v_msk(ji  ,jj-1,jk) )
240                  !
241                  IF( ocean_volume_crs_t(ji,jj,jk) .NE. 0._wp ) hdivb_crs(ji,jj,jk) = ( z2dcrsu + z2dcrsv ) / (facvol_t(ji,jj,jk)*ocean_volume_crs_t(ji,jj,jk) )
242               ENDIF
243            ENDDO
244         ENDDO
245      ENDDO
246      CALL crs_lbc_lnk( hdivn_crs, 'T', 1.0 )
247      !
248      CALL iom_put( "hdiv", hdivn_crs ) 
249
250
251      !  W-velocity
252      IF( ln_crs_wn ) THEN
253         CALL crs_dom_ope( wn, 'SUM', 'W', tmask, wn_crs, p_e12=e1e2t, p_surf_crs=e1e2w_msk, psgn=1.0 )
254      ELSE
255        wn_crs(:,:,jpk) = 0._wp
256        DO jk = jpkm1, 1, -1
257           wn_crs(:,:,jk) = e1e2w_msk(:,:,jk+1)*wn_crs(:,:,jk+1) - hdivn_crs(:,:,jk)
258           WHERE( e1e2w_msk(:,:,jk) .NE. 0._wp )  wn_crs(:,:,jk) =  wn_crs(:,:,jk) /e1e2w_msk(:,:,jk) 
259        ENDDO
260      ENDIF
261
262      CALL iom_put( "woce", wn_crs  )   ! vertical velocity
263      !  free memory
264
265      !  avt, avs
266      SELECT CASE ( nn_crs_kz )
267         CASE ( 0 )
268            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
269         CASE ( 1 )
270            CALL crs_dom_ope( avt, 'MAX', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
271         CASE ( 2 )
272            CALL crs_dom_ope( avt, 'MIN', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
273         CASE ( 3 )
274            CALL crs_dom_ope( avt, 'LOGVOL', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, p_mask_crs=tmask_crs, psgn=1.0 )
275         CASE ( 4 )
276            CALL crs_dom_ope( avt, 'MED', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
277         CASE ( 5 )
278            CALL crs_dom_ope( en , 'VOL', 'W', tmask, en_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
279            CALL crs_dom_ope( taum , 'SUM', 'T', tmask, taum_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
280            CALL crs_dom_ope( rn2(:,:,:), 'VOL', 'W', tmask, rn2_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
281            IF( kt==nit000 )CALL tke_avn_ini_crs
282            CALL tke_avn_crs
283            CALL zdf_evd_crs(kt)
284         CASE ( 6 )
285
286            avte_crs(:,:,:,:) = 0._wp
287            ztmp(:,:,:)=1.
288
289            zt(:,:,:) = 0._wp
290            zs(:,:,:) = 0._wp
291            DO jk=2,jpk 
292               WHERE( fse3w(:,:,jk) .NE. 0._wp) zs(:,:,jk) = ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) / fse3w(:,:,jk)
293               zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
294            ENDDO
295            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
296            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
297            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
298            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
299            zt_crs=tmask_crs*zt_crs
300            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
301            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,1) = zt_crs / zs_crs
302            zmin=MINVAL(avte_crs(:,:,:,1));zmax=MAXVAL(avte_crs(:,:,:,1));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
303
304            zt(:,:,:) = 0._wp
305            zs(:,:,:) = 0._wp
306            DO jk=2,jpk 
307               WHERE( fse3w(:,:,jk) .NE. 0._wp) zs(:,:,jk) = ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) / fse3w(:,:,jk)
308               zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
309            ENDDO
310            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
311            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
312            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax 
313            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax 
314            zt_crs=tmask_crs*zt_crs
315            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax
316            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,2) = zt_crs / zs_crs
317            zmin=MINVAL(avte_crs(:,:,:,2));zmax=MAXVAL(avte_crs(:,:,:,2));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax
318
319            zt(:,:,:) = 0._wp
320            zs(:,:,:) = 0._wp
321            DO jk=2,jpk
322                WHERE( fse3w(:,:,jk) .NE. 0._wp ) zs(:,:,jk)= ( rn_a0 * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) +  &
323                                                                  &   rn_b0 * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) )/ fse3w(:,:,jk) 
324                zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)                                                           
325            ENDDO
326            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
327            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
328            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax 
329            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax 
330            zt_crs=tmask_crs*zt_crs
331            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax
332            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,3) = zt_crs / zs_crs
333            zmin=MINVAL(avte_crs(:,:,:,3));zmax=MAXVAL(avte_crs(:,:,:,3));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax
334
335            zt(:,:,:) = 0._wp
336            zs(:,:,:) = 0._wp
337            DO jk=2,jpk
338                WHERE( fse3w(:,:,jk) .NE. 0._wp ) zs(:,:,jk)= ( rn_a0 * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) -  &
339                                                                  &   rn_b0 * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) )/ fse3w(:,:,jk) 
340                zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
341            ENDDO
342            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
343            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
344            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax 
345            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax 
346            zt_crs=tmask_crs*zt_crs
347            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax
348            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,4) = zt_crs / zs_crs
349            zmin=MINVAL(avte_crs(:,:,:,4));zmax=MAXVAL(avte_crs(:,:,:,4));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax
350
351            CALL iom_put( "avte_crs1", avte_crs(:,:,:,1) )   !  Kz
352            CALL iom_put( "avte_crs2", avte_crs(:,:,:,2) )   !  Kz
353            CALL iom_put( "avte_crs3", avte_crs(:,:,:,3) )   !  Kz
354            CALL iom_put( "avte_crs4", avte_crs(:,:,:,4) )   !  Kz
355!---------------------
356            CALL crs_dom_ope( avt, 'MED', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=zfse3w, p_mask_crs=tmask_crs, psgn=1.0 )
357!?            zmin=MINVAL(zs_crs*tmask_crs);zmax=MAXVAL(zs_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"logvol zs_crs*tmask ",zmin,zmax ; call flush(numout)
358            CALL iom_put( "zs_crs", zs_crs )   !  Kzlogvol
359!--------------------- ok
360
361            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, zmax_crs, p_e12=e1e2t, p_e3=zfse3w,  psgn=1.0 )
362            WRITE(narea+200,*)"zmax_crs ",SHAPE(zmax_crs) ; call flush(narea+200)
363            CALL iom_put( "zmax_crs", zmax_crs )   !  Kzlogvol
364            zmin=MINVAL(zmax_crs*tmask_crs);zmax=MAXVAL(zmax_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"vol zmax_crs*tmask ",zmin,zmax ; call flush(numout)
365!-------------------------nok
366            avt_crs=zs_crs
367
368
369            zerr0=0.01
370
371            icnt1=0
372            icnt2=0
373
374            zt_crs(:,:,:)=0._wp
375            zerr_crs(:,:,:)=0._wp
376            DO ji=1,jpi_crs 
377            DO jj=1,jpj_crs 
378            DO jk=1,jpk
379
380 
381!--------------
382               zwgt(1:4)=0._wp
383               DO jm=1,4 ; IF( avte_crs(ji,jj,jk,jm)  .GE. 0._wp .AND.  avte_crs(ji,jj,jk,jm)  .LE. zmax_crs(ji,jj,jk)  ) zwgt(jm) = 1._wp ; ENDDO
384!--------------
385               IF( SUM(zwgt(1:4)) .NE. 0._wp ) THEN 
386                  zmean = SUM( zwgt(1:4)*avte_crs(ji,jj,jk,1:4)) / SUM(zwgt(1:4) )
387                  zerr  = SQRT(SUM( zwgt(1:4)*(avte_crs(ji,jj,jk,1:4)-zmean)*(avte_crs(ji,jj,jk,1:4)-zmean) ) / SUM(zwgt(1:4) ) )
388               ELSE
389                  zmean=0._wp
390                 zerr=1.e10
391               ENDIF
392!--------------
393
394               zerr_crs(ji,jj,jk)=zerr
395
396               IF( zerr .LE. zerr0 .AND. zmean .GT. 0._wp )zt_crs(ji,jj,jk)=zmean
397               IF( zerr .LE. zerr0 .AND. zmean .GT. 0._wp )avt_crs(ji,jj,jk)=zmean
398
399               IF( tmask_crs(ji,jj,jk) == 1 ) icnt1=icnt1+1
400               IF( tmask_crs(ji,jj,jk) == 1 .AND.  zerr .LE. zerr0 .AND. zmean .GT. 0._wp ) icnt2=icnt2+1
401
402               !IF( ieee_is_nan(  zt_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANMEANEFF ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
403               !IF( ieee_is_nan(  zs_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANLOG ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
404               !IF( ieee_is_nan( avt_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANAVT ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
405            ENDDO
406            ENDDO
407            ENDDO
408            zmin=MINVAL(avt_crs);zmax=MAXVAL(avt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avt_crs ",zmin,zmax  ; call flush(numout)
409            zmin=MINVAL(avt_crs*tmask_crs);zmax=MAXVAL(avt_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avt_crs*tmask ",zmin,zmax  ; call flush(numout)
410
411            CALL mpp_sum(icnt1)
412            CALL mpp_sum(icnt2)
413            IF(lwp)WRITE(numout,*)"TOTO",kt,icnt1,icnt2
414            CALL iom_put( "zt_crs", zt_crs )   !  Kz
415            CALL iom_put( "zerr_crs", zerr_crs )   !  Kz
416
417      END SELECT
418      !
419      CALL iom_put( "avt", avt_crs )   !  Kz
420     
421      !deja dasn step CALL zdf_mxl_crs(kt)
422
423 
424      !  sbc fields 
425
426      CALL crs_dom_ope( sshb , 'VOL', 'T', tmask, sshb_crs , p_e12=e1e2t, p_e3=zfse3t         , psgn=1.0 ) 
427      CALL crs_dom_ope( sshn , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=zfse3t         , psgn=1.0 ) 
428      CALL crs_dom_ope( ssha , 'VOL', 'T', tmask, ssha_crs , p_e12=e1e2t, p_e3=zfse3t         , psgn=1.0 ) 
429      CALL crs_dom_ope( utau , 'SUM', 'U', umask, utau_crs , p_e12=e2u  , p_surf_crs=e2u_crs  , psgn=1.0 )
430      CALL crs_dom_ope( vtau , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v  , p_surf_crs=e1v_crs  , psgn=1.0 )
431      CALL crs_dom_ope( wndm , 'SUM', 'T', tmask, wndm_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
432      CALL crs_dom_ope( rnf  , 'MAX', 'T', tmask, rnf_crs                                     , psgn=1.0 )
433      CALL crs_dom_ope( qsr  , 'SUM', 'T', tmask, qsr_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
434      CALL crs_dom_ope( emp_b, 'SUM', 'T', tmask, emp_b_crs, p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
435      CALL crs_dom_ope( emp  , 'SUM', 'T', tmask, emp_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
436      CALL crs_dom_ope( fmmflx,'SUM', 'T', tmask, fmmflx_crs, p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
437      CALL crs_dom_ope( sfx  , 'SUM', 'T', tmask, sfx_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
438      CALL crs_dom_ope( fr_i , 'SUM', 'T', tmask, fr_i_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
439
440      z2d=REAL(nmln,wp)
441      CALL crs_dom_ope( z2d , 'MAX', 'T', tmask, z2d_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
442      nmln_crs=INT(z2d_crs) 
443      nmln_crs=MAX(nlb10,nmln_crs)   
444
445      CALL iom_put( "ssh"      , sshn_crs )   ! ssh output
446      CALL iom_put( "utau"     , utau_crs )   ! i-tau output
447      CALL iom_put( "vtau"     , vtau_crs )   ! j-tau output
448      CALL iom_put( "wspd"     , wndm_crs )   ! wind speed output
449      CALL iom_put( "runoffs"  , rnf_crs  )   ! runoff output
450      CALL iom_put( "qsr"      , qsr_crs  )   ! qsr output
451      CALL iom_put( "empmr"    , emp_crs  )   ! water flux output
452      CALL iom_put( "saltflx"  , sfx_crs  )   ! salt flux output
453      CALL iom_put( "ice_cover", fr_i_crs )   ! ice cover output
454
455      !  free memory
456      CALL wrk_dealloc( jpi, jpj, jpk, zfse3t, zfse3w )
457      CALL wrk_dealloc( jpi, jpj, jpk, zfse3u, zfse3v )
458      CALL wrk_dealloc( jpi, jpj, jpk, zt, zs, ztmp   )
459      CALL wrk_dealloc( jpi, jpj, z2d                 )
460      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs, zerr_crs,zmax_crs )
461      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, ztmp_crs )
462      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, 4, avte_crs )
463      CALL wrk_dealloc( jpi_crs, jpj_crs, z2d_crs     )
464      !
465      CALL iom_swap( "nemo" )     ! return back on high-resolution grid
466      !
467      IF( nn_timing == 1 )   CALL timing_stop('crs_fld')
468      !
469   END SUBROUTINE crs_fld
470
471   !!======================================================================
472END MODULE crsfld
Note: See TracBrowser for help on using the repository browser.