Ignore:
Timestamp:
2019-12-11T12:02:38+01:00 (11 months ago)
Author:
agn
Message:

updated trunk to v 11653

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/dia25h.F90

    r10641 r12178  
    5555      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics 
    5656      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 ) 
    57 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp ) 
     57901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_dia25h in reference namelist' ) 
    5858      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics 
    5959      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 ) 
    60 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp ) 
     60902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist' ) 
    6161      IF(lwm) WRITE ( numond, nam_dia25h ) 
    6262 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diacfl.F90

    r10824 r12178  
    2929   REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number  
    3030 
    31 !!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc ! 
    32 !!gm          I don't understand why. 
    33    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
    34 !!gm end 
    35  
    3631   PUBLIC   dia_cfl       ! routine called by step.F90 
    3732   PUBLIC   dia_cfl_init  ! routine called by nemogcm 
     
    5550      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    5651      ! 
    57       INTEGER                ::   ji, jj, jk                            ! dummy loop indices 
    58       REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars 
    59       INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace 
    60 !!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
     52      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices 
     53      REAL(wp)                         ::   z2dt, zCu_max, zCv_max, zCw_max  ! local scalars 
     54      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace 
     55      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace 
    6156      !!---------------------------------------------------------------------- 
    6257      ! 
     
    7166      DO jk = 1, jpk       ! calculate Courant numbers 
    7267         DO jj = 1, jpj 
    73             DO ji = 1, fs_jpim1   ! vector opt. 
     68            DO ji = 1, jpi 
    7469               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
    7570               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
     
    111106      !                    ! write out to file 
    112107      IF( lwp ) THEN 
    113          WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
     108         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
    114109         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) 
    115110         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) 
     
    172167      rCw_max = 0._wp 
    173168      ! 
    174 !!gm required to work 
    175       ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) ) 
    176 !!gm end 
    177       !       
    178169   END SUBROUTINE dia_cfl_init 
    179170 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diadct.F90

    r10425 r12178  
    1111   !!            3.4  ! 09/2011 (C Bricaud) 
    1212   !!---------------------------------------------------------------------- 
    13 #if defined key_diadct 
    14    !!---------------------------------------------------------------------- 
    15    !!   'key_diadct' : 
    16    !!---------------------------------------------------------------------- 
     13   !! does not work with agrif 
     14#if ! defined key_agrif 
    1715   !!---------------------------------------------------------------------- 
    1816   !!   dia_dct      :  Compute the transport through a sec. 
     
    4240 
    4341   PUBLIC   dia_dct      ! routine called by step.F90 
    44    PUBLIC   dia_dct_init ! routine called by opa.F90 
    45    PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90  
    46    PRIVATE  readsec 
    47    PRIVATE  removepoints 
    48    PRIVATE  transport 
    49    PRIVATE  dia_dct_wri 
    50  
    51    LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag 
    52  
    53    INTEGER :: nn_dct        ! Frequency of computation 
    54    INTEGER :: nn_dctwri     ! Frequency of output 
    55    INTEGER :: nn_secdebug   ! Number of the section to debug 
     42   PUBLIC   dia_dct_init ! routine called by nemogcm.F90 
     43 
     44   !                         !!** namelist variables ** 
     45   LOGICAL, PUBLIC ::   ln_diadct     !: Calculate transport thru a section or not 
     46   INTEGER         ::   nn_dct        !  Frequency of computation 
     47   INTEGER         ::   nn_dctwri     !  Frequency of output 
     48   INTEGER         ::   nn_secdebug   !  Number of the section to debug 
    5649    
    5750   INTEGER, PARAMETER :: nb_class_max  = 10 
     
    10497CONTAINS 
    10598  
    106   INTEGER FUNCTION diadct_alloc()  
    107      !!----------------------------------------------------------------------  
    108      !!                   ***  FUNCTION diadct_alloc  ***  
    109      !!----------------------------------------------------------------------  
    110      INTEGER :: ierr(2)  
    111      !!----------------------------------------------------------------------  
    112  
    113      ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) )  
    114      ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) )  
    115  
    116      diadct_alloc = MAXVAL( ierr )  
    117      IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )  
    118   
    119   END FUNCTION diadct_alloc  
    120  
     99   INTEGER FUNCTION diadct_alloc()  
     100      !!----------------------------------------------------------------------  
     101      !!                   ***  FUNCTION diadct_alloc  ***  
     102      !!----------------------------------------------------------------------  
     103 
     104      ALLOCATE( transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), & 
     105         &      transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=diadct_alloc )  
     106 
     107      CALL mpp_sum( 'diadct', diadct_alloc )  
     108      IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )  
     109 
     110   END FUNCTION diadct_alloc 
    121111 
    122112   SUBROUTINE dia_dct_init 
     
    130120      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    131121      !! 
    132       NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug 
     122      NAMELIST/nam_diadct/ln_diadct, nn_dct, nn_dctwri, nn_secdebug 
    133123      !!--------------------------------------------------------------------- 
    134124 
    135      REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections 
    136      READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901) 
    137 901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp ) 
    138  
    139      REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections 
    140      READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 ) 
    141 902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp ) 
    142      IF(lwm) WRITE ( numond, namdct ) 
     125     REWIND( numnam_ref )              ! Namelist nam_diadct in reference namelist : Diagnostic: transport through sections 
     126     READ  ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901) 
     127901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' ) 
     128 
     129     REWIND( numnam_cfg )              ! Namelist nam_diadct in configuration namelist : Diagnostic: transport through sections 
     130     READ  ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 ) 
     131902  IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' ) 
     132     IF(lwm) WRITE ( numond, nam_diadct ) 
    143133 
    144134     IF( lwp ) THEN 
     
    146136        WRITE(numout,*) "diadct_init: compute transports through sections " 
    147137        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
    148         WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
    149         WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     138        WRITE(numout,*) "       Calculate transport thru sections: ln_diadct = ", ln_diadct 
     139        WRITE(numout,*) "       Frequency of computation:          nn_dct    = ", nn_dct 
     140        WRITE(numout,*) "       Frequency of write:                nn_dctwri = ", nn_dctwri 
    150141 
    151142        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN 
     
    155146        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug 
    156147        ENDIF 
    157  
     148     ENDIF 
     149 
     150     IF( ln_diadct ) THEN 
     151        ! control 
    158152        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  & 
    159           &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' ) 
    160  
     153           &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' ) 
     154 
     155        ! allocate dia_dct arrays 
     156        IF( diadct_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' ) 
     157 
     158        !Read section_ijglobal.diadct 
     159        CALL readsec 
     160 
     161        !open output file 
     162        IF( lwm ) THEN 
     163           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     164           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     165           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     166        ENDIF 
     167 
     168        ! Initialise arrays to zero  
     169        transports_3d(:,:,:,:)=0.0  
     170        transports_2d(:,:,:)  =0.0  
     171        ! 
    161172     ENDIF 
    162  
    163      !Read section_ijglobal.diadct 
    164      CALL readsec 
    165  
    166      !open output file 
    167      IF( lwm ) THEN 
    168         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    169         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    170         CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    171      ENDIF 
    172  
    173      ! Initialise arrays to zero  
    174      transports_3d(:,:,:,:)=0.0  
    175      transports_2d(:,:,:)  =0.0  
    176173     ! 
    177174  END SUBROUTINE dia_dct_init 
     
    12411238#else 
    12421239   !!---------------------------------------------------------------------- 
    1243    !!   Default option :                                       Dummy module 
     1240   !!   Dummy module                                              
    12441241   !!---------------------------------------------------------------------- 
    1245    LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag 
    1246    PUBLIC  
    1247    !! $Id$ 
     1242   LOGICAL, PUBLIC ::   ln_diadct = .FALSE. 
    12481243CONTAINS 
    1249  
    1250    SUBROUTINE dia_dct_init          ! Dummy routine 
     1244   SUBROUTINE dia_dct_init 
    12511245      IMPLICIT NONE 
    1252       WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?' 
    12531246   END SUBROUTINE dia_dct_init 
    1254  
    1255    SUBROUTINE dia_dct( kt )         ! Dummy routine 
     1247   SUBROUTINE dia_dct( kt ) 
    12561248      IMPLICIT NONE 
    1257       INTEGER, INTENT( in ) :: kt   ! ocean time-step index 
    1258       WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt 
     1249      INTEGER, INTENT(in) ::   kt 
    12591250   END SUBROUTINE dia_dct 
     1251   ! 
    12601252#endif 
    12611253 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaharm.F90

    r10835 r12178  
    55   !!====================================================================== 
    66   !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code 
    7    !!---------------------------------------------------------------------- 
    8 #if defined key_diaharm 
    9    !!---------------------------------------------------------------------- 
    10    !!   'key_diaharm' 
    117   !!---------------------------------------------------------------------- 
    128   USE oce             ! ocean dynamics and tracers variables 
     
    2622   IMPLICIT NONE 
    2723   PRIVATE 
    28  
    29    LOGICAL, PUBLIC, PARAMETER :: lk_diaharm  = .TRUE. 
    3024    
    3125   INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo 
     
    3327 
    3428   !                         !!** namelist variables ** 
    35    INTEGER ::   nit000_han    ! First time step used for harmonic analysis 
    36    INTEGER ::   nitend_han    ! Last time step used for harmonic analysis 
    37    INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis 
    38    INTEGER ::   nb_ana        ! Number of harmonics to analyse 
     29   LOGICAL, PUBLIC ::   ln_diaharm    ! Choose tidal harmonic output or not 
     30   INTEGER         ::   nit000_han    ! First time step used for harmonic analysis 
     31   INTEGER         ::   nitend_han    ! Last time step used for harmonic analysis 
     32   INTEGER         ::   nstep_han     ! Time step frequency for harmonic analysis 
     33   INTEGER         ::   nb_ana        ! Number of harmonics to analyse 
    3934 
    4035   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
     
    5348   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...) 
    5449 
    55    PUBLIC   dia_harm   ! routine called by step.F90 
     50   PUBLIC   dia_harm        ! routine called by step.F90 
     51   PUBLIC   dia_harm_init   ! routine called by nemogcm.F90 
    5652 
    5753   !!---------------------------------------------------------------------- 
     
    7167      !! 
    7268      !!-------------------------------------------------------------------- 
    73       INTEGER :: jh, nhan, jk, ji 
     69      INTEGER ::   jh, nhan, ji 
    7470      INTEGER ::   ios                 ! Local integer output status for namelist read 
    7571 
    76       NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname 
     72      NAMELIST/nam_diaharm/ ln_diaharm, nit000_han, nitend_han, nstep_han, tname 
    7773      !!---------------------------------------------------------------------- 
    7874 
     
    8379      ENDIF 
    8480      ! 
    85       IF( .NOT. ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis') 
    86       ! 
    87       CALL tide_init_Wave 
    88       ! 
    8981      REWIND( numnam_ref )              ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis 
    9082      READ  ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901) 
    91 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_diaharm in reference namelist', lwp ) 
     83901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_diaharm in reference namelist' ) 
    9284      REWIND( numnam_cfg )              ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis 
    9385      READ  ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 ) 
    94 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist', lwp ) 
     86902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist' ) 
    9587      IF(lwm) WRITE ( numond, nam_diaharm ) 
    9688      ! 
    9789      IF(lwp) THEN 
    98          WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han 
    99          WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han 
    100          WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han 
     90         WRITE(numout,*) 'Tidal diagnostics = ', ln_diaharm 
     91         WRITE(numout,*) '   First time step used for analysis:         nit000_han= ', nit000_han 
     92         WRITE(numout,*) '   Last  time step used for analysis:         nitend_han= ', nitend_han 
     93         WRITE(numout,*) '   Time step frequency for harmonic analysis: nstep_han = ', nstep_han 
    10194      ENDIF 
    10295 
    103       ! Basic checks on harmonic analysis time window: 
    104       ! ---------------------------------------------- 
    105       IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   & 
    106          &                                       ' restart capability not implemented' ) 
    107       IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   & 
    108          &                                       'restart capability not implemented' ) 
    109  
    110       IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   & 
    111          &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' ) 
    112  
    113       nb_ana = 0 
    114       DO jk=1,jpmax_harmo 
    115          DO ji=1,jpmax_harmo 
    116             IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN 
    117                nb_ana=nb_ana+1 
    118             ENDIF 
    119          END DO 
    120       END DO 
    121       ! 
    122       IF(lwp) THEN 
    123          WRITE(numout,*) '        Namelist nam_diaharm' 
    124          WRITE(numout,*) '        nb_ana    = ', nb_ana 
    125          CALL flush(numout) 
     96      IF( ln_diaharm .AND. .NOT.ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis') 
     97 
     98      IF( ln_diaharm ) THEN 
     99 
     100         CALL tide_init_Wave 
     101         ! 
     102         ! Basic checks on harmonic analysis time window: 
     103         ! ---------------------------------------------- 
     104         IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   & 
     105            &                                       ' restart capability not implemented' ) 
     106         IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   & 
     107            &                                       'restart capability not implemented' ) 
     108 
     109         IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   & 
     110            &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' ) 
     111         ! 
     112         nb_ana = 0 
     113         DO jh=1,jpmax_harmo 
     114            DO ji=1,jpmax_harmo 
     115               IF(TRIM(tname(jh)) == Wave(ji)%cname_tide) THEN 
     116                  nb_ana=nb_ana+1 
     117               ENDIF 
     118            END DO 
     119         END DO 
     120         ! 
     121         IF(lwp) THEN 
     122            WRITE(numout,*) '        Namelist nam_diaharm' 
     123            WRITE(numout,*) '        nb_ana    = ', nb_ana 
     124            CALL flush(numout) 
     125         ENDIF 
     126         ! 
     127         IF (nb_ana > jpmax_harmo) THEN 
     128            WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo' 
     129            WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo 
     130            CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 ) 
     131         ENDIF 
     132 
     133         ALLOCATE(name    (nb_ana)) 
     134         DO jh=1,nb_ana 
     135            DO ji=1,jpmax_harmo 
     136               IF (TRIM(tname(jh)) ==  Wave(ji)%cname_tide) THEN 
     137                  name(jh) = ji 
     138                  EXIT 
     139               END IF 
     140            END DO 
     141         END DO 
     142 
     143         ! Initialize frequency array: 
     144         ! --------------------------- 
     145         ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) ) 
     146 
     147         CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana ) 
     148 
     149         IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency ' 
     150 
     151         DO jh = 1, nb_ana 
     152            IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh) 
     153         END DO 
     154 
     155         ! Initialize temporary arrays: 
     156         ! ---------------------------- 
     157         ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
     158         ana_temp(:,:,:,:) = 0._wp 
     159 
    126160      ENDIF 
    127       ! 
    128       IF (nb_ana > jpmax_harmo) THEN 
    129          WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo' 
    130          WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo 
    131          CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 ) 
    132       ENDIF 
    133  
    134       ALLOCATE(name    (nb_ana)) 
    135       DO jk=1,nb_ana 
    136        DO ji=1,jpmax_harmo 
    137           IF (TRIM(tname(jk)) ==  Wave(ji)%cname_tide) THEN 
    138              name(jk) = ji 
    139              EXIT 
    140           END IF 
    141        END DO 
    142       END DO 
    143  
    144       ! Initialize frequency array: 
    145       ! --------------------------- 
    146       ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) ) 
    147  
    148       CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana ) 
    149  
    150       IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency ' 
    151  
    152       DO jh = 1, nb_ana 
    153         IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh) 
    154       END DO 
    155  
    156       ! Initialize temporary arrays: 
    157       ! ---------------------------- 
    158       ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    159       ana_temp(:,:,:,:) = 0._wp 
    160161 
    161162   END SUBROUTINE dia_harm_init 
     
    177178      !!-------------------------------------------------------------------- 
    178179      IF( ln_timing )   CALL timing_start('dia_harm') 
    179       ! 
    180       IF( kt == nit000 )   CALL dia_harm_init 
    181180      ! 
    182181      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN 
     
    422421      INTEGER, INTENT(in) ::   init  
    423422      ! 
    424       INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
     423      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd 
    425424      REAL(wp)                        :: zval1, zval2, zx1 
    426425      REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2 
     
    434433         ztmp3(:,:) = 0._wp 
    435434         ! 
    436          DO jk1_sd = 1, nsparse 
    437             DO jk2_sd = 1, nsparse 
    438                nisparse(jk2_sd) = nisparse(jk2_sd) 
    439                njsparse(jk2_sd) = njsparse(jk2_sd) 
    440                IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
    441                   ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
    442                      &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
     435         DO jh1_sd = 1, nsparse 
     436            DO jh2_sd = 1, nsparse 
     437               nisparse(jh2_sd) = nisparse(jh2_sd) 
     438               njsparse(jh2_sd) = njsparse(jh2_sd) 
     439               IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN 
     440                  ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  & 
     441                     &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd) 
    443442               ENDIF 
    444443            END DO 
     
    515514   END SUBROUTINE SUR_DETERMINE 
    516515 
    517 #else 
    518    !!---------------------------------------------------------------------- 
    519    !!   Default case :   Empty module 
    520    !!---------------------------------------------------------------------- 
    521    LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE. 
    522 CONTAINS 
    523    SUBROUTINE dia_harm ( kt )     ! Empty routine 
    524       INTEGER, INTENT( IN ) :: kt   
    525       WRITE(*,*) 'dia_harm: you should not have seen this print' 
    526    END SUBROUTINE dia_harm 
    527 #endif 
    528  
    529516   !!====================================================================== 
    530517END MODULE diaharm 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahsb.F90

    r10425 r12178  
    362362      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist 
    363363      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901) 
    364 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp ) 
     364901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist' ) 
    365365      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist 
    366366      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 ) 
    367 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp ) 
     367902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist' ) 
    368368      IF(lwm) WRITE( numond, namhsb ) 
    369369 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahth.F90

    r11141 r12178  
    55   !!====================================================================== 
    66   !! History :  OPA  !  1994-09  (J.-P. Boulanger)  Original code 
    7    !!                 !  1996-11  (E. Guilyardi)  OPA8 
     7   !!                 !  1996-11  (E. Guilyardi)  OPA8  
    88   !!                 !  1997-08  (G. Madec)  optimization 
    9    !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content 
     9   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content  
    1010   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
    1111   !!            3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag 
     12   !!---------------------------------------------------------------------- 
     13#if   defined key_diahth 
     14   !!---------------------------------------------------------------------- 
     15   !!   'key_diahth' :                              thermocline depth diag. 
    1216   !!---------------------------------------------------------------------- 
    1317   !!   dia_hth      : Compute varius diagnostics associated with the mixed layer 
     
    2024   USE lib_mpp         ! MPP library 
    2125   USE iom             ! I/O library 
     26   USE timing          ! preformance summary 
    2227 
    2328   IMPLICIT NONE 
     
    2530 
    2631   PUBLIC   dia_hth       ! routine called by step.F90 
    27    PUBLIC   dia_hth_init  ! routine called by nemogcm.F90 
    28  
    29    LOGICAL, PUBLIC ::   ll_diahth   !: Compute further diagnostics of ML and thermocline depth 
     32   PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90 
     33 
     34   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag 
     35    
     36   ! note: following variables should move to local variables once iom_put is always used  
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m] 
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m] 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m] 
     40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W] 
    3041 
    3142   !!---------------------------------------------------------------------- 
     
    3647CONTAINS 
    3748 
     49   FUNCTION dia_hth_alloc() 
     50      !!--------------------------------------------------------------------- 
     51      INTEGER :: dia_hth_alloc 
     52      !!--------------------------------------------------------------------- 
     53      ! 
     54      ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc ) 
     55      ! 
     56      CALL mpp_sum ( 'diahth', dia_hth_alloc ) 
     57      IF(dia_hth_alloc /= 0)   CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' ) 
     58      ! 
     59   END FUNCTION dia_hth_alloc 
     60 
    3861 
    3962   SUBROUTINE dia_hth( kt ) 
    40      !!--------------------------------------------------------------------- 
    41      !!                  ***  ROUTINE dia_hth  *** 
    42      !! 
    43      !! ** Purpose : Computes 
    44      !!      the mixing layer depth (turbocline): avt = 5.e-4 
    45      !!      the depth of strongest vertical temperature gradient 
    46      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
    47      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 
    48      !!      the top of the thermochine: tn = tn(10m) - ztem2 
    49      !!      the pycnocline depth with density criteria equivalent to a temperature variation 
    50      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
    51      !!      the barrier layer thickness 
    52      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
    53      !!      the depth of the 20 degree isotherm (linear interpolation) 
    54      !!      the depth of the 28 degree isotherm (linear interpolation) 
    55      !!      the heat content of first 300 m 
    56      !! 
    57      !! ** Method : 
    58      !!------------------------------------------------------------------- 
    59      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    60      !! 
    61      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
    62      INTEGER                          ::   iid, ilevel           ! temporary integers 
    63      INTEGER, DIMENSION(jpi,jpj) ::   ik20, ik28  ! levels 
    64      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
    65      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
    66      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
    67      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
    68      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
    69      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
    70      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
    71      REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
    72      REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2 
    73      REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3 
    74      REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
    75      REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion 
    76      REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion 
    77      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
    78      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
    79      REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz 
    80      REAL(wp), DIMENSION(jpi,jpj) ::   zthick     ! vertical integration thickness 
    81      REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
    82    ! note: following variables should move to local variables once iom_put is always used 
    83      REAL(wp), DIMENSION(jpi,jpj) ::   zhth    !: depth of the max vertical temperature gradient [m] 
    84      REAL(wp), DIMENSION(jpi,jpj) ::   zhd20   !: depth of 20 C isotherm                         [m] 
    85      REAL(wp), DIMENSION(jpi,jpj) ::   zhd28   !: depth of 28 C isotherm                         [m] 
    86      REAL(wp), DIMENSION(jpi,jpj) ::   zhtc3   !: heat content of first 300 m                    [W] 
    87  
    88      IF (iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1")) THEN 
    89         ! ------------------------------------------------------------- ! 
    90         ! thermocline depth: strongest vertical gradient of temperature ! 
    91         ! turbocline depth (mixing layer depth): avt = zavt5            ! 
    92         ! MLD: rho = rho(1) + zrho3                                     ! 
    93         ! MLD: rho = rho(1) + zrho1                                     ! 
    94         ! ------------------------------------------------------------- ! 
    95         zmaxdzT(:,:) = 0._wp 
    96         IF( nla10 > 1 ) THEN 
    97            DO jj = 1, jpj 
    98               DO ji = 1, jpi 
    99                  zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
    100                  zrho0_3(ji,jj) = zztmp 
    101                  zrho0_1(ji,jj) = zztmp 
    102                  zhth(ji,jj) = zztmp 
    103               END DO 
    104            END DO 
    105         ELSE IF (iom_use("mlddzt")) THEN 
    106            DO jj = 1, jpj 
    107               DO ji = 1, jpi 
    108                  zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
    109                  zhth(ji,jj) = zztmp 
    110               END DO 
    111            END DO 
    112         ELSE 
    113            zhth(:,:) = 0._wp 
    114  
    115         ENDIF 
    116  
    117         DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
    118            DO jj = 1, jpj 
    119               DO ji = 1, jpi 
    120                  ! 
    121                  zzdep = gdepw_n(ji,jj,jk) 
    122                  zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
    123                  zzdep = zzdep * tmask(ji,jj,1) 
    124  
    125                  IF( zztmp > zmaxdzT(ji,jj) ) THEN 
    126                     zmaxdzT(ji,jj) = zztmp   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz 
    127                  ENDIF 
    128  
    129                  IF( nla10 > 1 ) THEN 
    130                     zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
    131                     IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
    132                     IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
    133                  ENDIF 
    134  
    135               END DO 
    136            END DO 
     63      !!--------------------------------------------------------------------- 
     64      !!                  ***  ROUTINE dia_hth  *** 
     65      !! 
     66      !! ** Purpose : Computes 
     67      !!      the mixing layer depth (turbocline): avt = 5.e-4 
     68      !!      the depth of strongest vertical temperature gradient 
     69      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
     70      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2        
     71      !!      the top of the thermochine: tn = tn(10m) - ztem2  
     72      !!      the pycnocline depth with density criteria equivalent to a temperature variation  
     73      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)  
     74      !!      the barrier layer thickness 
     75      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
     76      !!      the depth of the 20 degree isotherm (linear interpolation) 
     77      !!      the depth of the 28 degree isotherm (linear interpolation) 
     78      !!      the heat content of first 300 m 
     79      !! 
     80      !! ** Method :  
     81      !!------------------------------------------------------------------- 
     82      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     83      !! 
     84      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
     85      INTEGER                          ::   iid, ilevel           ! temporary integers 
     86      INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels 
     87      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
     88      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
     89      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
     90      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
     91      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
     92      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
     93      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
     94      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2  
     95      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2      
     96      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3       
     97      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
     98      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion 
     99      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion 
     100      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
     101      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
     102      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz 
     103      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness  
     104      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
     105      !!---------------------------------------------------------------------- 
     106      IF( ln_timing )   CALL timing_start('dia_hth') 
     107 
     108      IF( kt == nit000 ) THEN 
     109         !                                      ! allocate dia_hth array 
     110         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 
     111 
     112         IF(.NOT. ALLOCATED(ik20) ) THEN 
     113            ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 
     114               &      zabs2(jpi,jpj),   & 
     115               &      ztm2(jpi,jpj),    & 
     116               &      zrho10_3(jpi,jpj),& 
     117               &      zpycn(jpi,jpj),   & 
     118               &      ztinv(jpi,jpj),   & 
     119               &      zdepinv(jpi,jpj), & 
     120               &      zrho0_3(jpi,jpj), & 
     121               &      zrho0_1(jpi,jpj), & 
     122               &      zmaxdzT(jpi,jpj), & 
     123               &      zthick(jpi,jpj),  & 
     124               &      zdelr(jpi,jpj), STAT=ji) 
     125            CALL mpp_sum('diahth', ji) 
     126            IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' ) 
     127         END IF 
     128 
     129         IF(lwp) WRITE(numout,*) 
     130         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 
     131         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     132         IF(lwp) WRITE(numout,*) 
     133      ENDIF 
     134 
     135      ! initialization 
     136      ztinv  (:,:) = 0._wp   
     137      zdepinv(:,:) = 0._wp   
     138      zmaxdzT(:,:) = 0._wp   
     139      DO jj = 1, jpj 
     140         DO ji = 1, jpi 
     141            zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
     142            hth     (ji,jj) = zztmp 
     143            zabs2   (ji,jj) = zztmp 
     144            ztm2    (ji,jj) = zztmp 
     145            zrho10_3(ji,jj) = zztmp 
     146            zpycn   (ji,jj) = zztmp 
    137147        END DO 
    138  
    139         IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline 
    140         IF( nla10 > 1 ) THEN 
    141            IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.03 
    142            IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.01 
    143         ENDIF 
    144      ENDIF 
    145  
    146      IF (iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  & 
    147           & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti")) THEN 
    148         DO jj = 1, jpj 
    149            DO ji = 1, jpi 
    150               zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
    151               zabs2   (ji,jj) = zztmp 
    152               ztm2    (ji,jj) = zztmp 
    153               zrho10_3(ji,jj) = zztmp 
    154               zpycn   (ji,jj) = zztmp 
    155            END DO 
    156         END DO 
    157         ztinv  (:,:) = 0._wp 
    158         zdepinv(:,:) = 0._wp 
    159  
    160         IF (iom_use("pycndep")) THEN 
    161            ! Preliminary computation 
    162            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    163            DO jj = 1, jpj 
    164               DO ji = 1, jpi 
    165                  IF( tmask(ji,jj,nla10) == 1. ) THEN 
    166                     zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             & 
    167                          &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
    168                          &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
    169                     zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             & 
    170                          &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
    171                     zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
    172                     zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
    173                     zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
    174                     zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
    175                  ELSE 
    176                     zdelr(ji,jj) = 0._wp 
    177                  ENDIF 
    178               END DO 
    179            END DO 
    180         ELSE 
    181            zdelr(:,:) = 0._wp 
    182         ENDIF 
    183  
    184         ! ------------------------------------------------------------- ! 
    185         ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
    186         ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
    187         ! MLD: rho = rho10m + zrho3                                     ! 
    188         ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
    189         ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
    190         ! depth of temperature inversion                                ! 
    191         ! ------------------------------------------------------------- ! 
    192         DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
    193            DO jj = 1, jpj 
    194               DO ji = 1, jpi 
    195                  ! 
    196                  zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
    197                  ! 
    198                  zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
    199                  IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
    200                  IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
    201                  zztmp = -zztmp                                          ! delta T(10m) 
    202                  IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
    203                     ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
    204                  ENDIF 
    205  
    206                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
    207                  IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
    208                  IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
    209                  ! 
    210               END DO 
    211            END DO 
    212         END DO 
    213  
    214         IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1)        )   ! MLD abs(delta t) - 0.2 
    215         IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1)         )   ! T(10) - 0.2 
    216         IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1)     )   ! MLD delta rho(10m) = 0.03 
    217         IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1)        )   ! MLD delta rho equi. delta T(10m) = 0.2 
    218         IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref) 
    219         IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref) 
    220      ENDIF 
    221  
    222      IF(iom_use("20d") .OR. iom_use("28d")) THEN 
    223         ! ----------------------------------- ! 
    224         ! search deepest level above 20C/28C  ! 
    225         ! ----------------------------------- ! 
    226         ik20(:,:) = 1 
    227         ik28(:,:) = 1 
    228         DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
    229            DO jj = 1, jpj 
    230               DO ji = 1, jpi 
    231                  zztmp = tsn(ji,jj,jk,jp_tem) 
    232                  IF( zztmp >= 20. )   ik20(ji,jj) = jk 
    233                  IF( zztmp >= 28. )   ik28(ji,jj) = jk 
    234               END DO 
    235            END DO 
    236         END DO 
    237  
    238         ! --------------------------- ! 
    239         !  Depth of 20C/28C isotherm  ! 
    240         ! --------------------------- ! 
    241         DO jj = 1, jpj 
    242            DO ji = 1, jpi 
    243               ! 
    244               zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
    245               ! 
    246               iid = ik20(ji,jj) 
    247               IF( iid /= 1 ) THEN 
    248                  zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    249                       &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    250                       &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
    251                       &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    252                  zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    253               ELSE 
    254                  zhd20(ji,jj) = 0._wp 
    255               ENDIF 
    256               ! 
    257               iid = ik28(ji,jj) 
    258               IF( iid /= 1 ) THEN 
    259                  zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    260                       &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    261                       &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   & 
    262                       &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    263                  zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
    264               ELSE 
    265                  zhd28(ji,jj) = 0._wp 
    266               ENDIF 
    267  
    268            END DO 
    269         END DO 
    270         CALL iom_put( "20d", zhd20 )   ! depth of the 20 isotherm 
    271         CALL iom_put( "28d", zhd28 )   ! depth of the 28 isotherm 
    272      ENDIF 
    273  
    274      ! ----------------------------- ! 
    275      !  Heat content of first 300 m  ! 
    276      ! ----------------------------- ! 
    277      IF (iom_use("hc300")) THEN 
    278         ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 
    279         ilevel   = 0 
    280         zthick_0 = 0._wp 
    281         DO jk = 1, jpkm1 
    282            zthick_0 = zthick_0 + e3t_1d(jk) 
    283            IF( zthick_0 < 300. )   ilevel = jk 
    284         END DO 
    285         ! surface boundary condition 
    286         IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 
    287         ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp 
    288         ENDIF 
    289         ! integration down to ilevel 
    290         DO jk = 1, ilevel 
    291            zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 
    292            zhtc3  (:,:) = zhtc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 
    293         END DO 
    294         ! deepest layer 
    295         zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m 
    296         DO jj = 1, jpj 
    297            DO ji = 1, jpi 
    298               zhtc3(ji,jj) = zhtc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  & 
    299                    &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 
    300            END DO 
    301         END DO 
    302         ! from temperature to heat contain 
    303         zcoef = rau0 * rcp 
    304         zhtc3(:,:) = zcoef * zhtc3(:,:) 
    305         CALL iom_put( "hc300", zhtc3*tmask(:,:,1) )      ! first 300m heat content 
    306      ENDIF 
    307      ! 
     148      END DO 
     149      IF( nla10 > 1 ) THEN  
     150         DO jj = 1, jpj 
     151            DO ji = 1, jpi 
     152               zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
     153               zrho0_3(ji,jj) = zztmp 
     154               zrho0_1(ji,jj) = zztmp 
     155            END DO 
     156         END DO 
     157      ENDIF 
     158       
     159      ! Preliminary computation 
     160      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
     161      DO jj = 1, jpj 
     162         DO ji = 1, jpi 
     163            IF( tmask(ji,jj,nla10) == 1. ) THEN 
     164               zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             & 
     165                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
     166                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
     167               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             & 
     168                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
     169               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
     170               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
     171               zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
     172               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
     173            ELSE 
     174               zdelr(ji,jj) = 0._wp 
     175            ENDIF 
     176         END DO 
     177      END DO 
     178 
     179      ! ------------------------------------------------------------- ! 
     180      ! thermocline depth: strongest vertical gradient of temperature ! 
     181      ! turbocline depth (mixing layer depth): avt = zavt5            ! 
     182      ! MLD: rho = rho(1) + zrho3                                     ! 
     183      ! MLD: rho = rho(1) + zrho1                                     ! 
     184      ! ------------------------------------------------------------- ! 
     185      DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
     186         DO jj = 1, jpj 
     187            DO ji = 1, jpi 
     188               ! 
     189               zzdep = gdepw_n(ji,jj,jk) 
     190               zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
     191               zzdep = zzdep * tmask(ji,jj,1) 
     192 
     193               IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
     194                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
     195               ENDIF 
     196                
     197               IF( nla10 > 1 ) THEN  
     198                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
     199                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
     200                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
     201               ENDIF 
     202 
     203            END DO 
     204         END DO 
     205      END DO 
     206       
     207      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline 
     208      IF( nla10 > 1 ) THEN  
     209         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03 
     210         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01 
     211      ENDIF 
     212 
     213      ! ------------------------------------------------------------- ! 
     214      ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
     215      ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
     216      ! MLD: rho = rho10m + zrho3                                     ! 
     217      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
     218      ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
     219      ! depth of temperature inversion                                ! 
     220      ! ------------------------------------------------------------- ! 
     221      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
     222         DO jj = 1, jpj 
     223            DO ji = 1, jpi 
     224               ! 
     225               zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
     226               ! 
     227               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
     228               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
     229               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
     230               zztmp = -zztmp                                          ! delta T(10m) 
     231               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
     232                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
     233               ENDIF 
     234 
     235               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
     236               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
     237               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
     238               ! 
     239            END DO 
     240         END DO 
     241      END DO 
     242 
     243      CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2 
     244      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2 
     245      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03 
     246      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2 
     247      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)  
     248      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)  
     249 
     250 
     251      ! ----------------------------------- ! 
     252      ! search deepest level above 20C/28C  ! 
     253      ! ----------------------------------- ! 
     254      ik20(:,:) = 1 
     255      ik28(:,:) = 1 
     256      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
     257         DO jj = 1, jpj 
     258            DO ji = 1, jpi 
     259               zztmp = tsn(ji,jj,jk,jp_tem) 
     260               IF( zztmp >= 20. )   ik20(ji,jj) = jk 
     261               IF( zztmp >= 28. )   ik28(ji,jj) = jk 
     262            END DO 
     263         END DO 
     264      END DO 
     265 
     266      ! --------------------------- ! 
     267      !  Depth of 20C/28C isotherm  ! 
     268      ! --------------------------- ! 
     269      DO jj = 1, jpj 
     270         DO ji = 1, jpi 
     271            ! 
     272            zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
     273            ! 
     274            iid = ik20(ji,jj) 
     275            IF( iid /= 1 ) THEN  
     276               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
     277                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     278                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
     279                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
     280               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
     281            ELSE  
     282               hd20(ji,jj) = 0._wp 
     283            ENDIF 
     284            ! 
     285            iid = ik28(ji,jj) 
     286            IF( iid /= 1 ) THEN  
     287               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
     288                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     289                  &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   & 
     290                  &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
     291               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
     292            ELSE  
     293               hd28(ji,jj) = 0._wp 
     294            ENDIF 
     295 
     296         END DO 
     297      END DO 
     298      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm 
     299      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm 
     300 
     301      ! ----------------------------- ! 
     302      !  Heat content of first 300 m  ! 
     303      ! ----------------------------- ! 
     304 
     305      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 
     306      ilevel   = 0 
     307      zthick_0 = 0._wp 
     308      DO jk = 1, jpkm1                       
     309         zthick_0 = zthick_0 + e3t_1d(jk) 
     310         IF( zthick_0 < 300. )   ilevel = jk 
     311      END DO 
     312      ! surface boundary condition 
     313      IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
     314      ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                    
     315      ENDIF 
     316      ! integration down to ilevel 
     317      DO jk = 1, ilevel 
     318         zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 
     319         htc3  (:,:) = htc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 
     320      END DO 
     321      ! deepest layer 
     322      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m 
     323      DO jj = 1, jpj 
     324         DO ji = 1, jpi 
     325            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  & 
     326               &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 
     327         END DO 
     328      END DO 
     329      ! from temperature to heat contain 
     330      zcoef = rau0 * rcp 
     331      htc3(:,:) = zcoef * htc3(:,:) 
     332      CALL iom_put( "hc300", htc3 )      ! first 300m heat content 
     333      ! 
     334      IF( ln_timing )   CALL timing_stop('dia_hth') 
     335      ! 
    308336   END SUBROUTINE dia_hth 
    309337 
    310  
    311    SUBROUTINE dia_hth_init 
    312       !!--------------------------------------------------------------------------- 
    313       !!                  ***  ROUTINE dia_hth_init  *** 
    314       !! 
    315       !! ** Purpose: Initialization for ML and thermocline depths 
    316       !! 
    317       !! ** Action : If any upper ocean diagnostic required by xml file, set in dia_hth 
    318       !!--------------------------------------------------------------------------- 
    319        ! 
    320       IF(lwp) THEN 
    321          WRITE(numout,*) 
    322          WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics' 
    323          WRITE(numout,*) '~~~~~~~~~~~~ ' 
    324       ENDIF 
    325       ll_diahth = iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1") .OR.  & 
    326            & iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  & 
    327            & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti").OR. & 
    328            & iom_use("20d") .OR. iom_use("28d") .OR. iom_use("hc300") 
    329       IF(lwp) THEN 
    330          WRITE(numout,*) '   output upper ocean diagnostics (T) or not (F)       ll_diahth = ', ll_diahth 
    331       ENDIF 
    332       ! 
    333    END SUBROUTINE dia_hth_init 
     338#else 
     339   !!---------------------------------------------------------------------- 
     340   !!   Default option :                                       Empty module 
     341   !!---------------------------------------------------------------------- 
     342   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag 
     343CONTAINS 
     344   SUBROUTINE dia_hth( kt )         ! Empty routine 
     345      IMPLICIT NONE 
     346      INTEGER, INTENT( in ) :: kt 
     347      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 
     348   END SUBROUTINE dia_hth 
     349#endif 
     350 
     351   !!====================================================================== 
    334352END MODULE diahth 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaptr.F90

    r10425 r12178  
    393393      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport 
    394394      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901) 
    395 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp ) 
     395901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist' ) 
    396396 
    397397      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport 
    398398      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 ) 
    399 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp ) 
     399902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist' ) 
    400400      IF(lwm) WRITE ( numond, namptr ) 
    401401 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diatmb.F90

    r10499 r12178  
    4343      REWIND( numnam_ref )              ! Read Namelist nam_diatmb in reference namelist : TMB diagnostics 
    4444      READ  ( numnam_ref, nam_diatmb, IOSTAT=ios, ERR= 901 ) 
    45 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist', lwp ) 
     45901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist' ) 
    4646  
    4747      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics 
    4848      READ  ( numnam_cfg, nam_diatmb, IOSTAT = ios, ERR = 902 ) 
    49 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist', lwp ) 
     49902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist' ) 
    5050      IF(lwm) WRITE ( numond, nam_diatmb ) 
    5151 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90

    r11143 r12178  
    4343   USE zdfdrg         ! ocean vertical physics: top/bottom friction 
    4444   USE zdfmxl         ! mixed layer 
    45    USE zdfosm         ! mixed layer 
    4645   ! 
    4746   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    211210      ENDIF 
    212211 
     212      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
     213      ! 
    213214      CALL iom_put( "woce", wn )                   ! vertical velocity 
    214215      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
     
    221222         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    222223      ENDIF 
     224      ! 
     225      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output 
    223226 
    224227      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
     
    427430      !!      define all the NETCDF files and fields 
    428431      !!      At each time step call histdef to compute the mean if ncessary 
    429       !!      Each nwrite time step, output the instantaneous or mean fields 
     432      !!      Each nn_write time step, output the instantaneous or mean fields 
    430433      !!---------------------------------------------------------------------- 
    431434      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     
    443446      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
    444447      !!---------------------------------------------------------------------- 
    445       !  
    446       IF( ln_timing )   CALL timing_start('dia_wri') 
    447448      ! 
    448449      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==! 
     
    451452      ENDIF 
    452453      ! 
     454      IF( nn_write == -1 )   RETURN   ! we will never do any output 
     455      !  
     456      IF( ln_timing )   CALL timing_start('dia_wri') 
     457      ! 
    453458      ! 0. Initialisation 
    454459      ! ----------------- 
     
    460465      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes) 
    461466#if defined key_diainstant 
    462       zsto = nwrite * rdt 
     467      zsto = nn_write * rdt 
    463468      clop = "inst("//TRIM(clop)//")" 
    464469#else 
     
    466471      clop = "ave("//TRIM(clop)//")" 
    467472#endif 
    468       zout = nwrite * rdt 
     473      zout = nn_write * rdt 
    469474      zmax = ( nitend - nit000 + 1 ) * rdt 
    470475 
     
    497502         ! WRITE root name in date.file for use by postpro 
    498503         IF(lwp) THEN 
    499             CALL dia_nam( clhstnam, nwrite,' ' ) 
     504            CALL dia_nam( clhstnam, nn_write,' ' ) 
    500505            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    501506            WRITE(inum,*) clhstnam 
     
    505510         ! Define the T grid FILE ( nid_T ) 
    506511 
    507          CALL dia_nam( clhstnam, nwrite, 'grid_T' ) 
     512         CALL dia_nam( clhstnam, nn_write, 'grid_T' ) 
    508513         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
    509514         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     
    541546         ! Define the U grid FILE ( nid_U ) 
    542547 
    543          CALL dia_nam( clhstnam, nwrite, 'grid_U' ) 
     548         CALL dia_nam( clhstnam, nn_write, 'grid_U' ) 
    544549         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
    545550         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu 
     
    554559         ! Define the V grid FILE ( nid_V ) 
    555560 
    556          CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename 
     561         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename 
    557562         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam 
    558563         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv 
     
    567572         ! Define the W grid FILE ( nid_W ) 
    568573 
    569          CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename 
     574         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename 
    570575         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam 
    571576         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     
    658663         ENDIF 
    659664 
    660          IF( .NOT. ln_cpl ) THEN 
     665         IF( ln_ssr ) THEN 
    661666            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    662667               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    666671               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    667672         ENDIF 
    668  
    669          IF( ln_cpl .AND. nn_ice <= 1 ) THEN 
    670             CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    671                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    672             CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp 
    673                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    674             CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn 
    675                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    676          ENDIF 
    677           
     673        
    678674         clmx ="l_max(only(x))"    ! max index on a period 
    679675!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX  
     
    751747      ! donne le nombre d'elements, et ndex la liste des indices a sortir 
    752748 
    753       IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN  
     749      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN  
    754750         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step' 
    755751         WRITE(numout,*) '~~~~~~ ' 
     
    815811      ENDIF 
    816812 
    817       IF( .NOT. ln_cpl ) THEN 
     813      IF( ln_ssr ) THEN 
    818814         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    819815         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    820          IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
    821          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    822       ENDIF 
    823       IF( ln_cpl .AND. nn_ice <= 1 ) THEN 
    824          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    825          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    826          IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
     816         zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
    827817         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    828818      ENDIF 
     
    843833      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    844834 
    845       CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current 
     835      IF( ln_zad_Aimp ) THEN 
     836         CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current 
     837      ELSE 
     838         CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current 
     839      ENDIF 
    846840      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
    847841      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
     
    904898      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity 
    905899      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity 
    906       CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity 
     900      IF( ln_zad_Aimp ) THEN 
     901         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity 
     902      ELSE 
     903         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity 
     904      ENDIF 
    907905      IF( ALLOCATED(ahtu) ) THEN 
    908906         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
     
    928926         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity 
    929927      ENDIF 
    930  
    931       IF( ln_zdfosm ) THEN 
    932          CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)   )    ! now boundary-layer depth 
    933          CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)    )    ! now mixed-layer depth 
    934          CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask       )    ! w-level diffusion 
    935          CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask       )    ! now w-level viscosity 
    936          CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask       )    ! non-local t forcing 
    937          CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask       )    ! non-local s forcing 
    938          CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask       )    ! non-local u forcing 
    939          CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask       )    ! non-local v forcing 
    940       ENDIF 
    941      !     ! CALL histwrite( id_i, "zla", kt, zla*tmask(:,:,1)   , jpi*jpj, idex)         ! now Langmuir # 
    942      !     ! CALL histwrite( id_i, "zvstr", kt, zvstr*tmask(:,:,1)   , jpi*jpj, idex)     ! now mixed velocity scale 
    943      !     ! CALL histwrite( id_i, "zustar", kt, zustar*tmask(:,:,1)   , jpi*jpj, idex)   ! now friction velocity scale 
    944      !     ! CALL histwrite( id_i, "zwstrl", kt, zwstrl*tmask(:,:,1)   , jpi*jpj, idex)   ! now Langmuir velocity scale 
    945      !     ! CALL histwrite( id_i, "zwstrc", kt, zwstrc*tmask(:,:,1)   , jpi*jpj, idex)   ! now convective velocity scale 
    946      !     ! CALL histwrite( id_i, "zwb_ent", kt, zwb_ent*tmask(:,:,1)   , jpi*jpj, idex) ! now upward turb buoyancy entrainment flux 
    947      !     ! CALL histwrite( id_i, "zdb_bl", kt, zdb_bl*tmask(:,:,1)   , jpi*jpj, idex)   ! now db at ml base 
    948928  
    949929#if defined key_si3 
Note: See TracChangeset for help on using the changeset viewer.