Changeset 888 for trunk/NEMO/LIM_SRC_2
- Timestamp:
- 2008-04-11T19:05:03+02:00 (16 years ago)
- Location:
- trunk/NEMO/LIM_SRC_2
- Files:
-
- 1 added
- 2 deleted
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_2/dom_ice_2.F90
r823 r888 1 1 MODULE dom_ice_2 2 #if defined key_lim23 2 !!====================================================================== 4 3 !! *** MODULE dom_ice *** … … 7 6 !! History : 2.0 ! 03-08 (C. Ethe) Free form and module 8 7 !!---------------------------------------------------------------------- 9 8 #if defined key_lim2 10 9 !!---------------------------------------------------------------------- 11 10 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 12 !! $ Header$11 !! $ Id: $ 13 12 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 14 13 !!---------------------------------------------------------------------- -
trunk/NEMO/LIM_SRC_2/ice_2.F90
r823 r888 4 4 !! Sea Ice physics: diagnostics variables of ice defined in memory 5 5 !!===================================================================== 6 !! History : 2.0 ! 03-08 (C. Ethe) F90: Free form and module 7 !!---------------------------------------------------------------------- 6 8 #if defined key_lim2 7 9 !!---------------------------------------------------------------------- 8 10 !! 'key_lim2' : LIM 2.0 sea-ice model 9 11 !!---------------------------------------------------------------------- 10 !! History : 11 !! 2.0 ! 03-08 (C. Ethe) F90: Free form and module 12 !!---------------------------------------------------------------------- 13 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 14 !! $Header$ 15 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 12 !! LIM 2.0, UCL-LOCEAN-IPSL (2006) 13 !! $ Id: $ 14 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 16 15 !!---------------------------------------------------------------------- 17 16 !! * Modules used … … 21 20 PRIVATE 22 21 23 !! * Share Module variables 24 INTEGER , PUBLIC :: & !!: ** ice-dynamic namelist (namicedyn) ** 25 nbiter = 1 , & !: number of sub-time steps for relaxation 26 nbitdr = 250 !: maximum number of iterations for relaxation 22 !!* ice-dynamic namelist (namicedyn) * 23 INTEGER , PUBLIC :: nbiter = 1 !: number of sub-time steps for relaxation 24 INTEGER , PUBLIC :: nbitdr = 250 !: maximum number of iterations for relaxation 25 REAL(wp), PUBLIC :: epsd = 1.0e-20 !: tolerance parameter for dynamic 26 REAL(wp), PUBLIC :: alpha = 0.5 !: coefficient for semi-implicit coriolis 27 REAL(wp), PUBLIC :: dm = 0.6e+03 !: diffusion constant for dynamics 28 REAL(wp), PUBLIC :: om = 0.5 !: relaxation constant 29 REAL(wp), PUBLIC :: resl = 5.0e-05 !: maximum value for the residual of relaxation 30 REAL(wp), PUBLIC :: cw = 5.0e-03 !: drag coefficient for oceanic stress 31 REAL(wp), PUBLIC :: angvg = 0.e0 !: turning angle for oceanic stress 32 REAL(wp), PUBLIC :: pstar = 1.0e+04 !: first bulk-rheology parameter 33 REAL(wp), PUBLIC :: c_rhg = 20.e0 !: second bulk-rhelogy parameter 34 REAL(wp), PUBLIC :: etamn = 0.e+07 !: minimun value for viscosity 35 REAL(wp), PUBLIC :: creepl = 2.e-08 !: creep limit 36 REAL(wp), PUBLIC :: ecc = 2.e0 !: eccentricity of the elliptical yield curve 37 REAL(wp), PUBLIC :: ahi0 = 350.e0 !: sea-ice hor. eddy diffusivity coeff. (m2/s) 27 38 28 REAL(wp), PUBLIC :: & !!: ** ice-dynamic namelist (namicedyn) ** 29 epsd = 1.0e-20, & !: tolerance parameter for dynamic 30 alpha = 0.5 , & !: coefficient for semi-implicit coriolis 31 dm = 0.6e+03, & !: diffusion constant for dynamics 32 om = 0.5 , & !: relaxation constant 33 resl = 5.0e-05, & !: maximum value for the residual of relaxation 34 cw = 5.0e-03, & !: drag coefficient for oceanic stress 35 angvg = 0.e0 , & !: turning angle for oceanic stress 36 pstar = 1.0e+04, & !: first bulk-rheology parameter 37 c_rhg = 20.e0 , & !: second bulk-rhelogy parameter 38 etamn = 0.e+07, & !: minimun value for viscosity 39 creepl = 2.e-08, & !: creep limit 40 ecc = 2.e0 , & !: eccentricity of the elliptical yield curve 41 ahi0 = 350.e0 !: sea-ice hor. eddy diffusivity coeff. (m2/s) 39 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( ecc * ecc ) 40 REAL(wp), PUBLIC :: rhoco !: = rau0 * cw 41 REAL(wp), PUBLIC :: sangvg, cangvg !: sin and cos of the turning angle for ocean stress 42 REAL(wp), PUBLIC :: pstarh !: pstar / 2.0 42 43 43 REAL(wp), PUBLIC :: & !: 44 usecc2 , & !: = 1.0 / ( ecc * ecc ) 45 rhoco , & !: = rau0 * cw 46 sangvg, cangvg , & !: sin and cos of the turning angle for ocean stress 47 pstarh !: pstar / 2.0 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ahiu , ahiv !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 45 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: pahu , pahv !: ice hor. eddy diffusivity coef. at ocean U- and V-points 46 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hsnm , hicm !: mean snow and ice thicknesses 47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ust2s !: friction velocity 48 48 49 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 50 u_oce, v_oce, & !: surface ocean velocity used in ice dynamics 51 ahiu , ahiv , & !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 52 pahu , pahv , & !: ice hor. eddy diffusivity coef. at ocean U- and V-points 53 hsnm , hicm , & !: mean snow and ice thicknesses 54 ust2s !: friction velocity 49 !!* diagnostic quantities 50 !! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: firic !: IR flux over the ice (only used for outputs) 51 !! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fcsic !: Sensible heat flux over the ice (only used for outputs) 52 !! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fleic !: Latent heat flux over the ice (only used for outputs) 53 !! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qlatic !: latent flux (only used for outputs) 54 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rdvosif !: Variation of volume at surface (only used for outputs) 55 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rdvobif !: Variation of ice volume at the bottom ice (only used for outputs) 56 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fdvolif !: Total variation of ice volume (only used for outputs) 57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rdvonif !: Lateral Variation of ice volume (only used for outputs) 55 58 56 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 57 sst_ini, & !: sst read from a file for ice model initialization 58 sss_ini !: sss read from a file for ice model initialization 59 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sist !: Sea-Ice Surface Temperature (Kelvin ??? degree ??? I don't know) 60 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: tfu !: Freezing/Melting point temperature of sea water at SSS 61 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hicif !: Ice thickness 62 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hsnif !: Snow thickness 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hicifp !: Ice production/melting 64 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: frld !: Leads fraction = 1-a/totalarea 65 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: phicif !: ice thickness at previous time 66 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: pfrld !: Leads fraction at previous time 67 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qstoif !: Energy stored in the brine pockets 68 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fbif !: Heat flux at the ice base 69 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rdmsnif !: Variation of snow mass 70 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rdmicif !: Variation of ice mass 71 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qldif !: heat balance of the lead (or of the open ocean) 72 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qcmif !: Energy needed to bring the ocean surface layer until its freezing 73 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fdtcn !: net downward heat flux from the ice to the ocean 74 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qdtcn !: energy from the ice to the ocean point (at a factor 2) 75 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: thcm !: part of the solar energy used in the lead heat budget 76 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fstric !: Solar flux transmitted trough the ice 77 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ffltbif !: Array linked with the max heat contained in brine pockets (?) 78 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fscmbq !: Linked with the solar flux below the ice (?) 79 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: fsbbq !: Also linked with the solar flux below the ice (?) 80 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: qfvbq !: Array used to store energy in case of toral lateral ablation (?) 81 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: dmgwi !: Variation of the mass of snow ice 59 82 60 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 61 firic , & !: IR flux over the ice (only used for outputs) 62 fcsic , & !: Sensible heat flux over the ice (only used for outputs) 63 fleic , & !: Latent heat flux over the ice (only used for outputs) 64 qlatic , & !: latent flux 65 rdvosif, & !: Variation of volume at surface (only used for outputs) 66 rdvobif, & !: Variation of ice volume at the bottom ice (only used for outputs) 67 fdvolif, & !: Total variation of ice volume (only used for outputs) 68 rdvonif, & !: Lateral Variation of ice volume (only used for outputs) 69 sist , & !: Sea-Ice Surface Temperature (Kelvin ??? degree ??? I don't know) 70 tfu , & !: Melting point temperature of sea water 71 hsnif , & !: Snow thickness 72 hicif , & !: Ice thickness 73 hicifp , & !: Ice production/melting 74 frld , & !: Leads fraction = 1-a/totalarea 75 phicif , & !: ice thickness at previous time 76 pfrld , & !: Leads fraction at previous time 77 qstoif , & !: Energy stored in the brine pockets 78 fbif , & !: Heat flux at the ice base 79 rdmsnif, & !: Variation of snow mass 80 rdmicif, & !: Variation of ice mass 81 qldif , & !: heat balance of the lead (or of the open ocean) 82 qcmif , & !: Energy needed to bring the ocean surface layer until its freezing 83 fdtcn , & !: net downward heat flux from the ice to the ocean 84 qdtcn , & !: energy from the ice to the ocean 85 ! ! point (at a factor 2) 86 thcm , & !: part of the solar energy used in the lead heat budget 87 fstric , & !: Solar flux transmitted trough the ice 88 ffltbif, & !: Array linked with the max heat contained in brine pockets (?) 89 fscmbq , & !: Linked with the solar flux below the ice (?) 90 fsbbq , & !: Also linked with the solar flux below the ice (?) 91 qfvbq , & !: Array used to store energy in case of toral lateral ablation (?) 92 dmgwi !: Variation of the mass of snow ice 83 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: albege !: Albedo of the snow or ice (only for outputs) 84 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: albecn !: Albedo of the ocean (only for outputs) 85 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: tauc !: Cloud optical depth 93 86 94 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 95 albege , & !: Albedo of the snow or ice (only for outputs) 96 albecn , & !: Albedo of the ocean (only for outputs) 97 tauc , & !: Cloud optical depth 98 sdvt !: u*^2/(Stress/density) 87 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ui_ice, vi_ice !: two components of the ice velocity at I-point (m/s) 88 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ui_oce, vi_oce !: two components of the ocean velocity at I-point (m/s) 99 89 90 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax) :: scal0 !: ??? 91 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) :: tbif !: Temperature inside the ice/snow layer 100 92 101 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 102 u_ice, v_ice, & !: two components of the ice velocity (m/s) 103 tio_u, tio_v !: two components of the ice-ocean stress (N/m2) 93 !! REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) :: reslum !: Relative absorption of solar radiation in each ocean level 104 94 105 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax) :: & !: 106 scal0 !: ??? 107 108 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) :: & !: 109 tbif !: Temperature inside the ice/snow layer 110 111 REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) :: & !: 112 reslum !: Relative absorption of solar radiation in each ocean level 113 114 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 115 sxice, syice, sxxice, syyice, sxyice, & !: moments for advection 116 sxsn, sysn, sxxsn, syysn, sxysn, & !: 117 sxa, sya, sxxa, syya, sxya, & !: 118 sxc0, syc0, sxxc0, syyc0, sxyc0, & !: 119 sxc1, syc1, sxxc1, syyc1, sxyc1, & !: 120 sxc2, syc2, sxxc2, syyc2, sxyc2, & !: 121 sxst, syst, sxxst, syyst, sxyst !: 95 !!* moment used in the advection scheme 96 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sxice, syice, sxxice, syyice, sxyice !: for ice volume 97 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sxsn, sysn, sxxsn, syysn, sxysn !: for snow volume 98 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sxa, sya, sxxa, syya, sxya !: for ice cover area 99 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sxc0, syc0, sxxc0, syyc0, sxyc0 !: for heat content of snow 100 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sxc1, syc1, sxxc1, syyc1, sxyc1 !: for heat content of 1st ice layer 101 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sxc2, syc2, sxxc2, syyc2, sxyc2 !: for heat content of 2nd ice layer 102 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sxst, syst, sxxst, syyst, sxyst !: for heat content of brine pockets 122 103 123 104 #else -
trunk/NEMO/LIM_SRC_2/iceini_2.F90
r823 r888 17 17 USE dom_oce 18 18 USE dom_ice_2 19 USE in_out_manager20 19 USE ice_oce ! ice variables 21 USE flx_oce 20 USE sbc_oce ! surface boundary condition: ocean 21 USE sbc_ice ! surface boundary condition: ice 22 22 USE phycst ! Define parameters for the routines 23 23 USE ocfzpt … … 27 27 USE limrst_2 28 28 USE ini1d ! initialization of the 1D configuration 29 USE in_out_manager 29 30 30 31 IMPLICIT NONE … … 40 41 !!---------------------------------------------------------------------- 41 42 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 42 !! $ Header$43 !! $ Id: $ 43 44 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 44 45 !!---------------------------------------------------------------------- … … 62 63 63 64 ! Louvain la Neuve Ice model 64 IF( nacc == 1 ) THEN 65 dtsd2 = nfice * rdtmin * 0.5 66 rdt_ice = nfice * rdtmin 67 ELSE 68 dtsd2 = nfice * rdt * 0.5 69 rdt_ice = nfice * rdt 70 ENDIF 65 dtsd2 = nn_fsbc * rdttra(1) * 0.5 66 rdt_ice = nn_fsbc * rdttra(1) 71 67 72 68 CALL lim_msh_2 ! ice mesh initialization -
trunk/NEMO/LIM_SRC_2/limadv_2.F90
r823 r888 33 33 !!---------------------------------------------------------------------- 34 34 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 35 !! $ Header$35 !! $ Id: $ 36 36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 37 37 !!---------------------------------------------------------------------- -
trunk/NEMO/LIM_SRC_2/limdia_2.F90
r823 r888 19 19 USE par_ice_2 ! ice parameters 20 20 USE ice_oce ! ice variables 21 USE sbc_oce ! surface boundary condition variables 21 22 USE daymod ! 22 23 USE dom_ice_2 ! … … 28 29 PRIVATE 29 30 30 PUBLIC lim_dia_2 ! called by ice_step31 PUBLIC lim_dia_2 ! called by sbc_ice_lim_2 31 32 INTEGER, PUBLIC :: ntmoy = 1 , & !: instantaneous values of ice evolution or averaging ntmoy 32 33 & ninfo = 1 !: frequency of ouputs on file ice_evolu in case of averaging … … 58 59 !!---------------------------------------------------------------------- 59 60 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 60 !! $ Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdia.F90,v 1.9 2007/06/29 17:03:12 opalod Exp $61 !! $ Id: $ 61 62 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 62 63 !!---------------------------------------------------------------------- … … 87 88 88 89 nv = 1 89 vinfor(nv) = REAL( kt + n fice- 1 )90 vinfor(nv) = REAL( kt + nn_fsbc - 1 ) 90 91 nv = nv + 1 91 92 vinfor(nv) = nyear … … 107 108 zicevol = zarea * hicif(ji,jj) 108 109 zsnwvol = zarea * hsnif(ji,jj) 109 zicespd = zicevol * ( u _ice(ji,jj) * u_ice(ji,jj) &110 & + v _ice(ji,jj) * v_ice(ji,jj) )110 zicespd = zicevol * ( ui_ice(ji,jj) * ui_ice(ji,jj) & 111 & + vi_ice(ji,jj) * vi_ice(ji,jj) ) 111 112 vinfor(nv+ 1) = vinfor(nv+ 1) + zarea 112 113 vinfor(nv+ 3) = vinfor(nv+ 3) + zextent15 … … 133 134 zicevol = zarea * hicif(ji,jj) 134 135 zsnwvol = zarea * hsnif(ji,jj) 135 zicespd = zicevol * ( u _ice(ji,jj) * u_ice(ji,jj) &136 & + v _ice(ji,jj) * v_ice(ji,jj) )136 zicespd = zicevol * ( ui_ice(ji,jj) * ui_ice(ji,jj) & 137 & + vi_ice(ji,jj) * vi_ice(ji,jj) ) 137 138 vinfor(nv+ 1) = vinfor(nv+ 1) + zarea 138 139 vinfor(nv+ 3) = vinfor(nv+ 3) + zextent15 … … 154 155 155 156 ! oututs on file ice_evolu 156 IF( MOD( kt + n fice- 1, ninfo ) == 0 ) THEN157 IF( MOD( kt + nn_fsbc - 1, ninfo ) == 0 ) THEN 157 158 WRITE(numevo_ice,fmtw) ( titvar(jv), vinfom(jv)/naveg, jv = 1, nvinfo ) 158 159 naveg = 0 … … 227 228 228 229 ! Definition et Ecriture de l'entete : nombre d'enregistrements 229 ndeb = ( nit000 - 1 + n fice- 1 ) / ninfo230 IF( nit000 - 1 + n fice== 1 ) ndeb = -1231 232 nferme = ( nitend + n fice - 1 ) / ninfo ! nit000 - 1 + nfice- 1 + nitend - nit000 + 1230 ndeb = ( nit000 - 1 + nn_fsbc - 1 ) / ninfo 231 IF( nit000 - 1 + nn_fsbc == 1 ) ndeb = -1 232 233 nferme = ( nitend + nn_fsbc - 1 ) / ninfo ! nit000 - 1 + nn_fsbc - 1 + nitend - nit000 + 1 233 234 ntot = nferme - ndeb 234 235 ndeb = ninfo * ( 1 + ndeb ) -
trunk/NEMO/LIM_SRC_2/limdmp_2.F90
r823 r888 40 40 !!---------------------------------------------------------------------- 41 41 !! LIM 2.0 , UCL-LOCEAN-IPSL (2006) 42 !! $ Header$42 !! $ Id: $ 43 43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 44 !!---------------------------------------------------------------------- -
trunk/NEMO/LIM_SRC_2/limdyn_2.F90
r823 r888 4 4 !! Sea-Ice dynamics : 5 5 !!====================================================================== 6 !! History : 1.0 ! 01-04 (LIM) Original code 7 !! 2.0 ! 02-08 (C. Ethe, G. Madec) F90, mpp 8 !! 2.0 ! 03-08 (C. Ethe) add lim_dyn_init 9 !! 2.0 ! 06-07 (G. Madec) Surface module 10 !!--------------------------------------------------------------------- 6 11 #if defined key_lim2 7 12 !!---------------------------------------------------------------------- … … 11 16 !! lim_dyn_init_2 : initialization and namelist read 12 17 !!---------------------------------------------------------------------- 13 !! * Modules used 14 USE phycst 15 USE in_out_manager ! I/O manager 16 USE dom_ice_2 17 USE dom_oce ! ocean space and time domain 18 USE ice_2 19 USE ice_oce 20 USE iceini_2 21 USE limistate_2 22 USE limrhg_2 ! ice rheology 23 USE lbclnk 24 USE lib_mpp 25 USE prtctl ! Print control 18 USE dom_oce ! ocean space and time domain 19 USE sbc_oce ! 20 USE phycst ! 21 USE ice_2 ! 22 USE ice_oce ! 23 USE dom_ice_2 ! 24 USE iceini_2 ! 25 USE limistate_2 ! 26 USE limrhg_2 ! ice rheology 27 28 USE lbclnk ! 29 USE lib_mpp ! 30 USE in_out_manager ! I/O manager 31 USE prtctl ! Print control 26 32 27 33 IMPLICIT NONE 28 34 PRIVATE 29 35 30 !! * Accessibility 31 PUBLIC lim_dyn_2 ! routine called by ice_step 36 PUBLIC lim_dyn_2 ! routine called by sbc_ice_lim 32 37 33 38 !! * Module variables 34 39 REAL(wp) :: rone = 1.e0 ! constant value 35 40 36 !!---------------------------------------------------------------------- 37 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 38 !! $Header$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 41 # include "vectopt_loop_substitute.h90" 42 !!---------------------------------------------------------------------- 43 !! LIM 2.0, UCL-LOCEAN-IPSL (2006) 44 !! $ Id: $ 45 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 46 !!---------------------------------------------------------------------- 41 47 … … 46 52 !! *** ROUTINE lim_dyn_2 *** 47 53 !! 48 !! ** Purpose : compute ice velocity and ocean-ice stress54 !! ** Purpose : compute ice velocity and ocean-ice friction velocity 49 55 !! 50 56 !! ** Method : … … 52 58 !! ** Action : - Initialisation 53 59 !! - Call of the dynamic routine for each hemisphere 54 !! - computation of the stress at the ocean surface60 !! - computation of the friction velocity at the sea-ice base 55 61 !! - treatment of the case if no ice dynamic 56 !! History :57 !! 1.0 ! 01-04 (LIM) Original code58 !! 2.0 ! 02-08 (C. Ethe, G. Madec) F90, mpp59 62 !!--------------------------------------------------------------------- 60 63 INTEGER, INTENT(in) :: kt ! number of iteration 61 62 INTEGER :: ji, jj ! dummy loop indices 63 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 64 REAL(wp) :: & 65 ztairx, ztairy, & ! tempory scalars 66 zsang , zmod, & 67 ztglx , ztgly , & 68 zt11, zt12, zt21, zt22 , & 69 zustm, zsfrld, zsfrldm4, & 70 zu_ice, zv_ice, ztair2 71 REAL(wp),DIMENSION(jpj) :: & 72 zind, & ! i-averaged indicator of sea-ice 73 zmsk ! i-averaged of tmask 64 !! 65 INTEGER :: ji, jj ! dummy loop indices 66 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 67 REAL(wp) :: zcoef ! temporary scalar 68 REAL(wp), DIMENSION(jpj) :: zind ! i-averaged indicator of sea-ice 69 REAL(wp), DIMENSION(jpj) :: zmsk ! i-averaged of tmask 70 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io ! ice-ocean velocity 74 71 !!--------------------------------------------------------------------- 75 72 76 IF( kt == nit000 73 IF( kt == nit000 ) CALL lim_dyn_init_2 ! Initialization (first time-step only) 77 74 78 IF 79 75 IF( ln_limdyn ) THEN 76 ! 80 77 ! Mean ice and snow thicknesses. 81 78 hsnm(:,:) = ( 1.0 - frld(:,:) ) * hsnif(:,:) 82 79 hicm(:,:) = ( 1.0 - frld(:,:) ) * hicif(:,:) 83 84 u_oce(:,:) = u_io(:,:) * tmu(:,:) 85 v_oce(:,:) = v_io(:,:) * tmu(:,:) 86 87 ! ! Rheology (ice dynamics) 88 ! ! ======== 80 ! 81 ! ! Rheology (ice dynamics) 82 ! ! ======== 89 83 90 84 ! Define the j-limits where ice rheology is computed … … 94 88 i_j1 = 1 95 89 i_jpj = jpj 96 IF(ln_ctl) THEN 97 CALL prt_ctl_info('lim_dyn : i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj) 98 ENDIF 90 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 99 91 CALL lim_rhg_2( i_j1, i_jpj ) 100 92 ! 101 93 ELSE ! optimization of the computational area 102 94 ! 103 95 DO jj = 1, jpj 104 96 zind(jj) = SUM( frld (:,jj ) ) ! = FLOAT(jpj) if ocean everywhere on a j-line 105 97 zmsk(jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 106 !!i write(numout,*) narea, 'limdyn' , jj, zind(jj), zmsk(jj) 107 END DO 108 98 END DO 99 ! 109 100 IF( l_jeq ) THEN ! local domain include both hemisphere 110 101 ! ! Rheology is computed in each hemisphere … … 118 109 i_j1 = MAX( 1, i_j1-1 ) 119 110 IF(ln_ctl) WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 120 111 ! 121 112 CALL lim_rhg_2( i_j1, i_jpj ) 122 113 ! 123 114 ! Southern hemisphere 124 115 i_j1 = 1 … … 129 120 i_jpj = MIN( jpj, i_jpj+2 ) 130 121 IF(ln_ctl) WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 131 122 ! 132 123 CALL lim_rhg_2( i_j1, i_jpj ) 133 124 ! 134 125 ELSE ! local domain extends over one hemisphere only 135 126 ! ! Rheology is computed only over the ice cover … … 148 139 149 140 IF(ln_ctl) WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 150 141 ! 151 142 CALL lim_rhg_2( i_j1, i_jpj ) 152 143 ! 153 144 ENDIF 154 145 ! 155 146 ENDIF 156 147 157 IF(ln_ctl) THEN 158 CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :') 159 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 160 ENDIF 161 162 ! ! Ice-Ocean stress 163 ! ! ================ 148 IF(ln_ctl) CALL prt_ctl(tab2d_1=ui_ice , clinfo1=' lim_dyn : ui_ice :', tab2d_2=vi_ice , clinfo2=' vi_ice :') 149 150 ! computation of friction velocity 151 ! -------------------------------- 152 ! ice-ocean velocity at U & V-points (ui_ice vi_ice at I-point ; ssu_m, ssv_m at U- & V-points) 153 154 DO jj = 1, jpjm1 155 DO ji = 1, fs_jpim1 ! vector opt. 156 zu_io(ji,jj) = 0.5 * ( ui_ice(ji+1,jj+1) + ui_ice(ji+1,jj ) ) - ssu_m(ji,jj) 157 zv_io(ji,jj) = 0.5 * ( vi_ice(ji+1,jj+1) + vi_ice(ji ,jj+1) ) - ssv_m(ji,jj) 158 END DO 159 END DO 160 ! frictional velocity at T-point 164 161 DO jj = 2, jpjm1 165 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 166 DO ji = 2, jpim1 167 ! computation of wind stress over ocean in X and Y direction 168 #if defined key_coupled && defined key_lim_cp1 169 ztairx = frld(ji-1,jj ) * gtaux(ji-1,jj ) + frld(ji,jj ) * gtaux(ji,jj ) & 170 & + frld(ji-1,jj-1) * gtaux(ji-1,jj-1) + frld(ji,jj-1) * gtaux(ji,jj-1) 171 172 ztairy = frld(ji-1,jj ) * gtauy(ji-1,jj ) + frld(ji,jj ) * gtauy(ji,jj ) & 173 & + frld(ji-1,jj-1) * gtauy(ji-1,jj-1) + frld(ji,jj-1) * gtauy(ji,jj-1) 174 #else 175 zsfrld = frld(ji,jj) + frld(ji-1,jj) + frld(ji-1,jj-1) + frld(ji,jj-1) 176 ztairx = zsfrld * gtaux(ji,jj) 177 ztairy = zsfrld * gtauy(ji,jj) 178 #endif 179 zsfrldm4 = 4 - frld(ji,jj) - frld(ji-1,jj) - frld(ji-1,jj-1) - frld(ji,jj-1) 180 zu_ice = u_ice(ji,jj) - u_oce(ji,jj) 181 zv_ice = v_ice(ji,jj) - v_oce(ji,jj) 182 zmod = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 183 ztglx = zsfrldm4 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 184 ztgly = zsfrldm4 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 185 186 tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 4 * rau0 ) 187 tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 4 * rau0 ) 162 DO ji = fs_2, fs_jpim1 ! vector opt. 163 ust2s(ji,jj) = 0.5 * cw & 164 & * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 165 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tms(ji,jj) 188 166 END DO 189 167 END DO 190 191 ! computation of friction velocity 168 ! 169 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 170 ! 171 zcoef = SQRT( 0.5 ) / rau0 192 172 DO jj = 2, jpjm1 193 DO ji = 2, jpim1 194 195 zu_ice = u_ice(ji-1,jj-1) - u_oce(ji-1,jj-1) 196 zv_ice = v_ice(ji-1,jj-1) - v_oce(ji-1,jj-1) 197 zt11 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 198 199 zu_ice = u_ice(ji-1,jj) - u_oce(ji-1,jj) 200 zv_ice = v_ice(ji-1,jj) - v_oce(ji-1,jj) 201 zt12 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 202 203 zu_ice = u_ice(ji,jj-1) - u_oce(ji,jj-1) 204 zv_ice = v_ice(ji,jj-1) - v_oce(ji,jj-1) 205 zt21 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 206 207 zu_ice = u_ice(ji,jj) - u_oce(ji,jj) 208 zv_ice = v_ice(ji,jj) - v_oce(ji,jj) 209 zt22 = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 210 211 ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj) 212 213 zustm = ( 1 - frld(ji,jj) ) * 0.25 * ( zt11 + zt12 + zt21 + zt22 ) & 214 & + frld(ji,jj) * SQRT( ztair2 ) 215 216 ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj) 173 DO ji = fs_2, fs_jpim1 ! vector opt. 174 ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 175 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) 217 176 END DO 218 177 END DO 219 220 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 221 222 DO jj = 2, jpjm1 223 DO ji = 2, jpim1 224 #if defined key_coupled && defined key_lim_cp1 225 tio_u(ji,jj) = - ( gtaux(ji ,jj ) + gtaux(ji-1,jj ) & 226 & + gtaux(ji-1,jj-1) + gtaux(ji ,jj-1) ) / ( 4 * rau0 ) 227 228 tio_v(ji,jj) = - ( gtauy(ji ,jj ) + gtauy(ji-1,jj ) & 229 & + gtauy(ji-1,jj-1) + gtauy(ji ,jj-1) ) / ( 4 * rau0 ) 230 #else 231 tio_u(ji,jj) = - gtaux(ji,jj) / rau0 232 tio_v(ji,jj) = - gtauy(ji,jj) / rau0 233 #endif 234 ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj) 235 zustm = SQRT( ztair2 ) 236 237 ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj) 238 END DO 239 END DO 240 178 ! 241 179 ENDIF 242 180 ! 243 181 CALL lbc_lnk( ust2s, 'T', 1. ) ! T-point 244 CALL lbc_lnk( tio_u, 'I', -1. ) ! I-point (i.e. ice U-V point) 245 CALL lbc_lnk( tio_v, 'I', -1. ) ! I-point (i.e. ice U-V point) 246 247 IF(ln_ctl) THEN 248 CALL prt_ctl(tab2d_1=tio_u , clinfo1=' lim_dyn : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v :') 249 CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn : ust2s :') 250 ENDIF 182 ! 183 IF(ln_ctl) CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn : ust2s :') 251 184 252 185 END SUBROUTINE lim_dyn_2 … … 257 190 !! *** ROUTINE lim_dyn_init_2 *** 258 191 !! 259 !! ** Purpose : Physical constants and parameters linked to the ice260 !! dynamics261 !! 262 !! ** Method : Read the namicedyn namelist and check the ice-dynamic263 !! parameter values called at the first timestep (nit000)192 !! ** Purpose : Physical constants and parameters linked to the ice 193 !! dynamics 194 !! 195 !! ** Method : Read the namicedyn namelist and check the ice-dynamic 196 !! parameter values 264 197 !! 265 198 !! ** input : Namelist namicedyn 266 !!267 !! history :268 !! 8.5 ! 03-08 (C. Ethe) original code269 199 !!------------------------------------------------------------------- 270 200 NAMELIST/namicedyn/ epsd, alpha, & … … 273 203 !!------------------------------------------------------------------- 274 204 275 ! Define the initial parameters 276 ! ------------------------- 277 278 ! Read Namelist namicedyn 279 REWIND ( numnam_ice ) 205 REWIND ( numnam_ice ) ! Read Namelist namicedyn 280 206 READ ( numnam_ice , namicedyn ) 281 IF(lwp) THEN 207 208 IF(lwp) THEN ! Control print 282 209 WRITE(numout,*) 283 210 WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics ' … … 291 218 WRITE(numout,*) ' maximum value for the residual of relaxation resl = ', resl 292 219 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw 293 WRITE(numout,*) ' turning angle for oceanic stress angvg = ', angvg 220 WRITE(numout,*) ' turning angle for oceanic stress angvg = ', angvg, ' degrees' 294 221 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar 295 222 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg … … 303 230 usecc2 = 1.0 / ( ecc * ecc ) 304 231 rhoco = rau0 * cw 305 angvg = angvg * rad 232 angvg = angvg * rad ! convert angvg from degree to radian 306 233 sangvg = SIN( angvg ) 307 234 cangvg = COS( angvg ) 308 235 pstarh = pstar / 2.0 309 sdvt(:,:) = 0.e0 310 311 ! Diffusion coefficients. 312 ahiu(:,:) = ahi0 * umask(:,:,1) 236 ! 237 ahiu(:,:) = ahi0 * umask(:,:,1) ! Ice eddy Diffusivity coefficients. 313 238 ahiv(:,:) = ahi0 * vmask(:,:,1) 314 239 ! 315 240 END SUBROUTINE lim_dyn_init_2 316 241 -
trunk/NEMO/LIM_SRC_2/limhdf_2.F90
r823 r888 6 6 #if defined key_lim2 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim2' iLIM 2.0 sea-ice model8 !! 'key_lim2' LIM 2.0 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_hdf_2 : diffusion trend on sea-ice variable … … 34 34 !!---------------------------------------------------------------------- 35 35 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 36 !! $ Header$36 !! $ Id: $ 37 37 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 38 38 !!---------------------------------------------------------------------- -
trunk/NEMO/LIM_SRC_2/limistate_2.F90
r823 r888 4 4 !! Initialisation of diagnostics ice variables 5 5 !!====================================================================== 6 !! History : 2.0 ! 01-04 (C. Ethe, G. Madec) Original code 6 !! History : 1.0 ! 01-04 (C. Ethe, G. Madec) Original code 7 !! 2.0 ! 03-08 (G. Madec) add lim_istate_init 7 8 !! ! 04-04 (S. Theetten) initialization from a file 8 9 !! ! 06-07 (S. Masson) IOM to read the restart 10 !! ! 07-10 (G. Madec) surface module 9 11 !!-------------------------------------------------------------------- 10 12 #if defined key_lim2 … … 18 20 USE phycst 19 21 USE ocfzpt 20 USE oce ! dynamics and tracers variables !!gm used???21 USE dom_oce !!gm used???22 22 USE par_ice_2 ! ice parameters 23 23 USE ice_oce ! ice variables 24 24 USE dom_ice_2 25 25 USE lbclnk 26 USE oce 26 27 USE ice_2 27 28 USE iom … … 47 48 !!---------------------------------------------------------------------- 48 49 !! LIM 2.0, UCL-LOCEAN-IPSL (2006) 49 !! $ Header$50 !! $ Id: $ 50 51 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 51 52 !!---------------------------------------------------------------------- … … 66 67 REAL(wp), DIMENSION(jpi,jpj) :: ztn ! workspace 67 68 !-------------------------------------------------------------------- 68 69 CALL lim_istate_init_2 ! reading the initials parameters of the ice 70 71 !-- Initialisation of sst,sss,u,v do i=1,jpi 72 u_io(:,:) = 0.e0 ! ice velocity in x direction 73 v_io(:,:) = 0.e0 ! ice velocity in y direction 74 75 IF( ln_limini ) THEN ! 76 77 ! Initialisation at tn if no ice or sst_ini if ice 78 ! Idem for salinity 79 80 !--- Criterion for presence (zidto=1.) or absence (zidto=0.) of ice 81 DO jj = 1 , jpj 82 DO ji = 1 , jpi 83 84 zidto = MAX(zzero, - SIGN(1.,frld(ji,jj) - 1.)) 85 86 sst_io(ji,jj) = ( nfice - 1 ) * (zidto * sst_ini(ji,jj) + & ! use the ocean initial values 87 & (1.0 - zidto ) * ( tn(ji,jj,1) + rt0 )) ! tricky trick *(nfice-1) ! 88 sss_io(ji,jj) = ( nfice - 1 ) * (zidto * sss_ini(ji,jj) + & 89 & (1.0 - zidto ) * sn(ji,jj,1) ) 90 91 ! to avoid the the melting of ice, several layers (mixed layer) should be 92 ! set to sst_ini (sss_ini) if there is ice 93 ! example for one layer 94 ! tn(ji,jj,1) = zidto * ( sst_ini(ji,jj) - rt0 ) + (1.0 - zidto ) * tn(ji,jj,1) 95 ! sn(ji,jj,1) = zidto * sss_ini(ji,jj) + (1.0 - zidto ) * sn(ji,jj,1) 96 ! tb(ji,jj,1) = tn(ji,jj,1) 97 ! sb(ji,jj,1) = sn(ji,jj,1) 98 END DO 99 END DO 100 101 102 ! tfu: Melting point of sea water 103 tfu(:,:) = ztf 104 105 tfu(:,:) = ABS ( rt0 - 0.0575 * sss_ini(:,:) & 106 & + 1.710523e-03 * sss_ini(:,:) * SQRT( sss_ini(:,:) ) & 107 & - 2.154996e-04 * sss_ini(:,:) * sss_ini(:,:) ) 108 ELSE ! 109 69 70 CALL lim_istate_init_2 ! reading the initials parameters of the ice 71 72 IF( .NOT. ln_limini ) THEN 110 73 111 74 ! Initialisation at tn or -2 if ice … … 116 79 END DO 117 80 END DO 118 119 u_io (:,:) = 0.e0 120 v_io (:,:) = 0.e0 121 sst_io(:,:) = ( nfice - 1 ) * ( tn(:,:,1) + rt0 ) ! use the ocean initial values 122 sss_io(:,:) = ( nfice - 1 ) * sn(:,:,1) ! tricky trick *(nfice-1) ! 123 124 ! reference salinity 34psu 81 82 ! tfu: Melting point of sea water [Kelvin] 125 83 zs0 = 34.e0 126 ztf = ABS ( rt0 - 0.0575 * zs0 & 127 & + 1.710523e-03 * zs0 * SQRT( zs0 ) & 128 & - 2.154996e-04 * zs0 *zs0 ) 129 130 ! tfu: Melting point of sea water 131 tfu(:,:) = ztf 84 ztf = rt0 + ( - 0.0575 + 1.710523e-3 * SQRT( zs0 ) - 2.154996e-4 * zs0 ) * zs0 85 tfu(:,:) = ztf 132 86 133 87 DO jj = 1, jpj … … 152 106 tbif (:,:,2) = tfu(:,:) 153 107 tbif (:,:,3) = tfu(:,:) 154 108 155 109 ENDIF 110 156 111 fsbbq (:,:) = 0.e0 157 112 qstoif(:,:) = 0.e0 158 u _ice(:,:) = 0.e0159 v _ice(:,:) = 0.e0113 ui_ice(:,:) = 0.e0 114 vi_ice(:,:) = 0.e0 160 115 # if defined key_coupled 161 116 albege(:,:) = 0.8 * tms(:,:) … … 191 146 192 147 CALL lbc_lnk( hsnif, 'T', 1. ) 193 CALL lbc_lnk( sist , 'T', 1. )148 CALL lbc_lnk( sist , 'T', 1. , pval = rt0 ) ! set rt0 on closed boundary (required by bulk formulation) 194 149 DO jk = 1, jplayersp1 195 150 CALL lbc_lnk(tbif(:,:,jk), 'T', 1. ) … … 197 152 CALL lbc_lnk( fsbbq , 'T', 1. ) 198 153 CALL lbc_lnk( qstoif , 'T', 1. ) 199 CALL lbc_lnk( sss_io , 'T', 1. ) 200 ! 154 201 155 END SUBROUTINE lim_istate_2 202 156 … … 209 163 !! 210 164 !! ** Method : Read the namiceini namelist and check the parameter 211 !! values called at the first timestep (nit000) 212 !! or 213 !! Read 7 variables from a previous restart file 214 !! sst, sst, hicif, hsnif, frld, ts & tbif 165 !! values called at the first timestep (nit000) 215 166 !! 216 167 !! ** input : Namelist namiceini … … 222 173 & hnins, hgins, alins 223 174 !!------------------------------------------------------------------- 224 225 ! Read Namelist namiceini 226 REWIND ( numnam_ice ) 175 ! 176 REWIND ( numnam_ice ) ! Read Namelist namiceini 227 177 READ ( numnam_ice , namiceini ) 228 229 IF(.NOT. ln_limini) THEN 230 IF(lwp) THEN 231 WRITE(numout,*) 232 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 233 WRITE(numout,*) '~~~~~~~~~~~~~~~' 234 WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest 235 WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn 236 WRITE(numout,*) ' initial ice thickness in the north hginn = ', hginn 237 WRITE(numout,*) ' initial leads area in the north alinn = ', alinn 238 WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins 239 WRITE(numout,*) ' initial ice thickness in the south hgins = ', hgins 240 WRITE(numout,*) ' initial leads area in the south alins = ', alins 241 ENDIF 178 ! 179 IF(lwp) THEN 180 WRITE(numout,*) 181 WRITE(numout,*) 'lim_istate_init_2 : ice parameters inititialisation ' 182 WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 183 WRITE(numout,*) ' threshold water temp. for initial sea-ice ttest = ', ttest 184 WRITE(numout,*) ' initial snow thickness in the north hninn = ', hninn 185 WRITE(numout,*) ' initial ice thickness in the north hginn = ', hginn 186 WRITE(numout,*) ' initial leads area in the north alinn = ', alinn 187 WRITE(numout,*) ' initial snow thickness in the south hnins = ', hnins 188 WRITE(numout,*) ' initial ice thickness in the south hgins = ', hgins 189 WRITE(numout,*) ' initial leads area in the south alins = ', alins 190 WRITE(numout,*) ' Ice state initialization using input file ln_limini = ', ln_limini 242 191 ENDIF 243 192 244 193 IF( ln_limini ) THEN ! Ice initialization using input file 245 194 ! 246 195 CALL iom_open( 'Ice_initialization.nc', inum_ice ) 247 196 ! 248 197 IF( inum_ice > 0 ) THEN 249 IF(lwp) THEN 250 WRITE(numout,*) ' ' 251 WRITE(numout,*) 'lim_istate_init : ice state initialization with : Ice_initialization.nc' 252 WRITE(numout,*) '~~~~~~~~~~~~~~~' 253 WRITE(numout,*) ' Ice state initialization using input file ln_limini = ', ln_limini 254 WRITE(numout,*) ' ' 255 ENDIF 198 IF(lwp) WRITE(numout,*) 199 IF(lwp) WRITE(numout,*) ' ice state initialization with : Ice_initialization.nc' 256 200 257 CALL iom_get( inum_ice, jpdom_data, 'sst' , sst_ini(:,:) ) 258 CALL iom_get( inum_ice, jpdom_data, 'sss' , sss_ini(:,:) ) 259 CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif (:,:) ) 260 CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif (:,:) ) 261 CALL iom_get( inum_ice, jpdom_data, 'frld' , frld (:,:) ) 262 CALL iom_get( inum_ice, jpdom_data, 'ts' , sist (:,:) ) 201 CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif ) 202 CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif ) 203 CALL iom_get( inum_ice, jpdom_data, 'frld' , frld ) 204 CALL iom_get( inum_ice, jpdom_data, 'ts' , sist ) 263 205 CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif(1:nlci,1:nlcj,:), & 264 206 & kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) ) … … 268 210 269 211 CALL iom_close( inum_ice) 270 212 ! 271 213 ENDIF 272 214 ENDIF 273 ! 215 ! 274 216 END SUBROUTINE lim_istate_init_2 275 217 -
trunk/NEMO/LIM_SRC_2/limmsh_2.F90
r823 r888 25 25 !!---------------------------------------------------------------------- 26 26 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 27 !! $ Header$27 !! $ Id: $ 28 28 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 29 29 !!---------------------------------------------------------------------- -
trunk/NEMO/LIM_SRC_2/limrhg_2.F90
r823 r888 4 4 !! Ice rheology : performs sea ice rheology 5 5 !!====================================================================== 6 !! History : 0.0 ! 93-12 (M.A. Morales Maqueda.) Original code 7 !! 1.0 ! 94-12 (H. Goosse) 8 !! 2.0 ! 03-12 (C. Ethe, G. Madec) F90, mpp 9 !! " " ! 06-08 (G. Madec) surface module, ice-stress at I-point 10 !! " " ! 09-09 (G. Madec) Huge verctor optimisation 11 !!---------------------------------------------------------------------- 6 12 #if defined key_lim2 7 13 !!---------------------------------------------------------------------- 8 14 !! 'key_lim2' LIM 2.0 sea-ice model 9 15 !!---------------------------------------------------------------------- 16 !!---------------------------------------------------------------------- 10 17 !! lim_rhg_2 : computes ice velocities 11 18 !!---------------------------------------------------------------------- 12 !! * Modules used13 USE phycst14 USE par_oce15 USE ice_oce !ice variables16 USE dom_ice_217 USE ice_2 18 USE lbclnk 19 USE lib_mpp 20 USE in_out_manager 21 USE prtctl 19 USE par_oce ! ocean parameter 20 USE ice_oce ! ice variables 21 USE sbc_ice ! surface boundary condition: ice variables 22 USE dom_ice_2 ! domaine: ice variables 23 USE phycst ! physical constant 24 USE ice_2 ! ice variables 25 USE lbclnk ! lateral boundary condition 26 USE lib_mpp ! MPP library 27 USE in_out_manager ! I/O manager 28 USE prtctl ! Print control 22 29 23 30 IMPLICIT NONE 24 31 PRIVATE 25 32 26 !! * Routine accessibility27 PUBLIC lim_rhg_2 ! routine called by lim_dyn_2 28 29 !! * Module variables30 REAL(wp) :: & ! constant values 31 rzero = 0.e0 , &32 rone = 1.e0 33 !!---------------------------------------------------------------------- 34 !! LIM 2.0, UCL-LOCEAN-IPSL (200 5)35 !! $ Header$36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt33 PUBLIC lim_rhg_2 ! routine called by lim_dyn 34 35 REAL(wp) :: rzero = 0.e0 ! constant value: zero 36 REAL(wp) :: rone = 1.e0 ! and one 37 38 !! * Substitutions 39 # include "vectopt_loop_substitute.h90" 40 !!---------------------------------------------------------------------- 41 !! LIM 2.0, UCL-LOCEAN-IPSL (2006) 42 !! $ Id: $ 43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 37 44 !!---------------------------------------------------------------------- 38 45 … … 48 55 !! viscous-plastic law including shear strength and a bulk rheology. 49 56 !! 50 !! ** Action : - compute u_ice, v_ice the sea-ice velocity 57 !! ** Action : - compute ui_ice, vi_ice the sea-ice velocity defined 58 !! at I-point 59 !!------------------------------------------------------------------- 60 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation 61 INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation 51 62 !! 52 !! History : 53 !! 0.0 ! 93-12 (M.A. Morales Maqueda.) Original code 54 !! 1.0 ! 94-12 (H. Goosse) 55 !! 2.0 ! 03-12 (C. Ethe, G. Madec) F90, mpp 63 INTEGER :: ji, jj ! dummy loop indices 64 INTEGER :: iter, jter ! temporary integers 65 CHARACTER (len=50) :: charout 66 REAL(wp) :: ze11 , ze12 , ze22 , ze21 ! temporary scalars 67 REAL(wp) :: zt11 , zt12 , zt21 , zt22 ! " " 68 REAL(wp) :: zvis11, zvis21, zvis12, zvis22 ! " " 69 REAL(wp) :: zgphsx, ztagnx, zunw, zur, zusw ! " " 70 REAL(wp) :: zgphsy, ztagny, zvnw, zvr ! " " 71 REAL(wp) :: zresm, za, zac, zmod 72 REAL(wp) :: zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal 73 REAL(wp) :: ztrace2, zdeter, zdelta, zmask, zdgp, zdgi, zdiag 74 REAL(wp) :: za1, zb1, zc1, zd1 75 REAL(wp) :: za2, zb2, zc2, zd2, zden 76 REAL(wp) :: zs11_11, zs11_12, zs11_21, zs11_22 77 REAL(wp) :: zs12_11, zs12_12, zs12_21, zs12_22 78 REAL(wp) :: zs21_11, zs21_12, zs21_21, zs21_22 79 REAL(wp) :: zs22_11, zs22_12, zs22_21, zs22_22 80 REAL(wp), DIMENSION(jpi, jpj ) :: zfrld, zmass, zcorl 81 REAL(wp), DIMENSION(jpi, jpj ) :: za1ct, za2ct, zresr 82 REAL(wp), DIMENSION(jpi, jpj ) :: zc1u, zc1v, zc2u, zc2v 83 REAL(wp), DIMENSION(jpi, jpj ) :: zsang 84 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zu0, zv0 85 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zu_n, zv_n 86 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zu_a, zv_a 87 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zviszeta, zviseta 88 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zzfrld, zztms 89 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zi1, zi2, zmasst, zpresh 90 56 91 !!------------------------------------------------------------------- 57 ! * Arguments 58 INTEGER, INTENT(in) :: k_j1 , & ! southern j-index for ice computation 59 & k_jpj ! northern j-index for ice computation 60 61 ! * Local variables 62 INTEGER :: ji, jj ! dummy loop indices 63 64 INTEGER :: & 65 iim1, ijm1, iip1 , ijp1 , & ! temporary integers 66 iter, jter ! " " 67 68 CHARACTER (len=50) :: charout 69 70 REAL(wp) :: & 71 ze11 , ze12 , ze22 , ze21 , & ! temporary scalars 72 zt11 , zt12 , zt21 , zt22 , & ! " " 73 zvis11, zvis21, zvis12, zvis22, & ! " " 74 zgphsx, ztagnx, zusw , & ! " " 75 zgphsy, ztagny ! " " 76 REAL(wp) :: & 77 zresm, zunw, zvnw, zur, zvr, zmod, za, zac, & 78 zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal, & 79 ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag 80 REAL(wp),DIMENSION(jpi,jpj) :: & 81 zpresh, zfrld, zmass, zcorl, & 82 zu0, zv0, zviszeta, zviseta, & 83 zc1u, zc1v, zc2u, zc2v, za1ct, za2ct, za1, za2, zb1, zb2, & 84 zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr 85 REAL(wp),DIMENSION(jpi,jpj,2,2) :: & 86 zs11, zs12, zs22, zs21 87 !!------------------------------------------------------------------- 92 93 !!bug 94 !! ui_oce(:,:) = 0.e0 95 !! vi_oce(:,:) = 0.e0 96 !! write(*,*) 'rhg min, max u & v', maxval(ui_oce), minval(ui_oce), maxval(vi_oce), minval(vi_oce) 97 !!bug 88 98 89 99 ! Store initial velocities 90 ! ------------------------ 91 zu0(:,:) = u_ice(:,:) 92 zv0(:,:) = v_ice(:,:) 100 ! ---------------- 101 zztms(:,0 ) = 0.e0 ; zzfrld(:,0 ) = 0.e0 102 zztms(:,jpj+1) = 0.e0 ; zzfrld(:,jpj+1) = 0.e0 103 zu0(:,0 ) = 0.e0 ; zv0(:,0 ) = 0.e0 104 zu0(:,jpj+1) = 0.e0 ; zv0(:,jpj+1) = 0.e0 105 zztms(:,1:jpj) = tms(:,:) ; zzfrld(:,1:jpj) = frld(:,:) 106 zu0(:,1:jpj) = ui_ice(:,:) ; zv0(:,1:jpj) = vi_ice(:,:) 107 108 zu_a(:,:) = zu0(:,:) ; zv_a(:,:) = zv0(:,:) 109 zu_n(:,:) = zu0(:,:) ; zv_n(:,:) = zv0(:,:) 110 111 !i 112 zi1 (:,:) = 0.e0 113 zi2 (:,:) = 0.e0 114 zpresh(:,:) = 0.e0 115 zmasst(:,:) = 0.e0 116 !i 117 !!gm violant 118 zfrld(:,:) =0.e0 119 zcorl(:,:) =0.e0 120 zmass(:,:) =0.e0 121 za1ct(:,:) =0.e0 122 za2ct(:,:) =0.e0 123 !!gm end 124 125 zviszeta(:,:) = 0.e0 126 zviseta (:,:) = 0.e0 127 128 !i zviszeta(:,0 ) = 0.e0 ; zviseta(:,0 ) = 0.e0 129 !i zviszeta(:,jpj ) = 0.e0 ; zviseta(:,jpj ) = 0.e0 130 !i zviszeta(:,jpj+1) = 0.e0 ; zviseta(:,jpj+1) = 0.e0 131 93 132 94 133 ! Ice mass, ice strength, and wind stress at the center | … … 96 135 !------------------------------------------------------------------- 97 136 137 !CDIR NOVERRCHK 98 138 DO jj = k_j1 , k_jpj-1 139 !CDIR NOVERRCHK 99 140 DO ji = 1 , jpi 100 za1(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 141 ! only the sinus changes its sign with the hemisphere 142 zsang(ji,jj) = SIGN( 1.e0, fcor(ji,jj) ) * sangvg ! only the sinus changes its sign with the hemisphere 143 ! 144 zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 101 145 zpresh(ji,jj) = tms(ji,jj) * pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) ) 102 #if defined key_lim_cp1 && defined key_coupled 103 zb1(ji,jj) = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) ) 104 zb2(ji,jj) = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) ) 105 #else 106 zb1(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 107 zb2(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 108 #endif 146 !!gm :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld) 147 zi1(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 148 zi2(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 109 149 END DO 110 150 END DO … … 117 157 118 158 DO jj = k_j1+1, k_jpj-1 119 DO ji = 2, jpi120 zstms = tms(ji,jj ) * wght(ji,jj,2,2) +tms(ji-1,jj ) * wght(ji,jj,1,2) &121 & + tms(ji,jj-1) * wght(ji,jj,2,1) +tms(ji-1,jj-1) * wght(ji,jj,1,1)159 DO ji = fs_2, jpi 160 zstms = zztms(ji,jj ) * wght(ji,jj,2,2) + zztms(ji-1,jj ) * wght(ji,jj,1,2) & 161 & + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1) 122 162 zusw = 1.0 / MAX( zstms, epsd ) 123 163 124 zt11 = tms(ji ,jj ) *frld(ji ,jj )125 zt12 = tms(ji-1,jj ) *frld(ji-1,jj )126 zt21 = tms(ji ,jj-1) *frld(ji ,jj-1)127 zt22 = tms(ji-1,jj-1) *frld(ji-1,jj-1)164 zt11 = zztms(ji ,jj ) * zzfrld(ji ,jj ) 165 zt12 = zztms(ji-1,jj ) * zzfrld(ji-1,jj ) 166 zt21 = zztms(ji ,jj-1) * zzfrld(ji ,jj-1) 167 zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1) 128 168 129 169 ! Leads area. … … 131 171 & + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw 132 172 133 ! Mass and coriolis coeff. 134 zmass(ji,jj) = ( z a1(ji,jj ) * wght(ji,jj,2,2) + za1(ji-1,jj ) * wght(ji,jj,1,2) &135 & + z a1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw173 ! Mass and coriolis coeff. at I-point 174 zmass(ji,jj) = ( zmasst(ji,jj ) * wght(ji,jj,2,2) + zmasst(ji-1,jj ) * wght(ji,jj,1,2) & 175 & + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 136 176 zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj) 137 177 138 178 ! Wind stress. 139 #if defined key_lim_cp1 && defined key_coupled 140 ztagnx = ( zb1(ji,jj ) * wght(ji,jj,2,2) + zb1(ji-1,jj ) * wght(ji,jj,1,2) & 141 & + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 142 ztagny = ( zb2(ji,jj ) * wght(ji,jj,2,2) + zb2(ji-1,jj ) * wght(ji,jj,1,2) & 143 & + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 144 #else 145 ztagnx = ( zb1(ji,jj ) * wght(ji,jj,2,2) + zb1(ji-1,jj ) * wght(ji,jj,1,2) & 146 & + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj) 147 ztagny = ( zb2(ji,jj ) * wght(ji,jj,2,2) + zb2(ji-1,jj ) * wght(ji,jj,1,2) & 148 & + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj) 149 #endif 179 ! always provide stress at I-point (ocean F-point) 180 ztagnx = ( zi1(ji,jj ) * wght(ji,jj,2,2) + zi1(ji-1,jj ) * wght(ji,jj,1,2) & 181 & + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utaui_ice(ji,jj) 182 ztagny = ( zi2(ji,jj ) * wght(ji,jj,2,2) + zi2(ji-1,jj ) * wght(ji,jj,1,2) & 183 & + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtaui_ice(ji,jj) 150 184 151 185 ! Gradient of ice strength … … 161 195 162 196 ! Computation of the velocity field taking into account the ice-ice interaction. 163 ! Terms that are independent of the velocity field.164 za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v _oce(ji,jj) - zgphsx165 za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u _oce(ji,jj) - zgphsy197 ! Terms that are independent of the ice velocity field. 198 za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * vi_oce(ji,jj) - zgphsx 199 za2ct(ji,jj) = ztagny + zcorl(ji,jj) * ui_oce(ji,jj) - zgphsy 166 200 END DO 167 201 END DO 168 169 !! inutile!!170 !!?? CALL lbc_lnk( za1ct, 'I', -1. )171 !!?? CALL lbc_lnk( za2ct, 'I', -1. )172 202 173 203 … … 182 212 ! Computation of free drift field for free slip boundary conditions. 183 213 184 DO jj = k_j1, k_jpj-1 185 DO ji = 1, jpim1 186 !- Rate of strain tensor. 187 zt11 = akappa(ji,jj,1,1) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) - u_ice(ji,jj ) - u_ice(ji ,jj+1) ) & 188 & + akappa(ji,jj,1,2) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) + v_ice(ji,jj ) + v_ice(ji ,jj+1) ) 189 zt12 = - akappa(ji,jj,2,2) * ( u_ice(ji ,jj) + u_ice(ji+1,jj ) - u_ice(ji,jj+1) - u_ice(ji+1,jj+1) ) & 190 & - akappa(ji,jj,2,1) * ( v_ice(ji ,jj) + v_ice(ji+1,jj ) + v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) 191 zt22 = - akappa(ji,jj,2,2) * ( v_ice(ji ,jj) + v_ice(ji+1,jj ) - v_ice(ji,jj+1) - v_ice(ji+1,jj+1) ) & 192 & + akappa(ji,jj,2,1) * ( u_ice(ji ,jj) + u_ice(ji+1,jj ) + u_ice(ji,jj+1) + u_ice(ji+1,jj+1) ) 193 zt21 = akappa(ji,jj,1,1) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) - v_ice(ji,jj ) - v_ice(ji ,jj+1) ) & 194 & - akappa(ji,jj,1,2) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) + u_ice(ji,jj ) + u_ice(ji ,jj+1) ) 195 196 !- Rate of strain tensor. 197 zdgp = zt11 + zt22 198 zdgi = zt12 + zt21 199 ztrace2 = zdgp * zdgp 200 zdeter = zt11 * zt22 - 0.25 * zdgi * zdgi 201 202 ! Creep limit depends on the size of the grid. 203 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2), creepl) 204 205 !- Computation of viscosities. 206 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 207 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 208 END DO 209 END DO 210 !!?? CALL lbc_lnk( zviszeta, 'I', -1. ) ! or T point??? semble reellement inutile 211 !!?? CALL lbc_lnk( zviseta , 'I', -1. ) 212 213 214 !- Determination of zc1u, zc2u, zc1v and zc2v. 215 DO jj = k_j1+1, k_jpj-1 216 DO ji = 2, jpim1 217 ze11 = akappa(ji-1,jj-1,1,1) 218 ze12 = +akappa(ji-1,jj-1,2,2) 219 ze22 = akappa(ji-1,jj-1,2,1) 220 ze21 = -akappa(ji-1,jj-1,1,2) 221 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 222 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 223 zvis12 = zviseta (ji-1,jj-1) + dm 224 zvis21 = zviseta (ji-1,jj-1) 225 226 zdiag = zvis22 * ( ze11 + ze22 ) 227 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 228 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 229 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 230 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 231 232 ze11 = -akappa(ji,jj-1,1,1) 233 ze12 = +akappa(ji,jj-1,2,2) 234 ze22 = akappa(ji,jj-1,2,1) 235 ze21 = -akappa(ji,jj-1,1,2) 236 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 237 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 238 zvis12 = zviseta (ji,jj-1) + dm 239 zvis21 = zviseta (ji,jj-1) 240 241 zdiag = zvis22 * ( ze11 + ze22 ) 242 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 243 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 244 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 245 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 246 247 ze11 = akappa(ji-1,jj,1,1) 248 ze12 = -akappa(ji-1,jj,2,2) 249 ze22 = akappa(ji-1,jj,2,1) 250 ze21 = -akappa(ji-1,jj,1,2) 251 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 252 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 253 zvis12 = zviseta (ji-1,jj) + dm 254 zvis21 = zviseta (ji-1,jj) 255 256 zdiag = zvis22 * ( ze11 + ze22 ) 257 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 258 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 259 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 260 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 261 262 ze11 = -akappa(ji,jj,1,1) 263 ze12 = -akappa(ji,jj,2,2) 264 ze22 = akappa(ji,jj,2,1) 265 ze21 = -akappa(ji,jj,1,2) 266 zvis11 = 2.0 * zviseta (ji,jj) + dm 267 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 268 zvis12 = zviseta (ji,jj) + dm 269 zvis21 = zviseta (ji,jj) 270 271 zdiag = zvis22 * ( ze11 + ze22 ) 272 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 273 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 274 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 275 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 276 END DO 277 END DO 278 279 DO jj = k_j1+1, k_jpj-1 280 DO ji = 2, jpim1 281 zc1u(ji,jj) = & 282 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 283 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 284 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 285 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 286 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 287 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 288 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 289 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 290 291 zc2u(ji,jj) = & 292 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 293 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 294 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 295 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 296 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 297 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 298 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 299 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 300 END DO 301 END DO 302 303 DO jj = k_j1+1, k_jpj-1 304 DO ji = 2, jpim1 305 ! zc1v , zc2v. 306 ze11 = akappa(ji-1,jj-1,1,2) 307 ze12 = -akappa(ji-1,jj-1,2,1) 308 ze22 = +akappa(ji-1,jj-1,2,2) 309 ze21 = akappa(ji-1,jj-1,1,1) 310 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 311 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 312 zvis12 = zviseta (ji-1,jj-1) + dm 313 zvis21 = zviseta (ji-1,jj-1) 314 315 zdiag = zvis22 * ( ze11 + ze22 ) 316 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 317 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 318 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 319 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 214 !CDIR NOVERRCHK 215 DO jj = k_j1, k_jpj-1 216 !CDIR NOVERRCHK 217 DO ji = 1, fs_jpim1 218 !- Rate of strain tensor. 219 zt11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji,jj ) - zu_a(ji ,jj+1) ) & 220 & + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji,jj ) + zv_a(ji ,jj+1) ) 221 zt12 = - akappa(ji,jj,2,2) * ( zu_a(ji ,jj) + zu_a(ji+1,jj ) - zu_a(ji,jj+1) - zu_a(ji+1,jj+1) ) & 222 & - akappa(ji,jj,2,1) * ( zv_a(ji ,jj) + zv_a(ji+1,jj ) + zv_a(ji,jj+1) + zv_a(ji+1,jj+1) ) 223 zt22 = - akappa(ji,jj,2,2) * ( zv_a(ji ,jj) + zv_a(ji+1,jj ) - zv_a(ji,jj+1) - zv_a(ji+1,jj+1) ) & 224 & + akappa(ji,jj,2,1) * ( zu_a(ji ,jj) + zu_a(ji+1,jj ) + zu_a(ji,jj+1) + zu_a(ji+1,jj+1) ) 225 zt21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji,jj ) - zv_a(ji ,jj+1) ) & 226 & - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji,jj ) + zu_a(ji ,jj+1) ) 227 228 !- Rate of strain tensor. 229 zdgp = zt11 + zt22 230 zdgi = zt12 + zt21 231 ztrace2 = zdgp * zdgp 232 zdeter = zt11 * zt22 - 0.25 * zdgi * zdgi 233 234 ! Creep limit depends on the size of the grid. 235 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ), creepl) 236 237 !- Computation of viscosities. 238 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 239 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 240 END DO 241 END DO 242 243 !- Determination of zc1u, zc2u, zc1v and zc2v. 244 DO jj = k_j1+1, k_jpj-1 245 DO ji = fs_2, fs_jpim1 246 !* zc1u , zc2v 247 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 248 zvis12 = zviseta (ji-1,jj-1) + dm 249 zvis21 = zviseta (ji-1,jj-1) 250 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 251 zdiag = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) ) 252 zs11_11 = zvis11 * akappa(ji-1,jj-1,1,1) + zdiag 253 zs12_11 = zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2) 254 zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2) 255 zs22_11 = zvis11 * akappa(ji-1,jj-1,2,1) + zdiag 256 257 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 258 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 259 zvis12 = zviseta (ji,jj-1) + dm 260 zvis21 = zviseta (ji,jj-1) 261 zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) ) 262 zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag 263 zs12_21 = zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2) 264 zs22_21 = zvis11 * akappa(ji,jj-1,2,1) + zdiag 265 zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2) 266 267 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 268 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 269 zvis12 = zviseta (ji-1,jj) + dm 270 zvis21 = zviseta (ji-1,jj) 271 zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) ) 272 zs11_12 = zvis11 * akappa(ji-1,jj,1,1) + zdiag 273 zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2) 274 zs22_12 = zvis11 * akappa(ji-1,jj,2,1) + zdiag 275 zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2) 276 277 zvis11 = 2.0 * zviseta (ji,jj) + dm 278 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 279 zvis12 = zviseta (ji,jj) + dm 280 zvis21 = zviseta (ji,jj) 281 zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) ) 282 zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag 283 zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2) 284 zs22_22 = zvis11 * akappa(ji,jj,2,1) + zdiag 285 zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2) 286 287 zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 288 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 289 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 290 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 291 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 292 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 293 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 294 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 295 296 zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 297 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 298 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 299 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 300 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 301 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 302 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 303 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 304 305 !* zc1v , zc2v. 306 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 307 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 308 zvis12 = zviseta (ji-1,jj-1) + dm 309 zvis21 = zviseta (ji-1,jj-1) 310 zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) ) 311 zs11_11 = zvis11 * akappa(ji-1,jj-1,1,2) + zdiag 312 zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1) 313 zs22_11 = zvis11 * akappa(ji-1,jj-1,2,2) + zdiag 314 zs21_11 = zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1) 320 315 321 ze11 = akappa(ji,jj-1,1,2) 322 ze12 = -akappa(ji,jj-1,2,1) 323 ze22 = +akappa(ji,jj-1,2,2) 324 ze21 = -akappa(ji,jj-1,1,1) 325 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 326 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 327 zvis12 = zviseta (ji,jj-1) + dm 328 zvis21 = zviseta (ji,jj-1) 329 330 zdiag = zvis22 * ( ze11 + ze22 ) 331 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 332 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 333 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 334 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 335 336 ze11 = akappa(ji-1,jj,1,2) 337 ze12 = -akappa(ji-1,jj,2,1) 338 ze22 = -akappa(ji-1,jj,2,2) 339 ze21 = akappa(ji-1,jj,1,1) 340 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 341 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 342 zvis12 = zviseta (ji-1,jj) + dm 343 zvis21 = zviseta (ji-1,jj) 344 345 zdiag = zvis22 * ( ze11 + ze22 ) 346 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 347 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 348 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 349 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 350 351 ze11 = akappa(ji,jj,1,2) 352 ze12 = -akappa(ji,jj,2,1) 353 ze22 = -akappa(ji,jj,2,2) 354 ze21 = -akappa(ji,jj,1,1) 355 zvis11 = 2.0 * zviseta (ji,jj) + dm 356 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 357 zvis12 = zviseta (ji,jj) + dm 358 zvis21 = zviseta (ji,jj) 359 360 zdiag = zvis22 * ( ze11 + ze22 ) 361 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 362 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 363 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 364 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 365 366 END DO 367 END DO 368 369 DO jj = k_j1+1, k_jpj-1 370 DO ji = 2, jpim1 371 zc1v(ji,jj) = & 372 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 373 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 374 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 375 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 376 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 377 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 378 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 379 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 380 zc2v(ji,jj) = & 381 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 382 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 383 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 384 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 385 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 386 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 387 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 388 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 389 END DO 390 END DO 391 392 ! Relaxation. 393 394 iflag: DO jter = 1 , nbitdr 395 396 ! Store previous drift field. 397 DO jj = k_j1, k_jpj-1 398 zu_ice(:,jj) = u_ice(:,jj) 399 zv_ice(:,jj) = v_ice(:,jj) 316 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 317 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 318 zvis12 = zviseta (ji,jj-1) + dm 319 zvis21 = zviseta (ji,jj-1) 320 zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) ) 321 zs11_21 = zvis11 * akappa(ji,jj-1,1,2) + zdiag 322 zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1) 323 zs22_21 = zvis11 * akappa(ji,jj-1,2,2) + zdiag 324 zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1) 325 326 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 327 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 328 zvis12 = zviseta (ji-1,jj) + dm 329 zvis21 = zviseta (ji-1,jj) 330 zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) ) 331 zs11_12 = zvis11 * akappa(ji-1,jj,1,2) + zdiag 332 zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1) 333 zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag 334 zs21_12 = zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1) 335 336 zvis11 = 2.0 * zviseta (ji,jj) + dm 337 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 338 zvis12 = zviseta (ji,jj) + dm 339 zvis21 = zviseta (ji,jj) 340 zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) ) 341 zs11_22 = zvis11 * akappa(ji,jj,1,2) + zdiag 342 zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1) 343 zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag 344 zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1) 345 346 zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 347 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 348 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 349 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 350 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 351 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 352 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 353 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 354 355 zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 356 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 357 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 358 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 359 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 360 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 361 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 362 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 400 363 END DO 401 364 END DO 365 366 ! GAUSS-SEIDEL method 367 ! ! ================ ! 368 iflag: DO jter = 1 , nbitdr ! Relaxation ! 369 ! ! ================ ! 370 !CDIR NOVERRCHK 402 371 DO jj = k_j1+1, k_jpj-1 403 zsang = SIGN( 1.e0, fcor(1,jj) ) * sangvg ! only the sinus changes its sign with the hemisphere 404 DO ji = 2, jpim1 405 zur = u_ice(ji,jj) - u_oce(ji,jj) 406 zvr = v_ice(ji,jj) - v_oce(ji,jj) 407 zmod = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) ) 408 za = rhoco * zmod 409 zac = za * cangvg 410 zmpzas = alpha * zcorl(ji,jj) + za * zsang 372 !CDIR NOVERRCHK 373 DO ji = fs_2, fs_jpim1 374 ! 375 ze11 = akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj) 376 ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj) 377 ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj) 378 ze21 = akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj) 379 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 380 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 381 zvis12 = zviseta (ji,jj-1) + dm 382 zvis21 = zviseta (ji,jj-1) 383 zdiag = zvis22 * ( ze11 + ze22 ) 384 zs11_21 = zvis11 * ze11 + zdiag 385 zs12_21 = zvis12 * ze12 + zvis21 * ze21 386 zs22_21 = zvis11 * ze22 + zdiag 387 zs21_21 = zvis12 * ze21 + zvis21 * ze12 388 389 ze11 = akappa(ji-1,jj,1,1) * ( zu_a(ji ,jj+1) - zu_a(ji-1,jj+1) ) & 390 & + akappa(ji-1,jj,1,2) * ( zv_a(ji ,jj+1) + zv_a(ji-1,jj+1) ) 391 ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) & 392 & - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) 393 ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) & 394 & + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) 395 ze21 = akappa(ji-1,jj,1,1) * ( zv_a(ji ,jj+1) - zv_a(ji-1,jj+1) ) & 396 & - akappa(ji-1,jj,1,2) * ( zu_a(ji ,jj+1) + zu_a(ji-1,jj+1) ) 397 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 398 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 399 zvis12 = zviseta (ji-1,jj) + dm 400 zvis21 = zviseta (ji-1,jj) 401 zdiag = zvis22 * ( ze11 + ze22 ) 402 zs11_12 = zvis11 * ze11 + zdiag 403 zs12_12 = zvis12 * ze12 + zvis21 * ze21 404 zs22_12 = zvis11 * ze22 + zdiag 405 zs21_12 = zvis12 * ze21 + zvis21 * ze12 406 407 ze11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji ,jj+1) ) & 408 & + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji ,jj+1) ) 409 ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji ,jj+1) - zu_a(ji+1,jj+1) ) & 410 & - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji ,jj+1) + zv_a(ji+1,jj+1) ) 411 ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji ,jj+1) - zv_a(ji+1,jj+1) ) & 412 & + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji ,jj+1) + zu_a(ji+1,jj+1) ) 413 ze21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji ,jj+1) ) & 414 & - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji ,jj+1) ) 415 zvis11 = 2.0 * zviseta (ji,jj) + dm 416 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 417 zvis12 = zviseta (ji,jj) + dm 418 zvis21 = zviseta (ji,jj) 419 zdiag = zvis22 * ( ze11 + ze22 ) 420 zs11_22 = zvis11 * ze11 + zdiag 421 zs12_22 = zvis12 * ze12 + zvis21 * ze21 422 zs22_22 = zvis11 * ze22 + zdiag 423 zs21_22 = zvis12 * ze21 + zvis21 * ze12 424 425 ! 2nd part 426 ze11 = akappa(ji-1,jj-1,1,1) * ( zu_a(ji ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) ) & 427 & + akappa(ji-1,jj-1,1,2) * ( zv_a(ji ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) ) 428 ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) - zu_a(ji-1,jj) ) & 429 & - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) + zv_a(ji-1,jj) ) 430 ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) - zv_a(ji-1,jj) ) & 431 & + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) + zu_a(ji-1,jj) ) 432 ze21 = akappa(ji-1,jj-1,1,1) * ( zv_a(ji ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) ) & 433 & - akappa(ji-1,jj-1,1,2) * ( zu_a(ji ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) ) 434 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 435 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 436 zvis12 = zviseta (ji-1,jj-1) + dm 437 zvis21 = zviseta (ji-1,jj-1) 438 zdiag = zvis22 * ( ze11 + ze22 ) 439 zs11_11 = zvis11 * ze11 + zdiag 440 zs12_11 = zvis12 * ze12 + zvis21 * ze21 441 zs22_11 = zvis11 * ze22 + zdiag 442 zs21_11 = zvis12 * ze21 + zvis21 * ze12 443 444 ze11 = akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji ,jj-1) ) & 445 & + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji ,jj-1) ) 446 ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) & 447 & - akappa(ji,jj-1,2,1) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) 448 ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) & 449 & + akappa(ji,jj-1,2,1) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) 450 ze21 = akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji ,jj-1) ) & 451 & - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji ,jj-1) ) 452 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 453 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 454 zvis12 = zviseta (ji,jj-1) + dm 455 zvis21 = zviseta (ji,jj-1) 456 zdiag = zvis22 * ( ze11 + ze22 ) 457 zs11_21 = zs11_21 + zvis11 * ze11 + zdiag 458 zs12_21 = zs12_21 + zvis12 * ze12 + zvis21 * ze21 459 zs22_21 = zs22_21 + zvis11 * ze22 + zdiag 460 zs21_21 = zs21_21 + zvis12 * ze21 + zvis21 * ze12 461 462 ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj) 463 ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj) 464 ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj) 465 ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj) 466 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 467 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 468 zvis12 = zviseta (ji-1,jj) + dm 469 zvis21 = zviseta (ji-1,jj) 470 zdiag = zvis22 * ( ze11 + ze22 ) 471 zs11_12 = zs11_12 + zvis11 * ze11 + zdiag 472 zs12_12 = zs12_12 + zvis12 * ze12 + zvis21 * ze21 473 zs22_12 = zs22_12 + zvis11 * ze22 + zdiag 474 zs21_12 = zs21_12 + zvis12 * ze21 + zvis21 * ze12 475 476 zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 477 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 478 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 479 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 480 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 481 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 482 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 483 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 484 485 zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 486 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 487 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 488 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 489 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 490 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 491 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 492 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 493 494 zur = zu_a(ji,jj) - ui_oce(ji,jj) 495 zvr = zv_a(ji,jj) - vi_oce(ji,jj) 496 !!!! 497 zmod = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 498 za = rhoco * zmod 499 !!!! 500 !!gm chg resul za = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 501 zac = za * cangvg 502 zmpzas = alpha * zcorl(ji,jj) + za * zsang(ji,jj) 411 503 zmassdt = zusdtp * zmass(ji,jj) 412 504 zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj) 413 505 414 za1(ji,jj) = zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj) & 415 & + za * ( cangvg * u_oce(ji,jj) - zsang * v_oce(ji,jj) ) 416 417 za2(ji,jj) = zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj) & 418 & + za * ( cangvg * v_oce(ji,jj) + zsang * u_oce(ji,jj) ) 419 420 zb1(ji,jj) = zmassdt + zac - zc1u(ji,jj) 421 zb2(ji,jj) = zmpzas - zc2u(ji,jj) 422 zc1(ji,jj) = zmpzas + zc1v(ji,jj) 423 zc2(ji,jj) = zmassdt + zac - zc2v(ji,jj) 424 zdeter = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj) 425 zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 506 za1 = zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj) & 507 & + za * ( cangvg * ui_oce(ji,jj) - zsang(ji,jj) * vi_oce(ji,jj) ) 508 za2 = zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj) & 509 & + za * ( cangvg * vi_oce(ji,jj) + zsang(ji,jj) * ui_oce(ji,jj) ) 510 zb1 = zmassdt + zac - zc1u(ji,jj) 511 zb2 = zmpzas - zc2u(ji,jj) 512 zc1 = zmpzas + zc1v(ji,jj) 513 zc2 = zmassdt + zac - zc2v(ji,jj) 514 zdeter = zc1 * zb2 + zc2 * zb1 515 zden = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 516 zunw = ( ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden 517 zvnw = ( ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden 518 zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 519 520 zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask 521 zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask 426 522 END DO 427 523 END DO 428 524 429 ! The computation of ice interaction term is splitted into two parts 430 !------------------------------------------------------------------------- 431 432 ! Terms that do not involve already up-dated velocities. 433 434 DO jj = k_j1+1, k_jpj-1 435 DO ji = 2, jpim1 436 iim1 = ji 437 ijm1 = jj - 1 438 iip1 = ji + 1 439 ijp1 = jj 440 ze11 = akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1) 441 ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1) 442 ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1) 443 ze21 = akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1) 444 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 445 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 446 zvis12 = zviseta (iim1,ijm1) + dm 447 zvis21 = zviseta (iim1,ijm1) 448 zdiag = zvis22 * ( ze11 + ze22 ) 449 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 450 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 451 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 452 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 453 454 455 iim1 = ji - 1 456 ijm1 = jj 457 iip1 = ji 458 ijp1 = jj + 1 459 ze11 = akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) ) & 460 & + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 461 ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) & 462 & - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 463 ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) & 464 & + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 465 ze21 = akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) ) & 466 & - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 467 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 468 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 469 zvis12 = zviseta (iim1,ijm1) + dm 470 zvis21 = zviseta (iim1,ijm1) 471 zdiag = zvis22 * ( ze11 + ze22 ) 472 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 473 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 474 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 475 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 476 477 iim1 = ji 478 ijm1 = jj 479 iip1 = ji + 1 480 ijp1 = jj + 1 481 ze11 = akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) ) & 482 & + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 483 ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) ) & 484 & - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 485 ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) ) & 486 & + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 487 ze21 = akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) ) & 488 & - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 489 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 490 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 491 zvis12 = zviseta (iim1,ijm1) + dm 492 zvis21 = zviseta (iim1,ijm1) 493 494 zdiag = zvis22 * ( ze11 + ze22 ) 495 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 496 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 497 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 498 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 499 500 END DO 525 CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. ) 526 CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. ) 527 528 ! Test of Convergence 529 DO jj = k_j1+1 , k_jpj-1 530 zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) ) 501 531 END DO 502 503 ! Terms involving already up-dated velocities. 504 !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method; 505 ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method. 506 507 DO jj = k_j1+1, k_jpj-1 508 DO ji = 2, jpim1 509 iim1 = ji - 1 510 ijm1 = jj - 1 511 iip1 = ji 512 ijp1 = jj 513 ze11 = akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) ) & 514 & + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) ) 515 ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) ) & 516 & - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) ) 517 ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) ) & 518 & + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) ) 519 ze21 = akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) ) & 520 & - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) ) 521 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 522 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 523 zvis12 = zviseta (iim1,ijm1) + dm 524 zvis21 = zviseta (iim1,ijm1) 525 526 zdiag = zvis22 * ( ze11 + ze22 ) 527 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 528 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 529 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 530 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 531 532 #if defined key_agrif 533 END DO 534 END DO 535 536 DO jj = k_j1+1, k_jpj-1 537 DO ji = 2, jpim1 538 #endif 539 540 iim1 = ji 541 ijm1 = jj - 1 542 iip1 = ji + 1 543 ze11 = akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) ) & 544 & + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) ) 545 ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) ) & 546 & - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) ) 547 ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) ) & 548 & + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) ) 549 ze21 = akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) ) & 550 & - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) ) 551 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 552 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 553 zvis12 = zviseta (iim1,ijm1) + dm 554 zvis21 = zviseta (iim1,ijm1) 555 556 zdiag = zvis22 * ( ze11 + ze22 ) 557 zs11(ji,jj,2,1) = zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag 558 zs12(ji,jj,2,1) = zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21 559 zs22(ji,jj,2,1) = zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag 560 zs21(ji,jj,2,1) = zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12 561 562 563 iim1 = ji - 1 564 ijm1 = jj 565 ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1) 566 ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1) 567 ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1) 568 ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * zu_ice(iim1,ijm1) 569 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 570 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 571 zvis12 = zviseta (iim1,ijm1) + dm 572 zvis21 = zviseta (iim1,ijm1) 573 574 zdiag = zvis22 * ( ze11 + ze22 ) 575 zs11(ji,jj,1,2) = zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag 576 zs12(ji,jj,1,2) = zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21 577 zs22(ji,jj,1,2) = zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag 578 zs21(ji,jj,1,2) = zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12 579 580 #if defined key_agrif 581 END DO 582 END DO 583 584 DO jj = k_j1+1, k_jpj-1 585 DO ji = 2, jpim1 586 #endif 587 zd1(ji,jj) = & 588 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 589 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 590 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 591 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 592 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 593 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 594 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 595 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 596 zd2(ji,jj) = & 597 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 598 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 599 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 600 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 601 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 602 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 603 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 604 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 605 END DO 532 zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) ) 533 !!!! this should be faster on scalar processor 534 ! zresm = MAXVAL( MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ), & 535 ! & ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) ) ) 536 !!!! 537 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 538 539 DO jj = k_j1, k_jpj 540 zu_a(:,jj) = zu_n(:,jj) 541 zv_a(:,jj) = zv_n(:,jj) 606 542 END DO 607 543 608 DO jj = k_j1+1, k_jpj-1 609 DO ji = 2, jpim1 610 zunw = ( ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj) & 611 & + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj) 612 613 zvnw = ( ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj) & 614 & - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj) 615 616 zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 617 618 u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 619 v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 620 END DO 544 IF( zresm <= resl ) EXIT iflag 545 546 ! ! ================ ! 547 END DO iflag ! end Relaxation ! 548 ! ! ================ ! 549 550 IF( zindu == 0 ) THEN ! even iteration 551 DO jj = k_j1 , k_jpj-1 552 zu0(:,jj) = zu_a(:,jj) 553 zv0(:,jj) = zv_a(:,jj) 621 554 END DO 622 623 CALL lbc_lnk( u_ice, 'I', -1. ) 624 CALL lbc_lnk( v_ice, 'I', -1. ) 625 626 !--- 5.2.5.4. Convergence test. 627 DO jj = k_j1+1 , k_jpj-1 628 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 629 END DO 630 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 631 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 632 633 IF ( zresm <= resl) EXIT iflag 634 635 END DO iflag 636 637 zindu1 = 1.0 - zindu 638 DO jj = k_j1 , k_jpj-1 639 zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj) 640 zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj) 641 END DO 642 ! ! ==================== ! 555 ENDIF 556 ! ! ==================== ! 643 557 END DO ! end loop over iter ! 644 558 ! ! ==================== ! 559 560 ui_ice(:,:) = zu_a(:,1:jpj) 561 vi_ice(:,:) = zv_a(:,1:jpj) 645 562 646 563 IF(ln_ctl) THEN 647 564 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter 648 565 CALL prt_ctl_info(charout) 649 CALL prt_ctl(tab2d_1=u _ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')566 CALL prt_ctl(tab2d_1=ui_ice, clinfo1=' lim_rhg : ui_ice :', tab2d_2=vi_ice, clinfo2=' vi_ice :') 650 567 ENDIF 651 568 -
trunk/NEMO/LIM_SRC_2/limrst_2.F90
r823 r888 17 17 !!---------------------------------------------------------------------- 18 18 USE ice_2 19 USE dom_oce20 USE ice_oce ! ice variables19 USE sbc_oce 20 USE sbc_ice 21 21 USE daymod 22 22 … … 27 27 PRIVATE 28 28 29 PUBLIC lim_rst_opn_2 ! routine called by ??? module30 PUBLIC lim_rst_write_2 ! routine called by ??? module31 PUBLIC lim_rst_read_2 ! routine called by ??? module29 PUBLIC lim_rst_opn_2 ! routine called by sbcice_lim_2.F90 30 PUBLIC lim_rst_write_2 ! routine called by sbcice_lim_2.F90 31 PUBLIC lim_rst_read_2 ! routine called by iceini_2.F90 32 32 33 33 LOGICAL, PUBLIC :: lrst_ice !: logical to control the ice restart write … … 36 36 !!---------------------------------------------------------------------- 37 37 !! LIM 2.0, UCL-LOCEAN-IPSL (2006) 38 !! $ Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrst.F90,v 1.15 2007/06/29 14:54:06 opalod Exp $38 !! $ Id: $ 39 39 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 40 !!---------------------------------------------------------------------- … … 57 57 58 58 ! to get better performances with NetCDF format: 59 ! we open and define the ice restart file one ice time step before writing the data (-> at nitrst - 2*n fice+ 1)60 ! except if we write ice restart files every ice time step or if an ice restart file was writen at nitend - 2*n fice+ 161 IF( kt == nitrst - 2*n fice + 1 .OR. nstock == nfice .OR. ( kt == nitend - nfice+ 1 .AND. .NOT. lrst_ice ) ) THEN59 ! we open and define the ice restart file one ice time step before writing the data (-> at nitrst - 2*nn_fsbc + 1) 60 ! except if we write ice restart files every ice time step or if an ice restart file was writen at nitend - 2*nn_fsbc + 1 61 IF( kt == nitrst - 2*nn_fsbc + 1 .OR. nstock == nn_fsbc .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN 62 62 ! beware of the format used to write kt (default is i8.8, that should be large enough...) 63 63 IF( nitrst > 99999999 ) THEN ; WRITE(clkt, * ) nitrst … … 72 72 CASE DEFAULT ; WRITE(numout,*) ' open ice restart NetCDF file: '//clname 73 73 END SELECT 74 IF( kt == nitrst - 2*n fice+ 1 ) THEN75 WRITE(numout,*) ' kt = nitrst - 2*n fice+ 1 = ', kt,' date= ', ndastp76 ELSE ; WRITE(numout,*) ' kt = ' , kt,' date= ', ndastp74 IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN 75 WRITE(numout,*) ' kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 76 ELSE ; WRITE(numout,*) ' kt = ' , kt,' date= ', ndastp 77 77 ENDIF 78 78 ENDIF … … 90 90 !! ** purpose : output of sea-ice variable in a netcdf file 91 91 !!---------------------------------------------------------------------- 92 INTEGER, INTENT(in) :: kt 93 ! !94 INTEGER :: iter ! kt + nfice-195 !!---------------------------------------------------------------------- 96 97 iter = kt + n fice - 1 ! ice restarts are written at kt == nitrst - nfice+ 192 INTEGER, INTENT(in) :: kt ! number of iteration 93 ! 94 INTEGER :: iter ! kt + nn_fsbc -1 95 !!---------------------------------------------------------------------- 96 97 iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 98 98 99 99 IF( iter == nitrst ) THEN 100 100 IF(lwp) WRITE(numout,*) 101 101 IF(lwp) WRITE(numout,*) 'lim_rst_write_2 : write ice restart file kt =', kt 102 IF(lwp) WRITE(numout,*) '~~~~~~~ '102 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 103 103 ENDIF 104 104 … … 106 106 ! ------------------ 107 107 ! ! calendar control 108 CALL iom_rstput( iter, nitrst, numriw, 'nfice' , REAL( nfice, wp) ) ! time-step 109 CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter , wp) ) ! date 108 CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) ) 110 109 111 110 CALL iom_rstput( iter, nitrst, numriw, 'hicif' , hicif (:,:) ) ! prognostic variables … … 119 118 CALL iom_rstput( iter, nitrst, numriw, 'tbif2' , tbif (:,:,2) ) 120 119 CALL iom_rstput( iter, nitrst, numriw, 'tbif3' , tbif (:,:,3) ) 121 CALL iom_rstput( iter, nitrst, numriw, 'u_ice' , u_ice (:,:) ) 122 CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice (:,:) ) 123 CALL iom_rstput( iter, nitrst, numriw, 'gtaux' , gtaux (:,:) ) 124 CALL iom_rstput( iter, nitrst, numriw, 'gtauy' , gtauy (:,:) ) 120 CALL iom_rstput( iter, nitrst, numriw, 'ui_ice', ui_ice(:,:) ) 121 CALL iom_rstput( iter, nitrst, numriw, 'vi_ice', vi_ice(:,:) ) 125 122 CALL iom_rstput( iter, nitrst, numriw, 'qstoif', qstoif(:,:) ) 126 123 CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq (:,:) ) … … 175 172 !! ** purpose : read of sea-ice variable restart in a netcdf file 176 173 !!---------------------------------------------------------------------- 177 ! 178 REAL(wp) :: zfice, ziter 174 REAL(wp) :: ziter 179 175 !!---------------------------------------------------------------------- 180 176 … … 182 178 WRITE(numout,*) 183 179 WRITE(numout,*) 'lim_rst_read_2 : read ice NetCDF restart file' 184 WRITE(numout,*) '~~~~~~~~ '180 WRITE(numout,*) '~~~~~~~~~~~~~~' 185 181 ENDIF 186 182 187 183 CALL iom_open ( 'restart_ice_in', numrir, kiolib = jprstlib ) 188 184 189 CALL iom_get( numrir, 'nfice' , zfice ) 190 CALL iom_get( numrir, 'kt_ice', ziter ) 191 IF(lwp) WRITE(numout,*) ' read ice restart file at time step : ', ziter 185 CALL iom_get( numrir, 'kt_ice' , ziter ) 186 IF(lwp) WRITE(numout,*) ' read ice restart file at time step : ', INT( ziter ) 192 187 IF(lwp) WRITE(numout,*) ' in any case we force it to nit000 - 1 : ', nit000 - 1 193 188 … … 196 191 IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 ) & 197 192 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart', & 198 & ' verify the file or rerun with the value 0 for the', &199 & ' control of time parameter nrstdt' )200 IF( INT(zfice) /= nfice .AND. ABS( nrstdt ) == 1 ) &201 & CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nfice in ice restart', &202 193 & ' verify the file or rerun with the value 0 for the', & 203 194 & ' control of time parameter nrstdt' ) … … 213 204 CALL iom_get( numrir, jpdom_autoglo, 'tbif2' , tbif(:,:,2) ) 214 205 CALL iom_get( numrir, jpdom_autoglo, 'tbif3' , tbif(:,:,3) ) 215 CALL iom_get( numrir, jpdom_autoglo, 'u_ice' , u_ice ) 216 CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice ) 217 CALL iom_get( numrir, jpdom_autoglo, 'gtaux' , gtaux ) 218 CALL iom_get( numrir, jpdom_autoglo, 'gtauy' , gtauy ) 206 CALL iom_get( numrir, jpdom_autoglo, 'ui_ice', ui_ice ) 207 CALL iom_get( numrir, jpdom_autoglo, 'vi_ice', vi_ice ) 219 208 CALL iom_get( numrir, jpdom_autoglo, 'qstoif', qstoif ) 220 209 CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq ) -
trunk/NEMO/LIM_SRC_2/limtab_2.F90
r823 r888 21 21 !!---------------------------------------------------------------------- 22 22 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 23 !! $ Header$24 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt23 !! $ Id: $ 24 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 25 25 !!---------------------------------------------------------------------- 26 26 CONTAINS -
trunk/NEMO/LIM_SRC_2/limthd_2.F90
r823 r888 4 4 !! LIM thermo ice model : ice thermodynamic 5 5 !!====================================================================== 6 !! History : 1.0 ! 00-01 (LIM) 7 !! 2.0 ! 02-07 (C. Ethe, G. Madec) F90 8 !! 2.0 ! 03-08 (C. Ethe) add lim_thd_init 9 !!--------------------------------------------------------------------- 6 10 #if defined key_lim2 7 11 !!---------------------------------------------------------------------- … … 18 22 USE ice_2 ! LIM sea-ice variables 19 23 USE ice_oce ! sea-ice/ocean variables 20 USE flx_oce ! sea-ice/ocean forcings variables 24 USE sbc_oce ! 25 USE sbc_ice ! 21 26 USE thd_ice_2 ! LIM thermodynamic sea-ice variables 22 27 USE dom_ice_2 ! LIM sea-ice domain … … 30 35 PRIVATE 31 36 32 !! * Routine accessibility 33 PUBLIC lim_thd_2 ! called by lim_step_2 34 35 !! * Module variables 36 REAL(wp) :: & ! constant values 37 epsi20 = 1.e-20 , & 38 epsi16 = 1.e-16 , & 39 epsi04 = 1.e-04 , & 40 zzero = 0.e0 , & 41 zone = 1.e0 37 PUBLIC lim_thd_2 ! called by lim_step 38 39 REAL(wp) :: epsi20 = 1.e-20 , & ! constant values 40 & epsi16 = 1.e-16 , & 41 & epsi04 = 1.e-04 , & 42 & zzero = 0.e0 , & 43 & zone = 1.e0 42 44 43 45 !! * Substitutions … … 46 48 !!-------- ------------------------------------------------------------- 47 49 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 48 !! $ Header$49 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt50 !! $ Id: $ 51 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 50 52 !!---------------------------------------------------------------------- 51 53 … … 68 70 !! - back to the geographic grid 69 71 !! 70 !! ** References : 71 !! H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 72 !! 73 !! History : 74 !! 1.0 ! 00-01 (LIM) 75 !! 2.0 ! 02-07 (C. Ethe, G. Madec) F90 72 !! References : Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90 76 73 !!--------------------------------------------------------------------- 77 74 INTEGER, INTENT(in) :: kt ! number of iteration 78 75 !! 79 76 INTEGER :: ji, jj, & ! dummy loop indices 80 77 nbpb , & ! nb of icy pts for thermo. cal. … … 92 89 zfontn , & ! heat flux from snow thickness 93 90 zfntlat, zpareff ! test. the val. of lead heat budget 94 REAL(wp), DIMENSION(jpi,jpj) :: & 95 zhicifp , & ! ice thickness for outputs 96 zqlbsbq ! link with lead energy budget qldif 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 98 zmsk ! working array 91 REAL(wp), DIMENSION(jpi,jpj) :: zhicifp, & ! ice thickness for outputs 92 & zqlbsbq ! link with lead energy budget qldif 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmsk ! working array 99 94 !!------------------------------------------------------------------- 100 95 101 IF( kt == nit000 96 IF( kt == nit000 ) CALL lim_thd_init_2 ! Initialization (first time-step only) 102 97 103 98 !-------------------------------------------! … … 173 168 !-------------------------------------------------------------------------- 174 169 170 sst_m(:,:) = sst_m(:,:) + rt0 171 175 172 !CDIR NOVERRCHK 176 173 DO jj = 1, jpj … … 188 185 ! temperature and turbulent mixing (McPhee, 1992) 189 186 zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) ! friction velocity 190 fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_ io(ji,jj) - tfu(ji,jj) )187 fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) ) 191 188 qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice 192 189 193 190 ! partial computation of the lead energy budget (qldif) 194 191 zfontn = ( sprecip(ji,jj) / rhosn ) * xlsn ! energy for melting 195 zfnsol = qns r_oce(ji,jj) ! total non solar flux196 qldif(ji,jj) = tms(ji,jj) * ( qsr _oce(ji,jj) * ( 1.0 - thcm(ji,jj) ) &192 zfnsol = qns(ji,jj) ! total non solar flux over the ocean 193 qldif(ji,jj) = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) ) & 197 194 & + zfnsol + fdtcn(ji,jj) - zfontn & 198 195 & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & … … 206 203 207 204 ! energy needed to bring ocean surface layer until its freezing 208 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_ io(ji,jj) ) * ( 1 - zinda )205 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) 209 206 210 207 ! calculate oceanic heat flux. … … 216 213 END DO 217 214 215 sst_m(:,:) = sst_m(:,:) - rt0 218 216 219 217 ! Select icy points and fulfill arrays for the vectorial grid. … … 258 256 CALL tab_2d_1d_2( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0 , jpi, jpj, npb(1:nbpb) ) 259 257 CALL tab_2d_1d_2( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) ) 260 CALL tab_2d_1d_2( nbpb, qns r_ice_1d(1:nbpb) , qnsr_ice, jpi, jpj, npb(1:nbpb) )258 CALL tab_2d_1d_2( nbpb, qns_ice_1d (1:nbpb) , qns_ice , jpi, jpj, npb(1:nbpb) ) 261 259 #if ! defined key_coupled 262 260 CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb) , qla_ice , jpi, jpj, npb(1:nbpb) ) … … 404 402 CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif : ', tab2d_2=fsbbq , clinfo2=' fsbbq : ') 405 403 ENDIF 406 404 ! 407 405 END SUBROUTINE lim_thd_2 408 406 … … 419 417 !! 420 418 !! ** input : Namelist namicether 421 !!422 !! history :423 !! 8.5 ! 03-08 (C. Ethe) original code424 419 !!------------------------------------------------------------------- 425 420 NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax , & -
trunk/NEMO/LIM_SRC_2/limthd_lac_2.F90
r823 r888 30 30 !!---------------------------------------------------------------------- 31 31 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 32 !! $ Header$32 !! $ Id: $ 33 33 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 34 34 !!---------------------------------------------------------------------- -
trunk/NEMO/LIM_SRC_2/limthd_zdf_2.F90
r823 r888 4 4 !! thermodynamic growth and decay of the ice 5 5 !!====================================================================== 6 !! History : 1.0 ! 01-04 (LIM) Original code 7 !! 2.0 ! 02-08 (C. Ethe, G. Madec) F90 8 !!---------------------------------------------------------------------- 6 9 #if defined key_lim2 7 10 !!---------------------------------------------------------------------- 8 11 !! 'key_lim2' LIM 2.0 sea-ice model 12 !!---------------------------------------------------------------------- 9 13 !!---------------------------------------------------------------------- 10 14 !! lim_thd_zdf_2 : vertical accr./abl. and lateral ablation of sea ice … … 22 26 PRIVATE 23 27 24 !! * Routine accessibility 25 PUBLIC lim_thd_zdf_2 ! called by lim_thd_2 26 27 !! * Module variables 28 REAL(wp) :: & ! constant values 29 epsi20 = 1.e-20 , & 30 epsi13 = 1.e-13 , & 31 zzero = 0.e0 , & 32 zone = 1.e0 28 PUBLIC lim_thd_zdf_2 ! called by lim_thd_2 29 30 REAL(wp) :: epsi20 = 1.e-20 , & ! constant values 31 & epsi13 = 1.e-13 , & 32 & zzero = 0.e0 , & 33 & zone = 1.e0 33 34 !!---------------------------------------------------------------------- 34 35 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 35 !! $Header$ 36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 37 !!---------------------------------------------------------------------- 36 !! $ Id: $ 37 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 38 !!---------------------------------------------------------------------- 39 38 40 CONTAINS 39 41 … … 64 66 !! - Performs lateral ablation 65 67 !! 66 !! References : 67 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 68 !! Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268 68 !! References : Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 69 !! Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268 70 !!------------------------------------------------------------------ 71 INTEGER, INTENT(in) :: kideb ! Start point on which the the computation is applied 72 INTEGER, INTENT(in) :: kiut ! End point on which the the computation is applied 69 73 !! 70 !! History :71 !! original : 01-04 (LIM)72 !! addition : 02-08 (C. Ethe, G. Madec)73 !!------------------------------------------------------------------74 !! * Arguments75 INTEGER , INTENT (in) :: &76 kideb , & ! Start point on which the the computation is applied77 kiut ! End point on which the the computation is applied78 79 !! * Local variables80 74 INTEGER :: ji ! dummy loop indices 81 82 REAL(wp) , DIMENSION(jpij,2) :: & 83 zqcmlt ! energy due to surface( /1 ) and bottom melting( /2 ) 84 75 REAL(wp), DIMENSION(jpij,2) :: zqcmlt ! energy due to surface( /1 ) and bottom melting( /2 ) 85 76 REAL(wp), DIMENSION(jpij) :: & 86 77 ztsmlt & ! snow/ice surface melting temperature … … 97 88 , zts_old & ! previous surface temperature 98 89 , zidsn , z1midsn , zidsnic ! tempory variables 99 100 REAL(wp), DIMENSION(jpij) :: & 90 REAL(wp), DIMENSION(jpij) :: & 101 91 zfnet & ! net heat flux at the top surface( incl. conductive heat flux) 102 92 , zsprecip & ! snow accumulation … … 109 99 , zfrld_1d & ! new sea/ice fraction 110 100 , zep ! internal temperature of the 2nd layer of the snow/ice system 111 112 101 REAL(wp), DIMENSION(3) :: & 113 102 zplediag & ! principle diagonal, subdiag. and supdiag. of the … … 115 104 , zsupdiag & ! of the temperatures inside the snow-ice system 116 105 , zsmbr ! second member 117 118 106 REAL(wp) :: & 119 107 zhsu & ! thickness of surface layer … … 130 118 , zb2 , zd2 , zb3 , zd3 & 131 119 , ztint ! equivalent temperature at the snow-ice interface 132 133 120 REAL(wp) :: & 134 121 zexp & ! exponential function of the ice thickness … … 148 135 , zdtic & ! ice internal temp. increment 149 136 , zqnes ! conductive energy due to ice melting in the first ice layer 150 151 137 REAL(wp) :: & 152 138 ztbot & ! temperature at the bottom surface … … 162 148 , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables 163 149 , ztb2, ztb3 164 165 150 REAL(wp) :: & 166 151 zdrmh & ! change in snow/ice thick. after snow-ice formation … … 181 166 ! Computation of energies due to surface and bottom melting 182 167 !----------------------------------------------------------------------- 183 184 168 185 169 DO ji = kideb , kiut … … 201 185 END DO 202 186 203 204 187 !------------------------------------------- 205 188 ! 2. Calculate some intermediate variables. … … 265 248 ! - qstbif_1d, total energy stored in brine pockets (updating) 266 249 !------------------------------------------------------------------- 267 268 250 269 251 DO ji = kideb , kiut … … 288 270 END DO 289 271 290 291 272 !-------------------------------------------------------------------------------- 292 273 ! 4. Computation of the surface temperature : determined by considering the … … 333 314 !---computation of the energy balance function 334 315 zfts = - z1mi0 (ji) * qsr_ice_1d(ji) & ! net absorbed solar radiation 335 & - qns r_ice_1d(ji) & ! total non solar radiation336 & - zfcsu (ji) ! conductive heat flux from the surface316 & - qns_ice_1d(ji) & ! total non solar radiation 317 & - zfcsu (ji) ! conductive heat flux from the surface 337 318 !---computation of surface temperature increment 338 319 zdts = -zfts / zdfts … … 360 341 sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 361 342 #if ! defined key_coupled 362 qns r_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )343 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 363 344 qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 364 345 #endif … … 366 347 END DO 367 348 368 369 370 349 ! 5.2. Calculate available heat for surface ablation. 371 350 !--------------------------------------------------------------------- 372 351 373 352 DO ji = kideb, kiut 374 zfnet(ji) = qns r_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)353 zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji) 375 354 zfnet(ji) = MAX( zzero , zfnet(ji) ) 376 355 zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) ) … … 730 709 dvnbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) 731 710 dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn 732 ! case of natural freshwater flux 733 #if defined key_lim_fdd 734 rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic 711 !--- volume change of ice and snow (used for ocean-ice freshwater flux computation) 712 rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic 735 713 rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhsnnew - h_snow_1d(ji) ) * rhosn 736 714 737 #else738 rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( ( zhicnew - h_ice_1d(ji) ) * rhoic &739 & + ( zhsnnew - h_snow_1d(ji) ) * rhosn )740 #endif741 742 715 !--- Actualize new snow and ice thickness. 743 744 716 h_snow_1d(ji) = zhsnnew 745 h_ice_1d (ji) = zhicnew717 h_ice_1d (ji) = zhicnew 746 718 747 719 END DO … … 799 771 qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji) 800 772 frld_1d(ji) = zfrld_1d(ji) 801 802 END DO 803 773 ! 774 END DO 775 ! 804 776 END SUBROUTINE lim_thd_zdf_2 777 805 778 #else 806 !!====================================================================== 807 !! *** MODULE limthd_zdf_2 *** 808 !! no sea ice model 809 !!====================================================================== 779 !!---------------------------------------------------------------------- 780 !! Default Option NO sea-ice model 781 !!---------------------------------------------------------------------- 810 782 CONTAINS 811 783 SUBROUTINE lim_thd_zdf_2 ! Empty routine 812 784 END SUBROUTINE lim_thd_zdf_2 813 785 #endif 814 END MODULE limthd_zdf_2 786 787 !!====================================================================== 788 END MODULE limthd_zdf_2 -
trunk/NEMO/LIM_SRC_2/limtrp_2.F90
r823 r888 30 30 31 31 !! * Routine accessibility 32 PUBLIC lim_trp_2 ! called by ice_step_232 PUBLIC lim_trp_2 ! called by sbc_ice_lim_2 33 33 34 34 !! * Shared module variables … … 48 48 !!---------------------------------------------------------------------- 49 49 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 50 !! $ Header$50 !! $ Id: $ 51 51 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 52 52 !!---------------------------------------------------------------------- … … 112 112 DO jj = 1, jpjm1 113 113 DO ji = 1, jpim1 114 zui_u(ji,jj) = ( u _ice(ji+1,jj ) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj ) + tmu(ji+1,jj+1), zvbord ) )115 zvi_v(ji,jj) = ( v _ice(ji ,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji ,jj+1) + tmu(ji+1,jj+1), zvbord ) )114 zui_u(ji,jj) = ( ui_ice(ji+1,jj ) + ui_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj ) + tmu(ji+1,jj+1), zvbord ) ) 115 zvi_v(ji,jj) = ( vi_ice(ji ,jj+1) + vi_ice(ji+1,jj+1) ) / ( MAX( tmu(ji ,jj+1) + tmu(ji+1,jj+1), zvbord ) ) 116 116 END DO 117 117 END DO … … 128 128 IF (lk_mpp ) CALL mpp_max(zcfl) 129 129 130 IF ( zcfl > 0.5 .AND. lwp ) WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl130 IF ( zcfl > 0.5 .AND. lwp ) WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl 131 131 132 132 ! content of properties -
trunk/NEMO/LIM_SRC_2/limwri_2.F90
r823 r888 9 9 #if defined key_lim2 10 10 !!---------------------------------------------------------------------- 11 !! 'key_lim2' iLIM 2.0 sea-ice model11 !! 'key_lim2' LIM 2.0 sea-ice model 12 12 !!---------------------------------------------------------------------- 13 13 !!---------------------------------------------------------------------- … … 15 15 !! lim_wri_init_2 : initialization and namelist read 16 16 !!---------------------------------------------------------------------- 17 USE ioipsl18 USE dianam ! build name of file (routine)19 17 USE phycst 20 18 USE dom_oce 21 19 USE daymod 22 USE in_out_manager23 20 USE ice_oce ! ice variables 24 USE flx_oce 21 USE sbc_oce 22 USE sbc_ice 25 23 USE dom_ice_2 26 24 USE ice_2 25 27 26 USE lbclnk 27 USE dianam ! build name of file (routine) 28 USE in_out_manager 29 USE ioipsl 28 30 29 31 IMPLICIT NONE 30 32 PRIVATE 31 33 32 PUBLIC lim_wri_2 ! routine called by lim_step_2.F9034 PUBLIC lim_wri_2 ! routine called by sbc_ice_lim_2 33 35 34 36 INTEGER, PARAMETER :: jpnoumax = 40 ! maximum number of variable for ice output … … 49 51 zone = 1.e0 50 52 51 !!---------------------------------------------------------------------- 52 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 53 !! $Header$ 53 !! * Substitutions 54 # include "vectopt_loop_substitute.h90" 55 !!---------------------------------------------------------------------- 56 !! LIM 2.0, UCL-LOCEAN-IPSL (2006) 57 !! $ Id: $ 54 58 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 55 59 !!---------------------------------------------------------------------- … … 79 83 !!------------------------------------------------------------------- 80 84 INTEGER, INTENT(in) :: kt ! number of iteration 81 85 !! 82 86 INTEGER :: ji, jj, jf ! dummy loop indices 83 87 CHARACTER(len = 40) :: clhstnam, clop … … 90 94 91 95 ! !--------------------! 92 IF ( kt == nit000 ) THEN! Initialisation !96 IF( kt == nit000 ) THEN ! Initialisation ! 93 97 ! !--------------------! 94 98 CALL lim_wri_init_2 … … 97 101 !!Chris clop = "ave(only(x))" !ibug namelist parameter a ajouter 98 102 clop = "ave(x)" 99 zout = nwrite * rdt_ice / n fice103 zout = nwrite * rdt_ice / nn_fsbc 100 104 zsec = 0. 101 105 niter = 0 … … 110 114 111 115 DO jf = 1, noumef 112 IF 113 &, nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )116 IF( nc(jf) == 1 ) CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj & 117 & , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout ) 114 118 END DO 115 119 CALL histend( nice ) 116 120 ! 117 121 ENDIF 118 122 ! !--------------------! … … 120 124 ! !--------------------! 121 125 122 !!gm change the print below to have it only at output time step or when nitend =< 100 123 IF(lwp) THEN 124 WRITE(numout,*) 125 WRITE(numout,*) 'lim_wri_2 : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, kt + nfice - 1 126 WRITE(numout,*) '~~~~~~~~~ ' 127 ENDIF 128 129 !-- calculs des valeurs instantanees 126 !-- Store instantaneous values in zcmo 130 127 131 128 zcmo(:,:, 1:jpnoumax ) = 0.e0 132 129 DO jj = 2 , jpjm1 133 DO ji = 2 ,jpim1130 DO ji = fs_2 , fs_jpim1 134 131 zindh = MAX( zzero , SIGN( zone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) ) 135 132 zinda = MAX( zzero , SIGN( zone , ( 1.0 - frld(ji,jj) ) - 0.10 ) ) … … 142 139 zcmo(ji,jj,5) = sist (ji,jj) 143 140 zcmo(ji,jj,6) = fbif (ji,jj) 144 zcmo(ji,jj,7) = zindb * ( u _ice(ji,jj ) * tmu(ji,jj ) + u_ice(ji+1,jj ) * tmu(ji+1,jj ) &145 + u _ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &141 zcmo(ji,jj,7) = zindb * ( ui_ice(ji,jj ) * tmu(ji,jj ) + ui_ice(ji+1,jj ) * tmu(ji+1,jj ) & 142 + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 146 143 / ztmu 147 144 148 zcmo(ji,jj,8) = zindb * ( v _ice(ji,jj ) * tmu(ji,jj ) + v_ice(ji+1,jj ) * tmu(ji+1,jj ) &149 + v _ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &145 zcmo(ji,jj,8) = zindb * ( vi_ice(ji,jj ) * tmu(ji,jj ) + vi_ice(ji+1,jj ) * tmu(ji+1,jj ) & 146 + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 150 147 / ztmu 151 zcmo(ji,jj,9) = sst_io(ji,jj) 152 zcmo(ji,jj,10) = sss_io(ji,jj) 153 154 zcmo(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj) 155 zcmo(ji,jj,12) = fsolar (ji,jj) 156 zcmo(ji,jj,13) = fnsolar(ji,jj) 148 zcmo(ji,jj,9) = sst_m(ji,jj) 149 zcmo(ji,jj,10) = sss_m(ji,jj) 150 zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj) 151 zcmo(ji,jj,12) = qsr(ji,jj) 152 zcmo(ji,jj,13) = qns(ji,jj) 157 153 ! See thersf for the coefficient 158 zcmo(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce159 zcmo(ji,jj,15) = gtaux(ji,jj)160 zcmo(ji,jj,16) = gtauy(ji,jj)161 zcmo(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce(ji,jj)162 zcmo(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj)154 zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce !!gm ??? 155 zcmo(ji,jj,15) = utaui_ice(ji,jj) 156 zcmo(ji,jj,16) = vtaui_ice(ji,jj) 157 zcmo(ji,jj,17) = qsr_ice(ji,jj) 158 zcmo(ji,jj,18) = qns_ice(ji,jj) 163 159 zcmo(ji,jj,19) = sprecip(ji,jj) 164 160 END DO … … 175 171 END DO 176 172 177 IF ( jf == 7 .OR. jf == 8 .OR. jf == 15 .OR. jf == 16 ) THEN173 IF( jf == 7 .OR. jf == 8 .OR. jf == 15 .OR. jf == 16 ) THEN 178 174 CALL lbc_lnk( zfield, 'T', -1. ) 179 175 ELSE … … 181 177 ENDIF 182 178 183 IF ( nc(jf) == 1 )CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )179 IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 184 180 185 181 END DO 186 182 187 IF ( ( nfice * niter + nit000 - 1 ) >= nitend ) THEN 188 CALL histclo( nice ) 189 ENDIF 183 IF( ( nn_fsbc * niter + nit000 - 1 ) >= nitend ) CALL histclo( nice ) 190 184 ! 191 185 END SUBROUTINE lim_wri_2 … … 225 219 field_13, field_14, field_15, field_16, field_17, field_18, & 226 220 field_19 227 !!gm NAMELIST/namiceout/ noumef, & 228 !! zfield( 1), zfield( 2), zfield( 3), zfield( 4), zfield( 5), & 229 !! zfield( 6), zfield( 7), zfield( 8), zfield( 9), zfield(10), & 230 !! zfield(11), zfield(12), zfield(13), zfield(14), zfield(15), & 231 !!gm zfield(16), zfield(17), zfield(18), zfield(19) 232 !!------------------------------------------------------------------- 233 234 ! Read Namelist namicewri 235 REWIND ( numnam_ice ) 221 !!------------------------------------------------------------------- 222 223 REWIND ( numnam_ice ) ! Read Namelist namicewri 236 224 READ ( numnam_ice , namiceout ) 237 225 238 zfield( 1)= field_1239 zfield( 2)= field_2240 zfield( 3)= field_3241 zfield( 4)= field_4242 zfield( 5)= field_5243 zfield( 6)= field_6244 zfield( 7)= field_7245 zfield( 8)= field_8246 zfield( 9)= field_9226 zfield( 1) = field_1 227 zfield( 2) = field_2 228 zfield( 3) = field_3 229 zfield( 4) = field_4 230 zfield( 5) = field_5 231 zfield( 6) = field_6 232 zfield( 7) = field_7 233 zfield( 8) = field_8 234 zfield( 9) = field_9 247 235 zfield(10) = field_10 248 236 zfield(11) = field_11 … … 274 262 DO nf = 1 , noumef 275 263 WRITE(numout,*) ' ', titn(nf), ' ', nam(nf),' ', uni(nf),' ', nc(nf),' ', cmulti(nf), & 276 ' ', cadd(nf)264 & ' ', cadd(nf) 277 265 END DO 278 266 ENDIF -
trunk/NEMO/LIM_SRC_2/limwri_dimg_2.h90
r823 r888 2 2 !!---------------------------------------------------------------------- 3 3 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 4 !! $ Header$4 !! $ Id: $ 5 5 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 6 6 !!---------------------------------------------------------------------- … … 82 82 83 83 zsto = rdt_ice 84 zout = nwrite * rdt_ice / n fice84 zout = nwrite * rdt_ice / nn_fsbc 85 85 zsec = 0. 86 86 niter = 0 … … 106 106 zcmo(ji,jj,5) = sist (ji,jj) 107 107 zcmo(ji,jj,6) = fbif (ji,jj) 108 zcmo(ji,jj,7) = zindb * ( u _ice(ji,jj ) * tmu(ji,jj ) + u_ice(ji+1,jj ) * tmu(ji+1,jj ) &109 + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &108 zcmo(ji,jj,7) = zindb * ( ui_ice(ji,jj ) * tmu(ji,jj ) + ui_ice(ji+1,jj ) * tmu(ji+1,jj ) & 109 & + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 110 110 / ztmu 111 111 112 zcmo(ji,jj,8) = zindb * ( v _ice(ji,jj ) * tmu(ji,jj ) + v_ice(ji+1,jj ) * tmu(ji+1,jj ) &113 + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &112 zcmo(ji,jj,8) = zindb * ( vi_ice(ji,jj ) * tmu(ji,jj ) + vi_ice(ji+1,jj ) * tmu(ji+1,jj ) & 113 & + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 114 114 / ztmu 115 zcmo(ji,jj,9) = sst_ io(ji,jj)116 zcmo(ji,jj,10) = sss_ io(ji,jj)117 118 zcmo(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj)119 zcmo(ji,jj,12) = fsolar(ji,jj)120 zcmo(ji,jj,13) = fnsolar(ji,jj)115 zcmo(ji,jj,9) = sst_m(ji,jj) 116 zcmo(ji,jj,10) = sss_m(ji,jj) 117 118 zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj) 119 zcmo(ji,jj,12) = qsr(ji,jj) 120 zcmo(ji,jj,13) = qns(ji,jj) 121 121 ! See thersf for the coefficient 122 zcmo(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce123 zcmo(ji,jj,15) = gtaux(ji,jj)124 zcmo(ji,jj,16) = gtauy(ji,jj)125 zcmo(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce(ji,jj)126 zcmo(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj)122 zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 123 zcmo(ji,jj,15) = utaui_ice(ji,jj) 124 zcmo(ji,jj,16) = vtaui_ice(ji,jj) 125 zcmo(ji,jj,17) = qsr_ice(ji,jj) 126 zcmo(ji,jj,18) = qns_ice(ji,jj) 127 127 zcmo(ji,jj,19) = sprecip(ji,jj) 128 128 END DO … … 132 132 nmoyice = nmoyice + 1 133 133 ! compute mean value if it is time to write on file 134 IF ( MOD(kt+n fice-1-nit000+1,nwrite) == 0 ) THEN134 IF ( MOD(kt+nn_fsbc-1-nit000+1,nwrite) == 0 ) THEN 135 135 rcmoy(:,:,:) = rcmoy(:,:,:) / FLOAT(nmoyice) 136 136 #else 137 IF ( MOD(kt-n fice-1-nit000+1,nwrite) == 0 ) THEN137 IF ( MOD(kt-nn_fsbc-1-nit000+1,nwrite) == 0 ) THEN 138 138 ! case of instantaneaous output rcmoy(:,:, 1:jpnoumax ) = 0.e0 139 139 DO jj = 2 , jpjm1 … … 149 149 rcmoy(ji,jj,5) = sist (ji,jj) 150 150 rcmoy(ji,jj,6) = fbif (ji,jj) 151 rcmoy(ji,jj,7) = zindb * ( u _ice(ji,jj ) * tmu(ji,jj ) + u_ice(ji+1,jj ) * tmu(ji+1,jj ) &152 + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &151 rcmoy(ji,jj,7) = zindb * ( ui_ice(ji,jj ) * tmu(ji,jj ) + ui_ice(ji+1,jj ) * tmu(ji+1,jj ) & 152 & + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 153 153 / ztmu 154 154 155 rcmoy(ji,jj,8) = zindb * ( v _ice(ji,jj ) * tmu(ji,jj ) + v_ice(ji+1,jj ) * tmu(ji+1,jj ) &156 + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &155 rcmoy(ji,jj,8) = zindb * ( vi_ice(ji,jj ) * tmu(ji,jj ) + vi_ice(ji+1,jj ) * tmu(ji+1,jj ) & 156 & + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) & 157 157 / ztmu 158 rcmoy(ji,jj,9) = sst_ io(ji,jj)159 rcmoy(ji,jj,10) = sss_ io(ji,jj)160 161 rcmoy(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj)162 rcmoy(ji,jj,12) = fsolar(ji,jj)163 rcmoy(ji,jj,13) = fnsolar(ji,jj)158 rcmoy(ji,jj,9) = sst_m(ji,jj) 159 rcmoy(ji,jj,10) = sss_m(ji,jj) 160 161 rcmoy(ji,jj,11) = qns(ji,jj) + qsr(ji,jj) 162 rcmoy(ji,jj,12) = qsr(ji,jj) 163 rcmoy(ji,jj,13) = qns(ji,jj) 164 164 ! See thersf for the coefficient 165 rcmoy(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce166 rcmoy(ji,jj,15) = gtaux(ji,jj)167 rcmoy(ji,jj,16) = gtauy(ji,jj)168 rcmoy(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce(ji,jj)169 rcmoy(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj)165 rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 166 rcmoy(ji,jj,15) = utaui_ice(ji,jj) 167 rcmoy(ji,jj,16) = vtaui_ice(ji,jj) 168 rcmoy(ji,jj,17) = qsr_ice(ji,jj) 169 rcmoy(ji,jj,18) = qns_ice(ji,jj) 170 170 rcmoy(ji,jj,19) = sprecip(ji,jj) 171 171 END DO … … 201 201 rcmoy(:,:,:) = 0.0 202 202 nmoyice = 0 203 END IF ! MOD(kt+n fice-1-nit000+1, nwrite == 0 ) !203 END IF ! MOD(kt+nn_fsbc-1-nit000+1, nwrite == 0 ) ! 204 204 205 205 END SUBROUTINE lim_wri_2 -
trunk/NEMO/LIM_SRC_2/par_ice_2.F90
r823 r888 7 7 !!---------------------------------------------------------------------- 8 8 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 9 !! $ Header$9 !! $ Id: $ 10 10 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 11 11 !!---------------------------------------------------------------------- -
trunk/NEMO/LIM_SRC_2/thd_ice_2.F90
r823 r888 9 9 !!---------------------------------------------------------------------- 10 10 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 11 !! $ Header$11 !! $ Id: $ 12 12 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 13 13 !!---------------------------------------------------------------------- … … 57 57 fr1_i0_1d , & !: " " fr1_i0 58 58 fr2_i0_1d , & !: " " fr2_i0 59 qns r_ice_1d, & !: " " qns_ice59 qns_ice_1d , & !: " " qns_ice 60 60 qfvbq_1d , & !: " " qfvbq 61 61 sist_1d , & !: " " sist
Note: See TracChangeset
for help on using the changeset viewer.