source: LMDZ6/trunk/libf/phylmd/phyaqua_mod.F90 @ 5169

Last change on this file since 5169 was 5084, checked in by Laurent Fairhead, 12 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Id
File size: 31.8 KB
Line 
1!
2! $Id: phyaqua_mod.F90 5084 2024-07-19 16:40:44Z idelkadi $
3!
4MODULE phyaqua_mod
5  ! Routines complementaires pour la physique planetaire.
6  IMPLICIT NONE
7
8CONTAINS
9
10  SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
11
12    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13    ! Creation d'un etat initial et de conditions aux limites
14    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
15    ! du modele LMDZ dans sa version terrestre.
16    ! iflag_phys est un parametre qui controle
17    ! iflag_phys = N
18    ! de 100 a 199 : aqua planetes avec SST forcees
19    ! N-100 determine le type de SSTs
20    ! de 200 a 299 : terra planetes avec Ts calcule
21
22    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24    USE dimphy, ONLY: klon
25    USE geometry_mod, ONLY : latitude
26    USE surface_data, ONLY: type_ocean, ok_veget
27    USE pbl_surface_mod, ONLY: pbl_surface_init
28    USE fonte_neige_mod, ONLY: fonte_neige_init
29    USE phys_state_var_mod
30    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
31                                day_ini,day_end
32    USE indice_sol_mod
33    USE nrtype, ONLY: pi
34!    USE ioipsl
35    USE mod_phys_lmdz_para, ONLY: is_master
36    USE mod_phys_lmdz_transfert_para, ONLY: bcast
37    USE mod_grid_phy_lmdz
38    USE ioipsl_getin_p_mod, ONLY : getin_p
39    USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
40    IMPLICIT NONE
41
42    include "YOMCST.h"
43    include "clesphys.h"
44    include "dimsoil.h"
45
46    INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
47    ! IM ajout latfi, lonfi
48!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
49
50    INTEGER type_profil, type_aqua
51
52    ! Ajouts initialisation des surfaces
53    REAL :: run_off_lic_0(nlon)
54    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
55    REAL :: tsoil(nlon, nsoilmx, nbsrf)
56    REAL :: tslab(nlon), seaice(nlon)
57    REAL fder(nlon)
58
59
60
61    ! Arguments :
62    ! -----------
63
64    ! integer radpas
65    INTEGER it, unit, i, k, itap
66
67    REAL rugos, albedo
68    REAL tsurf
69    REAL time, timestep, day, day0
70    REAL qsol_f
71    REAL rugsrel(nlon)
72    LOGICAL alb_ocean
73
74    CHARACTER *80 ans, file_forctl, file_fordat, file_start
75    CHARACTER *100 file, var
76    CHARACTER *2 cnbl
77
78    REAL phy_nat(nlon, year_len)
79    REAL phy_alb(nlon, year_len)
80    REAL phy_sst(nlon, year_len)
81    REAL phy_bil(nlon, year_len)
82    REAL phy_rug(nlon, year_len)
83    REAL phy_ice(nlon, year_len)
84    REAL phy_fter(nlon, year_len)
85    REAL phy_foce(nlon, year_len)
86    REAL phy_fsic(nlon, year_len)
87    REAL phy_flic(nlon, year_len)
88
89    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
90!$OMP THREADPRIVATE(read_climoz)
91
92    ! -------------------------------------------------------------------------
93    ! declaration pour l'appel a phyredem
94    ! -------------------------------------------------------------------------
95
96    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
97    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
98    ! real pbl_tke(nlon,llm,nbsrf)
99    ! real rain_fall(nlon),snow_fall(nlon)
100    ! real solsw(nlon), sollw(nlon),radsol(nlon)
101    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
102    ! real ratqs(nlon,llm)
103    ! real clwcon(nlon,llm)
104
105    INTEGER longcles
106    PARAMETER (longcles=20)
107    REAL clesphy0(longcles)
108
109
110    ! -----------------------------------------------------------------------
111    ! dynamial tendencies :
112    ! ---------------------
113
114    INTEGER l, ierr, aslun
115
116    REAL paire
117
118    ! Local
119    CHARACTER (LEN=20) :: modname='phyaqua'
120    CHARACTER (LEN=80) :: abort_message
121
122
123    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
124    ! INITIALISATIONS
125    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
126
127    ! -----------------------------------------------------------------------
128    ! Initialisations  des constantes
129    ! -------------------------------
130
131    !IF (calend .EQ. "earth_360d") Then
132      year_len_phy = year_len
133    !END IF
134   
135    if (year_len.ne.360) then
136      write (*,*) year_len
137      call abort_physic("iniaqua", 'iniaqua: 360 day calendar is required !', 1)
138    endif
139
140    type_aqua = iflag_phys/100
141    type_profil = iflag_phys - type_aqua*100
142    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
143
144    IF (klon/=nlon) THEN
145      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
146      abort_message= 'probleme de dimensions dans iniaqua'
147      CALL abort_physic(modname,abort_message,1)
148    END IF
149    CALL phys_state_var_init(read_climoz)
150
151
152    read_climoz = 0
153    day0 = 217.
154    day = day0
155    it = 0
156    time = 0.
157
158    ! -----------------------------------------------------------------------
159    ! initialisations de la physique
160    ! -----------------------------------------------------------------------
161
162    day_ini = day_ref
163    day_end = day_ini + ndays
164
165    nbapp_rad = 24
166    CALL getin_p('nbapp_rad', nbapp_rad)
167
168    ! ---------------------------------------------------------------------
169    ! Creation des conditions aux limites:
170    ! ------------------------------------
171    ! Initialisations des constantes
172    ! Ajouter les manquants dans planete.def... (albedo etc)
173    co2_ppm = 348.
174    CALL getin_p('co2_ppm', co2_ppm)
175
176    solaire = 1365.
177    CALL getin_p('solaire', solaire)
178 
179    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
180    ! type_aqua
181    alb_ocean = .TRUE.
182    CALL getin_p('alb_ocean', alb_ocean)
183
184    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
185    WRITE (*, *) 'iniaqua: solaire=', solaire
186    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
187
188    radsol = 0.
189    qsol_f = 10.
190
191    ! Conditions aux limites:
192    ! -----------------------
193
194    qsol(:) = qsol_f
195    rugsrel = 0.0 ! (rugsrel = rugoro)
196    rugoro = 0.0
197    u_ancien = 0.0
198    v_ancien = 0.0
199    agesno = 50.0
200    ! Relief plat
201    zmea = 0.
202    zstd = 0.
203    zsig = 0.
204    zgam = 0.
205    zthe = 0.
206    zpic = 0.
207    zval = 0.
208
209    ! Une seule surface
210    pctsrf = 0.
211    IF (type_aqua==1) THEN
212      rugos = 1.E-4
213      albedo = 0.19
214      pctsrf(:, is_oce) = 1.
215    ELSE IF (type_aqua==2) THEN
216      rugos = 0.03
217      albedo = 0.1
218      pctsrf(:, is_ter) = 1.
219    END IF
220
221    CALL getin_p('rugos', rugos)
222
223    WRITE (*, *) 'iniaqua: rugos=', rugos
224    zmasq(:) = pctsrf(:, is_ter)
225
226    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
227    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
228
229    ! Si alb_ocean on calcule un albedo oceanique moyen
230    ! if (alb_ocean) then
231    ! Voir pourquoi on avait ca.
232    ! CALL ini_alb_oce(phy_alb)
233    ! else
234    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
235    ! endif !alb_ocean
236
237    DO i = 1, year_len
238      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
239      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
240
241      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
242      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
243      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
244      phy_ice(:, i) = 0.0 ! fraction de glace (?)
245      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
246      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
247      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
248      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
249    END DO
250    ! IM calcul profil sst
251    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
252
253    IF (grid_type==unstructured) THEN
254      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
255                             phy_fter, phy_foce, phy_flic, phy_fsic)
256    ELSE
257     
258       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
259                     phy_fter, phy_foce, phy_flic, phy_fsic)
260    ENDIF
261
262    ! ---------------------------------------------------------------------
263    ! Ecriture de l'etat initial:
264    ! ---------------------------
265
266
267    ! Ecriture etat initial physique
268
269    timestep = pdtphys
270    radpas = nint(rday/timestep/float(nbapp_rad))
271
272    DO i = 1, longcles
273      clesphy0(i) = 0.
274    END DO
275    clesphy0(1) = float(iflag_con)
276    clesphy0(2) = float(nbapp_rad)
277    ! IF( cycle_diurne  ) clesphy0(3) =  1.
278    clesphy0(3) = 1. ! cycle_diurne
279    clesphy0(4) = 1. ! soil_model
280    clesphy0(5) = 1. ! new_oliq
281    clesphy0(6) = 0. ! ok_orodr
282    clesphy0(7) = 0. ! ok_orolf
283    clesphy0(8) = 0. ! ok_limitvrai
284
285
286    ! =======================================================================
287    ! Profils initiaux
288    ! =======================================================================
289
290    ! On initialise les temperatures de surfaces comme les sst
291    DO i = 1, nlon
292      ftsol(i, :) = phy_sst(i, 1)
293      tsoil(i, :, :) = phy_sst(i, 1)
294      tslab(i) = phy_sst(i, 1)
295    END DO
296
297    falbe(:, :) = albedo
298    falblw(:, :) = albedo
299    rain_fall(:) = 0.
300    snow_fall(:) = 0.
301    solsw(:) = 0.
302    sollw(:) = 0.
303    radsol(:) = 0.
304
305    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
306    ! intialisation bidon mais pas grave
307    t_ancien(:, :) = 0.
308    q_ancien(:, :) = 0.
309    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310    rnebcon = 0.
311    ratqs = 0.
312    clwcon = 0.
313    pbl_tke = 1.E-8
314
315    ! variables supplementaires pour appel a plb_surface_init
316    fder(:) = 0.
317    seaice(:) = 0.
318    run_off_lic_0 = 0.
319    fevap = 0.
320
321
322    ! Initialisations necessaires avant phyredem
323    type_ocean = 'force'
324    CALL fonte_neige_init(run_off_lic_0)
325    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
326    snsrf(:, :) = 0. ! couverture de neige des sous surface
327    z0m(:, :) = rugos ! couverture de neige des sous surface
328    z0h=z0m
329
330
331    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
332
333    PRINT *, 'iniaqua: before phyredem'
334
335    pbl_tke(:,:,:) = 1.e-8
336    falb1 = albedo
337    falb2 = albedo
338    zmax0 = 0.
339    f0 = 0.
340    sig1 = 0.
341    w01 = 0.
342    wake_deltat = 0.
343    wake_deltaq = 0.
344    wake_s = 0.
345    wake_dens = 0.
346    wake_cstar = 0.
347    wake_pe = 0.
348    wake_fip = 0.
349    fm_therm = 0.
350    entr_therm = 0.
351    detr_therm = 0.
352    ale_bl = 0.
353    ale_bl_trig =0.
354    alp_bl =0.
355    treedrg(:,:,:)=0.
356
357    u10m = 0.
358    v10m = 0.
359
360    ql_ancien   = 0.
361    qs_ancien   = 0.
362    qbs_ancien  = 0.
363    u_ancien    = 0.
364    v_ancien    = 0.
365    prw_ancien  = 0.
366    prlw_ancien = 0.
367    prsw_ancien = 0. 
368    prbsw_ancien= 0.
369
370    ale_wake    = 0.
371    ale_bl_stat = 0. 
372
373
374!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
375!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
376!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
377    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
378    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
379
380!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
381!ym probably the uninitialized value was 0 for standard (regular grid) case
382    falb_dif(:,:,:)=0
383
384
385    CALL phyredem('startphy.nc')
386
387    PRINT *, 'iniaqua: after phyredem'
388    CALL phys_state_var_end
389
390    RETURN
391  END SUBROUTINE iniaqua
392
393
394  ! ====================================================================
395  ! ====================================================================
396  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
397    USE dimphy
398    IMPLICIT NONE
399    ! ====================================================================
400    ! =============================================================
401    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
402    ! Auteur : A. Campoy et F. Hourdin
403    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
404    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
405    ! connaissant la declinaison, la latitude et la longitude.
406
407    ! Dans cette version particuliere, on calcule le rayonnement
408    ! moyen sur l'année à chaque latitude.
409    ! angle zenithal calculé pour obtenir un
410    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
411    ! en moyenne annuelle.
412    ! Spécifique de la terre. Utilisé pour les aqua planetes.
413
414    ! Rque   : Different de la routine angle en ce sens que zenang
415    ! fournit des moyennes de pmu0 et non des valeurs
416    ! instantanees, du coup frac prend toutes les valeurs
417    ! entre 0 et 1.
418    ! Date   : premiere version le 13 decembre 1994
419    ! revu pour  GCM  le 30 septembre 1996
420    ! ===============================================================
421    ! longi----INPUT : la longitude vraie de la terre dans son plan
422    ! solaire a partir de l'equinoxe de printemps (degre)
423    ! gmtime---INPUT : temps universel en fraction de jour
424    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
425    ! lat------INPUT : latitude en degres
426    ! long-----INPUT : longitude en degres
427    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
428    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
429    ! ================================================================
430    include "YOMCST.h"
431    ! ================================================================
432    LOGICAL cycle_diurne
433    REAL gmtime
434    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
435    ! ================================================================
436    INTEGER i
437    REAL gmtime1, gmtime2
438    REAL pi_local
439
440
441    REAL rmu0m(klon), rmu0a(klon)
442
443
444    pi_local = 4.0*atan(1.0)
445
446    ! ================================================================
447    ! Calcul de l'angle zenithal moyen sur la journee
448    ! ================================================================
449
450    DO i = 1, klon
451      fract(i) = 1.
452      ! Calcule du flux moyen
453      IF (abs(rlat(i))<=28.75) THEN
454        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
455      ELSE IF (abs(rlat(i))<=43.75) THEN
456        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
457      ELSE IF (abs(rlat(i))<=71.25) THEN
458        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
459      ELSE
460        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
461      END IF
462    END DO
463
464    ! ================================================================
465    ! Avec ou sans cycle diurne
466    ! ================================================================
467
468    IF (cycle_diurne) THEN
469
470      ! On redecompose flux  au sommet suivant un cycle diurne idealise
471      ! identique a toutes les latitudes.
472
473      DO i = 1, klon
474        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
475        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
476          rmu0a(i)/sqrt(2.)
477      END DO
478
479      DO i = 1, klon
480        IF (rmu0(i)<=0.) THEN
481          rmu0(i) = 0.
482          fract(i) = 0.
483        ELSE
484          fract(i) = 1.
485        END IF
486      END DO
487
488      ! Affichage de l'angel zenitale
489      ! print*,'************************************'
490      ! print*,'************************************'
491      ! print*,'************************************'
492      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
493      ! print*,'rmu0m=',rmu0m(i)
494      ! print*,'rmu0a=',rmu0a(i)
495      ! print*,'rmu0=',rmu0(i)
496
497    ELSE
498
499      DO i = 1, klon
500        fract(i) = 0.5
501        rmu0(i) = rmu0m(i)*2.
502      END DO
503
504    END IF
505
506    RETURN
507  END SUBROUTINE zenang_an
508
509  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
510
511  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
512      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
513
514    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
515    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
516    USE lmdz_xios
517    IMPLICIT NONE
518
519    include "netcdf.inc"
520
521    INTEGER, INTENT (IN) :: klon
522    REAL, INTENT (IN) :: phy_nat(klon, 360)
523    REAL, INTENT (IN) :: phy_alb(klon, 360)
524    REAL, INTENT (IN) :: phy_sst(klon, 360)
525    REAL, INTENT (IN) :: phy_bil(klon, 360)
526    REAL, INTENT (IN) :: phy_rug(klon, 360)
527    REAL, INTENT (IN) :: phy_ice(klon, 360)
528    REAL, INTENT (IN) :: phy_fter(klon, 360)
529    REAL, INTENT (IN) :: phy_foce(klon, 360)
530    REAL, INTENT (IN) :: phy_flic(klon, 360)
531    REAL, INTENT (IN) :: phy_fsic(klon, 360)
532
533    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
534      ! on the whole physics grid
535 
536    IF (using_xios) THEN
537      PRINT *, 'writelim: Ecriture du fichier limit'
538
539      CALL gather_omp(phy_foce, phy_mpi)
540      IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
541
542      CALL gather_omp(phy_fsic, phy_mpi)
543      IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
544     
545      CALL gather_omp(phy_fter, phy_mpi)
546      IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
547     
548      CALL gather_omp(phy_flic, phy_mpi)
549      IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
550
551      CALL gather_omp(phy_sst, phy_mpi)
552      IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
553
554      CALL gather_omp(phy_bil, phy_mpi)
555      IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
556
557      CALL gather_omp(phy_alb, phy_mpi)
558      IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
559
560      CALL gather_omp(phy_rug, phy_mpi)
561      IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
562    ENDIF
563  END SUBROUTINE writelim_unstruct
564
565
566
567  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
568      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
569
570    USE mod_phys_lmdz_para, ONLY: is_master
571    USE mod_grid_phy_lmdz, ONLY: klon_glo
572    USE mod_phys_lmdz_transfert_para, ONLY: gather
573    USE phys_cal_mod, ONLY: year_len
574    use netcdf, only: nf90_def_var, nf90_double, nf90_float
575    IMPLICIT NONE
576    include "netcdf.inc"
577
578    INTEGER, INTENT (IN) :: klon
579    REAL, INTENT (IN) :: phy_nat(klon, year_len)
580    REAL, INTENT (IN) :: phy_alb(klon, year_len)
581    REAL, INTENT (IN) :: phy_sst(klon, year_len)
582    REAL, INTENT (IN) :: phy_bil(klon, year_len)
583    REAL, INTENT (IN) :: phy_rug(klon, year_len)
584    REAL, INTENT (IN) :: phy_ice(klon, year_len)
585    REAL, INTENT (IN) :: phy_fter(klon, year_len)
586    REAL, INTENT (IN) :: phy_foce(klon, year_len)
587    REAL, INTENT (IN) :: phy_flic(klon, year_len)
588    REAL, INTENT (IN) :: phy_fsic(klon, year_len)
589
590    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
591      ! on the whole physics grid
592    INTEGER :: k
593    INTEGER ierr
594    INTEGER dimfirst(3)
595    INTEGER dimlast(3)
596
597    INTEGER nid, ndim, ntim
598    INTEGER dims(2), debut(2), epais(2)
599    INTEGER id_tim
600    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
601    INTEGER id_fter, id_foce, id_fsic, id_flic
602
603    IF (is_master) THEN
604
605      PRINT *, 'writelim: Ecriture du fichier limit'
606
607      ierr = nf_create('limit.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
608
609      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
610        'Fichier conditions aux limites')
611      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
612      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
613      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
614
615      dims(1) = ndim
616      dims(2) = ntim
617
618#ifdef NC_DOUBLE
619      ierr = nf90_def_var(nid, 'TEMPS', nf90_double, [ntim], id_tim)
620#else
621      ierr = nf90_def_var(nid, 'TEMPS', nf90_float, [ntim], id_tim)
622#endif
623      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
624
625#ifdef NC_DOUBLE
626      ierr = nf90_def_var(nid, 'NAT', nf90_double, dims, id_nat)
627#else
628      ierr = nf90_def_var(nid, 'NAT', nf90_float, dims, id_nat)
629#endif
630      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
631        'Nature du sol (0,1,2,3)')
632
633#ifdef NC_DOUBLE
634      ierr = nf90_def_var(nid, 'SST', nf90_double, dims, id_sst)
635#else
636      ierr = nf90_def_var(nid, 'SST', nf90_float, dims, id_sst)
637#endif
638      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
639        'Temperature superficielle de la mer')
640
641#ifdef NC_DOUBLE
642      ierr = nf90_def_var(nid, 'BILS', nf90_double, dims, id_bils)
643#else
644      ierr = nf90_def_var(nid, 'BILS', nf90_float, dims, id_bils)
645#endif
646      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
647        'Reference flux de chaleur au sol')
648
649#ifdef NC_DOUBLE
650      ierr = nf90_def_var(nid, 'ALB', nf90_double, dims, id_alb)
651#else
652      ierr = nf90_def_var(nid, 'ALB', nf90_float, dims, id_alb)
653#endif
654      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
655
656#ifdef NC_DOUBLE
657      ierr = nf90_def_var(nid, 'RUG', nf90_double, dims, id_rug)
658#else
659      ierr = nf90_def_var(nid, 'RUG', nf90_float, dims, id_rug)
660#endif
661      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
662
663#ifdef NC_DOUBLE
664      ierr = nf90_def_var(nid, 'FTER', nf90_double, dims, id_fter)
665#else
666      ierr = nf90_def_var(nid, 'FTER', nf90_float, dims, id_fter)
667#endif
668      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
669#ifdef NC_DOUBLE
670      ierr = nf90_def_var(nid, 'FOCE', nf90_double, dims, id_foce)
671#else
672      ierr = nf90_def_var(nid, 'FOCE', nf90_float, dims, id_foce)
673#endif
674      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
675#ifdef NC_DOUBLE
676      ierr = nf90_def_var(nid, 'FSIC', nf90_double, dims, id_fsic)
677#else
678      ierr = nf90_def_var(nid, 'FSIC', nf90_float, dims, id_fsic)
679#endif
680      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
681#ifdef NC_DOUBLE
682      ierr = nf90_def_var(nid, 'FLIC', nf90_double, dims, id_flic)
683#else
684      ierr = nf90_def_var(nid, 'FLIC', nf90_float, dims, id_flic)
685#endif
686      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
687
688      ierr = nf_enddef(nid)
689      IF (ierr/=nf_noerr) THEN
690        WRITE (*, *) 'writelim error: failed to end define mode'
691        WRITE (*, *) nf_strerror(ierr)
692      END IF
693
694
695      ! write the 'times'
696      DO k = 1, year_len
697#ifdef NC_DOUBLE
698        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
699#else
700        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
701#endif
702        IF (ierr/=nf_noerr) THEN
703          WRITE (*, *) 'writelim error with temps(k),k=', k
704          WRITE (*, *) nf_strerror(ierr)
705        END IF
706      END DO
707
708    END IF ! of if (is_master)
709
710    ! write the fields, after having collected them on master
711
712    CALL gather(phy_nat, phy_glo)
713    IF (is_master) THEN
714#ifdef NC_DOUBLE
715      ierr = nf_put_var_double(nid, id_nat, phy_glo)
716#else
717      ierr = nf_put_var_real(nid, id_nat, phy_glo)
718#endif
719      IF (ierr/=nf_noerr) THEN
720        WRITE (*, *) 'writelim error with phy_nat'
721        WRITE (*, *) nf_strerror(ierr)
722      END IF
723    END IF
724
725    CALL gather(phy_sst, phy_glo)
726    IF (is_master) THEN
727#ifdef NC_DOUBLE
728      ierr = nf_put_var_double(nid, id_sst, phy_glo)
729#else
730      ierr = nf_put_var_real(nid, id_sst, phy_glo)
731#endif
732      IF (ierr/=nf_noerr) THEN
733        WRITE (*, *) 'writelim error with phy_sst'
734        WRITE (*, *) nf_strerror(ierr)
735      END IF
736    END IF
737
738    CALL gather(phy_bil, phy_glo)
739    IF (is_master) THEN
740#ifdef NC_DOUBLE
741      ierr = nf_put_var_double(nid, id_bils, phy_glo)
742#else
743      ierr = nf_put_var_real(nid, id_bils, phy_glo)
744#endif
745      IF (ierr/=nf_noerr) THEN
746        WRITE (*, *) 'writelim error with phy_bil'
747        WRITE (*, *) nf_strerror(ierr)
748      END IF
749    END IF
750
751    CALL gather(phy_alb, phy_glo)
752    IF (is_master) THEN
753#ifdef NC_DOUBLE
754      ierr = nf_put_var_double(nid, id_alb, phy_glo)
755#else
756      ierr = nf_put_var_real(nid, id_alb, phy_glo)
757#endif
758      IF (ierr/=nf_noerr) THEN
759        WRITE (*, *) 'writelim error with phy_alb'
760        WRITE (*, *) nf_strerror(ierr)
761      END IF
762    END IF
763
764    CALL gather(phy_rug, phy_glo)
765    IF (is_master) THEN
766#ifdef NC_DOUBLE
767      ierr = nf_put_var_double(nid, id_rug, phy_glo)
768#else
769      ierr = nf_put_var_real(nid, id_rug, phy_glo)
770#endif
771      IF (ierr/=nf_noerr) THEN
772        WRITE (*, *) 'writelim error with phy_rug'
773        WRITE (*, *) nf_strerror(ierr)
774      END IF
775    END IF
776
777    CALL gather(phy_fter, phy_glo)
778    IF (is_master) THEN
779#ifdef NC_DOUBLE
780      ierr = nf_put_var_double(nid, id_fter, phy_glo)
781#else
782      ierr = nf_put_var_real(nid, id_fter, phy_glo)
783#endif
784      IF (ierr/=nf_noerr) THEN
785        WRITE (*, *) 'writelim error with phy_fter'
786        WRITE (*, *) nf_strerror(ierr)
787      END IF
788    END IF
789
790    CALL gather(phy_foce, phy_glo)
791    IF (is_master) THEN
792#ifdef NC_DOUBLE
793      ierr = nf_put_var_double(nid, id_foce, phy_glo)
794#else
795      ierr = nf_put_var_real(nid, id_foce, phy_glo)
796#endif
797      IF (ierr/=nf_noerr) THEN
798        WRITE (*, *) 'writelim error with phy_foce'
799        WRITE (*, *) nf_strerror(ierr)
800      END IF
801    END IF
802
803    CALL gather(phy_fsic, phy_glo)
804    IF (is_master) THEN
805#ifdef NC_DOUBLE
806      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
807#else
808      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
809#endif
810      IF (ierr/=nf_noerr) THEN
811        WRITE (*, *) 'writelim error with phy_fsic'
812        WRITE (*, *) nf_strerror(ierr)
813      END IF
814    END IF
815
816    CALL gather(phy_flic, phy_glo)
817    IF (is_master) THEN
818#ifdef NC_DOUBLE
819      ierr = nf_put_var_double(nid, id_flic, phy_glo)
820#else
821      ierr = nf_put_var_real(nid, id_flic, phy_glo)
822#endif
823      IF (ierr/=nf_noerr) THEN
824        WRITE (*, *) 'writelim error with phy_flic'
825        WRITE (*, *) nf_strerror(ierr)
826      END IF
827    END IF
828
829    ! close file:
830    IF (is_master) THEN
831      ierr = nf_close(nid)
832    END IF
833
834  END SUBROUTINE writelim
835
836  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
837
838  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
839    USE dimphy
840    USE phys_cal_mod , ONLY: year_len
841    IMPLICIT NONE
842
843    INTEGER nlon, type_profil, i, k, j
844    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
845    INTEGER imn, imx, amn, amx, kmn, kmx
846    INTEGER p, pplus, nlat_max
847    PARAMETER (nlat_max=72)
848    REAL x_anom_sst(nlat_max)
849    CHARACTER (LEN=20) :: modname='profil_sst'
850    CHARACTER (LEN=80) :: abort_message
851
852    IF (klon/=nlon) THEN
853       abort_message='probleme de dimensions dans profil_sst'
854       CALL abort_physic(modname,abort_message,1)
855    ENDIF
856    WRITE (*, *) ' profil_sst: type_profil=', type_profil
857    DO i = 1, year_len
858      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
859
860      ! Rajout fbrlmd
861
862      IF (type_profil==1) THEN
863        ! Méthode 1 "Control" faible plateau à l'Equateur
864        DO j = 1, klon
865          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
866          ! PI/3=1.047197551
867          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
868            phy_sst(j, i) = 273.
869          END IF
870        END DO
871      END IF
872      IF (type_profil==2) THEN
873        ! Méthode 2 "Flat" fort plateau à l'Equateur
874        DO j = 1, klon
875          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
876          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
877            phy_sst(j, i) = 273.
878          END IF
879        END DO
880      END IF
881
882
883      IF (type_profil==3) THEN
884        ! Méthode 3 "Qobs" plateau réel à l'Equateur
885        DO j = 1, klon
886          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
887            rlatd(j))**4)
888          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
889            phy_sst(j, i) = 273.
890          END IF
891        END DO
892      END IF
893
894      IF (type_profil==4) THEN
895        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
896        DO j = 1, klon
897          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
898            rlatd(j))**4)
899          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
900            phy_sst(j, i) = 273.
901          END IF
902        END DO
903      END IF
904
905      IF (type_profil==5) THEN
906        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
907        DO j = 1, klon
908          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
909            *rlatd(j))**4)
910          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
911            phy_sst(j, i) = 273. + 2.
912          END IF
913
914        END DO
915      END IF
916
917      IF (type_profil==6) THEN
918        ! Méthode 6 "cst" valeur constante de SST
919        DO j = 1, klon
920          phy_sst(j, i) = 288.
921        END DO
922      END IF
923
924
925      IF (type_profil==7) THEN
926        ! Méthode 7 "cst" valeur constante de SST +2
927        DO j = 1, klon
928          phy_sst(j, i) = 288. + 2.
929        END DO
930      END IF
931
932      p = 0
933      IF (type_profil==8) THEN
934        ! Méthode 8 profil anomalies SST du modèle couplé AR4
935        DO j = 1, klon
936          IF (rlatd(j)==rlatd(j-1)) THEN
937            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
938              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
939            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
940              phy_sst(j, i) = 273. + x_anom_sst(pplus)
941            END IF
942          ELSE
943            p = p + 1
944            pplus = 73 - p
945            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
946              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
947            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
948              phy_sst(j, i) = 273. + x_anom_sst(pplus)
949            END IF
950            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
951          END IF
952        END DO
953      END IF
954
955      IF (type_profil==9) THEN
956        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
957        DO j = 1, klon
958          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
959            *rlatd(j))**4)
960          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
961            phy_sst(j, i) = 273. - 2.
962          END IF
963        END DO
964      END IF
965
966
967      IF (type_profil==10) THEN
968        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
969        DO j = 1, klon
970          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
971            *rlatd(j))**4)
972          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
973            phy_sst(j, i) = 273. + 4.
974          END IF
975        END DO
976      END IF
977
978      IF (type_profil==11) THEN
979        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
980        DO j = 1, klon
981          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
982            rlatd(j))**4)
983          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
984            phy_sst(j, i) = 273.
985          END IF
986        END DO
987      END IF
988
989      IF (type_profil==12) THEN
990        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
991        DO j = 1, klon
992          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
993            *rlatd(j))**4)
994          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
995            phy_sst(j, i) = 273. + 4.
996          END IF
997        END DO
998      END IF
999
1000      IF (type_profil==13) THEN
1001        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
1002        DO j = 1, klon
1003          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
1004            rlatd(j))**4)
1005          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1006            phy_sst(j, i) = 273.
1007          END IF
1008        END DO
1009      END IF
1010
1011      IF (type_profil==14) THEN
1012        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
1013        DO j = 1, klon
1014          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
1015            *rlatd(j))**4)
1016          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1017            phy_sst(j, i) = 273.
1018          END IF
1019        END DO
1020      END IF
1021
1022      if (type_profil.EQ.20) then
1023      print*,'Profile SST 20'
1024!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1025
1026      do j=1,klon
1027        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
1028      enddo
1029      endif
1030
1031      if (type_profil.EQ.21) then
1032      print*,'Profile SST 21'
1033!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1034      do j=1,klon
1035        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
1036      enddo
1037      endif
1038
1039
1040
1041    END DO
1042
1043    ! IM beg : verif profil SST: phy_sst
1044    amn = min(phy_sst(1,1), 1000.)
1045    amx = max(phy_sst(1,1), -1000.)
1046    imn = 1
1047    kmn = 1
1048    imx = 1
1049    kmx = 1
1050    DO k = 1, year_len
1051      DO i = 2, nlon
1052        IF (phy_sst(i,k)<amn) THEN
1053          amn = phy_sst(i, k)
1054          imn = i
1055          kmn = k
1056        END IF
1057        IF (phy_sst(i,k)>amx) THEN
1058          amx = phy_sst(i, k)
1059          imx = i
1060          kmx = k
1061        END IF
1062      END DO
1063    END DO
1064
1065    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1066    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1067    ! IM end : verif profil SST: phy_sst
1068
1069    RETURN
1070  END SUBROUTINE profil_sst
1071
1072END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.