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