MODULE limmsh_2 !!====================================================================== !! *** MODULE limmsh_2 *** !! LIM 2.0 ice model : definition of the ice mesh parameters !!====================================================================== #if defined key_lim2 !!---------------------------------------------------------------------- !! 'key_lim2' LIM 2.0sea-ice model !!---------------------------------------------------------------------- !! lim_msh_2 : definition of the ice mesh !!---------------------------------------------------------------------- !! * Modules used USE phycst USE dom_oce USE dom_ice_2 USE lbclnk USE in_out_manager IMPLICIT NONE PRIVATE !! * Accessibility PUBLIC lim_msh_2 ! routine called by ice_ini_2.F90 !!---------------------------------------------------------------------- !! LIM 2.0, UCL-LOCEAN-IPSL (2005) !! $Id: limmsh_2.F90 1694 2009-10-30 16:06:23Z smasson $ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_msh_2 !!------------------------------------------------------------------- !! *** ROUTINE lim_msh_2 *** !! !! ** Purpose : Definition of the charact. of the numerical grid !! !! ** Action : - Initialisation of some variables !! - Definition of some constants linked with the grid !! - Definition of the metric coef. for the sea/ice !! - Initialization of the ice masks (tmsk, umsk) !! !! ** Refer. : Deleersnijder et al. Ocean Modelling 100, 7-10 !! !! ** History : !! original : 01-04 (LIM) !! addition : 02-08 (C. Ethe, G. Madec) !!--------------------------------------------------------------------- !! * Local variables INTEGER :: ji, jj ! dummy loop indices REAL(wp), DIMENSION(jpi,jpj) :: & zd2d1 , zd1d2 ! Derivative of zh2 (resp. zh1) in the x direction ! ! (resp. y direction) (defined at the center) REAL(wp) :: & zh1p , zh2p , & ! Idem zh1, zh2 for the bottom left corner of the grid zd2d1p, zd1d2p , & ! Idem zd2d1, zd1d2 for the bottom left corner of the grid zusden, zusden2 ! temporary scalars !!--------------------------------------------------------------------- IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'lim_msh_2 : LIM 2.0 sea-ice model, mesh initialization' WRITE(numout,*) '~~~~~~~~~' ENDIF !---------------------------------------------------------- ! Initialization of local and some global (common) variables !------------------------------------------------------------------ njeq = INT( jpj / 2 ) !i bug mpp potentiel njeqm1 = njeq - 1 fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad ) ! coriolis factor !i DO jj = 1, jpj !i zmsk(jj) = SUM( tmask(:,jj,:) ) ! = 0 if land everywhere on a j-line !!ii write(numout,*) jj, zind(jj) !i END DO IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN ! local domain include both hemisphere l_jeq = .TRUE. njeq = 1 DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 ) njeq = njeq + 1 END DO IF(lwp ) WRITE(numout,*) ' the equator is inside the domain at about njeq = ', njeq ELSEIF( fcor(1,1) < 0.e0 ) THEN l_jeq = .FALSE. njeq = jpj IF(lwp ) WRITE(numout,*) ' the model domain is entirely in the southern hemisphere: njeq = ', njeq ELSE l_jeq = .FALSE. njeq = 2 IF(lwp ) WRITE(numout,*) ' the model domain is entirely in the northern hemisphere: njeq = ', njeq ENDIF njeqm1 = njeq - 1 ! For each grid, definition of geometric tables !------------------------------------------------------------------ !------------------- ! Conventions : ! !------------------- ! indices 1 \ 2 <-> localisation in the 2 direction x \ y ! 3rd indice <-> localisation on the mesh : ! 0 = Centre ; 1 = corner W x(i-1/2) ; 2 = corner S y(j-1/2) ; ! 3 = corner SW x(i-1/2),y(j-1/2) !------------------- !!ibug ??? akappa(:,:,:,:) = 0.e0 wght(:,:,:,:) = 0.e0 alambd(:,:,:,:,:,:) = 0.e0 tmu(:,:) = 0.e0 !!i ! metric coefficients for sea ice dynamic !---------------------------------------- ! ! akappa DO jj = 2, jpj zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1) END DO CALL lbc_lnk( zd1d2, 'T', -1. ) DO ji = 2, jpi zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:) END DO CALL lbc_lnk( zd2d1, 'T', -1. ) akappa(:,:,1,1) = 1.0 / ( 2.0 * e1t(:,:) ) akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) ) akappa(:,:,2,1) = zd2d1(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) ) akappa(:,:,2,2) = 1.0 / ( 2.0 * e2t(:,:) ) ! ! weights (wght) DO jj = 2, jpj DO ji = 2, jpi zusden = 1. / ( ( e1t(ji,jj) + e1t(ji-1,jj ) ) & & * ( e2t(ji,jj) + e2t(ji ,jj-1) ) ) wght(ji,jj,1,1) = zusden * e1t(ji ,jj) * e2t(ji,jj ) wght(ji,jj,1,2) = zusden * e1t(ji ,jj) * e2t(ji,jj-1) wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj ) wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1) END DO END DO CALL lbc_lnk( wght(:,:,1,1), 'I', 1. ) ! CAUTION: even with the lbc_lnk at ice U-V-point CALL lbc_lnk( wght(:,:,1,2), 'I', 1. ) ! the value of wght at jpj is wrong CALL lbc_lnk( wght(:,:,2,1), 'I', 1. ) ! but it is never used CALL lbc_lnk( wght(:,:,2,2), 'I', 1. ) ! Coefficients for divergence of the stress tensor !------------------------------------------------- DO jj = 2, jpj DO ji = 2, jpi ! NO vector opt. zh1p = e1t(ji ,jj ) * wght(ji,jj,2,2) & & + e1t(ji-1,jj ) * wght(ji,jj,1,2) & & + e1t(ji ,jj-1) * wght(ji,jj,2,1) & & + e1t(ji-1,jj-1) * wght(ji,jj,1,1) zh2p = e2t(ji ,jj ) * wght(ji,jj,2,2) & & + e2t(ji-1,jj ) * wght(ji,jj,1,2) & & + e2t(ji ,jj-1) * wght(ji,jj,2,1) & & + e2t(ji-1,jj-1) * wght(ji,jj,1,1) ! better written but change the last digit and thus solver in less than 100 timestep ! zh1p = e1t(ji-1,jj ) * wght(ji,jj,1,2) + e1t(ji,jj ) * wght(ji,jj,2,2) & ! & + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1) ! zh2p = e2t(ji-1,jj ) * wght(ji,jj,1,2) + e2t(ji,jj ) * wght(ji,jj,2,2) & ! & + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1) !!ibug =0 zusden = 1.0 / ( zh1p * zh2p * 4.e0 ) zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 ) zusden2 = zusden * 2.0 zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj ) - e1t(ji,jj-1) + e1t(ji ,jj) ) zd2d1p = zusden * 0.5 * ( e2t(ji ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj ) - e2t(ji-1,jj) ) alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji ,jj-1) alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji ,jj ) alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1) alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj ) alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji ,jj-1) alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji ,jj ) alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1) alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj ) alambd(ji,jj,1,2,2,1) = zd1d2p alambd(ji,jj,1,2,2,2) = zd1d2p alambd(ji,jj,1,2,1,1) = zd1d2p alambd(ji,jj,1,2,1,2) = zd1d2p alambd(ji,jj,2,1,2,1) = zd2d1p alambd(ji,jj,2,1,2,2) = zd2d1p alambd(ji,jj,2,1,1,1) = zd2d1p alambd(ji,jj,2,1,1,2) = zd2d1p END DO END DO CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. ) ! CAUTION: even with the lbc_lnk at ice U-V point CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. ) ! the value of wght at jpj is wrong CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. ) ! but it is never used CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. ) ! CAUTION: idem CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. ) ! CAUTION: idem CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. ) ! CAUTION: idem CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. ) ! CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. ) ! ! Initialization of ice masks !---------------------------- tms(:,:) = tmask(:,:,1) ! ice T-point : use surface tmask !i here we can use umask with a i and j shift of -1,-1 tmu(:,1) = 0.e0 tmu(1,:) = 0.e0 DO jj = 2, jpj ! ice U.V-point: computed from ice T-point mask DO ji = 2, jpim1 ! NO vector opt. tmu(ji,jj) = tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1) END DO END DO !--lateral boundary conditions CALL lbc_lnk( tmu(:,:), 'I', 1. ) ! unmasked and masked area of T-grid cell area(:,:) = e1t(:,:) * e2t(:,:) END SUBROUTINE lim_msh_2 #else !!---------------------------------------------------------------------- !! Default option Dummy Module NO LIM sea-ice model !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_msh_2 ! Dummy routine END SUBROUTINE lim_msh_2 #endif !!====================================================================== END MODULE limmsh_2