MODULE crsini !!====================================================================== !! *** MODULE crsini *** !! Manage the grid coarsening module initialization !!====================================================================== !! History 2012-05 (J. Simeon, G. Madec, C. Ethe) Original code !!---------------------------------------------------------------------- USE timing ! Timing USE par_oce ! For parameter jpi,jpj,jphgr_msh USE dom_oce ! For parameters in par_oce (jperio, lk_vvl) USE crs_dom ! Coarse grid domain USE phycst, ONLY: omega, rad ! physical constants ! USE wrk_nemo USE in_out_manager USE par_kind, ONLY: wp USE crs USE crsdomwri USE crslbclnk IMPLICIT NONE PRIVATE PUBLIC crs_init !! * Substitutions # include "domzgr_substitute.h90" CONTAINS SUBROUTINE crs_init !!------------------------------------------------------------------- !! *** SUBROUTINE crs_init !! ** Purpose : Initialization of the grid coarsening module !! 1. Read namelist !! X2. MOVE TO crs_dom.F90 Set the domain definitions for coarse grid !! a. Define the coarse grid starting/ending indices on parent grid !! Here is where the T-pivot or F-pivot grids are discerned !! b. TODO. Include option for north-centric or equator-centric binning. !! (centered only for odd reduction factors; even reduction bins bias equator to the south) !! 3. Mask and mesh creation. => calls to crsfun !! a. Use crsfun_mask to generate tmask,umask, vmask, fmask. !! b. Use crsfun_coordinates to get coordinates !! c. Use crsfun_UV to get horizontal scale factors !! d. Use crsfun_TW to get initial vertical scale factors !! 4. Volumes and weights jes.... TODO. Updates for vvl? Where to do this? crsstp.F90? !! a. Calculate initial coarse grid box volumes: pvol_T, pvol_W !! b. Calculate initial coarse grid surface-averaging weights !! c. Calculate initial coarse grid volume-averaging weights !! !! X5. MOVE TO crs_dom_wri.F90 Using iom_rstput output the initial meshmask. !! ?. Another set of "masks" to generate !! are the u- and v- surface areas for U- and V- area-weighted means. !! Need to put this somewhere in section 3? !! jes. What do to about the vvl? GM. could separate the weighting (denominator), so !! output C*dA or C*dV as summation not mran, then do mean (division) at moment of output. !! As is, crsfun takes into account vvl. !! Talked about pre-setting the surface array to avoid IF/ENDIFS and division. !! But have then to make that preset array here and elsewhere. !! that is called every timestep... !! !! - Read in pertinent data ? !!------------------------------------------------------------------- !! Local variables INTEGER :: ji,jj,jk,ijjgloT,ijis,ijie,ijjs,ijje ! dummy indices INTEGER :: ierr ! allocation error status REAL(wp) :: zrestx, zresty ! for determining odd or even reduction factor REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmbk REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfse3t, zfse3u, zfse3v, zfse3f REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfse3w, zfse3t_n, zfse3t_b LOGICAL :: llok NAMELIST/namcrs/ nn_factx, nn_facty, cn_binref, nn_fcrs, nn_msh_crs, cn_ocerstcrs !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('crs_init') ! IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'crs_init : Initializing the grid coarsening module ' ENDIF !--------------------------------------------------------- ! 1. Read Namelist file !--------------------------------------------------------- ! READ( numnam, namcrs ) ! Namelist namcrs : Grid-coarsening utility namelist IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'crs_init: Namelist namcrs ' WRITE(numout,*) ' nn_factx = ', nn_factx WRITE(numout,*) ' nn_facty = ', nn_facty WRITE(numout,*) ' cn_binref = ', cn_binref WRITE(numout,*) ' nn_fcrs = ', nn_fcrs WRITE(numout,*) ' nn_msh_crs = ', nn_msh_crs ENDIF rfactx_r = 1./nn_factx rfacty_r = 1./nn_facty !--------------------------------------------------------- ! 2. Define Global Dimensions of the coarsened grid !--------------------------------------------------------- ! 2.a. Define global domain indices jpiglo_crs = INT( (jpiglo - 2) / nn_factx ) + 2 jpjglo_crs = INT( (jpjglo - 2) / nn_facty ) + 2 ! the -2 removes j=1, j=jpj jpiglo_crsm1 = jpiglo_crs - 1 jpjglo_crsm1 = jpjglo_crs - 1 jpkm1 = jpk - 1 ! 2.b. Define local domain indices jpi_crs = ( jpiglo_crs-2*jpreci + (jpni-1) ) / jpni + 2*jpreci jpj_crs = ( jpjglo_crs-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj jpi_crsm1 = jpi_crs - 1 jpj_crsm1 = jpj_crs - 1 nperio_crs = jperio npolj_crs = npolj IF ( jpnij == 1 ) THEN jpnij_crs = jpnij narea_crs = narea nimpp_crs = nimpp njmpp_crs = njmpp ELSE WRITE(numout,*) 'crsini.F90. mpp not supported... Stopping' STOP ENDIF nlcj_crs = jpj_crs nlci_crs = jpi_crs nldi_crs = 1 nlei_crs = jpi_crs nlej_crs = jpj_crs nldj_crs = 1 IF (lwp) THEN WRITE(numout,*) WRITE(numout,*) 'crs_init : coarse grid dimensions' WRITE(numout,*) '~~~~~~~ coarse domain global j-dimension jpjglo_crs = ', jpjglo_crs WRITE(numout,*) '~~~~~~~ coarse domain global i-dimension jpiglo_crs = ', jpiglo_crs WRITE(numout,*) '~~~~~~~ coarse domain local i-dimension jpi_crs = ', jpi_crs WRITE(numout,*) '~~~~~~~ coarse domain local j-dimension jpj_crs = ', jpj_crs ENDIF mxbinctr = INT( nn_factx * 0.5 ) mybinctr = INT( nn_facty * 0.5 ) zrestx = MOD( nn_factx, 2 ) ! check if even- or odd- numbered reduction factor zresty = MOD( nn_facty, 2 ) IF ( zrestx == 0 ) THEN mxbinctr = mxbinctr - 1 ENDIF IF ( zresty == 0 ) THEN mybinctr = mybinctr - 1 IF ( jperio == 3 .OR. jperio == 4 ) nperio_crs = jperio + 2 IF ( jperio == 5 .OR. jperio == 6 ) nperio_crs = jperio - 2 IF ( npolj == 3 ) npolj_crs = 5 IF ( npolj == 5 ) npolj_crs = 3 ENDIF rfactxy = nn_factx * nn_facty !jes. TODO Need to deallocate these if ln_crs = F ierr = crs_dom_alloc() ! allocate most coarse grid arrays ! jes. TODO. Add the next two lines when mpp is done ! IF( lk_mpp ) CALL mpp_sum( ierr ) ! IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'nemo_alloc : unable to allocate standard ocean arrays' ) ! 2.c. Set up bins for coarse grid, horizontal only. mis_crs(:) = 0; mie_crs(:) = 0 mjs_crs(:) = 0; mje_crs(:) = 0 SELECT CASE ( cn_binref ) CASE ( 'NORTH' ) SELECT CASE ( nperio ) CASE ( 2 ) WRITE(numout,*) 'crs_init, jperio=2 not supported' CASE ( 3, 4 ) ! T-Pivot at North Fold DO ji = 2, jpiglo_crsm1 !cc ijie = (ji*nn_factx)-nn_factx+1 ijie = (ji*nn_factx)-nn_factx !cc ijis = ijie-nn_factx+1 IF ( ji == jpiglo_crsm1 ) THEN IF ( ((jpiglo-1)-ijie) <= nn_factx ) ijie = jpiglo-2 ! ijie = jpiglo-1 !cc ENDIF ! Handle first the northernmost bin IF ( nn_facty == 2 ) THEN ijjgloT=jpjglo-1 ELSE ijjgloT=jpjglo ENDIF DO jj = 2, jpjglo_crsm1 ! cc ijje = ijjgloT-nn_facty*(jj-2) ijje = ijjgloT-nn_facty*(jj-2) - 1 ijjs = ijje-nn_facty+1 IF ( ijjs <= nn_facty ) ijjs = 2 mis_crs(ji) = ijis mie_crs(ji) = ijie mjs_crs(jpjglo_crs-jj+1) = ijjs mje_crs(jpjglo_crs-jj+1) = ijje ENDDO ENDDO CASE ( 5, 6 ) ! F-pivot at North Fold DO ji = 2, jpiglo_crsm1 ijie = (ji*nn_factx)-nn_factx+1 ijis = ijie-nn_factx+1 IF ( ji == jpiglo_crsm1 ) THEN IF ( ((jpiglo-1)-ijie) <= nn_factx ) ijie = jpiglo-1 ENDIF ! Treat the northernmost bin separately. jj = 2 ijje = jpjglo-nn_facty*(jj-2) IF ( nn_facty == 3 ) THEN ijjs=ijje-1 ELSE ijjs=ijje-nn_facty+1 ENDIF mis_crs(ji) = ijis mie_crs(ji) = ijie mjs_crs(jpjglo_crs-jj+1) = ijjs mje_crs(jpjglo_crs-jj+1) = ijje ! Now bin the rest, any remainder at the south is lumped in the southern bin DO jj = 3, jpjglo_crsm1 ijje = jpjglo-nn_facty*(jj-2) ijjs = ijje-nn_facty+1 IF ( ijjs <= nn_facty ) ijjs = 2 mis_crs(ji) = ijis mie_crs(ji) = ijie mjs_crs(jpjglo_crs-jj+1) = ijjs mje_crs(jpjglo_crs-jj+1) = ijje ENDDO ENDDO CASE DEFAULT WRITE(numout,*) 'crs_init. Only jperio = 3, 4, 5, 6 supported' END SELECT CASE ( 'EQUAT' ) WRITE(numout,*) 'crs_init. Equator-centered bins option not yet available' END SELECT ! Pad the boundaries, do not know if it is necessary mis_crs(1) = 1 ; mis_crs(jpiglo_crs) = mie_crs(jpiglo_crs - 1) + 1 !cc mie_crs(1) = nn_factx ; mie_crs(jpiglo_crs) = jpiglo !cc mjs_crs(1) = 1 ; mjs_crs(jpjglo_crs) = mje_crs(jpjglo_crs - 1) + 1 mje_crs(1) = mjs_crs(2)-1; mje_crs(jpjglo_crs) = jpjglo ! WRITE(numout,*) 'crs_init. coarse grid bounds on parent grid' ! WRITE(numout,'(1x,a,62(1x,i3),/)') 'mis_crs=', mis_crs ! WRITE(numout,'(1x,a,62(1x,i3),/)') 'mie_crs=', mie_crs ! WRITE(numout,'(1x,a,51(1x,i3),/)') 'mjs_crs=', mjs_crs ! WRITE(numout,'(1x,a,51(1x,i3),/)') 'mje_crs=', mje_crs !--------------------------------------------------------- ! 3. Mask and Mesh !--------------------------------------------------------- ! Set up the masks and meshes ! 3.a. Get the masks CALL crsfun( p_ptmask=tmask, cd_type='T', p_cmask=tmask_crs ) CALL crsfun( p_ptmask=umask, cd_type='U', p_cmask=umask_crs ) CALL crsfun( p_ptmask=vmask, cd_type='V', p_cmask=vmask_crs ) CALL crsfun( p_ptmask=fmask, cd_type='F', p_cmask=fmask_crs ) ! CALL crsfun( p_ptmask=tmask, cd_type='T', p_pmask2d=rnfmsk, p_cmask2d=rnfmsk_crs ) ! rnfmsk_crs(:,:) = rnfmsk_crs(:,:) * tmask_crs(:,:,1) WRITE(numout,*) 'crsini. Finished masks' ! 3.b. Get the coordinates ! Odd-numbered reduction factor, center coordinate on T-cell ! Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner. ! IF ( zresty /= 0 .AND. zrestx /= 0 ) THEN CALL crsfun( gphit, glamt, 'T', gphit_crs, glamt_crs ) WRITE(numout,*) 'crsini. gphit_crs(15,15)', gphit_crs(15,15) WRITE(numout,*) 'crsini. glamt_crs(15,15)', glamt_crs(15,15) WRITE(numout,*) 'crsini. count 1' CALL crsfun( gphiu, glamu, 'U', gphiu_crs, glamu_crs ) !cc WRITE(numout,*) 'crsini. gphiu_crs(15,15)', gphiu_crs(15,15) !cc WRITE(numout,*) 'crsini. glamu_crs(15,15)', glamu_crs(15,15) !cc WRITE(numout,*) 'crsini. count 2' CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) !cc WRITE(numout,*) 'crsini. gphiv_crs(15,15)', gphiv_crs(15,15) !cc WRITE(numout,*) 'crsini. glamv_crs(15,15)', glamv_crs(15,15) !cc WRITE(numout,*) 'crsini. count 3' CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) !cc WRITE(numout,*) 'crsini. gphif_crs(15,15)', gphif_crs(15,15) !cc WRITE(numout,*) 'crsini. glamf_crs(15,15)', glamf_crs(15,15) !cc WRITE(numout,*) 'crsini. count 4' ELSEIF ( zresty /= 0 .AND. zrestx == 0 ) THEN CALL crsfun( p_pgphi=gphiu, p_pglam=glamu, cd_type='T', p_cgphi=gphit_crs, p_cglam=glamt_crs ) CALL crsfun( p_pgphi=gphiu, p_pglam=glamu, cd_type='U', p_cgphi=gphiu_crs, p_cglam=glamu_crs ) CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) ELSEIF ( zresty == 0 .AND. zrestx /= 0 ) THEN CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='T', p_cgphi=gphit_crs, p_cglam=glamt_crs ) CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='U', p_cgphi=gphiu_crs, p_cglam=glamu_crs ) CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) ELSE CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='T', p_cgphi=gphit_crs, p_cglam=glamt_crs ) CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='U', p_cgphi=gphiu_crs, p_cglam=glamu_crs ) CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) ENDIF WRITE(numout,*) 'crsini. Finished coordinates' ! 3.c. Get the horizontal mesh ! 3.c.1 Horizontal scale factors ! CALL crsfun( cd_type='T', cd_op='SUM', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield2d_1=e1t_crs, p_cfield2d_2=e2t_crs ) ! CALL crsfun( cd_type='U', cd_op='SUM', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_cfield2d_1=e1u_crs, p_cfield2d_2=e2u_crs ) ! CALL crsfun( cd_type='V', cd_op='SUM', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_cfield2d_1=e1v_crs, p_cfield2d_2=e2v_crs ) ! CALL crsfun( cd_type='F', cd_op='SUM', p_pmask=fmask, p_e1=e1f, p_e2=e2f, p_cfield2d_1=e1f_crs, p_cfield2d_2=e2f_crs ) CALL crsfun( cd_type='T', cd_op='POS', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield2d_1=e1t_crs, p_cfield2d_2=e2t_crs ) CALL crsfun( cd_type='U', cd_op='POS', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_cfield2d_1=e1u_crs, p_cfield2d_2=e2u_crs ) CALL crsfun( cd_type='V', cd_op='POS', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_cfield2d_1=e1v_crs, p_cfield2d_2=e2v_crs ) CALL crsfun( cd_type='F', cd_op='POS', p_pmask=fmask, p_e1=e1f, p_e2=e2f, p_cfield2d_1=e1f_crs, p_cfield2d_2=e2f_crs ) e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:) WRITE(numout,*) 'crsini. Finished horizontal scale factors' ! 3.c.2 Coriolis factor SELECT CASE( jphgr_msh ) ! type of horizontal mesh CASE ( 0, 1, 4 ) ! mesh on the sphere ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) ) CASE DEFAULT WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' END SELECT ! 3.d. Get the initial vertical mesh ! nav_lon(jpi,jpj), nav_lat(jpi,jpj) ! nav_lev(jpk), e3t_0(jpk), e3w_0(jpk), gdep[t,w]_0(jpk) -> these stay constant ! 2D: mbathy (k-levels) ! 3D: fse3[t,u,v,w,f] (scale factors), gdep[t,u,v,w] (depth in meters) ! jes. TODO. Could probably make optional e1e2t in crsfun_TW... ! 3.d.1 mbathy ( vertical k-levels of bathymetry ) mbathy_crs(:,:) = jpkm1 mbkt_crs(:,:) = 1 mbku_crs(:,:) = 1 mbkv_crs(:,:) = 1 DO jj = 1, jpj_crs DO ji = 1, jpi_crs jk = 0 DO WHILE( tmask_crs(ji,jj,jk+1) > 0.) jk = jk + 1 ENDDO mbathy_crs(ji,jj) = float( jk ) ENDDO ENDDO ALLOCATE( zmbk(jpi_crs,jpj_crs) ) zmbk(:,:) = 0.0 zmbk(:,:) = REAL( mbathy_crs(:,:), wp ) ; CALL crs_lbc_lnk('T',1.0,zmbk) ; mbathy_crs(:,:) = INT( zmbk(:,:) ) ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' crsini : mbkt is ocean bottom k-index of T-, U-, V- and W-levels ' IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' ! mbkt_crs(:,:) = MAX( mbathy_crs(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) ! ! bottom k-index of W-level = mbkt+1 DO jj = 1, jpj_crsm1 ! bottom k-index of u- (v-) level DO ji = 1, jpi_crsm1 mbku_crs(ji,jj) = MIN( mbkt_crs(ji+1,jj ) , mbkt_crs(ji,jj) ) mbkv_crs(ji,jj) = MIN( mbkt_crs(ji ,jj+1) , mbkt_crs(ji,jj) ) END DO END DO WRITE(numout,*) 'crsini. Set mbku, mkbv' ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk zmbk(:,:) = 1.e0 zmbk(:,:) = REAL( mbku_crs(:,:), wp ) ; CALL crs_lbc_lnk('U',1.0,zmbk) ; mbku_crs (:,:) = MAX( INT( zmbk(:,:) ), 1 ) zmbk(:,:) = REAL( mbkv_crs(:,:), wp ) ; CALL crs_lbc_lnk('V',1.0,zmbk) ; mbkv_crs (:,:) = MAX( INT( zmbk(:,:) ), 1 ) ! WRITE(numout,*) 'crs_init = finished section 3d.1 jpi=', jpi, 'jpj=',jpj, ' jpk=', jpk ! 3.d.2 Vertical scale factors ALLOCATE( zfse3t(jpi,jpj,jpk), zfse3u(jpi,jpj,jpk), zfse3v(jpi,jpj,jpk), zfse3f(jpi,jpj,jpk), & & zfse3w(jpi,jpj,jpk), zfse3t_n(jpi,jpj,jpk), zfse3t_b(jpi,jpj,jpk) ) zfse3t(:,:,:) = fse3t(:,:,:) zfse3u(:,:,:) = fse3u(:,:,:) zfse3v(:,:,:) = fse3v(:,:,:) zfse3f(:,:,:) = fse3f(:,:,:) zfse3w(:,:,:) = fse3w(:,:,:) !CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3t, p_cfield3d=e3t_crs ) !CALL crsfun( p_e1e2t=e1e2t, cd_type='W', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3w, p_cfield3d=e3w_crs ) !CALL crsfun( p_e1e2t=e1e2t, cd_type='U', cd_op='MIN', p_cmask=umask_crs, p_ptmask=umask, p_pfield3d_1=zfse3u, p_cfield3d=e3u_crs ) !CALL crsfun( p_e1e2t=e1e2t, cd_type='V', cd_op='MIN', p_cmask=vmask_crs, p_ptmask=vmask, p_pfield3d_1=zfse3v, p_cfield3d=e3v_crs ) !CALL crsfun( p_e1e2t=e1e2t, cd_type='F', cd_op='MIN', p_cmask=fmask_crs, p_ptmask=fmask, p_pfield3d_1=zfse3f, p_cfield3d=e3f_crs ) CALL crs_e3_max( p_e3=zfse3t, cd_type='T', p_mask=tmask, p_e3_crs=e3t_crs) CALL crs_e3_max( p_e3=zfse3w, cd_type='W', p_mask=tmask, p_e3_crs=e3w_crs) ! Reset 0 to e3t_0 or e3w_0 DO jk = 1, jpk DO ji = 1, jpi_crs DO jj = 1, jpj_crs IF ( e3t_crs(ji,jj,jk) == 0 ) e3t_crs(ji,jj,jk) = e3t_0(jk) IF ( e3w_crs(ji,jj,jk) == 0 ) e3w_crs(ji,jj,jk) = e3w_0(jk) IF ( e3u_crs(ji,jj,jk) == 0 ) e3u_crs(ji,jj,jk) = e3t_0(jk) IF ( e3v_crs(ji,jj,jk) == 0 ) e3v_crs(ji,jj,jk) = e3t_0(jk) IF ( e3f_crs(ji,jj,jk) == 0 ) e3f_crs(ji,jj,jk) = e3t_0(jk) ENDDO ENDDO ENDDO ! 3.d.3 Vertical depth (meters) CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, & & p_pfield3d_1=gdept, p_cfield3d=gdept_crs ) CALL crsfun( p_e1e2t=e1e2t, cd_type='W', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, & & p_pfield3d_1=gdepw, p_cfield3d=gdepw_crs ) ! 3.d.4 Surfaces CALL crs_surf(p_e1=e1t, p_e2=e2t ,p_e3=zfse3w, cd_type='W', p_mask=tmask, surf_crs=e1e2w, surf_msk_crs=e1e2w_msk) CALL crs_surf(p_e1=e1u, p_e2=e2u ,p_e3=zfse3u, cd_type='U', p_mask=umask, surf_crs=e2e3u, surf_msk_crs=e2e3u_msk) CALL crs_surf(p_e1=e1v, p_e2=e2v ,p_e3=zfse3v, cd_type='V', p_mask=vmask, surf_crs=e1e3v, surf_msk_crs=e1e3v_msk) !--------------------------------------------------------- ! 4. Coarse grid ocean volume and averaging weights !--------------------------------------------------------- ! 4.a. Ocean volume or area unmasked and masked !! ! jes. May not need ocean_volume_crs_t, ocean_volume_crs_w as calculated already in trc_init as cvol CALL crsfun( cd_type='T', cd_op='VOL', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3t, & & p_cfield3d_1=ocean_volume_crs_t, p_cfield3d_2=facvol_t ) r1_bt_crs(:,:,:) = 0._wp bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:)* facvol_t(:,:,:) WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp/bt_crs(:,:,:) CALL crsfun( cd_type='W', cd_op='VOL', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3w, & & p_cfield3d_1=ocean_volume_crs_w, p_cfield3d_2=facvol_w ) CALL crsfun( cd_type='U', cd_op='SUM', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_fse3=zfse3u, p_cfield3d_1=facsurfu ) CALL crsfun( cd_type='V', cd_op='SUM', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_fse3=zfse3v, p_cfield3d_1=facsurfv ) ! 4.b. V, U and W surface area weights CALL crsfun( cd_type='U', cd_op='WGT', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_fse3=zfse3u, p_cfield3d_1=crs_surfu_wgt ) CALL crsfun( cd_type='V', cd_op='WGT', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_fse3=zfse3v, p_cfield3d_1=crs_surfv_wgt ) CALL crsfun( cd_type='W', cd_op='WGT', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield3d_1=crs_surfw_wgt ) ! 4.c. T volume weights CALL crsfun( cd_type='T', cd_op='WGT', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3t, p_cfield3d_1=crs_volt_wgt ) !--------------------------------------------------------- ! 5. Write out coarse meshmask (see OPA_SRC/DOM/domwri.F90 for ideas later) !--------------------------------------------------------- IF ( nn_msh_crs > 0 ) CALL crs_dom_wri WRITE(numout,*) 'crs_init done' !--------------------------------------------------------- ! 7. Finish and clean-up !--------------------------------------------------------- DEALLOCATE( zmbk ) DEALLOCATE( zfse3t, zfse3u, zfse3v, zfse3f ) DEALLOCATE( zfse3w, zfse3t_n, zfse3t_b ) END SUBROUTINE crs_init !!====================================================================== END MODULE crsini