Changeset 2528 for trunk/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90
- Property svn:eol-style deleted
r1715 r2528 2 2 !!====================================================================== 3 3 !! *** MODULE limdmp_2 *** 4 !! Ice model : restoring Ice thickness and Fraction leads4 !! LIM-2 ice model : restoring Ice thickness and Fraction leads 5 5 !!====================================================================== 6 !! History : 2.0 ! 04-04 (S. Theetten) Original code 6 !! History : 2.0 ! 2004-04 (S. Theetten) Original code 7 !! 3.3 ! 2010-06 (J.-M. Molines) use of fldread 7 8 !!---------------------------------------------------------------------- 8 #if defined key_lim2 && defined key_tradmp9 #if defined key_lim2 9 10 !!---------------------------------------------------------------------- 10 !! 'key_lim2' AND LIM 2.0 sea-ice model 11 !! 'key_tradmp' Damping 12 !!---------------------------------------------------------------------- 11 !! 'key_lim2' LIM 2.0 sea-ice model 13 12 !!---------------------------------------------------------------------- 14 13 !! lim_dmp_2 : ice model damping 15 14 !!---------------------------------------------------------------------- 16 15 USE in_out_manager ! I/O manager 17 USE phycst ! physical constants 18 USE ice_2 19 USE tradmp 20 USE dom_oce 21 USE oce 22 USE iom 16 USE ice_2 ! ice variables 17 USE sbc_oce, ONLY : nn_fsbc ! for fldread 18 USE dom_oce ! for mi0; mi1 etc ... 19 USE fldread ! read input fields 23 20 24 21 IMPLICIT NONE 25 22 PRIVATE 26 23 27 PUBLIC lim_dmp_2 ! called by ice_step_2 24 PUBLIC lim_dmp_2 ! called by sbc_ice_lim2 25 26 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: resto_ice ! restoring coeff. on ICE [s-1] 27 28 INTEGER, PARAMETER :: jp_hicif = 1 , jp_frld = 2 29 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_icedmp ! structure of ice damping input 28 30 29 INTEGER :: nice1, nice2, & ! first and second record used30 & inumice_dmp ! logical unit for ice variables (damping)31 REAL(wp), DIMENSION(jpi,jpj) :: hicif_dta , & ! ice thickness at a given time32 & frld_dta ! fraction lead at a given time33 REAL(wp), DIMENSION(jpi,jpj,2) :: hicif_data , & ! ice thickness data at two consecutive times34 & frld_data ! fraction lead data at two consecutive times35 36 31 !! * Substitution 37 32 # include "vectopt_loop_substitute.h90" 38 33 !!---------------------------------------------------------------------- 39 !! LIM 2.0 , UCL-LOCEAN-IPSL (2006)34 !! NEMO/LIM 3.3 , UCL-NEMO-consortium (2010) 40 35 !! $Id$ 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)36 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 42 37 !!---------------------------------------------------------------------- 43 44 38 CONTAINS 45 39 46 SUBROUTINE lim_dmp_2( kt)40 SUBROUTINE lim_dmp_2( kt ) 47 41 !!------------------------------------------------------------------- 48 !! *** ROUTINE lim_dmp_2 ***42 !! *** ROUTINE lim_dmp_2 *** 49 43 !! 50 !! ** purpose : ice model damping : restoring ice thickness and 51 !! fraction leads 44 !! ** purpose : ice model damping : restoring ice thickness and fraction leads 52 45 !! 53 !! ** method : the key_tradmp must be used to compute resto(:,: ) coef.46 !! ** method : the key_tradmp must be used to compute resto(:,:,1) coef. 54 47 !!--------------------------------------------------------------------- 55 INTEGER, INTENT(in) :: kt 48 INTEGER, INTENT(in) :: kt ! ocean time-step 56 49 ! 57 INTEGER :: ji, jj ! dummy loop indices 50 INTEGER :: ji, jj ! dummy loop indices 51 REAL(wp) :: zfrld, zhice ! local scalars 58 52 !!--------------------------------------------------------------------- 59 53 ! 60 CALL dta_lim_2( kt ) 61 62 DO jj = 2, jpjm1 63 DO ji = fs_2, fs_jpim1 ! vector opt. 64 hicif(ji,jj) = hicif(ji,jj) - rdt_ice * resto(ji,jj,1) * ( hicif(ji,jj) - hicif_dta(ji,jj) ) 65 frld(ji,jj) = frld (ji,jj) - rdt_ice * resto(ji,jj,1) * ( frld(ji,jj) - frld_dta (ji,jj) ) 66 END DO 67 END DO 54 IF (kt == nit000) THEN 55 IF(lwp) WRITE(numout,*) 56 IF(lwp) WRITE(numout,*) 'lim_dmp_2 : Ice thickness and ice concentration restoring' 57 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 58 ! 59 ! ice_resto_init create resto_ice (in 1/s) for restoring ice parameters near open boundaries. 60 ! Double check this routine to verify if it corresponds to your config 61 CALL lim_dmp_init 62 ENDIF 63 ! 64 IF( ln_limdmp ) THEN ! ice restoring in this case 65 ! 66 CALL fld_read( kt, nn_fsbc, sf_icedmp ) 67 ! 68 !CDIR COLLAPSE 69 hicif(:,:) = MAX( 0._wp, & ! h >= 0 avoid spurious out of physical range 70 & hicif(:,:) - rdt_ice * resto_ice(:,:,1) * ( hicif(:,:) - sf_icedmp(jp_hicif)%fnow(:,:,1) ) ) 71 !CDIR COLLAPSE 72 hicif(:,:) = MAX( 0._wp, MIN( 1._wp, & ! 0<= frld<=1 values which blow the run up 73 & frld (:,:) - rdt_ice * resto_ice(:,:,1) * ( frld (:,:) - sf_icedmp(jp_frld )%fnow(:,:,1) ) ) ) 74 ! 75 ENDIF 68 76 ! 69 77 END SUBROUTINE lim_dmp_2 70 78 71 79 72 SUBROUTINE dta_lim_2( kt )80 SUBROUTINE lim_dmp_init 73 81 !!---------------------------------------------------------------------- 74 !! *** ROUTINE dta_lim_2***82 !! *** ROUTINE lim_dmp_init *** 75 83 !! 76 !! ** Purpose : Reads monthly ice thickness and fraction lead data 84 !! ** Purpose : Initialization for the ice thickness and concentration 85 !! restoring 86 !! restoring will be used. It is used to mimic ice open 87 !! boundaries. 77 88 !! 78 !! ** Method : Read on unit numicedt the interpolated ice variable 79 !! onto the model grid. 80 !! Data begin at january. 81 !! The value is centered at the middle of month. 82 !! In the opa model, kt=1 agree with january 1. 83 !! At each time step, a linear interpolation is applied between 84 !! two monthly values. 89 !! ** Method : ????? 85 90 !! 86 !! ** Action : define hicif_dta and frld_dta arrays at time-step kt91 !! ** Action : define resto_ice(:,:,1) 87 92 !!---------------------------------------------------------------------- 88 INTEGER, INTENT(in) :: kt ! ocean time-step 93 INTEGER :: ji, jj, jk ! dummy loop indices 94 INTEGER :: irelax, ierror ! error flag for allocation 89 95 ! 90 INTEGER :: imois, iman, i15 ! temporary integers 91 REAL(wp) :: zxy 96 REAL(wp) :: zdmpmax, zdmpmin, zfactor, zreltim ! temporary scalar 97 ! 98 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 99 TYPE(FLD_N), DIMENSION (2) :: sl_icedmp ! informations about the icedmp field to be read 100 TYPE(FLD_N) :: sn_hicif ! 101 TYPE(FLD_N) :: sn_frld ! 102 NAMELIST/namice_dmp/ cn_dir, ln_limdmp, sn_hicif, sn_frld 92 103 !!---------------------------------------------------------------------- 104 ! 105 ! 1) initialize fld read structure for input data 106 ! -------------------------------------------- 107 ln_limdmp = .false. !* set file information (default values) 108 cn_dir = './' 109 ! (NB: frequency positive => hours, negative => months) 110 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 111 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 112 sn_hicif = FLD_N( 'ice_damping ', -1 , 'hicif' , .true. , .true. , 'yearly' , '' , '' ) 113 sn_frld = FLD_N( 'ice_damping ', -1 , 'frld' , .true. , .true. , 'yearly' , '' , '' ) 93 114 94 ! 0. Initialization 95 ! ----------------- 96 iman = INT( raamo ) 97 !!! better but change the results i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 98 i15 = nday / 16 99 imois = nmonth + i15 - 1 100 IF( imois == 0 ) imois = iman 101 102 ! 1. First call kt=nit000: Initialization and Open 103 ! ----------------------- 104 IF( kt == nit000 ) THEN 105 nice1 = 0 106 IF(lwp) WRITE(numout,*) 107 IF(lwp) WRITE(numout,*) 'dtalim : Ice thickness and lead fraction monthly fields' 108 IF(lwp) WRITE(numout,*) '~~~~~~' 109 IF(lwp) WRITE(numout,*) ' NetCDF FORMAT' 110 IF(lwp) WRITE(numout,*) 111 ! open file 112 CALL iom_open( 'ice_damping_ATL4.nc', inumice_dmp ) 115 REWIND( numnam_ice ) !* read in namelist_ice namicedmp 116 READ ( numnam_ice, namice_dmp ) 117 ! 118 IF ( lwp ) THEN !* control print 119 WRITE (numout,*)' lim_dmp_init : lim_dmp initialization ' 120 WRITE (numout,*)' Namelist namicedmp read ' 121 WRITE (numout,*)' Ice restoring (T) or not (F) ln_limdmp =', ln_limdmp 122 WRITE (numout,*) 123 WRITE (numout,*)' CAUTION : here hard coded ice restoring along northern and southern boundaries' 124 WRITE (numout,*)' adapt the lim_dmp_init routine to your needs' 113 125 ENDIF 114 126 127 ! 2) initialise resto_ice ==> config dependant ! 128 ! -------------------- ++++++++++++++++ 129 ! 130 IF( ln_limdmp ) THEN !* ice restoring is used, follow initialization 131 ! 132 sl_icedmp ( jp_hicif ) = sn_hicif 133 sl_icedmp ( jp_frld ) = sn_frld 134 ALLOCATE ( sf_icedmp (2) , resto_ice(jpi,jpj,1), STAT=ierror ) 135 IF( ierror > 0 ) THEN 136 CALL ctl_stop( 'lim_dmp_init: unable to allocate sf_icedmp structure or resto_ice array' ) ; RETURN 137 ENDIF 138 ALLOCATE( sf_icedmp(jp_hicif)%fnow(jpi,jpj,1) , sf_icedmp(jp_hicif)%fdta(jpi,jpj,1,2) ) 139 ALLOCATE( sf_icedmp(jp_frld )%fnow(jpi,jpj,1) , sf_icedmp(jp_frld )%fdta(jpi,jpj,1,2) ) 140 ! ! fill sf_icedmp with sn_icedmp and control print 141 CALL fld_fill( sf_icedmp, sl_icedmp, cn_dir, 'lim_dmp_init', 'Ice restoring input data', 'namicedmp' ) 142 143 resto_ice(:,:,:) = 0._wp 144 ! Re-calculate the North and South boundary restoring term 145 ! because those boundaries may change with the prescribed zoom area. 146 ! 147 irelax = 16 ! width of buffer zone with respect to close boundary 148 zdmpmax = 10._wp ! max restoring time scale (days) (low restoring) 149 zdmpmin = rdt_ice / 86400._wp ! min restoring time scale (days) (high restoring) 150 ! ! days / grid-point 151 zfactor = ( zdmpmax - zdmpmin ) / REAL( irelax, wp ) 115 152 116 ! 2. Read monthly file 117 ! -------------------- 118 IF( ( kt == nit000 ) .OR. imois /= nice1 ) THEN 119 ! 120 ! Calendar computation 121 nice1 = imois ! first file record used 122 nice2 = nice1 + 1 ! last file record used 123 nice1 = MOD( nice1, iman ) 124 nice2 = MOD( nice2, iman ) 125 IF( nice1 == 0 ) nice1 = iman 126 IF( nice2 == 0 ) nice2 = iman 127 IF(lwp) WRITE(numout,*) 'first record file used nice1 ', nice1 128 IF(lwp) WRITE(numout,*) 'last record file used nice2 ', nice2 129 130 ! Read monthly ice thickness Levitus 131 CALL iom_get( inumice_dmp, jpdom_data, 'iicethic', hicif_data(:,:,1), nice1 ) 132 CALL iom_get( inumice_dmp, jpdom_data, 'iicethic', hicif_data(:,:,2), nice2 ) 133 134 ! Read monthly ice thickness Levitus 135 CALL iom_get( inumice_dmp, jpdom_data, 'ileadfra', frld_data(:,:,1), nice1 ) 136 CALL iom_get( inumice_dmp, jpdom_data, 'ileadfra', frld_data(:,:,2), nice2 ) 137 138 ! The fraction lead read in the file is in fact the 139 ! ice concentration which is 1 - the fraction lead 140 frld_data = 1 - frld_data 141 142 IF(lwp) THEN 143 WRITE(numout,*) 144 WRITE(numout,*) ' Ice thickness month ', nice1,' and ', nice2 145 WRITE(numout,*) 146 WRITE(numout,*) ' Ice thickness month = ', nice1 147 CALL prihre( hicif_data(1,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 148 WRITE(numout,*) 149 WRITE(numout,*) ' Fraction lead months ', nice1,' and ', nice2 150 WRITE(numout,*) 151 WRITE(numout,*) ' Fraction lead month = ', nice1 152 CALL prihre( frld_data(1,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 153 ENDIF 153 ! South boundary restoring term 154 ! REM: if there is no ice in the model and in the data, 155 ! no restoring even with non zero resto_ice 156 DO jj = mj0(jpjzoom - 1 + 1), mj1(jpjzoom -1 + irelax) 157 zreltim = zdmpmin + zfactor * ( mjg(jj) - jpjzoom + 1 ) 158 resto_ice(:,jj,:) = 1._wp / ( zreltim * 86400._wp ) 159 END DO 154 160 155 CALL FLUSH(numout) 156 161 ! North boundary restoring term 162 DO jj = mj0(jpjzoom -1 + jpjglo - irelax), mj1(jpjzoom - 1 + jpjglo) 163 zreltim = zdmpmin + zfactor * (jpjglo - ( mjg(jj) - jpjzoom + 1 )) 164 resto_ice(:,jj,:) = 1.e0 / ( zreltim * 86400 ) 165 END DO 157 166 ENDIF 158 159 ! 3. At every time step compute ice thickness and fraction lead data160 ! ------------------------------------------------------------------161 zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.162 hicif_dta(:,:) = (1.-zxy) * hicif_data(:,:,1) + zxy * hicif_data(:,:,2)163 frld_dta (:,:) = (1.-zxy) * frld_data (:,:,1) + zxy * frld_data (:,:,2)164 165 IF( kt == nitend ) CALL iom_close( inumice_dmp )166 167 ! 167 END SUBROUTINE dta_lim_2168 168 END SUBROUTINE lim_dmp_init 169 169 170 #else 170 171 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.