New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5168 for branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 – NEMO

Ignore:
Timestamp:
2015-03-25T10:01:32+01:00 (9 years ago)
Author:
pabouttier
Message:

Fixed several bugs described in Ticket #1360

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r3784 r5168  
    99   !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR 
    1010   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2 
    11    !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init  
     11   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init 
    1212   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization 
    1313   !!---------------------------------------------------------------------- 
     
    4343   IMPLICIT NONE 
    4444   PRIVATE 
    45     
     45 
     46   PUBLIC asm_inc_rea_nam  !: Read namelist for ASM 
    4647   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights 
    4748   PUBLIC   calc_date      !: Compute the calendar date YYYYMMDD on a given step 
     
    7071   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components 
    7172   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S 
    72    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components  
     73   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components 
    7374   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step 
    7475#if defined key_asminc 
     
    7879   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term 
    7980   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization 
    80    INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval  
     81   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval 
    8182   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval 
    82    !  
     83   ! 
    8384   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting 
    84    !                                 !: = 1   Linear hat-like, centred in middle of IAU interval  
     85   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval 
    8586   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix) 
    8687 
     
    100101CONTAINS 
    101102 
    102    SUBROUTINE asm_inc_init 
    103       !!---------------------------------------------------------------------- 
    104       !!                    ***  ROUTINE asm_inc_init  *** 
    105       !!           
    106       !! ** Purpose : Initialize the assimilation increment and IAU weights. 
    107       !! 
    108       !! ** Method  : Initialize the assimilation increment and IAU weights. 
    109       !! 
    110       !! ** Action  :  
    111       !!---------------------------------------------------------------------- 
    112       !! 
    113       !! 
    114       INTEGER :: ji,jj,jk 
    115       INTEGER :: jt 
    116       INTEGER :: imid 
    117       INTEGER :: inum 
    118       INTEGER :: iiauper         ! Number of time steps in the IAU period 
    119       INTEGER :: icycper         ! Number of time steps in the cycle 
    120       INTEGER :: iitend_date     ! Date YYYYMMDD of final time step 
    121       INTEGER :: iitbkg_date     ! Date YYYYMMDD of background time step for Jb term 
    122       INTEGER :: iitdin_date     ! Date YYYYMMDD of background time step for DI 
    123       INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step 
    124       INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step 
    125  
    126       REAL(wp) :: znorm        ! Normalization factor for IAU weights 
    127       REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights  
    128                                ! (should be equal to one) 
    129       REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid 
    130       REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid 
    131       REAL(wp) :: zdate_bkg    ! Date in background state file for DI 
    132       REAL(wp) :: zdate_inc    ! Time axis in increments file 
    133  
    134       REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 
    135       !! 
     103   SUBROUTINE asm_inc_rea_nam 
    136104      NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
    137105         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
     
    158126      salfixmin = -9999 
    159127      nitbkg    = 0 
    160       nitdin    = 0       
     128      nitdin    = 0 
    161129      nitiaustr = 1 
    162130      nitiaufin = 150      ! = 10 days with ORCA2 
     
    187155         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin 
    188156      ENDIF 
     157   END SUBROUTINE 
     158 
     159   SUBROUTINE asm_inc_init 
     160      !!---------------------------------------------------------------------- 
     161      !!                    ***  ROUTINE asm_inc_init  *** 
     162      !! 
     163      !! ** Purpose : Initialize the assimilation increment and IAU weights. 
     164      !! 
     165      !! ** Method  : Initialize the assimilation increment and IAU weights. 
     166      !! 
     167      !! ** Action  : 
     168      !!---------------------------------------------------------------------- 
     169      !! 
     170      !! 
     171      INTEGER :: ji,jj,jk 
     172      INTEGER :: jt 
     173      INTEGER :: imid 
     174      INTEGER :: inum 
     175      INTEGER :: iiauper         ! Number of time steps in the IAU period 
     176      INTEGER :: icycper         ! Number of time steps in the cycle 
     177      INTEGER :: iitend_date     ! Date YYYYMMDD of final time step 
     178      INTEGER :: iitbkg_date     ! Date YYYYMMDD of background time step for Jb term 
     179      INTEGER :: iitdin_date     ! Date YYYYMMDD of background time step for DI 
     180      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step 
     181      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step 
     182 
     183      REAL(wp) :: znorm        ! Normalization factor for IAU weights 
     184      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights 
     185                               ! (should be equal to one) 
     186      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid 
     187      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid 
     188      REAL(wp) :: zdate_bkg    ! Date in background state file for DI 
     189      REAL(wp) :: zdate_inc    ! Time axis in increments file 
     190 
     191      REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 
     192      !! 
    189193 
    190194      nitbkg_r    = nitbkg    + nit000 - 1  ! Background time referenced to nit000 
     
    293297 
    294298            !--------------------------------------------------------- 
    295             ! Constant IAU forcing  
     299            ! Constant IAU forcing 
    296300            !--------------------------------------------------------- 
    297301 
     
    303307 
    304308            !--------------------------------------------------------- 
    305             ! Linear hat-like, centred in middle of IAU interval  
     309            ! Linear hat-like, centred in middle of IAU interval 
    306310            !--------------------------------------------------------- 
    307311 
     
    309313            znorm = 0.0 
    310314            IF ( MOD( iiauper, 2 ) == 0 ) THEN  ! Even number of time steps in IAU interval 
    311                imid = iiauper / 2  
     315               imid = iiauper / 2 
    312316               DO jt = 1, imid 
    313317                  znorm = znorm + REAL( jt ) 
     
    315319               znorm = 2.0 * znorm 
    316320            ELSE                               ! Odd number of time steps in IAU interval 
    317                imid = ( iiauper + 1 ) / 2         
     321               imid = ( iiauper + 1 ) / 2 
    318322               DO jt = 1, imid - 1 
    319323                  znorm = znorm + REAL( jt ) 
     
    342346             DO jt = 1, icycper 
    343347                ztotwgt = ztotwgt + wgtiau(jt) 
    344                 WRITE(numout,*) '         ', jt, '       ', wgtiau(jt)  
    345              END DO    
     348                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
     349             END DO 
    346350             WRITE(numout,*) '         ===================================' 
    347351             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt 
    348352             WRITE(numout,*) '         ===================================' 
    349353          ENDIF 
    350           
     354 
    351355      ENDIF 
    352356 
     
    381385         CALL iom_open( c_asminc, inum ) 
    382386 
    383          CALL iom_get( inum, 'time', zdate_inc )  
     387         CALL iom_get( inum, 'time', zdate_inc ) 
    384388 
    385389         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb ) 
     
    389393 
    390394         IF(lwp) THEN 
    391             WRITE(numout,*)  
     395            WRITE(numout,*) 
    392396            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', & 
    393397               &            ' between dates ', NINT( z_inc_dateb ),' and ',  & 
     
    405409            &                ' not agree with Direct Initialization time' ) 
    406410 
    407          IF ( ln_trainc ) THEN    
     411         IF ( ln_trainc ) THEN 
    408412            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
    409413            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
     
    417421         ENDIF 
    418422 
    419          IF ( ln_dyninc ) THEN    
    420             CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )               
    421             CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )               
     423         IF ( ln_dyninc ) THEN 
     424            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 ) 
     425            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 ) 
    422426            ! Apply the masks 
    423427            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:) 
     
    428432            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0 
    429433         ENDIF 
    430          
     434 
    431435         IF ( ln_sshinc ) THEN 
    432436            CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 ) 
     
    448452 
    449453         CALL iom_close( inum ) 
    450   
     454 
    451455      ENDIF 
    452456 
     
    457461      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 
    458462 
    459          CALL wrk_alloc(jpi,jpj,hdiv)  
     463         CALL wrk_alloc(jpi,jpj,hdiv) 
    460464 
    461465         DO  jt = 1, nn_divdmp 
     
    482486                     u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj)   & 
    483487                                                                        - e1t(ji  ,jj)*e2t(ji  ,jj) * hdiv(ji  ,jj) ) & 
    484                                                                       / e1u(ji,jj) * umask(ji,jj,jk)  
     488                                                                      / e1u(ji,jj) * umask(ji,jj,jk) 
    485489                     v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1)   & 
    486490                                                                        - e1t(ji,jj  )*e2t(ji,jj  ) * hdiv(ji,jj  ) ) & 
    487                                                                       / e2v(ji,jj) * vmask(ji,jj,jk)  
     491                                                                      / e2v(ji,jj) * vmask(ji,jj,jk) 
    488492                  END DO 
    489493               END DO 
     
    493497         END DO 
    494498 
    495          CALL wrk_dealloc(jpi,jpj,hdiv)  
     499         CALL wrk_dealloc(jpi,jpj,hdiv) 
    496500 
    497501      ENDIF 
     
    523527         CALL iom_open( c_asmdin, inum ) 
    524528 
    525          CALL iom_get( inum, 'rdastp', zdate_bkg )  
    526          
     529         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
     530 
    527531         IF(lwp) THEN 
    528             WRITE(numout,*)  
     532            WRITE(numout,*) 
    529533            WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', & 
    530534               &  NINT( zdate_bkg ) 
     
    536540            &                ' not agree with Direct Initialization time' ) 
    537541 
    538          IF ( ln_trainc ) THEN    
     542         IF ( ln_trainc ) THEN 
    539543            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) 
    540544            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg ) 
     
    543547         ENDIF 
    544548 
    545          IF ( ln_dyninc ) THEN    
     549         IF ( ln_dyninc ) THEN 
    546550            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) 
    547551            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg ) 
     
    549553            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:) 
    550554         ENDIF 
    551          
     555 
    552556         IF ( ln_sshinc ) THEN 
    553557            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg ) 
     
    565569      !!---------------------------------------------------------------------- 
    566570      !!                    ***  ROUTINE calc_date  *** 
    567       !!           
     571      !! 
    568572      !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step. 
    569573      !! 
    570574      !! ** Method  : Compute the calendar date YYYYMMDD at a given time step. 
    571575      !! 
    572       !! ** Action  :  
     576      !! ** Action  : 
    573577      !!---------------------------------------------------------------------- 
    574578      INTEGER, INTENT(IN) :: kit000  ! Initial time step 
     
    595599      iyea0 =   kdate0 / 10000 
    596600      imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100 
    597       iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )  
     601      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 ) 
    598602 
    599603      ! Check that kt >= kit000 - 1 
     
    610614      ! Compute the number of days from the initial date 
    611615      idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. ) 
    612     
     616 
    613617      iday    = iday0 
    614618      imon    = imon0 
     
    628632            iyea = iyea + 1 
    629633            CALL calc_month_len( iyea, imonth_len )  ! update month lengths 
    630          ENDIF                  
     634         ENDIF 
    631635         idaycnt = idaycnt + 1 
    632636      END DO 
     
    640644      !!---------------------------------------------------------------------- 
    641645      !!                    ***  ROUTINE calc_month_len  *** 
    642       !!           
     646      !! 
    643647      !! ** Purpose : Compute the number of days in a months given a year. 
    644648      !! 
    645       !! ** Method  :  
     649      !! ** Method  : 
    646650      !!---------------------------------------------------------------------- 
    647651      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year 
     
    650654      ! 
    651655      ! length of the month of the current year (from nleapy, read in namelist) 
    652       IF ( nleapy < 2 ) THEN  
     656      IF ( nleapy < 2 ) THEN 
    653657         imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) 
    654658         IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years 
     
    667671      !!---------------------------------------------------------------------- 
    668672      !!                    ***  ROUTINE tra_asm_inc  *** 
    669       !!           
     673      !! 
    670674      !! ** Purpose : Apply the tracer (T and S) assimilation increments 
    671675      !! 
    672676      !! ** Method  : Direct initialization or Incremental Analysis Updating 
    673677      !! 
    674       !! ** Action  :  
     678      !! ** Action  : 
    675679      !!---------------------------------------------------------------------- 
    676680      INTEGER, INTENT(IN) :: kt               ! Current time step 
     
    682686      !!---------------------------------------------------------------------- 
    683687 
    684       ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)  
    685       ! used to prevent the applied increments taking the temperature below the local freezing point  
    686  
    687 #if defined key_cice  
     688      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 
     689      ! used to prevent the applied increments taking the temperature below the local freezing point 
     690 
     691#if defined key_cice 
    688692        fzptnz(:,:,:) = -1.8_wp 
    689 #else  
     693#else 
    690694        DO jk = 1, jpk 
    691695           DO jj = 1, jpj 
    692696              DO ji = 1, jpk 
    693                  fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) )                   &  
    694                                                   - 2.154996e-4_wp *       tsn(ji,jj,jk,jp_sal)   ) * tsn(ji,jj,jk,jp_sal)  &  
    695                                                   - 7.53e-4_wp * fsdepw(ji,jj,jk)       ! (pressure in dbar)  
     697                 fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) )                   & 
     698                                                  - 2.154996e-4_wp *       tsn(ji,jj,jk,jp_sal)   ) * tsn(ji,jj,jk,jp_sal)  & 
     699                                                  - 7.53e-4_wp * fsdepw(ji,jj,jk)       ! (pressure in dbar) 
    696700              END DO 
    697701           END DO 
    698702        END DO 
    699 #endif  
     703#endif 
    700704 
    701705      IF ( ln_asmiau ) THEN 
     
    711715 
    712716            IF(lwp) THEN 
    713                WRITE(numout,*)  
     717               WRITE(numout,*) 
    714718               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', & 
    715719                  &  kt,' with IAU weight = ', wgtiau(it) 
     
    722726                  ! Do not apply negative increments if the temperature will fall below freezing 
    723727                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 
    724                      &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )  
    725                      tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt   
     728                     &   t_bkg(:,:,jk) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
     729                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
    726730                  END WHERE 
    727731               ELSE 
    728                   tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt   
     732                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
    729733               ENDIF 
    730734               IF (ln_salfix) THEN 
     
    732736                  ! minimum value salfixmin 
    733737                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 
    734                      &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )  
     738                     &  s_bkg(:,:,jk) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
    735739                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 
    736740                  END WHERE 
     
    753757         ! Direct Initialization 
    754758         !-------------------------------------------------------------------- 
    755              
     759 
    756760         IF ( kt == nitdin_r ) THEN 
    757761 
     
    762766               ! Do not apply negative increments if the temperature will fall below freezing 
    763767               WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 
    764                   &   tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )  
    765                   tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     768                  &   tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
     769                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 
    766770               END WHERE 
    767771            ELSE 
    768                tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     772               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 
    769773            ENDIF 
    770774            IF (ln_salfix) THEN 
     
    772776               ! minimum value salfixmin 
    773777               WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 
    774                   &   tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )  
    775                   tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
     778                  &   tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
     779                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 
    776780               END WHERE 
    777781            ELSE 
    778                tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
     782               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 
    779783            ENDIF 
    780784 
     
    782786 
    783787            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities 
    784           
     788 
    785789            IF( ln_zps .AND. .NOT. lk_c1d ) & 
    786790               &  CALL zps_hde( nit000, jpts, tsb, &  ! Partial steps: before horizontal derivative 
     
    797801            DEALLOCATE( s_bkg    ) 
    798802         ENDIF 
    799          !   
     803         ! 
    800804      ENDIF 
    801805      ! Perhaps the following call should be in step 
     
    808812      !!---------------------------------------------------------------------- 
    809813      !!                    ***  ROUTINE dyn_asm_inc  *** 
    810       !!           
     814      !! 
    811815      !! ** Purpose : Apply the dynamics (u and v) assimilation increments. 
    812816      !! 
    813817      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    814818      !! 
    815       !! ** Action  :  
     819      !! ** Action  : 
    816820      !!---------------------------------------------------------------------- 
    817821      INTEGER, INTENT(IN) :: kt   ! Current time step 
     
    834838 
    835839            IF(lwp) THEN 
    836                WRITE(numout,*)  
     840               WRITE(numout,*) 
    837841               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', & 
    838842                  &  kt,' with IAU weight = ', wgtiau(it) 
     
    845849               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 
    846850            END DO 
    847             
     851 
    848852            IF ( kt == nitiaufin_r ) THEN 
    849853               DEALLOCATE( u_bkginc ) 
     
    853857         ENDIF 
    854858 
    855       ELSEIF ( ln_asmdin ) THEN  
     859      ELSEIF ( ln_asmdin ) THEN 
    856860 
    857861         !-------------------------------------------------------------------- 
    858862         ! Direct Initialization 
    859863         !-------------------------------------------------------------------- 
    860           
     864 
    861865         IF ( kt == nitdin_r ) THEN 
    862866 
     
    865869            ! Initialize the now fields with the background + increment 
    866870            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 
    867             vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)   
     871            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
    868872 
    869873            ub(:,:,:) = un(:,:,:)         ! Update before fields 
    870874            vb(:,:,:) = vn(:,:,:) 
    871   
     875 
    872876            DEALLOCATE( u_bkg    ) 
    873877            DEALLOCATE( v_bkg    ) 
     
    884888      !!---------------------------------------------------------------------- 
    885889      !!                    ***  ROUTINE ssh_asm_inc  *** 
    886       !!           
     890      !! 
    887891      !! ** Purpose : Apply the sea surface height assimilation increment. 
    888892      !! 
    889893      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    890894      !! 
    891       !! ** Action  :  
     895      !! ** Action  : 
    892896      !!---------------------------------------------------------------------- 
    893897      INTEGER, INTENT(IN) :: kt   ! Current time step 
     
    910914 
    911915            IF(lwp) THEN 
    912                WRITE(numout,*)  
     916               WRITE(numout,*) 
    913917               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', & 
    914918                  &  kt,' with IAU weight = ', wgtiau(it) 
     
    938942 
    939943            ! Initialize the now fields the background + increment 
    940             sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   
     944            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 
    941945 
    942946            ! Update before fields 
    943             sshb(:,:) = sshn(:,:)          
     947            sshb(:,:) = sshn(:,:) 
    944948 
    945949            IF( lk_vvl ) THEN 
     
    961965      !!---------------------------------------------------------------------- 
    962966      !!                    ***  ROUTINE seaice_asm_inc  *** 
    963       !!           
     967      !! 
    964968      !! ** Purpose : Apply the sea ice assimilation increment. 
    965969      !! 
    966970      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    967971      !! 
    968       !! ** Action  :  
     972      !! ** Action  : 
    969973      !! 
    970974      !! History : 
     
    9981002 
    9991003            it = kt - nit000 + 1 
    1000             zincwgt = wgtiau(it)      ! IAU weight for the current time step  
     1004            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
    10011005            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 
    10021006 
    10031007            IF(lwp) THEN 
    1004                WRITE(numout,*)  
     1008               WRITE(numout,*) 
    10051009               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 
    10061010                  &  kt,' with IAU weight = ', wgtiau(it) 
     
    10211025            ! Nudge sea ice depth to bring it up to a required minimum depth 
    10221026 
    1023             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
    1024                zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
     1027            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
     1028               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 
    10251029            ELSEWHERE 
    10261030               zhicifinc(:,:) = 0.0_wp 
     
    10291033! nudge ice depth 
    10301034            hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 
    1031             phicif(:,:)=phicif(:,:) + zhicifinc(:,:)        
     1035            phicif(:,:)=phicif(:,:) + zhicifinc(:,:) 
    10321036 
    10331037! seaice salinity balancing (to add) 
     
    10711075            zofrld(:,:)=frld(:,:) 
    10721076            zohicif(:,:)=hicif(:,:) 
    1073   
     1077 
    10741078            ! Initialize the now fields the background + increment 
    10751079 
    10761080            frld(:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 
    1077             pfrld(:,:) = frld(:,:)  
     1081            pfrld(:,:) = frld(:,:) 
    10781082            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction 
    10791083 
     
    10821086            ! Nudge sea ice depth to bring it up to a required minimum depth 
    10831087 
    1084             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
    1085                zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
     1088            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
     1089               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 
    10861090            ELSEWHERE 
    10871091               zhicifinc(:,:) = 0.0_wp 
     
    10901094! nudge ice depth 
    10911095            hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 
    1092             phicif(:,:)=phicif(:,:)        
     1096            phicif(:,:)=phicif(:,:) 
    10931097 
    10941098! seaice salinity balancing (to add) 
    1095    
    1096 #endif 
    1097   
     1099 
     1100#endif 
     1101 
    10981102#if defined key_cice 
    10991103 
     
    11101114#if defined key_cice 
    11111115 
    1112 ! Zero ice increment tendency into CICE  
     1116! Zero ice increment tendency into CICE 
    11131117            ndaice_da(:,:) = 0.0_wp 
    11141118 
    11151119#endif 
    1116           
     1120 
    11171121         ENDIF 
    11181122 
    11191123!#if defined key_lim2 || defined key_cice 
    11201124! 
    1121 !            IF (ln_seaicebal ) THEN        
     1125!            IF (ln_seaicebal ) THEN 
    11221126!             !! balancing salinity increments 
    11231127!             !! simple case from limflx.F90 (doesn't include a mass flux) 
     
    11311135! 
    11321136!             DO jj = 1, jpj 
    1133 !               DO ji = 1, jpi  
     1137!               DO ji = 1, jpi 
    11341138!           ! calculate change in ice and snow mass per unit area 
    11351139!           ! positive values imply adding salt to the ocean (results from ice formation) 
     
    11421146! 
    11431147!           ! prevent small mld 
    1144 !           ! less than 10m can cause salinity instability  
     1148!           ! less than 10m can cause salinity instability 
    11451149!                 IF (mld < 10) mld=10 
    11461150! 
    1147 !           ! set to bottom of a level  
     1151!           ! set to bottom of a level 
    11481152!                 DO jk = jpk-1, 2, -1 
    1149 !                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN  
     1153!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 
    11501154!                     mld=gdepw(ji,jj,jk+1) 
    11511155!                     jkmax=jk 
     
    11541158! 
    11551159!            ! avoid applying salinity balancing in shallow water or on land 
    1156 !            !  
     1160!            ! 
    11571161! 
    11581162!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 
     
    11651169! 
    11661170!           ! put increments in for levels in the mixed layer 
    1167 !           ! but prevent salinity below a threshold value  
    1168 ! 
    1169 !                   DO jk = 1, jkmax               
    1170 ! 
    1171 !                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN  
     1171!           ! but prevent salinity below a threshold value 
     1172! 
     1173!                   DO jk = 1, jkmax 
     1174! 
     1175!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 
    11721176!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 
    11731177!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 
     
    11801184!      ! 
    11811185!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean 
    1182 !      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt   
    1183 !      !!                
    1184 !      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)  
     1186!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
     1187!      !! 
     1188!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d) 
    11851189!      !!                                                     ! E-P (kg m-2 s-2) 
    11861190!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2) 
Note: See TracChangeset for help on using the changeset viewer.