! $Id: iniacademic.F90 5186 2024-09-11 16:03:07Z abarral $ SUBROUTINE iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0) USE lmdz_filtreg, ONLY: inifilr USE lmdz_infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName USE control_mod, ONLY: day_step, planet_type USE exner_hyb_m, ONLY: exner_hyb USE exner_milieu_m, ONLY: exner_milieu USE IOIPSL, ONLY: getin USE lmdz_write_field USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm USE logic_mod, ONLY: iflag_phys, read_start USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner USE temps_mod, ONLY: annee_ref, day_ini, day_ref USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0 USE lmdz_readTracFiles, ONLY: addPhase USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_noerr, nf90_inq_varid, nf90_close, nf90_get_var USE lmdz_ran1, ONLY: ran1 USE lmdz_iniprint, ONLY: lunout, prt_level USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4 USE lmdz_comgeom ! Author: Frederic Hourdin original: 15/01/93 ! The forcing defined here is from Held and Suarez, 1994, Bulletin ! of the American Meteorological Society, 75, 1825. USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet USE lmdz_check_isotopes, ONLY: check_isotopes_seq IMPLICIT NONE ! Declararations: ! --------------- ! Arguments: ! ---------- REAL, INTENT(OUT) :: time_0 ! fields REAL, INTENT(OUT) :: vcov(ip1jm, llm) ! meridional covariant wind REAL, INTENT(OUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind REAL, INTENT(OUT) :: teta(ip1jmp1, llm) ! potential temperature (K) REAL, INTENT(OUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers (.../kg_of_air) REAL, INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa) REAL, INTENT(OUT) :: masse(ip1jmp1, llm) ! air mass in grid cell (kg) REAL, INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential ! Local: ! ------ REAL p (ip1jmp1, llmp1) ! pression aux interfac.des couches REAL pks(ip1jmp1) ! exner au sol REAL pk(ip1jmp1, llm) ! exner au milieu des couches REAL phi(ip1jmp1, llm) ! geopotentiel REAL ddsin, zsig, tetapv, w_pv ! variables auxiliaires REAL tetastrat ! potential temperature in the stratosphere, in K REAL tetajl(jjp1, llm) INTEGER i, j, l, lsup, ij, iq, iName, iPhase, iqParent INTEGER :: nid_relief, varid, ierr REAL, DIMENSION(iip1, jjp1) :: relief REAL teta0, ttp, delt_y, delt_z, eps ! Constantes pour profil de T REAL k_f, k_c_a, k_c_s ! Constantes de rappel LOGICAL ok_geost ! Initialisation vent geost. ou nul LOGICAL ok_pv ! Polar Vortex REAL phi_pv, dphi_pv, gam_pv, tetanoise ! Constantes pour polar vortex REAL zz INTEGER idum REAL zdtvr, tnat, alpha_ideal LOGICAL, PARAMETER :: tnat1 = .TRUE. CHARACTER(LEN = *), parameter :: modname = "iniacademic" CHARACTER(LEN = 80) :: abort_message ! Sanity check: verify that options selected by user are not incompatible IF ((iflag_phys==1).AND. .NOT. read_start) THEN WRITE(lunout, *) trim(modname), " error: if read_start is set to ", & " false then iflag_phys should not be 1" WRITE(lunout, *) "You most likely want an aquaplanet initialisation", & " (iflag_phys >= 100)" CALL abort_gcm(modname, "incompatible iflag_phys==1 and read_start==.FALSE.", 1) ENDIF !----------------------------------------------------------------------- ! 1. Initializations for Earth-like case ! -------------------------------------- ! initialize planet radius, rotation rate,... CALL conf_planete time_0 = 0. day_ref = 1 annee_ref = 0 im = iim jm = jjm day_ini = 1 dtvr = daysec / REAL(day_step) zdtvr = dtvr etot0 = 0. ptot0 = 0. ztot0 = 0. stot0 = 0. ang0 = 0. IF (llm == 1) THEN ! specific initializations for the shallow water case kappa = 1 ENDIF CALL iniconst CALL inigeom CALL inifilr !------------------------------------------------------------------ ! Initialize pressure and mass field if read_start=.FALSE. !------------------------------------------------------------------ IF (.NOT. read_start) THEN !------------------------------------------------------------------ ! Lecture eventuelle d'un fichier de relief interpollee sur la grille ! du modele. ! On suppose que le fichier relief_in.nc est stoké sur une grille ! iim*jjp1 ! Facile a créer à partir de la commande ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc !------------------------------------------------------------------ relief = 0. ierr = nf90_open ('relief_in.nc', nf90_nowrite, nid_relief) IF (ierr==nf90_noerr) THEN ierr = nf90_inq_varid(nid_relief, 'RELIEF', varid) IF (ierr==nf90_noerr) THEN ierr = nf90_get_var(nid_relief, varid, relief(1:iim, 1:jjp1)) relief(iip1, :) = relief(1, :) else CALL abort_gcm ('iniacademic', 'variable RELIEF pas la', 1) endif endif ierr = nf90_close (nid_relief) !------------------------------------------------------------------ ! Initialisation du geopotentiel au sol et de la pression !------------------------------------------------------------------ PRINT*, 'relief=', minval(relief), maxval(relief), 'g=', g DO j = 1, jjp1 DO i = 1, iip1 phis((j - 1) * iip1 + i) = g * relief(i, j) enddo enddo PRINT*, 'phis=', minval(phis), maxval(phis), 'g=', g ! ground geopotential !phis(:)=0. ps(:) = preff CALL pression (ip1jmp1, ap, bp, ps, p) IF (pressure_exner) THEN CALL exner_hyb(ip1jmp1, ps, p, pks, pk) else CALL exner_milieu(ip1jmp1, ps, p, pks, pk) endif CALL massdair(p, masse) ENDIF IF (llm == 1) THEN ! initialize fields for the shallow water case, if required IF (.NOT.read_start) THEN phis(:) = 0. q(:, :, :) = 0 CALL sw_case_williamson91_6(vcov, ucov, teta, masse, ps) endif ENDIF academic_case: if (iflag_phys >= 2) THEN ! initializations ! 1. local parameters ! by convention, winter is in the southern hemisphere ! Geostrophic wind or no wind? ok_geost = .TRUE. CALL getin('ok_geost', ok_geost) ! Constants for Newtonian relaxation and friction k_f = 1. !friction CALL getin('k_j', k_f) k_f = 1. / (daysec * k_f) k_c_s = 4. !cooling surface CALL getin('k_c_s', k_c_s) k_c_s = 1. / (daysec * k_c_s) k_c_a = 40. !cooling free atm CALL getin('k_c_a', k_c_a) k_c_a = 1. / (daysec * k_c_a) ! Constants for Teta equilibrium profile teta0 = 315. ! mean Teta (S.H. 315K) CALL getin('teta0', teta0) ttp = 200. ! Tropopause temperature (S.H. 200K) CALL getin('ttp', ttp) eps = 0. ! Deviation to N-S symmetry(~0-20K) CALL getin('eps', eps) delt_y = 60. ! Merid Temp. Gradient (S.H. 60K) CALL getin('delt_y', delt_y) delt_z = 10. ! Vertical Gradient (S.H. 10K) CALL getin('delt_z', delt_z) ! Polar vortex ok_pv = .FALSE. CALL getin('ok_pv', ok_pv) phi_pv = -50. ! Latitude of edge of vortex CALL getin('phi_pv', phi_pv) phi_pv = phi_pv * pi / 180. dphi_pv = 5. ! Width of the edge CALL getin('dphi_pv', dphi_pv) dphi_pv = dphi_pv * pi / 180. gam_pv = 4. ! -dT/dz vortex (in K/km) CALL getin('gam_pv', gam_pv) tetanoise = 0.005 CALL getin('tetanoise', tetanoise) ! 2. Initialize fields towards which to relax ! Friction knewt_g = k_c_a DO l = 1, llm zsig = presnivs(l) / preff knewt_t(l) = (k_c_s - k_c_a) * MAX(0., (zsig - 0.7) / 0.3) kfrict(l) = k_f * MAX(0., (zsig - 0.7) / 0.3) ENDDO DO j = 1, jjp1 clat4((j - 1) * iip1 + 1:j * iip1) = cos(rlatu(j))**4 ENDDO ! Potential temperature DO l = 1, llm zsig = presnivs(l) / preff tetastrat = ttp * zsig**(-kappa) tetapv = tetastrat IF ((ok_pv).AND.(zsig<0.1)) THEN tetapv = tetastrat * (zsig * 10.)**(kappa * cpp * gam_pv / 1000. / g) ENDIF DO j = 1, jjp1 ! Troposphere ddsin = sin(rlatu(j)) tetajl(j, l) = teta0 - delt_y * ddsin * ddsin + eps * ddsin & - delt_z * (1. - ddsin * ddsin) * log(zsig) IF (planet_type=="giant") THEN tetajl(j, l) = teta0 + (delt_y * & ((sin(rlatu(j) * 3.14159 * eps + 0.0001))**2) & / ((rlatu(j) * 3.14159 * eps + 0.0001)**2)) & - delt_z * log(zsig) endif ! Profil stratospherique isotherme (+vortex) w_pv = (1. - tanh((rlatu(j) - phi_pv) / dphi_pv)) / 2. tetastrat = tetastrat * (1. - w_pv) + tetapv * w_pv tetajl(j, l) = MAX(tetajl(j, l), tetastrat) ENDDO ENDDO ! CALL writefield('theta_eq',tetajl) DO l = 1, llm DO j = 1, jjp1 DO i = 1, iip1 ij = (j - 1) * iip1 + i tetarappel(ij, l) = tetajl(j, l) enddo enddo enddo ! 3. Initialize fields (if necessary) IF (.NOT. read_start) THEN ! bulk initialization of temperature IF (iflag_phys>10000) THEN ! Particular case to impose a constant temperature T0=0.01*iflag_physx teta(:, :) = 0.01 * iflag_phys / (pk(:, :) / cpp) ELSE teta(:, :) = tetarappel(:, :) ENDIF ! geopotential CALL geopot(ip1jmp1, teta, pk, pks, phis, phi) DO l = 1, llm PRINT*, 'presnivs,play,l', presnivs(l), (pk(1, l) / cpp)**(1. / kappa) * preff !pks(ij) = (cpp/preff) * ps(ij) !pk(ij,1) = .5*pks(ij) ! pk = cpp * (p/preff)^kappa ENDDO ! winds IF (ok_geost) THEN CALL ugeostr(phi, ucov) else ucov(:, :) = 0. endif vcov(:, :) = 0. ! bulk initialization of tracers IF (planet_type=="earth") THEN ! Earth: first two tracers will be water DO iq = 1, nqtot q(:, :, iq) = 0. IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:, :, iq) = 1.e-10 IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:, :, iq) = 1.e-15 ! CRisi: init des isotopes ! distill de Rayleigh très simplifiée iName = tracers(iq)%iso_iName IF (niso <= 0 .OR. iName <= 0) CYCLE iPhase = tracers(iq)%iso_iPhase iqParent = tracers(iq)%iqParent IF(tracers(iq)%iso_iZone == 0) THEN IF (tnat1) THEN tnat = 1.0 alpha_ideal = 1.0 WRITE(*, *) 'Attention dans iniacademic: alpha_ideal=1' else IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) endif q(:, :, iq) = q(:, :, iqParent) * tnat * (q(:, :, iqParent) / 30.e-3)**(alpha_ideal - 1.) ELSE !IF(tracers(iq)%iso_iZone == 0) THEN IF(tracers(iq)%iso_iZone == 1) THEN ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1. ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs. q(:, :, iq) = q(:, :, iqIsoPha(iName, iPhase)) else !IF(tracers(iq)%iso_iZone == 1) THEN q(:, :, iq) = 0. endif !IF(tracers(iq)%iso_iZone == 1) THEN END IF !IF(tracers(iq)%iso_iZone == 0) THEN enddo else q(:, :, :) = 0 endif ! of if (planet_type=="earth") CALL check_isotopes_seq(q, ip1jmp1, 'iniacademic_loc') ! add random perturbation to temperature idum = -1 zz = ran1(idum) idum = 0 DO l = 1, llm DO ij = iip2, ip1jm teta(ij, l) = teta(ij, l) * (1. + tetanoise * ran1(idum)) enddo enddo ! maintain periodicity in longitude DO l = 1, llm DO ij = 1, ip1jmp1, iip1 teta(ij + iim, l) = teta(ij, l) enddo enddo ENDIF ! of IF (.NOT. read_start) ENDIF academic_case END SUBROUTINE iniacademic