Changeset 2715 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
- Timestamp:
- 2011-03-30T17:58:35+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r2528 r2715 4 4 !! LIM ice model : horizontal diffusion of sea-ice quantities 5 5 !!====================================================================== 6 !! History : LIM ! 2000-01 (LIM) Original code 7 !! - ! 2001-05 (G. Madec, R. Hordoir) opa norm 8 !! 1.0 ! 2002-08 (C. Ethe) F90, free form 9 !!---------------------------------------------------------------------- 6 10 #if defined key_lim3 7 11 !!---------------------------------------------------------------------- … … 10 14 !! lim_hdf : diffusion trend on sea-ice variable 11 15 !!---------------------------------------------------------------------- 12 !! * Modules used 13 USE dom_oce 14 USE in_out_manager 15 USE ice 16 USE lbclnk 17 USE lib_mpp 18 USE prtctl ! Print control 16 USE dom_oce ! ocean domain 17 USE ice ! LIM-3: ice variables 18 USE lbclnk ! lateral boundary condition - MPP exchanges 19 USE lib_mpp ! MPP library 20 USE prtctl ! Print control 21 USE in_out_manager ! I/O manager 19 22 20 23 IMPLICIT NONE 21 24 PRIVATE 22 25 23 !! * Routine accessibility 24 PUBLIC lim_hdf ! called by lim_tra 26 PUBLIC lim_hdf ! called by lim_tra 25 27 26 !! * Module variables 27 LOGICAL :: linit = .TRUE. ! ??? 28 LOGICAL :: linit = .TRUE. ! initialization flag (set to flase after the 1st call) 28 29 REAL(wp) :: epsi04 = 1e-04 ! constant 29 REAL(wp), DIMENSION(jpi,jpj) :: zfact ! ???30 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: efact ! metric coefficient 30 31 31 32 !! * Substitution 32 33 # include "vectopt_loop_substitute.h90" 33 34 !!---------------------------------------------------------------------- 34 !! NEMO/LIM3 3.3, UCL - NEMO Consortium (2010)35 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 35 36 !! $Id$ 36 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)37 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 37 38 !!---------------------------------------------------------------------- 38 39 39 CONTAINS 40 40 … … 43 43 !! *** ROUTINE lim_hdf *** 44 44 !! 45 !! ** purpose : Compute and add the diffusive trend on sea-ice 46 !! variables 45 !! ** purpose : Compute and add the diffusive trend on sea-ice variables 47 46 !! 48 47 !! ** method : Second order diffusive operator evaluated using a 49 !! Cranck-Nicholson time Scheme.48 !! Cranck-Nicholson time Scheme. 50 49 !! 51 50 !! ** Action : update ptab with the diffusive contribution 52 !!53 !! History :54 !! ! 00-01 (LIM) Original code55 !! ! 01-05 (G. Madec, R. Hordoir) opa norm56 !! ! 02-08 (C. Ethe) F90, free form57 51 !!------------------------------------------------------------------- 58 ! * Arguments 59 REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: & 60 ptab ! Field on which the diffusion is applied 61 REAL(wp), DIMENSION(jpi,jpj) :: & 62 ptab0 ! ??? 63 64 ! * Local variables 65 INTEGER :: ji, jj ! dummy loop indices 66 INTEGER :: & 67 its, iter ! temporary integers 68 CHARACTER (len=55) :: charout 69 REAL(wp) :: & 70 zalfa, zrlxint, zconv, zeps ! temporary scalars 71 REAL(wp), DIMENSION(jpi,jpj) :: & 72 zrlx, zflu, zflv, & ! temporary workspaces 73 zdiv0, zdiv ! " " 52 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 53 USE wrk_nemo, ONLY: zflu => wrk_2d_11, zdiv => wrk_2d_13, zrlx => wrk_2d_15 54 USE wrk_nemo, ONLY: zflv => wrk_2d_12, zdiv0 => wrk_2d_14, ztab0 => wrk_2d_16 55 ! 56 REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: ptab ! Field on which the diffusion is applied 57 ! 58 INTEGER :: ji, jj ! dummy loop indices 59 INTEGER :: its, iter, ierr ! local integers 60 REAL(wp) :: zalfa, zrlxint, zconv, zeps ! local scalars 61 CHARACTER(lc) :: charout ! local character 74 62 !!------------------------------------------------------------------- 75 76 ! Initialisation 77 ! --------------- 78 ! Time integration parameters 79 zalfa = 0.5 ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 80 its = 100 ! Maximum number of iteration 81 zeps = 2. * epsi04 82 83 ! Arrays initialization 84 ptab0 (:, : ) = ptab(:,:) 85 !bug zflu (:,jpj) = 0.e0 86 !bug zflv (:,jpj) = 0.e0 87 zdiv0(:, 1 ) = 0.e0 88 zdiv0(:,jpj) = 0.e0 89 IF( .NOT.lk_vopt_loop ) THEN 90 zflu (jpi,:) = 0.e0 91 zflv (jpi,:) = 0.e0 92 zdiv0(1, :) = 0.e0 93 zdiv0(jpi,:) = 0.e0 63 64 IF( wrk_in_use(2, 11,12,13,14,15,16) ) THEN 65 CALL ctl_stop( 'lim_hdf: requested workspace arrays unavailable' ) ; RETURN 94 66 ENDIF 95 67 96 ! Metric coefficient (compute at the first call and saved in 97 IF( linit ) THEN 68 ! !== Initialisation ==! 69 ! 70 IF( linit ) THEN ! Metric coefficient (compute at the first call and saved in efact) 71 ALLOCATE( efact(jpi,jpj) , STAT=ierr ) 72 IF( lk_mpp ) CALL mpp_sum( ierr ) 73 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' ) 98 74 DO jj = 2, jpjm1 99 75 DO ji = fs_2 , fs_jpim1 ! vector opt. 100 zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj ) + e1v(ji,jj) + e1v(ji,jj-1) ) & 101 & / ( e1t(ji,jj) * e2t(ji,jj) ) 76 efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 102 77 END DO 103 78 END DO 104 79 linit = .FALSE. 105 80 ENDIF 81 ! ! Time integration parameters 82 zalfa = 0.5_wp ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 83 its = 100 ! Maximum number of iteration 84 zeps = 2._wp * epsi04 85 ! 86 ztab0(:, : ) = ptab(:,:) ! Arrays initialization 87 zdiv0(:, 1 ) = 0._wp 88 zdiv0(:,jpj) = 0._wp 89 IF( .NOT.lk_vopt_loop ) THEN 90 zflu (jpi,:) = 0._wp 91 zflv (jpi,:) = 0._wp 92 zdiv0(1, :) = 0._wp 93 zdiv0(jpi,:) = 0._wp 94 ENDIF 106 95 107 108 ! Sub-time step loop 109 zconv = 1.e0 96 zconv = 1._wp !== horizontal diffusion using a Crant-Nicholson scheme ==! 110 97 iter = 0 111 112 ! !=================== 113 DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) ) ! Sub-time step loop 114 ! !=================== 115 ! incrementation of the sub-time step number 116 iter = iter + 1 117 118 ! diffusive fluxes in U- and V- direction 119 DO jj = 1, jpjm1 98 ! 99 DO WHILE( zconv > zeps .AND. iter <= its ) ! Sub-time step loop 100 ! 101 iter = iter + 1 ! incrementation of the sub-time step number 102 ! 103 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 120 104 DO ji = 1 , fs_jpim1 ! vector opt. 121 105 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) … … 123 107 END DO 124 108 END DO 125 126 ! diffusive trend : divergence of the fluxes 127 DO jj= 2, jpjm1 109 ! 110 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 128 111 DO ji = fs_2 , fs_jpim1 ! vector opt. 129 112 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & … … 131 114 END DO 132 115 END DO 133 134 ! save the first evaluation of the diffusive trend in zdiv0 135 IF( iter == 1 ) zdiv0(:,:) = zdiv(:,:) 136 137 ! XXXX iterative evaluation????? 138 DO jj = 2, jpjm1 116 ! 117 IF( iter == 1 ) zdiv0(:,:) = zdiv(:,:) ! save the 1st evaluation of the diffusive trend in zdiv0 118 ! 119 DO jj = 2, jpjm1 ! iterative evaluation 139 120 DO ji = fs_2 , fs_jpim1 ! vector opt. 140 zrlxint = ( ptab0(ji,jj) &141 & + rdt_ice * ( zalfa * ( zdiv(ji,jj) + zfact(ji,jj) * ptab(ji,jj) ) &121 zrlxint = ( ztab0(ji,jj) & 122 & + rdt_ice * ( zalfa * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) ) & 142 123 & + ( 1.0 - zalfa ) * zdiv0(ji,jj) ) ) & 143 & / ( 1.0 + zalfa * rdt_ice * zfact(ji,jj) )124 & / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 144 125 zrlx(ji,jj) = ptab(ji,jj) + om * ( zrlxint - ptab(ji,jj) ) 145 126 END DO 146 127 END DO 147 148 ! lateral boundary condition on zrlx 149 CALL lbc_lnk( zrlx, 'T', 1. ) 150 151 ! convergence test 152 zconv = 0.e0 128 CALL lbc_lnk( zrlx, 'T', 1. ) ! lateral boundary condition 129 ! 130 zconv = 0._wp ! convergence test 153 131 DO jj = 2, jpjm1 154 132 DO ji = fs_2, fs_jpim1 … … 156 134 END DO 157 135 END DO 158 IF( lk_mpp ) CALL mpp_max( zconv ) ! max over the global domain 159 160 DO jj = 1, jpj 161 DO ji = 1 , jpi 162 ptab(ji,jj) = zrlx(ji,jj) 163 END DO 164 END DO 165 166 ! !========================== 136 IF( lk_mpp ) CALL mpp_max( zconv ) ! max over the global domain 137 ! 138 ptab(:,:) = zrlx(:,:) 139 ! 167 140 END DO ! end of sub-time step loop 168 ! !==========================169 141 170 142 IF(ln_ctl) THEN 171 zrlx(:,:) = ptab(:,:) - ptab0(:,:)143 zrlx(:,:) = ptab(:,:) - ztab0(:,:) 172 144 WRITE(charout,FMT="(' lim_hdf : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter 173 CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout)145 CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout ) 174 146 ENDIF 175 176 147 ! 148 IF( wrk_not_released(2, 11,12,13,14,15,16) ) CALL ctl_stop('lim_hdf: failed to release workspace arrays') 149 ! 177 150 END SUBROUTINE lim_hdf 178 151
Note: See TracChangeset
for help on using the changeset viewer.