Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5034 r5600  
    2222   USE thd_ice        ! LIM thermodynamics 
    2323   USE dom_ice        ! LIM domain 
    24    USE par_ice        ! LIM parameters 
    2524   USE ice            ! LIM variables 
    2625   USE limtab         ! LIM 2D <==> 1D 
     
    3231   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3332   USE limthd_ent 
     33   USE limvar 
    3434 
    3535   IMPLICIT NONE 
     
    106106      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i 
    107107      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i 
    108       REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d  ! 1-D version of oa_i 
    109108      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i 
    110109 
     
    112111 
    113112      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity 
     113 
     114      REAL(wp) :: zcai = 1.4e-3_wp 
    114115      !!-----------------------------------------------------------------------! 
    115116 
     
    117118      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    118119      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
    119       CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
    120       CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d ) 
     120      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 
     121      CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d ) 
    121122      CALL wrk_alloc( jpi,jpj, zvrel ) 
    122123 
     124      CALL lim_var_agg(1) 
     125      CALL lim_var_glo2eqv 
    123126      !------------------------------------------------------------------------------| 
    124127      ! 2) Convert units for ice internal energy 
     
    129132               DO ji = 1, jpi 
    130133                  !Energy of melting q(S,T) [J.m-3] 
    131                   rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes 
    132                   e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) & 
    133                       &   / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp ) 
    134                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 
     134                  rswitch          = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )   !0 if no ice 
     135                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp ) 
    135136               END DO 
    136137            END DO 
     
    155156 
    156157      ! Default new ice thickness  
    157       hicol(:,:) = hiccrit 
    158  
    159       IF( fraz_swi == 1 ) THEN 
     158      hicol(:,:) = rn_hnewice 
     159 
     160      IF( ln_frazil ) THEN 
    160161 
    161162         !-------------------- 
     
    166167         zhicrit = 0.04 ! frazil ice thickness 
    167168         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 
    168          zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag) 
     169         zsqcd   = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 
    169170         zgamafr = 0.03 
    170171 
     
    176177                  !------------- 
    177178                  ! C-grid wind stress components 
    178                   ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   & 
    179                      &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp 
    180                   ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   & 
    181                      &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp 
     179                  ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   & 
     180                     &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp 
     181                  ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   & 
     182                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp 
    182183                  ! Square root of wind stress 
    183184                  ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 
     
    195196                  ! C-grid ice velocity 
    196197                  rswitch = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  ) 
    197                   zvgx    = rswitch * ( u_ice(ji-1,jj  ) * tmu(ji-1,jj  )  + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp 
    198                   zvgy    = rswitch * ( v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)  + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp 
     198                  zvgx    = rswitch * ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
     199                  zvgy    = rswitch * ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
    199200 
    200201                  !----------------------------------- 
     
    222223                  iterate_frazil = .true. 
    223224 
    224                   DO WHILE ( iter .LT. 100 .AND. iterate_frazil )  
     225                  DO WHILE ( iter < 100 .AND. iterate_frazil )  
    225226                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 
    226227                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 
     
    266267      ! debug point to follow 
    267268      jiindex_1d = 0 
    268       IF( ln_nicep ) THEN 
    269          DO ji = mi0(jiindx), mi1(jiindx) 
    270             DO jj = mj0(jjindx), mj1(jjindx) 
     269      IF( ln_icectl ) THEN 
     270         DO ji = mi0(iiceprt), mi1(iiceprt) 
     271            DO jj = mj0(jiceprt), mj1(jiceprt) 
    271272               IF ( qlead(ji,jj)  <  0._wp ) THEN 
    272273                  jiindex_1d = (jj - 1) * jpi + ji 
     
    276277      ENDIF 
    277278    
    278       IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
     279      IF( ln_icectl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
    279280 
    280281      !------------------------------ 
     
    290291            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    291292            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    292             CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    293293            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    294294            DO jk = 1, nlay_i 
    295295               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 
    296             END DO ! jk 
    297          END DO ! jl 
     296            END DO 
     297         END DO 
    298298 
    299299         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) ) 
     
    320320         !---------------------- 
    321321         DO ji = 1, nbpac 
    322             zh_newice(ji) = hiccrit 
    323          END DO 
    324          IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 
     322            zh_newice(ji) = rn_hnewice 
     323         END DO 
     324         IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 
    325325 
    326326         !---------------------- 
    327327         ! Salinity of new ice  
    328328         !---------------------- 
    329          SELECT CASE ( num_sal ) 
     329         SELECT CASE ( nn_icesal ) 
    330330         CASE ( 1 )                    ! Sice = constant  
    331             zs_newice(1:nbpac) = bulk_sal 
     331            zs_newice(1:nbpac) = rn_icesal 
    332332         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
    333333            DO ji = 1, nbpac 
    334334               ii =   MOD( npac(ji) - 1 , jpi ) + 1 
    335335               ij =      ( npac(ji) - 1 ) / jpi + 1 
    336                zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  ) 
     336               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij)  ) 
    337337            END DO 
    338338         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
     
    345345         ! We assume that new ice is formed at the seawater freezing point 
    346346         DO ji = 1, nbpac 
    347             ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K) 
     347            ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K) 
    348348            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                             & 
    349                &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) )   & 
    350                &                       - rcp  *         ( ztmelts - rtt )  ) 
    351          END DO ! ji 
     349               &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
     350               &                       - rcp  *         ( ztmelts - rt0 )  ) 
     351         END DO 
    352352 
    353353         !---------------- 
     
    356356         DO ji = 1, nbpac 
    357357            zo_newice(ji) = 0._wp 
    358          END DO ! ji 
     358         END DO 
    359359 
    360360         !------------------- 
     
    363363         DO ji = 1, nbpac 
    364364 
    365             zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg] 
    366  
    367             zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg] 
     365            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg] 
     366 
     367            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg] 
    368368                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)  
    369369                                                                    
     
    372372            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)  
    373373                                                                   ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point    
    374             zv_newice(ji) = - zfmdt / rhoic 
     374            zv_newice(ji) = - zfmdt * r1_rhoic 
    375375 
    376376            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux   
     
    387387            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    388388            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 
    389             zfrazb        = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 
     389            zfrazb        = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 
    390390            zv_frazb(ji)  =         zfrazb   * zv_newice(ji) 
    391391            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
     
    409409         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    410410         DO ji = 1, nbpac 
    411             IF ( za_newice(ji) >  ( amax - zat_i_1d(ji) ) ) THEN 
    412                zda_res(ji)   = za_newice(ji) - ( amax - zat_i_1d(ji) ) 
     411            IF ( za_newice(ji) >  ( rn_amax - zat_i_1d(ji) ) ) THEN 
     412               zda_res(ji)   = za_newice(ji) - ( rn_amax - zat_i_1d(ji) ) 
    413413               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    414414               za_newice(ji) = za_newice(ji) - zda_res  (ji) 
     
    459459            DO jk = 1, nlay_i 
    460460               DO ji = 1, nbpac 
    461                   h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i ) 
     461                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 
    462462                  qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 
    463463               END DO 
     
    478478         ENDDO 
    479479 
    480          !------------ 
    481          ! Update age  
    482          !------------ 
    483          DO jl = 1, jpl 
    484             DO ji = 1, nbpac 
    485                rswitch          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes 
    486                zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch    
    487             END DO  
    488          END DO    
    489  
    490480         !----------------- 
    491481         ! Update salinity 
     
    504494            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 
    505495            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 
    506             CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 
    507496            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 
    508497            DO jk = 1, nlay_i 
     
    525514            DO jj = 1, jpj 
    526515               DO ji = 1, jpi 
    527                   ! heat content in Joules 
    528                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac )  
     516                  ! heat content in J/m2 
     517                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i  
    529518               END DO 
    530519            END DO 
     
    536525      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    537526      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
    538       CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
    539       CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 
     527      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 
     528      CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d ) 
    540529      CALL wrk_dealloc( jpi,jpj, zvrel ) 
    541530      ! 
Note: See TracChangeset for help on using the changeset viewer.